(git:1f9fd2c)
Loading...
Searching...
No Matches
molden_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Functions handling the MOLDEN format. Split from mode_selective.
10!> \author Teodoro Laino, 03.2009
11! **************************************************************************************************
12MODULE molden_utils
13 USE admm_types, ONLY: admm_type
14 USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
15 admm_uncorrect_for_eigenvalues
16 USE atomic_kind_types, ONLY: get_atomic_kind
17 USE basis_set_types, ONLY: get_gto_basis_set,&
18 gto_basis_set_type
19 USE cell_types, ONLY: cell_type
20 USE cp_array_utils, ONLY: cp_1d_r_p_type
21 USE cp_control_types, ONLY: dft_control_type
22 USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
23 dbcsr_type
24 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
25 USE cp_fm_types, ONLY: cp_fm_get_info,&
26 cp_fm_get_submatrix,&
27 cp_fm_type
28 USE cp_log_handling, ONLY: cp_get_default_logger,&
29 cp_logger_type
30 USE cp_output_handling, ONLY: cp_p_file,&
31 cp_print_key_finished_output,&
32 cp_print_key_should_output,&
33 cp_print_key_unit_nr
34 USE input_constants, ONLY: gto_cartesian,&
35 gto_spherical
36 USE input_section_types, ONLY: section_vals_type,&
37 section_vals_val_get
38 USE kinds, ONLY: dp
39 USE mathconstants, ONLY: pi
40 USE orbital_pointers, ONLY: nco,&
41 nso
42 USE orbital_transformation_matrices, ONLY: orbtramat
43 USE particle_types, ONLY: particle_type
44 USE periodic_table, ONLY: get_ptable_info
45 USE physcon, ONLY: angstrom,&
46 massunit
47 USE qs_environment_types, ONLY: get_qs_env,&
48 qs_environment_type
49 USE qs_kind_types, ONLY: get_qs_kind,&
50 get_qs_kind_set,&
51 qs_kind_type
52 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
53 USE qs_mo_types, ONLY: get_mo_set,&
54 mo_set_type
55#include "./base/base_uses.f90"
56
57 IMPLICIT NONE
58
59 PRIVATE
60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molden_utils'
61 LOGICAL, PARAMETER :: debug_this_module = .false.
62
63 INTEGER, PARAMETER :: molden_lmax = 4
64 INTEGER, PARAMETER :: molden_ncomax = (molden_lmax + 1)*(molden_lmax + 2)/2 ! 15
65
66 PUBLIC :: write_vibrations_molden, write_mos_molden
67
68CONTAINS
69
70! **************************************************************************************************
71!> \brief Write out the MOs in molden format for visualisation
72!> \param mos the set of MOs (both spins, if UKS)
73!> \param qs_kind_set for basis set info
74!> \param particle_set particles data structure, for positions and kinds
75!> \param print_section input section containing relevant print key
76!> \param cell ...
77!> \param unoccupied_orbs optional: unoccupied orbital coefficients from make_lumo_gpw
78!> \param unoccupied_evals optional: unoccupied orbital eigenvalues
79!> \param qs_env ...
80!> \param calc_energies ...
81!> \author MattW, IainB
82! **************************************************************************************************
83 SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section, cell, &
84 unoccupied_orbs, unoccupied_evals, qs_env, calc_energies)
85 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
86 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
87 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
88 TYPE(section_vals_type), POINTER :: print_section
89 TYPE(cell_type), OPTIONAL, POINTER :: cell
90 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
91 OPTIONAL :: unoccupied_orbs
92 TYPE(cp_1d_r_p_type), DIMENSION(:), INTENT(IN), &
93 OPTIONAL :: unoccupied_evals
94 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
95 LOGICAL, INTENT(IN), OPTIONAL :: calc_energies
96
97 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_mos_molden'
98 CHARACTER(LEN=molden_lmax+1), PARAMETER :: angmom = "spdfg"
99
100 CHARACTER(LEN=15) :: fmtstr1, fmtstr2
101 CHARACTER(LEN=2) :: element_symbol
102 INTEGER :: gto_kind, handle, i, iatom, icgf, icol, ikind, ipgf, irow, irow_in, iset, isgf, &
103 ishell, ispin, iw, lshell, ncgf, ncol_global, ndigits, nrow_global, nset, nsgf, numos, &
104 unit_choice, z
105 INTEGER, DIMENSION(:), POINTER :: npgf, nshell
106 INTEGER, DIMENSION(:, :), POINTER :: l
107 INTEGER, DIMENSION(molden_ncomax, 0:molden_lmax) :: orbmap
108 LOGICAL :: do_calc_energies, ghost_atom, &
109 mark_ghost, print_warn, write_cell, &
110 write_pseudo
111 REAL(kind=dp) :: expzet, prefac, scale_factor, zeff
112 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
113 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
114 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
115 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
116 TYPE(admm_type), POINTER :: admm_env
117 TYPE(cp_fm_type), POINTER :: mo_coeff
118 TYPE(cp_logger_type), POINTER :: logger
119 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks
120 TYPE(dbcsr_type), POINTER :: matrix_ks, mo_coeff_deriv
121 TYPE(dft_control_type), POINTER :: dft_control
122 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
123
124 CALL timeset(routinen, handle)
125
126 logger => cp_get_default_logger()
127 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, ""), cp_p_file)) THEN
128
129 iw = cp_print_key_unit_nr(logger, print_section, "", &
130 extension=".molden", file_status='REPLACE')
131
132 print_warn = .true.
133
134 CALL section_vals_val_get(print_section, "UNIT", i_val=unit_choice)
135 IF (unit_choice == 2) THEN
136 scale_factor = angstrom
137 ELSE
138 scale_factor = 1.0_dp
139 END IF
140
141 CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
142 ndigits = min(max(3, ndigits), 30)
143 WRITE (unit=fmtstr1, fmt='("(I6,1X,ES",I0,".",I0,")")') ndigits + 7, ndigits
144 WRITE (unit=fmtstr2, fmt='("((T51,2F",I0,".",I0,"))")') ndigits + 10, ndigits
145
146 CALL section_vals_val_get(print_section, "GTO_KIND", i_val=gto_kind)
147 CALL section_vals_val_get(print_section, "WRITE_CELL", l_val=write_cell)
148 CALL section_vals_val_get(print_section, "WRITE_PSEUDO", l_val=write_pseudo)
149 CALL section_vals_val_get(print_section, "MARK_GHOST", l_val=mark_ghost)
150
151 IF (mos(1)%use_mo_coeff_b) THEN
152 ! we are using the dbcsr mo_coeff
153 ! we copy it to the fm anyway
154 DO ispin = 1, SIZE(mos)
155 IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
156 cpassert(.false.)
157 END IF
158 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
159 mos(ispin)%mo_coeff) !fm->dbcsr
160 END DO
161 END IF
162
163 IF (iw > 0) THEN
164 WRITE (iw, '(T2,A)') "[Molden Format]"
165 IF (write_cell) THEN
166 IF (unit_choice == 2) THEN
167 WRITE (iw, '(T2,A)') "[Cell] Angs"
168 ELSE
169 WRITE (iw, '(T2,A)') "[Cell] AU"
170 END IF
171 WRITE (iw, '(T2,3(F12.6,3X))') &
172 cell%hmat(1, 1)*scale_factor, cell%hmat(2, 1)*scale_factor, cell%hmat(3, 1)*scale_factor
173 WRITE (iw, '(T2,3(F12.6,3X))') &
174 cell%hmat(1, 2)*scale_factor, cell%hmat(2, 2)*scale_factor, cell%hmat(3, 2)*scale_factor
175 WRITE (iw, '(T2,3(F12.6,3X))') &
176 cell%hmat(1, 3)*scale_factor, cell%hmat(2, 3)*scale_factor, cell%hmat(3, 3)*scale_factor
177 END IF
178 IF (unit_choice == 2) THEN
179 WRITE (iw, '(T2,A)') "[Atoms] Angs"
180 ELSE
181 WRITE (iw, '(T2,A)') "[Atoms] AU"
182 END IF
183 DO i = 1, SIZE(particle_set)
184 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
185 element_symbol=element_symbol)
186 CALL get_ptable_info(element_symbol, number=z)
187 IF (mark_ghost) THEN
188 CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_atom)
189 IF (ghost_atom) z = 0
190 END IF
191
192 WRITE (iw, '(T2,A2,I6,I6,3X,3(F12.6,3X))') &
193 element_symbol, i, z, particle_set(i)%r(:)*scale_factor
194 END DO
195 IF (write_pseudo) THEN
196 WRITE (iw, '(T2,A)') "[Pseudo]"
197 DO i = 1, SIZE(particle_set)
198 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
199 element_symbol=element_symbol)
200 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
201 WRITE (iw, '(T2,A2,I6,I6)') &
202 element_symbol, i, nint(zeff)
203 END DO
204 END IF
205
206 WRITE (iw, '(T2,A)') "[GTO]"
207
208 DO i = 1, SIZE(particle_set)
209 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
210 element_symbol=element_symbol)
211 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
212 IF (ASSOCIATED(orb_basis_set)) THEN
213 WRITE (iw, '(T2,I8,I8)') i, 0
214 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
215 nset=nset, &
216 npgf=npgf, &
217 nshell=nshell, &
218 l=l, &
219 zet=zet, &
220 gcc=gcc)
221
222 DO iset = 1, nset
223 DO ishell = 1, nshell(iset)
224 lshell = l(ishell, iset)
225 IF (lshell <= molden_lmax) THEN
226 WRITE (unit=iw, fmt='(T25,A2,4X,I4,4X,F4.2)') &
227 angmom(lshell + 1:lshell + 1), npgf(iset), 1.0_dp
228 ! MOLDEN expects the contraction coefficient of spherical NOT CARTESIAN NORMALISED
229 ! functions. So we undo the normalisation factors included in the gccs
230 ! Reverse engineered from basis_set_types, normalise_gcc_orb
231 prefac = 2_dp**lshell*(2/pi)**0.75_dp
232 expzet = 0.25_dp*(2*lshell + 3.0_dp)
233 WRITE (unit=iw, fmt=fmtstr2) &
234 (zet(ipgf, iset), gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet), &
235 ipgf=1, npgf(iset))
236 ELSE
237 IF (print_warn) THEN
238 CALL cp_warn(__location__, &
239 "MOLDEN format does not support Gaussian orbitals with l > 4.")
240 print_warn = .false.
241 END IF
242 END IF
243 END DO
244 END DO
245
246 WRITE (iw, '(A4)') " "
247
248 END IF
249
250 END DO
251
252 IF (gto_kind == gto_spherical) THEN
253 WRITE (iw, '(T2,A)') "[5D7F]"
254 WRITE (iw, '(T2,A)') "[9G]"
255 END IF
256
257 WRITE (iw, '(T2,A)') "[MO]"
258 END IF
259
260 !------------------------------------------------------------------------
261 ! convert from CP2K to MOLDEN format ordering
262 ! http://www.cmbi.ru.nl/molden/molden_format.html
263 !"The following order of D, F and G functions is expected:
264 !
265 ! 5D: D 0, D+1, D-1, D+2, D-2
266 ! 6D: xx, yy, zz, xy, xz, yz
267 !
268 ! 7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
269 ! 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
270 !
271 ! 9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
272 ! 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
273 ! xxyy xxzz yyzz xxyz yyxz zzxy
274 !"
275 ! CP2K has x in the outer (slower loop), so
276 ! xx, xy, xz, yy, yz,zz for l=2, for instance
277 !
278 ! iorb_cp2k = orbmap(iorb_molden, l), l = 0 .. 4
279 ! -----------------------------------------------------------------------
280 IF (iw > 0) THEN
281 IF (gto_kind == gto_cartesian) THEN
282 ! -----------------------------------------------------------------
283 ! Use cartesian (6D, 10F, 15G) representation.
284 ! This is only format VMD can process.
285 ! -----------------------------------------------------------------
286 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
287 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
288 1, 4, 6, 2, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
289 1, 7, 10, 4, 2, 3, 6, 9, 8, 5, 0, 0, 0, 0, 0, &
290 1, 11, 15, 2, 3, 7, 12, 10, 14, 4, 6, 13, 5, 8, 9], &
291 [molden_ncomax, molden_lmax + 1])
292 ELSE IF (gto_kind == gto_spherical) THEN
293 ! -----------------------------------------------------------------
294 ! Use spherical (5D, 7F, 9G) representation.
295 ! -----------------------------------------------------------------
296 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
297 3, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
298 3, 4, 2, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
299 4, 5, 3, 6, 2, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
300 5, 6, 4, 7, 3, 8, 2, 9, 1, 0, 0, 0, 0, 0, 0], &
301 [molden_ncomax, molden_lmax + 1])
302 END IF
303 END IF
304
305 DO ispin = 1, SIZE(mos)
306 do_calc_energies = .false.
307 IF (PRESENT(calc_energies)) do_calc_energies = calc_energies
308
309 IF (PRESENT(qs_env) .AND. do_calc_energies) THEN
310 CALL get_qs_env(qs_env, matrix_ks=ks, dft_control=dft_control)
311
312 matrix_ks => ks(ispin)%matrix
313
314 ! With ADMM, we have to modify the Kohn-Sham matrix
315 IF (dft_control%do_admm) THEN
316 CALL get_qs_env(qs_env, admm_env=admm_env)
317 CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks)
318 END IF
319
320 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues)
321
322 IF (ASSOCIATED(qs_env%mo_derivs)) THEN
323 mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
324 ELSE
325 mo_coeff_deriv => null()
326 END IF
327
328 ! Update the eigenvalues of the occupied orbitals
329 CALL calculate_subspace_eigenvalues(orbitals=mo_coeff, &
330 ks_matrix=matrix_ks, &
331 evals_arg=mo_eigenvalues, &
332 co_rotate_dbcsr=mo_coeff_deriv)
333
334 ! With ADMM, we have to undo the modification of the Kohn-Sham matrix
335 IF (dft_control%do_admm) THEN
336 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks)
337 END IF
338 END IF
339
340 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
341 nrow_global=nrow_global, &
342 ncol_global=ncol_global)
343 ALLOCATE (smatrix(nrow_global, ncol_global))
344 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
345
346 IF (iw > 0) THEN
347 IF (gto_kind == gto_cartesian) THEN
348 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
349
350 ALLOCATE (cmatrix(ncgf, ncgf))
351
352 cmatrix = 0.0_dp
353
354 ! Transform spherical MOs to Cartesian MOs
355
356 icgf = 1
357 isgf = 1
358 DO iatom = 1, SIZE(particle_set)
359 NULLIFY (orb_basis_set)
360 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
361 CALL get_qs_kind(qs_kind_set(ikind), &
362 basis_set=orb_basis_set)
363 IF (ASSOCIATED(orb_basis_set)) THEN
364 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
365 nset=nset, &
366 nshell=nshell, &
367 l=l)
368 DO iset = 1, nset
369 DO ishell = 1, nshell(iset)
370 lshell = l(ishell, iset)
371 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, nso(lshell), 1.0_dp, &
372 orbtramat(lshell)%c2s, nso(lshell), &
373 smatrix(isgf, 1), nsgf, 0.0_dp, &
374 cmatrix(icgf, 1), ncgf)
375 icgf = icgf + nco(lshell)
376 isgf = isgf + nso(lshell)
377 END DO
378 END DO
379 END IF
380 END DO ! iatom
381 END IF
382
383 DO icol = 1, mos(ispin)%nmo
384 ! index of the first basis function for the given atom, set, and shell
385 irow = 1
386
387 ! index of the first basis function in MOLDEN file.
388 ! Due to limitation of the MOLDEN format, basis functions with l > molden_lmax
389 ! cannot be exported, so we need to renumber atomic orbitals
390 irow_in = 1
391
392 WRITE (iw, '(A,ES20.10)') 'Ene=', mos(ispin)%eigenvalues(icol)
393 IF (ispin < 2) THEN
394 WRITE (iw, '(A)') 'Spin= Alpha'
395 ELSE
396 WRITE (iw, '(A)') 'Spin= Beta'
397 END IF
398 WRITE (iw, '(A,F12.7)') 'Occup=', mos(ispin)%occupation_numbers(icol)
399
400 DO iatom = 1, SIZE(particle_set)
401 NULLIFY (orb_basis_set)
402 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
403 element_symbol=element_symbol, kind_number=ikind)
404 CALL get_qs_kind(qs_kind_set(ikind), &
405 basis_set=orb_basis_set)
406 IF (ASSOCIATED(orb_basis_set)) THEN
407 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
408 nset=nset, &
409 nshell=nshell, &
410 l=l)
411
412 IF (gto_kind == gto_cartesian) THEN
413 ! ----------------------------------------------
414 ! Use cartesian (6D, 10F, 15G) representation.
415 ! ----------------------------------------------
416 icgf = 1
417 DO iset = 1, nset
418 DO ishell = 1, nshell(iset)
419 lshell = l(ishell, iset)
420
421 IF (lshell <= molden_lmax) THEN
422 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
423 cmatrix(irow:irow + nco(lshell) - 1, icol))
424 irow_in = irow_in + nco(lshell)
425 END IF
426
427 irow = irow + nco(lshell)
428 END DO ! ishell
429 END DO
430
431 ELSE IF (gto_kind == gto_spherical) THEN
432 ! ----------------------------------------------
433 ! Use spherical (5D, 7F, 9G) representation.
434 ! ----------------------------------------------
435 DO iset = 1, nset
436 DO ishell = 1, nshell(iset)
437 lshell = l(ishell, iset)
438
439 IF (lshell <= molden_lmax) THEN
440 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
441 smatrix(irow:irow + nso(lshell) - 1, icol))
442 irow_in = irow_in + nso(lshell)
443 END IF
444
445 irow = irow + nso(lshell)
446 END DO
447 END DO
448 END IF
449
450 END IF
451 END DO ! iatom
452 END DO
453 END IF
454
455 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
456 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
457 END DO
458
459 ! Write unoccupied (virtual) orbitals if provided; only used with OT
460 IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN
461 DO ispin = 1, SIZE(unoccupied_orbs)
462 CALL cp_fm_get_info(unoccupied_orbs(ispin), &
463 nrow_global=nrow_global, &
464 ncol_global=numos)
465 ALLOCATE (smatrix(nrow_global, numos))
466 CALL cp_fm_get_submatrix(unoccupied_orbs(ispin), smatrix)
467
468 IF (iw > 0) THEN
469 IF (gto_kind == gto_cartesian) THEN
470 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
471 ALLOCATE (cmatrix(ncgf, numos))
472 cmatrix = 0.0_dp
473
474 icgf = 1
475 isgf = 1
476 DO iatom = 1, SIZE(particle_set)
477 NULLIFY (orb_basis_set)
478 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
479 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
480 IF (ASSOCIATED(orb_basis_set)) THEN
481 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
482 nset=nset, nshell=nshell, l=l)
483 DO iset = 1, nset
484 DO ishell = 1, nshell(iset)
485 lshell = l(ishell, iset)
486 CALL dgemm("T", "N", nco(lshell), numos, nso(lshell), 1.0_dp, &
487 orbtramat(lshell)%c2s, nso(lshell), &
488 smatrix(isgf, 1), nsgf, 0.0_dp, &
489 cmatrix(icgf, 1), ncgf)
490 icgf = icgf + nco(lshell)
491 isgf = isgf + nso(lshell)
492 END DO
493 END DO
494 END IF
495 END DO
496 END IF
497
498 DO icol = 1, numos
499 irow = 1
500 irow_in = 1
501
502 WRITE (iw, '(A,ES20.10)') 'Ene=', unoccupied_evals(ispin)%array(icol)
503 IF (ispin < 2) THEN
504 WRITE (iw, '(A)') 'Spin= Alpha'
505 ELSE
506 WRITE (iw, '(A)') 'Spin= Beta'
507 END IF
508 WRITE (iw, '(A,F12.7)') 'Occup=', 0.0_dp
509
510 DO iatom = 1, SIZE(particle_set)
511 NULLIFY (orb_basis_set)
512 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
513 element_symbol=element_symbol, kind_number=ikind)
514 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
515 IF (ASSOCIATED(orb_basis_set)) THEN
516 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
517 nset=nset, nshell=nshell, l=l)
518
519 IF (gto_kind == gto_cartesian) THEN
520 icgf = 1
521 DO iset = 1, nset
522 DO ishell = 1, nshell(iset)
523 lshell = l(ishell, iset)
524 IF (lshell <= molden_lmax) THEN
525 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
526 cmatrix(irow:irow + nco(lshell) - 1, icol))
527 irow_in = irow_in + nco(lshell)
528 END IF
529 irow = irow + nco(lshell)
530 END DO
531 END DO
532 ELSE IF (gto_kind == gto_spherical) THEN
533 DO iset = 1, nset
534 DO ishell = 1, nshell(iset)
535 lshell = l(ishell, iset)
536 IF (lshell <= molden_lmax) THEN
537 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
538 smatrix(irow:irow + nso(lshell) - 1, icol))
539 irow_in = irow_in + nso(lshell)
540 END IF
541 irow = irow + nso(lshell)
542 END DO
543 END DO
544 END IF
545
546 END IF
547 END DO ! iatom
548 END DO ! icol
549 END IF
550
551 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
552 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
553 END DO ! ispin
554 END IF
555
556 CALL cp_print_key_finished_output(iw, logger, print_section, "")
557
558 END IF
559
560 CALL timestop(handle)
561
562 END SUBROUTINE write_mos_molden
563
564! **************************************************************************************************
565!> \brief Output MO coefficients formatted correctly for MOLDEN, omitting those <= 1E(-digits)
566!> \param iw output file unit
567!> \param fmtstr1 format string
568!> \param ndigits number of significant digits in MO coefficients
569!> \param irow_in index of the first atomic orbital: mo_coeff(orbmap(1))
570!> \param orbmap array to map Gaussian functions from MOLDEN to CP2K ordering
571!> \param mo_coeff MO coefficients
572! **************************************************************************************************
573 SUBROUTINE print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap, mo_coeff)
574 INTEGER, INTENT(in) :: iw
575 CHARACTER(LEN=*), INTENT(in) :: fmtstr1
576 INTEGER, INTENT(in) :: ndigits, irow_in
577 INTEGER, DIMENSION(molden_ncomax), INTENT(in) :: orbmap
578 REAL(kind=dp), DIMENSION(:), INTENT(in) :: mo_coeff
579
580 INTEGER :: orbital
581
582 DO orbital = 1, molden_ncomax
583 IF (orbmap(orbital) /= 0) THEN
584 IF (abs(mo_coeff(orbmap(orbital))) >= 10.0_dp**(-ndigits)) THEN
585 WRITE (iw, fmtstr1) irow_in + orbital - 1, mo_coeff(orbmap(orbital))
586 END IF
587 END IF
588 END DO
589
590 END SUBROUTINE print_coeffs
591
592! **************************************************************************************************
593!> \brief writes the output for vibrational analysis in MOLDEN format
594!> \param input ...
595!> \param particles ...
596!> \param freq ...
597!> \param eigen_vec ...
598!> \param intensities ...
599!> \param calc_intens ...
600!> \param dump_only_positive ...
601!> \param logger ...
602!> \param list array of mobile atom indices
603!> \author Florian Schiffmann 11.2007
604! **************************************************************************************************
605 SUBROUTINE write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, &
606 dump_only_positive, logger, list)
607
608 TYPE(section_vals_type), POINTER :: input
609 TYPE(particle_type), DIMENSION(:), POINTER :: particles
610 REAL(kind=dp), DIMENSION(:) :: freq
611 REAL(kind=dp), DIMENSION(:, :) :: eigen_vec
612 REAL(kind=dp), DIMENSION(:), POINTER :: intensities
613 LOGICAL, INTENT(in) :: calc_intens, dump_only_positive
614 TYPE(cp_logger_type), POINTER :: logger
615 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: list
616
617 CHARACTER(len=*), PARAMETER :: routinen = 'write_vibrations_molden'
618
619 CHARACTER(LEN=2) :: element_symbol
620 INTEGER :: handle, i, iw, j, k, l, z
621 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
622 REAL(kind=dp) :: fint
623
624 CALL timeset(routinen, handle)
625
626 iw = cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
627 extension=".mol", file_status='REPLACE')
628
629 IF (iw > 0) THEN
630 cpassert(mod(SIZE(eigen_vec, 1), 3) == 0)
631 cpassert(SIZE(freq, 1) == SIZE(eigen_vec, 2))
632 ALLOCATE (my_list(SIZE(particles)))
633 ! Either we have a list of the subset of mobile atoms,
634 ! Or the eigenvectors must span the full space (all atoms)
635 IF (PRESENT(list)) THEN
636 my_list(:) = 0
637 DO i = 1, SIZE(list)
638 my_list(list(i)) = i
639 END DO
640 ELSE
641 cpassert(SIZE(particles) == SIZE(eigen_vec, 1)/3)
642 DO i = 1, SIZE(my_list)
643 my_list(i) = i
644 END DO
645 END IF
646 WRITE (iw, '(T2,A)') "[Molden Format]"
647 WRITE (iw, '(T2,A)') "[Atoms] AU"
648 DO i = 1, SIZE(particles)
649 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
650 element_symbol=element_symbol)
651 CALL get_ptable_info(element_symbol, number=z)
652
653 WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
654 element_symbol, i, z, particles(i)%r(:)
655
656 END DO
657 WRITE (iw, '(T2,A)') "[FREQ]"
658 DO i = 1, SIZE(freq, 1)
659 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(T5,F12.6)') freq(i)
660 END DO
661 WRITE (iw, '(T2,A)') "[FR-COORD]"
662 DO i = 1, SIZE(particles)
663 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
664 element_symbol=element_symbol)
665 WRITE (iw, '(T2,A2,3X,3(F12.6,3X))') &
666 element_symbol, particles(i)%r(:)
667 END DO
668 WRITE (iw, '(T2,A)') "[FR-NORM-COORD]"
669 l = 0
670 DO i = 1, SIZE(eigen_vec, 2)
671 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) THEN
672 l = l + 1
673 WRITE (iw, '(T2,A,1X,I6)') "vibration", l
674 DO j = 1, SIZE(particles)
675 IF (my_list(j) /= 0) THEN
676 k = (my_list(j) - 1)*3
677 WRITE (iw, '(T2,3(F12.6,3X))') eigen_vec(k + 1, i), eigen_vec(k + 2, i), eigen_vec(k + 3, i)
678 ELSE
679 WRITE (iw, '(T2,3(F12.6,3X))') 0.0_dp, 0.0_dp, 0.0_dp
680 END IF
681 END DO
682 END IF
683 END DO
684 IF (calc_intens) THEN
685 fint = massunit
686 ! intensity units are a.u./amu
687 WRITE (iw, '(T2,A)') "[INT]"
688 DO i = 1, SIZE(intensities)
689 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(3X,F18.6)') fint*intensities(i)**2
690 END DO
691 END IF
692 DEALLOCATE (my_list)
693 END IF
694 CALL cp_print_key_finished_output(iw, logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
695
696 CALL timestop(handle)
697
698 END SUBROUTINE write_vibrations_molden
699
700END MODULE molden_utils
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
qs_mo_methods::calculate_subspace_eigenvalues
Definition qs_mo_methods.F:68
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_utils
Contains methods used in the context of density fitting.
Definition admm_utils.F:15
admm_utils::admm_uncorrect_for_eigenvalues
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:127
admm_utils::admm_correct_for_eigenvalues
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:53
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1037
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:924
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1067
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::gto_cartesian
integer, parameter, public gto_cartesian
Definition input_constants.F:1328
input_constants::gto_spherical
integer, parameter, public gto_spherical
Definition input_constants.F:1328
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
molden_utils
Functions handling the MOLDEN format. Split from mode_selective.
Definition molden_utils.F:12
molden_utils::write_mos_molden
subroutine, public write_mos_molden(mos, qs_kind_set, particle_set, print_section, cell, unoccupied_orbs, unoccupied_evals, qs_env, calc_energies)
Write out the MOs in molden format for visualisation.
Definition molden_utils.F:85
molden_utils::write_vibrations_molden
subroutine, public write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, dump_only_positive, logger, list)
writes the output for vibrational analysis in MOLDEN format
Definition molden_utils.F:607
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
physcon::massunit
real(kind=dp), parameter, public massunit
Definition physcon.F:141
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:913
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition qs_mo_methods.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
cp_control_types::dft_control_type
Definition cp_control_types.F:751
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
orbital
Orbital angular momentum.
Definition grid_common.h:133
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47