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