(git:c835e6b)
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 cpassert(ASSOCIATED(mos(ispin)%mo_coeff_b))
156 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
157 mos(ispin)%mo_coeff) !fm->dbcsr
158 END DO
159 END IF
160
161 IF (iw > 0) THEN
162 WRITE (iw, '(T2,A)') "[Molden Format]"
163 IF (write_cell) THEN
164 IF (unit_choice == 2) THEN
165 WRITE (iw, '(T2,A)') "[Cell] Angs"
166 ELSE
167 WRITE (iw, '(T2,A)') "[Cell] AU"
168 END IF
169 WRITE (iw, '(T2,3(F12.6,3X))') &
170 cell%hmat(1, 1)*scale_factor, cell%hmat(2, 1)*scale_factor, cell%hmat(3, 1)*scale_factor
171 WRITE (iw, '(T2,3(F12.6,3X))') &
172 cell%hmat(1, 2)*scale_factor, cell%hmat(2, 2)*scale_factor, cell%hmat(3, 2)*scale_factor
173 WRITE (iw, '(T2,3(F12.6,3X))') &
174 cell%hmat(1, 3)*scale_factor, cell%hmat(2, 3)*scale_factor, cell%hmat(3, 3)*scale_factor
175 END IF
176 IF (unit_choice == 2) THEN
177 WRITE (iw, '(T2,A)') "[Atoms] Angs"
178 ELSE
179 WRITE (iw, '(T2,A)') "[Atoms] AU"
180 END IF
181 DO i = 1, SIZE(particle_set)
182 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
183 element_symbol=element_symbol)
184 CALL get_ptable_info(element_symbol, number=z)
185 IF (mark_ghost) THEN
186 CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_atom)
187 IF (ghost_atom) z = 0
188 END IF
189
190 WRITE (iw, '(T2,A2,I6,I6,3X,3(F12.6,3X))') &
191 element_symbol, i, z, particle_set(i)%r(:)*scale_factor
192 END DO
193 IF (write_pseudo) THEN
194 WRITE (iw, '(T2,A)') "[Pseudo]"
195 DO i = 1, SIZE(particle_set)
196 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
197 element_symbol=element_symbol)
198 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
199 WRITE (iw, '(T2,A2,I6,I6)') &
200 element_symbol, i, nint(zeff)
201 END DO
202 END IF
203
204 WRITE (iw, '(T2,A)') "[GTO]"
205
206 DO i = 1, SIZE(particle_set)
207 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
208 element_symbol=element_symbol)
209 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
210 IF (ASSOCIATED(orb_basis_set)) THEN
211 WRITE (iw, '(T2,I8,I8)') i, 0
212 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
213 nset=nset, &
214 npgf=npgf, &
215 nshell=nshell, &
216 l=l, &
217 zet=zet, &
218 gcc=gcc)
219
220 DO iset = 1, nset
221 DO ishell = 1, nshell(iset)
222 lshell = l(ishell, iset)
223 IF (lshell <= molden_lmax) THEN
224 WRITE (unit=iw, fmt='(T25,A2,4X,I4,4X,F4.2)') &
225 angmom(lshell + 1:lshell + 1), npgf(iset), 1.0_dp
226 ! MOLDEN expects the contraction coefficient of spherical NOT CARTESIAN NORMALISED
227 ! functions. So we undo the normalisation factors included in the gccs
228 ! Reverse engineered from basis_set_types, normalise_gcc_orb
229 prefac = 2_dp**lshell*(2/pi)**0.75_dp
230 expzet = 0.25_dp*(2*lshell + 3.0_dp)
231 WRITE (unit=iw, fmt=fmtstr2) &
232 (zet(ipgf, iset), gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet), &
233 ipgf=1, npgf(iset))
234 ELSE
235 IF (print_warn) THEN
236 CALL cp_warn(__location__, &
237 "MOLDEN format does not support Gaussian orbitals with l > 4.")
238 print_warn = .false.
239 END IF
240 END IF
241 END DO
242 END DO
243
244 WRITE (iw, '(A4)') " "
245
246 END IF
247
248 END DO
249
250 IF (gto_kind == gto_spherical) THEN
251 WRITE (iw, '(T2,A)') "[5D7F]"
252 WRITE (iw, '(T2,A)') "[9G]"
253 END IF
254
255 WRITE (iw, '(T2,A)') "[MO]"
256 END IF
257
258 !------------------------------------------------------------------------
259 ! convert from CP2K to MOLDEN format ordering
260 ! http://www.cmbi.ru.nl/molden/molden_format.html
261 !"The following order of D, F and G functions is expected:
262 !
263 ! 5D: D 0, D+1, D-1, D+2, D-2
264 ! 6D: xx, yy, zz, xy, xz, yz
265 !
266 ! 7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
267 ! 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
268 !
269 ! 9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
270 ! 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
271 ! xxyy xxzz yyzz xxyz yyxz zzxy
272 !"
273 ! CP2K has x in the outer (slower loop), so
274 ! xx, xy, xz, yy, yz,zz for l=2, for instance
275 !
276 ! iorb_cp2k = orbmap(iorb_molden, l), l = 0 .. 4
277 ! -----------------------------------------------------------------------
278 IF (iw > 0) THEN
279 IF (gto_kind == gto_cartesian) THEN
280 ! -----------------------------------------------------------------
281 ! Use cartesian (6D, 10F, 15G) representation.
282 ! This is only format VMD can process.
283 ! -----------------------------------------------------------------
284 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
285 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
286 1, 4, 6, 2, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
287 1, 7, 10, 4, 2, 3, 6, 9, 8, 5, 0, 0, 0, 0, 0, &
288 1, 11, 15, 2, 3, 7, 12, 10, 14, 4, 6, 13, 5, 8, 9], &
289 [molden_ncomax, molden_lmax + 1])
290 ELSE IF (gto_kind == gto_spherical) THEN
291 ! -----------------------------------------------------------------
292 ! Use spherical (5D, 7F, 9G) representation.
293 ! -----------------------------------------------------------------
294 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
295 3, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
296 3, 4, 2, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
297 4, 5, 3, 6, 2, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
298 5, 6, 4, 7, 3, 8, 2, 9, 1, 0, 0, 0, 0, 0, 0], &
299 [molden_ncomax, molden_lmax + 1])
300 END IF
301 END IF
302
303 DO ispin = 1, SIZE(mos)
304 do_calc_energies = .false.
305 IF (PRESENT(calc_energies)) do_calc_energies = calc_energies
306
307 IF (PRESENT(qs_env) .AND. do_calc_energies) THEN
308 CALL get_qs_env(qs_env, matrix_ks=ks, dft_control=dft_control)
309
310 matrix_ks => ks(ispin)%matrix
311
312 ! With ADMM, we have to modify the Kohn-Sham matrix
313 IF (dft_control%do_admm) THEN
314 CALL get_qs_env(qs_env, admm_env=admm_env)
315 CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks)
316 END IF
317
318 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues)
319
320 IF (ASSOCIATED(qs_env%mo_derivs)) THEN
321 mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
322 ELSE
323 mo_coeff_deriv => null()
324 END IF
325
326 ! Update the eigenvalues of the occupied orbitals
327 CALL calculate_subspace_eigenvalues(orbitals=mo_coeff, &
328 ks_matrix=matrix_ks, &
329 evals_arg=mo_eigenvalues, &
330 co_rotate_dbcsr=mo_coeff_deriv)
331
332 ! With ADMM, we have to undo the modification of the Kohn-Sham matrix
333 IF (dft_control%do_admm) THEN
334 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks)
335 END IF
336 END IF
337
338 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
339 nrow_global=nrow_global, &
340 ncol_global=ncol_global)
341 ALLOCATE (smatrix(nrow_global, ncol_global))
342 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
343
344 IF (iw > 0) THEN
345 IF (gto_kind == gto_cartesian) THEN
346 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
347
348 ALLOCATE (cmatrix(ncgf, ncgf))
349
350 cmatrix = 0.0_dp
351
352 ! Transform spherical MOs to Cartesian MOs
353
354 icgf = 1
355 isgf = 1
356 DO iatom = 1, SIZE(particle_set)
357 NULLIFY (orb_basis_set)
358 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
359 CALL get_qs_kind(qs_kind_set(ikind), &
360 basis_set=orb_basis_set)
361 IF (ASSOCIATED(orb_basis_set)) THEN
362 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
363 nset=nset, &
364 nshell=nshell, &
365 l=l)
366 DO iset = 1, nset
367 DO ishell = 1, nshell(iset)
368 lshell = l(ishell, iset)
369 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, nso(lshell), 1.0_dp, &
370 orbtramat(lshell)%c2s, nso(lshell), &
371 smatrix(isgf, 1), nsgf, 0.0_dp, &
372 cmatrix(icgf, 1), ncgf)
373 icgf = icgf + nco(lshell)
374 isgf = isgf + nso(lshell)
375 END DO
376 END DO
377 END IF
378 END DO ! iatom
379 END IF
380
381 DO icol = 1, mos(ispin)%nmo
382 ! index of the first basis function for the given atom, set, and shell
383 irow = 1
384
385 ! index of the first basis function in MOLDEN file.
386 ! Due to limitation of the MOLDEN format, basis functions with l > molden_lmax
387 ! cannot be exported, so we need to renumber atomic orbitals
388 irow_in = 1
389
390 WRITE (iw, '(A,ES20.10)') 'Ene=', mos(ispin)%eigenvalues(icol)
391 IF (ispin < 2) THEN
392 WRITE (iw, '(A)') 'Spin= Alpha'
393 ELSE
394 WRITE (iw, '(A)') 'Spin= Beta'
395 END IF
396 WRITE (iw, '(A,F12.7)') 'Occup=', mos(ispin)%occupation_numbers(icol)
397
398 DO iatom = 1, SIZE(particle_set)
399 NULLIFY (orb_basis_set)
400 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
401 element_symbol=element_symbol, kind_number=ikind)
402 CALL get_qs_kind(qs_kind_set(ikind), &
403 basis_set=orb_basis_set)
404 IF (ASSOCIATED(orb_basis_set)) THEN
405 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
406 nset=nset, &
407 nshell=nshell, &
408 l=l)
409
410 IF (gto_kind == gto_cartesian) THEN
411 ! ----------------------------------------------
412 ! Use cartesian (6D, 10F, 15G) representation.
413 ! ----------------------------------------------
414 icgf = 1
415 DO iset = 1, nset
416 DO ishell = 1, nshell(iset)
417 lshell = l(ishell, iset)
418
419 IF (lshell <= molden_lmax) THEN
420 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
421 cmatrix(irow:irow + nco(lshell) - 1, icol))
422 irow_in = irow_in + nco(lshell)
423 END IF
424
425 irow = irow + nco(lshell)
426 END DO ! ishell
427 END DO
428
429 ELSE IF (gto_kind == gto_spherical) THEN
430 ! ----------------------------------------------
431 ! Use spherical (5D, 7F, 9G) representation.
432 ! ----------------------------------------------
433 DO iset = 1, nset
434 DO ishell = 1, nshell(iset)
435 lshell = l(ishell, iset)
436
437 IF (lshell <= molden_lmax) THEN
438 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
439 smatrix(irow:irow + nso(lshell) - 1, icol))
440 irow_in = irow_in + nso(lshell)
441 END IF
442
443 irow = irow + nso(lshell)
444 END DO
445 END DO
446 END IF
447
448 END IF
449 END DO ! iatom
450 END DO
451 END IF
452
453 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
454 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
455 END DO
456
457 ! Write unoccupied (virtual) orbitals if provided; only used with OT
458 IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN
459 DO ispin = 1, SIZE(unoccupied_orbs)
460 CALL cp_fm_get_info(unoccupied_orbs(ispin), &
461 nrow_global=nrow_global, &
462 ncol_global=numos)
463 ALLOCATE (smatrix(nrow_global, numos))
464 CALL cp_fm_get_submatrix(unoccupied_orbs(ispin), smatrix)
465
466 IF (iw > 0) THEN
467 IF (gto_kind == gto_cartesian) THEN
468 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
469 ALLOCATE (cmatrix(ncgf, numos))
470 cmatrix = 0.0_dp
471
472 icgf = 1
473 isgf = 1
474 DO iatom = 1, SIZE(particle_set)
475 NULLIFY (orb_basis_set)
476 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
477 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
478 IF (ASSOCIATED(orb_basis_set)) THEN
479 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
480 nset=nset, nshell=nshell, l=l)
481 DO iset = 1, nset
482 DO ishell = 1, nshell(iset)
483 lshell = l(ishell, iset)
484 CALL dgemm("T", "N", nco(lshell), numos, nso(lshell), 1.0_dp, &
485 orbtramat(lshell)%c2s, nso(lshell), &
486 smatrix(isgf, 1), nsgf, 0.0_dp, &
487 cmatrix(icgf, 1), ncgf)
488 icgf = icgf + nco(lshell)
489 isgf = isgf + nso(lshell)
490 END DO
491 END DO
492 END IF
493 END DO
494 END IF
495
496 DO icol = 1, numos
497 irow = 1
498 irow_in = 1
499
500 WRITE (iw, '(A,ES20.10)') 'Ene=', unoccupied_evals(ispin)%array(icol)
501 IF (ispin < 2) THEN
502 WRITE (iw, '(A)') 'Spin= Alpha'
503 ELSE
504 WRITE (iw, '(A)') 'Spin= Beta'
505 END IF
506 WRITE (iw, '(A,F12.7)') 'Occup=', 0.0_dp
507
508 DO iatom = 1, SIZE(particle_set)
509 NULLIFY (orb_basis_set)
510 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
511 element_symbol=element_symbol, kind_number=ikind)
512 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
513 IF (ASSOCIATED(orb_basis_set)) THEN
514 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
515 nset=nset, nshell=nshell, l=l)
516
517 IF (gto_kind == gto_cartesian) THEN
518 icgf = 1
519 DO iset = 1, nset
520 DO ishell = 1, nshell(iset)
521 lshell = l(ishell, iset)
522 IF (lshell <= molden_lmax) THEN
523 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
524 cmatrix(irow:irow + nco(lshell) - 1, icol))
525 irow_in = irow_in + nco(lshell)
526 END IF
527 irow = irow + nco(lshell)
528 END DO
529 END DO
530 ELSE IF (gto_kind == gto_spherical) THEN
531 DO iset = 1, nset
532 DO ishell = 1, nshell(iset)
533 lshell = l(ishell, iset)
534 IF (lshell <= molden_lmax) THEN
535 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
536 smatrix(irow:irow + nso(lshell) - 1, icol))
537 irow_in = irow_in + nso(lshell)
538 END IF
539 irow = irow + nso(lshell)
540 END DO
541 END DO
542 END IF
543
544 END IF
545 END DO ! iatom
546 END DO ! icol
547 END IF
548
549 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
550 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
551 END DO ! ispin
552 END IF
553
554 CALL cp_print_key_finished_output(iw, logger, print_section, "")
555
556 END IF
557
558 CALL timestop(handle)
559
560 END SUBROUTINE write_mos_molden
561
562! **************************************************************************************************
563!> \brief Output MO coefficients formatted correctly for MOLDEN, omitting those <= 1E(-digits)
564!> \param iw output file unit
565!> \param fmtstr1 format string
566!> \param ndigits number of significant digits in MO coefficients
567!> \param irow_in index of the first atomic orbital: mo_coeff(orbmap(1))
568!> \param orbmap array to map Gaussian functions from MOLDEN to CP2K ordering
569!> \param mo_coeff MO coefficients
570! **************************************************************************************************
571 SUBROUTINE print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap, mo_coeff)
572 INTEGER, INTENT(in) :: iw
573 CHARACTER(LEN=*), INTENT(in) :: fmtstr1
574 INTEGER, INTENT(in) :: ndigits, irow_in
575 INTEGER, DIMENSION(molden_ncomax), INTENT(in) :: orbmap
576 REAL(kind=dp), DIMENSION(:), INTENT(in) :: mo_coeff
577
578 INTEGER :: orbital
579
580 DO orbital = 1, molden_ncomax
581 IF (orbmap(orbital) /= 0) THEN
582 IF (abs(mo_coeff(orbmap(orbital))) >= 10.0_dp**(-ndigits)) THEN
583 WRITE (iw, fmtstr1) irow_in + orbital - 1, mo_coeff(orbmap(orbital))
584 END IF
585 END IF
586 END DO
587
588 END SUBROUTINE print_coeffs
589
590! **************************************************************************************************
591!> \brief writes the output for vibrational analysis in MOLDEN format
592!> \param input ...
593!> \param particles ...
594!> \param freq ...
595!> \param eigen_vec ...
596!> \param intensities ...
597!> \param calc_intens ...
598!> \param dump_only_positive ...
599!> \param logger ...
600!> \param list array of mobile atom indices
601!> \author Florian Schiffmann 11.2007
602! **************************************************************************************************
603 SUBROUTINE write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, &
604 dump_only_positive, logger, list)
605
606 TYPE(section_vals_type), POINTER :: input
607 TYPE(particle_type), DIMENSION(:), POINTER :: particles
608 REAL(kind=dp), DIMENSION(:) :: freq
609 REAL(kind=dp), DIMENSION(:, :) :: eigen_vec
610 REAL(kind=dp), DIMENSION(:), POINTER :: intensities
611 LOGICAL, INTENT(in) :: calc_intens, dump_only_positive
612 TYPE(cp_logger_type), POINTER :: logger
613 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: list
614
615 CHARACTER(len=*), PARAMETER :: routinen = 'write_vibrations_molden'
616
617 CHARACTER(LEN=2) :: element_symbol
618 INTEGER :: handle, i, iw, j, k, l, z
619 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
620 REAL(kind=dp) :: fint
621
622 CALL timeset(routinen, handle)
623
624 iw = cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
625 extension=".mol", file_status='REPLACE')
626
627 IF (iw > 0) THEN
628 cpassert(mod(SIZE(eigen_vec, 1), 3) == 0)
629 cpassert(SIZE(freq, 1) == SIZE(eigen_vec, 2))
630 ALLOCATE (my_list(SIZE(particles)))
631 ! Either we have a list of the subset of mobile atoms,
632 ! Or the eigenvectors must span the full space (all atoms)
633 IF (PRESENT(list)) THEN
634 my_list(:) = 0
635 DO i = 1, SIZE(list)
636 my_list(list(i)) = i
637 END DO
638 ELSE
639 cpassert(SIZE(particles) == SIZE(eigen_vec, 1)/3)
640 DO i = 1, SIZE(my_list)
641 my_list(i) = i
642 END DO
643 END IF
644 WRITE (iw, '(T2,A)') "[Molden Format]"
645 WRITE (iw, '(T2,A)') "[Atoms] AU"
646 DO i = 1, SIZE(particles)
647 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
648 element_symbol=element_symbol)
649 CALL get_ptable_info(element_symbol, number=z)
650
651 WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
652 element_symbol, i, z, particles(i)%r(:)
653
654 END DO
655 WRITE (iw, '(T2,A)') "[FREQ]"
656 DO i = 1, SIZE(freq, 1)
657 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(T5,F12.6)') freq(i)
658 END DO
659 WRITE (iw, '(T2,A)') "[FR-COORD]"
660 DO i = 1, SIZE(particles)
661 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
662 element_symbol=element_symbol)
663 WRITE (iw, '(T2,A2,3X,3(F12.6,3X))') &
664 element_symbol, particles(i)%r(:)
665 END DO
666 WRITE (iw, '(T2,A)') "[FR-NORM-COORD]"
667 l = 0
668 DO i = 1, SIZE(eigen_vec, 2)
669 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) THEN
670 l = l + 1
671 WRITE (iw, '(T2,A,1X,I6)') "vibration", l
672 DO j = 1, SIZE(particles)
673 IF (my_list(j) /= 0) THEN
674 k = (my_list(j) - 1)*3
675 WRITE (iw, '(T2,3(F12.6,3X))') eigen_vec(k + 1, i), eigen_vec(k + 2, i), eigen_vec(k + 3, i)
676 ELSE
677 WRITE (iw, '(T2,3(F12.6,3X))') 0.0_dp, 0.0_dp, 0.0_dp
678 END IF
679 END DO
680 END IF
681 END DO
682 IF (calc_intens) THEN
683 fint = massunit
684 ! intensity units are a.u./amu
685 WRITE (iw, '(T2,A)') "[INT]"
686 DO i = 1, SIZE(intensities)
687 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(3X,F18.6)') fint*intensities(i)**2
688 END DO
689 END IF
690 DEALLOCATE (my_list)
691 END IF
692 CALL cp_print_key_finished_output(iw, logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
693
694 CALL timestop(handle)
695
696 END SUBROUTINE write_vibrations_molden
697
698END 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:605
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