(git:0341268)
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! **************************************************************************************************
13 USE admm_types, ONLY: admm_type
19 USE cell_types, ONLY: cell_type
22 USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
25 USE cp_fm_types, ONLY: cp_fm_get_info,&
30 USE cp_output_handling, ONLY: cp_p_file,&
38 USE kinds, ONLY: dp
39 USE mathconstants, ONLY: pi
40 USE orbital_pointers, ONLY: nco,&
41 nso
45 USE physcon, ONLY: angstrom,&
49 USE qs_kind_types, ONLY: get_qs_kind,&
53 USE qs_mo_types, ONLY: get_mo_set,&
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
67
68CONTAINS
69
70! **************************************************************************************************
71!> \brief Write the CP2K [Cell] extension to a MOLDEN file
72!> \param iw output unit
73!> \param cell simulation cell
74!> \param unit_choice 1 for atomic units, 2 for Angstrom
75! **************************************************************************************************
76 SUBROUTINE write_cell_molden(iw, cell, unit_choice)
77 INTEGER, INTENT(IN) :: iw
78 TYPE(cell_type), INTENT(IN) :: cell
79 INTEGER, INTENT(IN) :: unit_choice
80
81 REAL(KIND=dp) :: scale_factor
82
83 IF (unit_choice == 2) THEN
84 scale_factor = angstrom
85 WRITE (iw, '(T2,A)') "[Cell] Angs"
86 ELSE
87 scale_factor = 1.0_dp
88 WRITE (iw, '(T2,A)') "[Cell] AU"
89 END IF
90 WRITE (iw, '(T2,3(F12.6,3X))') &
91 cell%hmat(1, 1)*scale_factor, cell%hmat(2, 1)*scale_factor, cell%hmat(3, 1)*scale_factor
92 WRITE (iw, '(T2,3(F12.6,3X))') &
93 cell%hmat(1, 2)*scale_factor, cell%hmat(2, 2)*scale_factor, cell%hmat(3, 2)*scale_factor
94 WRITE (iw, '(T2,3(F12.6,3X))') &
95 cell%hmat(1, 3)*scale_factor, cell%hmat(2, 3)*scale_factor, cell%hmat(3, 3)*scale_factor
96 END SUBROUTINE write_cell_molden
97
98! **************************************************************************************************
99!> \brief Write out the MOs in molden format for visualisation
100!> \param mos the set of MOs (both spins, if UKS)
101!> \param qs_kind_set for basis set info
102!> \param particle_set particles data structure, for positions and kinds
103!> \param print_section input section containing relevant print key
104!> \param cell ...
105!> \param unoccupied_orbs optional: unoccupied orbital coefficients from make_lumo_gpw
106!> \param unoccupied_evals optional: unoccupied orbital eigenvalues
107!> \param qs_env ...
108!> \param calc_energies ...
109!> \author MattW, IainB
110! **************************************************************************************************
111 SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section, cell, &
112 unoccupied_orbs, unoccupied_evals, qs_env, calc_energies)
113 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
114 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
115 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
116 TYPE(section_vals_type), POINTER :: print_section
117 TYPE(cell_type), OPTIONAL, POINTER :: cell
118 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
119 OPTIONAL :: unoccupied_orbs
120 TYPE(cp_1d_r_p_type), DIMENSION(:), INTENT(IN), &
121 OPTIONAL :: unoccupied_evals
122 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
123 LOGICAL, INTENT(IN), OPTIONAL :: calc_energies
124
125 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_mos_molden'
126 CHARACTER(LEN=molden_lmax+1), PARAMETER :: angmom = "spdfg"
127
128 CHARACTER(LEN=15) :: fmtstr1, fmtstr2
129 CHARACTER(LEN=2) :: element_symbol
130 INTEGER :: gto_kind, handle, i, iatom, icgf, icol, ikind, ipgf, irow, irow_in, iset, isgf, &
131 ishell, ispin, iw, lshell, ncgf, ncol_global, ndigits, nrow_global, nset, nsgf, numos, &
132 unit_choice, z
133 INTEGER, DIMENSION(:), POINTER :: npgf, nshell
134 INTEGER, DIMENSION(:, :), POINTER :: l
135 INTEGER, DIMENSION(molden_ncomax, 0:molden_lmax) :: orbmap
136 LOGICAL :: do_calc_energies, ghost_atom, &
137 mark_ghost, print_warn, write_cell, &
138 write_pseudo
139 REAL(kind=dp) :: expzet, prefac, scale_factor, zeff
140 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
141 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
142 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
143 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
144 TYPE(admm_type), POINTER :: admm_env
145 TYPE(cp_fm_type), POINTER :: mo_coeff
146 TYPE(cp_logger_type), POINTER :: logger
147 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks
148 TYPE(dbcsr_type), POINTER :: matrix_ks, mo_coeff_deriv
149 TYPE(dft_control_type), POINTER :: dft_control
150 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
151
152 CALL timeset(routinen, handle)
153
154 logger => cp_get_default_logger()
155 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, ""), cp_p_file)) THEN
156
157 iw = cp_print_key_unit_nr(logger, print_section, "", &
158 extension=".molden", file_status='REPLACE')
159
160 print_warn = .true.
161
162 CALL section_vals_val_get(print_section, "UNIT", i_val=unit_choice)
163 IF (unit_choice == 2) THEN
164 scale_factor = angstrom
165 ELSE
166 scale_factor = 1.0_dp
167 END IF
168
169 CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
170 ndigits = min(max(3, ndigits), 30)
171 WRITE (unit=fmtstr1, fmt='("(I6,1X,ES",I0,".",I0,")")') ndigits + 7, ndigits
172 WRITE (unit=fmtstr2, fmt='("((T51,2F",I0,".",I0,"))")') ndigits + 10, ndigits
173
174 CALL section_vals_val_get(print_section, "GTO_KIND", i_val=gto_kind)
175 CALL section_vals_val_get(print_section, "WRITE_CELL", l_val=write_cell)
176 CALL section_vals_val_get(print_section, "WRITE_PSEUDO", l_val=write_pseudo)
177 CALL section_vals_val_get(print_section, "MARK_GHOST", l_val=mark_ghost)
178
179 IF (mos(1)%use_mo_coeff_b) THEN
180 ! we are using the dbcsr mo_coeff
181 ! we copy it to the fm anyway
182 DO ispin = 1, SIZE(mos)
183 cpassert(ASSOCIATED(mos(ispin)%mo_coeff_b))
184 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
185 mos(ispin)%mo_coeff) !fm->dbcsr
186 END DO
187 END IF
188
189 IF (iw > 0) THEN
190 WRITE (iw, '(T2,A)') "[Molden Format]"
191 IF (write_cell) THEN
192 cpassert(PRESENT(cell))
193 cpassert(ASSOCIATED(cell))
194 CALL write_cell_molden(iw, cell, unit_choice)
195 END IF
196 IF (unit_choice == 2) THEN
197 WRITE (iw, '(T2,A)') "[Atoms] Angs"
198 ELSE
199 WRITE (iw, '(T2,A)') "[Atoms] AU"
200 END IF
201 DO i = 1, SIZE(particle_set)
202 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
203 element_symbol=element_symbol)
204 CALL get_ptable_info(element_symbol, number=z)
205 IF (mark_ghost) THEN
206 CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_atom)
207 IF (ghost_atom) z = 0
208 END IF
209
210 WRITE (iw, '(T2,A2,I6,I6,3X,3(F12.6,3X))') &
211 element_symbol, i, z, particle_set(i)%r(:)*scale_factor
212 END DO
213 IF (write_pseudo) THEN
214 WRITE (iw, '(T2,A)') "[Pseudo]"
215 DO i = 1, SIZE(particle_set)
216 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
217 element_symbol=element_symbol)
218 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
219 WRITE (iw, '(T2,A2,I6,I6)') &
220 element_symbol, i, nint(zeff)
221 END DO
222 END IF
223
224 WRITE (iw, '(T2,A)') "[GTO]"
225
226 DO i = 1, SIZE(particle_set)
227 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
228 element_symbol=element_symbol)
229 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
230 IF (ASSOCIATED(orb_basis_set)) THEN
231 WRITE (iw, '(T2,I8,I8)') i, 0
232 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
233 nset=nset, &
234 npgf=npgf, &
235 nshell=nshell, &
236 l=l, &
237 zet=zet, &
238 gcc=gcc)
239
240 DO iset = 1, nset
241 DO ishell = 1, nshell(iset)
242 lshell = l(ishell, iset)
243 IF (lshell <= molden_lmax) THEN
244 WRITE (unit=iw, fmt='(T25,A2,4X,I4,4X,F4.2)') &
245 angmom(lshell + 1:lshell + 1), npgf(iset), 1.0_dp
246 ! MOLDEN expects the contraction coefficient of spherical NOT CARTESIAN NORMALISED
247 ! functions. So we undo the normalisation factors included in the gccs
248 ! Reverse engineered from basis_set_types, normalise_gcc_orb
249 prefac = 2_dp**lshell*(2/pi)**0.75_dp
250 expzet = 0.25_dp*(2*lshell + 3.0_dp)
251 WRITE (unit=iw, fmt=fmtstr2) &
252 (zet(ipgf, iset), gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet), &
253 ipgf=1, npgf(iset))
254 ELSE
255 IF (print_warn) THEN
256 CALL cp_warn(__location__, &
257 "MOLDEN format does not support Gaussian orbitals with l > 4.")
258 print_warn = .false.
259 END IF
260 END IF
261 END DO
262 END DO
263
264 WRITE (iw, '(A4)') " "
265
266 END IF
267
268 END DO
269
270 IF (gto_kind == gto_spherical) THEN
271 WRITE (iw, '(T2,A)') "[5D7F]"
272 WRITE (iw, '(T2,A)') "[9G]"
273 END IF
274
275 WRITE (iw, '(T2,A)') "[MO]"
276 END IF
277
278 !------------------------------------------------------------------------
279 ! convert from CP2K to MOLDEN format ordering
280 ! http://www.cmbi.ru.nl/molden/molden_format.html
281 !"The following order of D, F and G functions is expected:
282 !
283 ! 5D: D 0, D+1, D-1, D+2, D-2
284 ! 6D: xx, yy, zz, xy, xz, yz
285 !
286 ! 7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
287 ! 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
288 !
289 ! 9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
290 ! 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
291 ! xxyy xxzz yyzz xxyz yyxz zzxy
292 !"
293 ! CP2K has x in the outer (slower loop), so
294 ! xx, xy, xz, yy, yz,zz for l=2, for instance
295 !
296 ! iorb_cp2k = orbmap(iorb_molden, l), l = 0 .. 4
297 ! -----------------------------------------------------------------------
298 IF (iw > 0) THEN
299 IF (gto_kind == gto_cartesian) THEN
300 ! -----------------------------------------------------------------
301 ! Use cartesian (6D, 10F, 15G) representation.
302 ! This is only format VMD can process.
303 ! -----------------------------------------------------------------
304 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
305 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
306 1, 4, 6, 2, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
307 1, 7, 10, 4, 2, 3, 6, 9, 8, 5, 0, 0, 0, 0, 0, &
308 1, 11, 15, 2, 3, 7, 12, 10, 14, 4, 6, 13, 5, 8, 9], &
309 [molden_ncomax, molden_lmax + 1])
310 ELSE IF (gto_kind == gto_spherical) THEN
311 ! -----------------------------------------------------------------
312 ! Use spherical (5D, 7F, 9G) representation.
313 ! -----------------------------------------------------------------
314 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
315 3, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
316 3, 4, 2, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
317 4, 5, 3, 6, 2, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
318 5, 6, 4, 7, 3, 8, 2, 9, 1, 0, 0, 0, 0, 0, 0], &
319 [molden_ncomax, molden_lmax + 1])
320 END IF
321 END IF
322
323 DO ispin = 1, SIZE(mos)
324 do_calc_energies = .false.
325 IF (PRESENT(calc_energies)) do_calc_energies = calc_energies
326
327 IF (PRESENT(qs_env) .AND. do_calc_energies) THEN
328 CALL get_qs_env(qs_env, matrix_ks=ks, dft_control=dft_control)
329
330 matrix_ks => ks(ispin)%matrix
331
332 ! With ADMM, we have to modify the Kohn-Sham matrix
333 IF (dft_control%do_admm) THEN
334 CALL get_qs_env(qs_env, admm_env=admm_env)
335 CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks)
336 END IF
337
338 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues)
339
340 IF (ASSOCIATED(qs_env%mo_derivs)) THEN
341 mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
342 ELSE
343 mo_coeff_deriv => null()
344 END IF
345
346 ! Update the eigenvalues of the occupied orbitals
347 CALL calculate_subspace_eigenvalues(orbitals=mo_coeff, &
348 ks_matrix=matrix_ks, &
349 evals_arg=mo_eigenvalues, &
350 co_rotate_dbcsr=mo_coeff_deriv)
351
352 ! With ADMM, we have to undo the modification of the Kohn-Sham matrix
353 IF (dft_control%do_admm) THEN
354 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks)
355 END IF
356 END IF
357
358 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
359 nrow_global=nrow_global, &
360 ncol_global=ncol_global)
361 ALLOCATE (smatrix(nrow_global, ncol_global))
362 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
363
364 IF (iw > 0) THEN
365 IF (gto_kind == gto_cartesian) THEN
366 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
367
368 ALLOCATE (cmatrix(ncgf, ncgf))
369
370 cmatrix = 0.0_dp
371
372 ! Transform spherical MOs to Cartesian MOs
373
374 icgf = 1
375 isgf = 1
376 DO iatom = 1, SIZE(particle_set)
377 NULLIFY (orb_basis_set)
378 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
379 CALL get_qs_kind(qs_kind_set(ikind), &
380 basis_set=orb_basis_set)
381 IF (ASSOCIATED(orb_basis_set)) THEN
382 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
383 nset=nset, &
384 nshell=nshell, &
385 l=l)
386 DO iset = 1, nset
387 DO ishell = 1, nshell(iset)
388 lshell = l(ishell, iset)
389 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, nso(lshell), 1.0_dp, &
390 orbtramat(lshell)%c2s, nso(lshell), &
391 smatrix(isgf, 1), nsgf, 0.0_dp, &
392 cmatrix(icgf, 1), ncgf)
393 icgf = icgf + nco(lshell)
394 isgf = isgf + nso(lshell)
395 END DO
396 END DO
397 END IF
398 END DO ! iatom
399 END IF
400
401 DO icol = 1, mos(ispin)%nmo
402 ! index of the first basis function for the given atom, set, and shell
403 irow = 1
404
405 ! index of the first basis function in MOLDEN file.
406 ! Due to limitation of the MOLDEN format, basis functions with l > molden_lmax
407 ! cannot be exported, so we need to renumber atomic orbitals
408 irow_in = 1
409
410 WRITE (iw, '(A,ES20.10)') 'Ene=', mos(ispin)%eigenvalues(icol)
411 IF (ispin < 2) THEN
412 WRITE (iw, '(A)') 'Spin= Alpha'
413 ELSE
414 WRITE (iw, '(A)') 'Spin= Beta'
415 END IF
416 WRITE (iw, '(A,F12.7)') 'Occup=', mos(ispin)%occupation_numbers(icol)
417
418 DO iatom = 1, SIZE(particle_set)
419 NULLIFY (orb_basis_set)
420 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
421 element_symbol=element_symbol, kind_number=ikind)
422 CALL get_qs_kind(qs_kind_set(ikind), &
423 basis_set=orb_basis_set)
424 IF (ASSOCIATED(orb_basis_set)) THEN
425 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
426 nset=nset, &
427 nshell=nshell, &
428 l=l)
429
430 IF (gto_kind == gto_cartesian) THEN
431 ! ----------------------------------------------
432 ! Use cartesian (6D, 10F, 15G) representation.
433 ! ----------------------------------------------
434 icgf = 1
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 cmatrix(irow:irow + nco(lshell) - 1, icol))
442 irow_in = irow_in + nco(lshell)
443 END IF
444
445 irow = irow + nco(lshell)
446 END DO ! ishell
447 END DO
448
449 ELSE IF (gto_kind == gto_spherical) THEN
450 ! ----------------------------------------------
451 ! Use spherical (5D, 7F, 9G) representation.
452 ! ----------------------------------------------
453 DO iset = 1, nset
454 DO ishell = 1, nshell(iset)
455 lshell = l(ishell, iset)
456
457 IF (lshell <= molden_lmax) THEN
458 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
459 smatrix(irow:irow + nso(lshell) - 1, icol))
460 irow_in = irow_in + nso(lshell)
461 END IF
462
463 irow = irow + nso(lshell)
464 END DO
465 END DO
466 END IF
467
468 END IF
469 END DO ! iatom
470 END DO
471 END IF
472
473 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
474 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
475 END DO
476
477 ! Write unoccupied (virtual) orbitals if provided; only used with OT
478 IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN
479 DO ispin = 1, SIZE(unoccupied_orbs)
480 CALL cp_fm_get_info(unoccupied_orbs(ispin), &
481 nrow_global=nrow_global, &
482 ncol_global=numos)
483 ALLOCATE (smatrix(nrow_global, numos))
484 CALL cp_fm_get_submatrix(unoccupied_orbs(ispin), smatrix)
485
486 IF (iw > 0) THEN
487 IF (gto_kind == gto_cartesian) THEN
488 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
489 ALLOCATE (cmatrix(ncgf, numos))
490 cmatrix = 0.0_dp
491
492 icgf = 1
493 isgf = 1
494 DO iatom = 1, SIZE(particle_set)
495 NULLIFY (orb_basis_set)
496 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
497 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
498 IF (ASSOCIATED(orb_basis_set)) THEN
499 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
500 nset=nset, nshell=nshell, l=l)
501 DO iset = 1, nset
502 DO ishell = 1, nshell(iset)
503 lshell = l(ishell, iset)
504 CALL dgemm("T", "N", nco(lshell), numos, nso(lshell), 1.0_dp, &
505 orbtramat(lshell)%c2s, nso(lshell), &
506 smatrix(isgf, 1), nsgf, 0.0_dp, &
507 cmatrix(icgf, 1), ncgf)
508 icgf = icgf + nco(lshell)
509 isgf = isgf + nso(lshell)
510 END DO
511 END DO
512 END IF
513 END DO
514 END IF
515
516 DO icol = 1, numos
517 irow = 1
518 irow_in = 1
519
520 WRITE (iw, '(A,ES20.10)') 'Ene=', unoccupied_evals(ispin)%array(icol)
521 IF (ispin < 2) THEN
522 WRITE (iw, '(A)') 'Spin= Alpha'
523 ELSE
524 WRITE (iw, '(A)') 'Spin= Beta'
525 END IF
526 WRITE (iw, '(A,F12.7)') 'Occup=', 0.0_dp
527
528 DO iatom = 1, SIZE(particle_set)
529 NULLIFY (orb_basis_set)
530 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
531 element_symbol=element_symbol, kind_number=ikind)
532 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
533 IF (ASSOCIATED(orb_basis_set)) THEN
534 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
535 nset=nset, nshell=nshell, l=l)
536
537 IF (gto_kind == gto_cartesian) THEN
538 icgf = 1
539 DO iset = 1, nset
540 DO ishell = 1, nshell(iset)
541 lshell = l(ishell, iset)
542 IF (lshell <= molden_lmax) THEN
543 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
544 cmatrix(irow:irow + nco(lshell) - 1, icol))
545 irow_in = irow_in + nco(lshell)
546 END IF
547 irow = irow + nco(lshell)
548 END DO
549 END DO
550 ELSE IF (gto_kind == gto_spherical) THEN
551 DO iset = 1, nset
552 DO ishell = 1, nshell(iset)
553 lshell = l(ishell, iset)
554 IF (lshell <= molden_lmax) THEN
555 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
556 smatrix(irow:irow + nso(lshell) - 1, icol))
557 irow_in = irow_in + nso(lshell)
558 END IF
559 irow = irow + nso(lshell)
560 END DO
561 END DO
562 END IF
563
564 END IF
565 END DO ! iatom
566 END DO ! icol
567 END IF
568
569 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
570 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
571 END DO ! ispin
572 END IF
573
574 CALL cp_print_key_finished_output(iw, logger, print_section, "")
575
576 END IF
577
578 CALL timestop(handle)
579
580 END SUBROUTINE write_mos_molden
581
582! **************************************************************************************************
583!> \brief Output MO coefficients formatted correctly for MOLDEN, omitting those <= 1E(-digits)
584!> \param iw output file unit
585!> \param fmtstr1 format string
586!> \param ndigits number of significant digits in MO coefficients
587!> \param irow_in index of the first atomic orbital: mo_coeff(orbmap(1))
588!> \param orbmap array to map Gaussian functions from MOLDEN to CP2K ordering
589!> \param mo_coeff MO coefficients
590! **************************************************************************************************
591 SUBROUTINE print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap, mo_coeff)
592 INTEGER, INTENT(in) :: iw
593 CHARACTER(LEN=*), INTENT(in) :: fmtstr1
594 INTEGER, INTENT(in) :: ndigits, irow_in
595 INTEGER, DIMENSION(molden_ncomax), INTENT(in) :: orbmap
596 REAL(kind=dp), DIMENSION(:), INTENT(in) :: mo_coeff
597
598 INTEGER :: orbital
599
600 DO orbital = 1, molden_ncomax
601 IF (orbmap(orbital) /= 0) THEN
602 IF (abs(mo_coeff(orbmap(orbital))) >= 10.0_dp**(-ndigits)) THEN
603 WRITE (iw, fmtstr1) irow_in + orbital - 1, mo_coeff(orbmap(orbital))
604 END IF
605 END IF
606 END DO
607
608 END SUBROUTINE print_coeffs
609
610! **************************************************************************************************
611!> \brief writes the output for vibrational analysis in MOLDEN format
612!> \param input ...
613!> \param particles ...
614!> \param freq ...
615!> \param eigen_vec ...
616!> \param intensities ...
617!> \param calc_intens ...
618!> \param dump_only_positive ...
619!> \param logger ...
620!> \param list array of mobile atom indices
621!> \param cell optional simulation cell for the CP2K [Cell] extension
622!> \author Florian Schiffmann 11.2007
623! **************************************************************************************************
624 SUBROUTINE write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, &
625 dump_only_positive, logger, list, cell)
626
627 TYPE(section_vals_type), POINTER :: input
628 TYPE(particle_type), DIMENSION(:), POINTER :: particles
629 REAL(kind=dp), DIMENSION(:) :: freq
630 REAL(kind=dp), DIMENSION(:, :) :: eigen_vec
631 REAL(kind=dp), DIMENSION(:), POINTER :: intensities
632 LOGICAL, INTENT(in) :: calc_intens, dump_only_positive
633 TYPE(cp_logger_type), POINTER :: logger
634 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: list
635 TYPE(cell_type), OPTIONAL, POINTER :: cell
636
637 CHARACTER(len=*), PARAMETER :: routinen = 'write_vibrations_molden'
638
639 CHARACTER(LEN=2) :: element_symbol
640 INTEGER :: handle, i, iw, j, k, l, z
641 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
642 LOGICAL :: write_cell
643 REAL(kind=dp) :: fint
644
645 CALL timeset(routinen, handle)
646
647 iw = cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
648 extension=".mol", file_status='REPLACE')
649
650 IF (iw > 0) THEN
651 cpassert(mod(SIZE(eigen_vec, 1), 3) == 0)
652 cpassert(SIZE(freq, 1) == SIZE(eigen_vec, 2))
653 ALLOCATE (my_list(SIZE(particles)))
654 ! Either we have a list of the subset of mobile atoms,
655 ! Or the eigenvectors must span the full space (all atoms)
656 IF (PRESENT(list)) THEN
657 my_list(:) = 0
658 DO i = 1, SIZE(list)
659 my_list(list(i)) = i
660 END DO
661 ELSE
662 cpassert(SIZE(particles) == SIZE(eigen_vec, 1)/3)
663 DO i = 1, SIZE(my_list)
664 my_list(i) = i
665 END DO
666 END IF
667 WRITE (iw, '(T2,A)') "[Molden Format]"
668 CALL section_vals_val_get(input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB%WRITE_CELL", &
669 l_val=write_cell)
670 IF (write_cell) THEN
671 cpassert(PRESENT(cell))
672 cpassert(ASSOCIATED(cell))
673 CALL write_cell_molden(iw, cell, 1)
674 END IF
675 WRITE (iw, '(T2,A)') "[Atoms] AU"
676 DO i = 1, SIZE(particles)
677 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
678 element_symbol=element_symbol)
679 CALL get_ptable_info(element_symbol, number=z)
680
681 WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
682 element_symbol, i, z, particles(i)%r(:)
683
684 END DO
685 WRITE (iw, '(T2,A)') "[FREQ]"
686 DO i = 1, SIZE(freq, 1)
687 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(T5,F12.6)') freq(i)
688 END DO
689 WRITE (iw, '(T2,A)') "[FR-COORD]"
690 DO i = 1, SIZE(particles)
691 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
692 element_symbol=element_symbol)
693 WRITE (iw, '(T2,A2,3X,3(F12.6,3X))') &
694 element_symbol, particles(i)%r(:)
695 END DO
696 WRITE (iw, '(T2,A)') "[FR-NORM-COORD]"
697 l = 0
698 DO i = 1, SIZE(eigen_vec, 2)
699 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) THEN
700 l = l + 1
701 WRITE (iw, '(T2,A,1X,I6)') "vibration", l
702 DO j = 1, SIZE(particles)
703 IF (my_list(j) /= 0) THEN
704 k = (my_list(j) - 1)*3
705 WRITE (iw, '(T2,3(F12.6,3X))') eigen_vec(k + 1, i), eigen_vec(k + 2, i), eigen_vec(k + 3, i)
706 ELSE
707 WRITE (iw, '(T2,3(F12.6,3X))') 0.0_dp, 0.0_dp, 0.0_dp
708 END IF
709 END DO
710 END IF
711 END DO
712 IF (calc_intens) THEN
713 fint = massunit
714 ! intensity units are a.u./amu
715 WRITE (iw, '(T2,A)') "[INT]"
716 DO i = 1, SIZE(intensities)
717 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(3X,F18.6)') fint*intensities(i)**2
718 END DO
719 END IF
720 DEALLOCATE (my_list)
721 END IF
722 CALL cp_print_key_finished_output(iw, logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
723
724 CALL timestop(handle)
725
726 END SUBROUTINE write_vibrations_molden
727
728END MODULE molden_utils
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.
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
Contains methods used in the context of density fitting.
Definition admm_utils.F:15
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:127
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:53
Define the atomic kind types and their sub types.
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.
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)
...
Handles all functions related to the CELL.
Definition cell_types.F:15
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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
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,...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
integer, parameter, public cp_p_file
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...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public gto_cartesian
integer, parameter, public gto_spherical
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Functions handling the MOLDEN format. Split from mode_selective.
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.
subroutine, public write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, dump_only_positive, logger, list, cell)
writes the output for vibrational analysis in MOLDEN format
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nco
integer, dimension(:), allocatable, public nso
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
type(orbtramat_type), dimension(:), pointer, public orbtramat
Define the data structure for the particle information.
Periodic Table related data definitions.
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 of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
real(kind=dp), parameter, public massunit
Definition physcon.F:141
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, matrix_vhxc, 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.
Define the quickstep kind type and their sub types.
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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, 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.
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, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
collects routines that perform operations directly related to MOs
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
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.
stores some data used in wavefunction fitting
Definition admm_types.F:120
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a pointer to a 1d array
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Orbital angular momentum.
Provides all information about a quickstep kind.