(git:ce05c02)
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 atomic_kind_types, ONLY: get_atomic_kind
14 USE basis_set_types, ONLY: get_gto_basis_set,&
15 gto_basis_set_type
16 USE cp_array_utils, ONLY: cp_1d_r_p_type
17 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
18 USE cp_fm_types, ONLY: cp_fm_get_info,&
19 cp_fm_get_submatrix,&
20 cp_fm_type
21 USE cp_log_handling, ONLY: cp_get_default_logger,&
22 cp_logger_type
23 USE cp_output_handling, ONLY: cp_p_file,&
24 cp_print_key_finished_output,&
25 cp_print_key_should_output,&
26 cp_print_key_unit_nr
27 USE input_constants, ONLY: gto_cartesian,&
28 gto_spherical
29 USE input_section_types, ONLY: section_vals_type,&
30 section_vals_val_get
31 USE kinds, ONLY: dp
32 USE mathconstants, ONLY: pi
33 USE orbital_pointers, ONLY: nco,&
34 nso
35 USE orbital_transformation_matrices, ONLY: orbtramat
36 USE particle_types, ONLY: particle_type
37 USE periodic_table, ONLY: get_ptable_info
38 USE physcon, ONLY: massunit
39 USE qs_kind_types, ONLY: get_qs_kind,&
40 get_qs_kind_set,&
41 qs_kind_type
42 USE qs_mo_types, ONLY: mo_set_type
43#include "./base/base_uses.f90"
44
45 IMPLICIT NONE
46
47 PRIVATE
48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molden_utils'
49 LOGICAL, PARAMETER :: debug_this_module = .false.
50
51 INTEGER, PARAMETER :: molden_lmax = 4
52 INTEGER, PARAMETER :: molden_ncomax = (molden_lmax + 1)*(molden_lmax + 2)/2 ! 15
53
54 PUBLIC :: write_vibrations_molden, write_mos_molden
55
56CONTAINS
57
58! **************************************************************************************************
59!> \brief Write out the MOs in molden format for visualisation
60!> \param mos the set of MOs (both spins, if UKS)
61!> \param qs_kind_set for basis set info
62!> \param particle_set particles data structure, for positions and kinds
63!> \param print_section input section containing relevant print key
64!> \param unoccupied_orbs optional: unoccupied orbital coefficients from make_lumo_gpw
65!> \param unoccupied_evals optional: unoccupied orbital eigenvalues
66!> \author MattW, IainB
67! **************************************************************************************************
68 SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section, &
69 unoccupied_orbs, unoccupied_evals)
70 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
71 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
72 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
73 TYPE(section_vals_type), POINTER :: print_section
74 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN), &
75 OPTIONAL :: unoccupied_orbs
76 TYPE(cp_1d_r_p_type), DIMENSION(:), INTENT(IN), &
77 OPTIONAL :: unoccupied_evals
78
79 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_mos_molden'
80 CHARACTER(LEN=molden_lmax+1), PARAMETER :: angmom = "spdfg"
81
82 CHARACTER(LEN=15) :: fmtstr1, fmtstr2
83 CHARACTER(LEN=2) :: element_symbol
84 INTEGER :: gto_kind, handle, i, iatom, icgf, icol, ikind, ipgf, irow, irow_in, iset, isgf, &
85 ishell, ispin, iw, lshell, ncgf, ncol_global, ndigits, nrow_global, nset, nsgf, numos, z
86 INTEGER, DIMENSION(:), POINTER :: npgf, nshell
87 INTEGER, DIMENSION(:, :), POINTER :: l
88 INTEGER, DIMENSION(molden_ncomax, 0:molden_lmax) :: orbmap
89 LOGICAL :: print_warn
90 REAL(kind=dp) :: expzet, prefac
91 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
92 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
93 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
94 TYPE(cp_logger_type), POINTER :: logger
95 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
96
97 CALL timeset(routinen, handle)
98
99 logger => cp_get_default_logger()
100 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, ""), cp_p_file)) THEN
101
102 iw = cp_print_key_unit_nr(logger, print_section, "", &
103 extension=".molden", file_status='REPLACE')
104
105 print_warn = .true.
106
107 CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
108 ndigits = min(max(3, ndigits), 30)
109 WRITE (unit=fmtstr1, fmt='("(I6,1X,ES",I0,".",I0,")")') ndigits + 7, ndigits
110 WRITE (unit=fmtstr2, fmt='("((T51,2F",I0,".",I0,"))")') ndigits + 10, ndigits
111
112 CALL section_vals_val_get(print_section, "GTO_KIND", i_val=gto_kind)
113
114 IF (mos(1)%use_mo_coeff_b) THEN
115 ! we are using the dbcsr mo_coeff
116 ! we copy it to the fm anyway
117 DO ispin = 1, SIZE(mos)
118 IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
119 cpassert(.false.)
120 END IF
121 CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
122 mos(ispin)%mo_coeff) !fm->dbcsr
123 END DO
124 END IF
125
126 IF (iw > 0) THEN
127 WRITE (iw, '(T2,A)') "[Molden Format]"
128 WRITE (iw, '(T2,A)') "[Atoms] AU"
129 DO i = 1, SIZE(particle_set)
130 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
131 element_symbol=element_symbol)
132 CALL get_ptable_info(element_symbol, number=z)
133
134 WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
135 element_symbol, i, z, particle_set(i)%r(:)
136 END DO
137
138 WRITE (iw, '(T2,A)') "[GTO]"
139
140 DO i = 1, SIZE(particle_set)
141 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
142 element_symbol=element_symbol)
143 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
144 IF (ASSOCIATED(orb_basis_set)) THEN
145 WRITE (iw, '(T2,I8,I8)') i, 0
146 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
147 nset=nset, &
148 npgf=npgf, &
149 nshell=nshell, &
150 l=l, &
151 zet=zet, &
152 gcc=gcc)
153
154 DO iset = 1, nset
155 DO ishell = 1, nshell(iset)
156 lshell = l(ishell, iset)
157 IF (lshell <= molden_lmax) THEN
158 WRITE (unit=iw, fmt='(T25,A2,4X,I4,4X,F4.2)') &
159 angmom(lshell + 1:lshell + 1), npgf(iset), 1.0_dp
160 ! MOLDEN expects the contraction coefficient of spherical NOT CARTESIAN NORMALISED
161 ! functions. So we undo the normalisation factors included in the gccs
162 ! Reverse engineered from basis_set_types, normalise_gcc_orb
163 prefac = 2_dp**lshell*(2/pi)**0.75_dp
164 expzet = 0.25_dp*(2*lshell + 3.0_dp)
165 WRITE (unit=iw, fmt=fmtstr2) &
166 (zet(ipgf, iset), gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet), &
167 ipgf=1, npgf(iset))
168 ELSE
169 IF (print_warn) THEN
170 CALL cp_warn(__location__, &
171 "MOLDEN format does not support Gaussian orbitals with l > 4.")
172 print_warn = .false.
173 END IF
174 END IF
175 END DO
176 END DO
177
178 WRITE (iw, '(A4)') " "
179
180 END IF
181
182 END DO
183
184 IF (gto_kind == gto_spherical) THEN
185 WRITE (iw, '(T2,A)') "[5D7F]"
186 WRITE (iw, '(T2,A)') "[9G]"
187 END IF
188
189 WRITE (iw, '(T2,A)') "[MO]"
190 END IF
191
192 !------------------------------------------------------------------------
193 ! convert from CP2K to MOLDEN format ordering
194 ! http://www.cmbi.ru.nl/molden/molden_format.html
195 !"The following order of D, F and G functions is expected:
196 !
197 ! 5D: D 0, D+1, D-1, D+2, D-2
198 ! 6D: xx, yy, zz, xy, xz, yz
199 !
200 ! 7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
201 ! 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
202 !
203 ! 9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
204 ! 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
205 ! xxyy xxzz yyzz xxyz yyxz zzxy
206 !"
207 ! CP2K has x in the outer (slower loop), so
208 ! xx, xy, xz, yy, yz,zz for l=2, for instance
209 !
210 ! iorb_cp2k = orbmap(iorb_molden, l), l = 0 .. 4
211 ! -----------------------------------------------------------------------
212 IF (iw > 0) THEN
213 IF (gto_kind == gto_cartesian) THEN
214 ! -----------------------------------------------------------------
215 ! Use cartesian (6D, 10F, 15G) representation.
216 ! This is only format VMD can process.
217 ! -----------------------------------------------------------------
218 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
219 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
220 1, 4, 6, 2, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
221 1, 7, 10, 4, 2, 3, 6, 9, 8, 5, 0, 0, 0, 0, 0, &
222 1, 11, 15, 2, 3, 7, 12, 10, 14, 4, 6, 13, 5, 8, 9], &
223 [molden_ncomax, molden_lmax + 1])
224 ELSE IF (gto_kind == gto_spherical) THEN
225 ! -----------------------------------------------------------------
226 ! Use spherical (5D, 7F, 9G) representation.
227 ! -----------------------------------------------------------------
228 orbmap = reshape([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
229 3, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
230 3, 4, 2, 5, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
231 4, 5, 3, 6, 2, 7, 1, 0, 0, 0, 0, 0, 0, 0, 0, &
232 5, 6, 4, 7, 3, 8, 2, 9, 1, 0, 0, 0, 0, 0, 0], &
233 [molden_ncomax, molden_lmax + 1])
234 END IF
235 END IF
236
237 DO ispin = 1, SIZE(mos)
238 CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
239 nrow_global=nrow_global, &
240 ncol_global=ncol_global)
241 ALLOCATE (smatrix(nrow_global, ncol_global))
242 CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
243
244 IF (iw > 0) THEN
245 IF (gto_kind == gto_cartesian) THEN
246 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
247
248 ALLOCATE (cmatrix(ncgf, ncgf))
249
250 cmatrix = 0.0_dp
251
252 ! Transform spherical MOs to Cartesian MOs
253
254 icgf = 1
255 isgf = 1
256 DO iatom = 1, SIZE(particle_set)
257 NULLIFY (orb_basis_set)
258 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
259 CALL get_qs_kind(qs_kind_set(ikind), &
260 basis_set=orb_basis_set)
261 IF (ASSOCIATED(orb_basis_set)) THEN
262 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
263 nset=nset, &
264 nshell=nshell, &
265 l=l)
266 DO iset = 1, nset
267 DO ishell = 1, nshell(iset)
268 lshell = l(ishell, iset)
269 CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, nso(lshell), 1.0_dp, &
270 orbtramat(lshell)%c2s, nso(lshell), &
271 smatrix(isgf, 1), nsgf, 0.0_dp, &
272 cmatrix(icgf, 1), ncgf)
273 icgf = icgf + nco(lshell)
274 isgf = isgf + nso(lshell)
275 END DO
276 END DO
277 END IF
278 END DO ! iatom
279 END IF
280
281 DO icol = 1, mos(ispin)%nmo
282 ! index of the first basis function for the given atom, set, and shell
283 irow = 1
284
285 ! index of the first basis function in MOLDEN file.
286 ! Due to limitation of the MOLDEN format, basis functions with l > molden_lmax
287 ! cannot be exported, so we need to renumber atomic orbitals
288 irow_in = 1
289
290 WRITE (iw, '(A,ES20.10)') 'Ene=', mos(ispin)%eigenvalues(icol)
291 IF (ispin < 2) THEN
292 WRITE (iw, '(A)') 'Spin= Alpha'
293 ELSE
294 WRITE (iw, '(A)') 'Spin= Beta'
295 END IF
296 WRITE (iw, '(A,F12.7)') 'Occup=', mos(ispin)%occupation_numbers(icol)
297
298 DO iatom = 1, SIZE(particle_set)
299 NULLIFY (orb_basis_set)
300 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
301 element_symbol=element_symbol, kind_number=ikind)
302 CALL get_qs_kind(qs_kind_set(ikind), &
303 basis_set=orb_basis_set)
304 IF (ASSOCIATED(orb_basis_set)) THEN
305 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
306 nset=nset, &
307 nshell=nshell, &
308 l=l)
309
310 IF (gto_kind == gto_cartesian) THEN
311 ! ----------------------------------------------
312 ! Use cartesian (6D, 10F, 15G) representation.
313 ! ----------------------------------------------
314 icgf = 1
315 DO iset = 1, nset
316 DO ishell = 1, nshell(iset)
317 lshell = l(ishell, iset)
318
319 IF (lshell <= molden_lmax) THEN
320 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
321 cmatrix(irow:irow + nco(lshell) - 1, icol))
322 irow_in = irow_in + nco(lshell)
323 END IF
324
325 irow = irow + nco(lshell)
326 END DO ! ishell
327 END DO
328
329 ELSE IF (gto_kind == gto_spherical) THEN
330 ! ----------------------------------------------
331 ! Use spherical (5D, 7F, 9G) representation.
332 ! ----------------------------------------------
333 DO iset = 1, nset
334 DO ishell = 1, nshell(iset)
335 lshell = l(ishell, iset)
336
337 IF (lshell <= molden_lmax) THEN
338 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
339 smatrix(irow:irow + nso(lshell) - 1, icol))
340 irow_in = irow_in + nso(lshell)
341 END IF
342
343 irow = irow + nso(lshell)
344 END DO
345 END DO
346 END IF
347
348 END IF
349 END DO ! iatom
350 END DO
351 END IF
352
353 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
354 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
355 END DO
356
357 ! Write unoccupied (virtual) orbitals if provided; only used with OT
358 IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN
359 DO ispin = 1, SIZE(unoccupied_orbs)
360 CALL cp_fm_get_info(unoccupied_orbs(ispin), &
361 nrow_global=nrow_global, &
362 ncol_global=numos)
363 ALLOCATE (smatrix(nrow_global, numos))
364 CALL cp_fm_get_submatrix(unoccupied_orbs(ispin), smatrix)
365
366 IF (iw > 0) THEN
367 IF (gto_kind == gto_cartesian) THEN
368 CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
369 ALLOCATE (cmatrix(ncgf, numos))
370 cmatrix = 0.0_dp
371
372 icgf = 1
373 isgf = 1
374 DO iatom = 1, SIZE(particle_set)
375 NULLIFY (orb_basis_set)
376 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
377 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
378 IF (ASSOCIATED(orb_basis_set)) THEN
379 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
380 nset=nset, nshell=nshell, l=l)
381 DO iset = 1, nset
382 DO ishell = 1, nshell(iset)
383 lshell = l(ishell, iset)
384 CALL dgemm("T", "N", nco(lshell), numos, nso(lshell), 1.0_dp, &
385 orbtramat(lshell)%c2s, nso(lshell), &
386 smatrix(isgf, 1), nsgf, 0.0_dp, &
387 cmatrix(icgf, 1), ncgf)
388 icgf = icgf + nco(lshell)
389 isgf = isgf + nso(lshell)
390 END DO
391 END DO
392 END IF
393 END DO
394 END IF
395
396 DO icol = 1, numos
397 irow = 1
398 irow_in = 1
399
400 WRITE (iw, '(A,ES20.10)') 'Ene=', unoccupied_evals(ispin)%array(icol)
401 IF (ispin < 2) THEN
402 WRITE (iw, '(A)') 'Spin= Alpha'
403 ELSE
404 WRITE (iw, '(A)') 'Spin= Beta'
405 END IF
406 WRITE (iw, '(A,F12.7)') 'Occup=', 0.0_dp
407
408 DO iatom = 1, SIZE(particle_set)
409 NULLIFY (orb_basis_set)
410 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
411 element_symbol=element_symbol, kind_number=ikind)
412 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
413 IF (ASSOCIATED(orb_basis_set)) THEN
414 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
415 nset=nset, nshell=nshell, l=l)
416
417 IF (gto_kind == gto_cartesian) THEN
418 icgf = 1
419 DO iset = 1, nset
420 DO ishell = 1, nshell(iset)
421 lshell = l(ishell, iset)
422 IF (lshell <= molden_lmax) THEN
423 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
424 cmatrix(irow:irow + nco(lshell) - 1, icol))
425 irow_in = irow_in + nco(lshell)
426 END IF
427 irow = irow + nco(lshell)
428 END DO
429 END DO
430 ELSE IF (gto_kind == gto_spherical) THEN
431 DO iset = 1, nset
432 DO ishell = 1, nshell(iset)
433 lshell = l(ishell, iset)
434 IF (lshell <= molden_lmax) THEN
435 CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
436 smatrix(irow:irow + nso(lshell) - 1, icol))
437 irow_in = irow_in + nso(lshell)
438 END IF
439 irow = irow + nso(lshell)
440 END DO
441 END DO
442 END IF
443
444 END IF
445 END DO ! iatom
446 END DO ! icol
447 END IF
448
449 IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
450 IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
451 END DO ! ispin
452 END IF
453
454 CALL cp_print_key_finished_output(iw, logger, print_section, "")
455
456 END IF
457
458 CALL timestop(handle)
459
460 END SUBROUTINE write_mos_molden
461
462! **************************************************************************************************
463!> \brief Output MO coefficients formatted correctly for MOLDEN, omitting those <= 1E(-digits)
464!> \param iw output file unit
465!> \param fmtstr1 format string
466!> \param ndigits number of significant digits in MO coefficients
467!> \param irow_in index of the first atomic orbital: mo_coeff(orbmap(1))
468!> \param orbmap array to map Gaussian functions from MOLDEN to CP2K ordering
469!> \param mo_coeff MO coefficients
470! **************************************************************************************************
471 SUBROUTINE print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap, mo_coeff)
472 INTEGER, INTENT(in) :: iw
473 CHARACTER(LEN=*), INTENT(in) :: fmtstr1
474 INTEGER, INTENT(in) :: ndigits, irow_in
475 INTEGER, DIMENSION(molden_ncomax), INTENT(in) :: orbmap
476 REAL(kind=dp), DIMENSION(:), INTENT(in) :: mo_coeff
477
478 INTEGER :: orbital
479
480 DO orbital = 1, molden_ncomax
481 IF (orbmap(orbital) /= 0) THEN
482 IF (abs(mo_coeff(orbmap(orbital))) >= 10.0_dp**(-ndigits)) THEN
483 WRITE (iw, fmtstr1) irow_in + orbital - 1, mo_coeff(orbmap(orbital))
484 END IF
485 END IF
486 END DO
487
488 END SUBROUTINE print_coeffs
489
490! **************************************************************************************************
491!> \brief writes the output for vibrational analysis in MOLDEN format
492!> \param input ...
493!> \param particles ...
494!> \param freq ...
495!> \param eigen_vec ...
496!> \param intensities ...
497!> \param calc_intens ...
498!> \param dump_only_positive ...
499!> \param logger ...
500!> \param list array of mobile atom indices
501!> \author Florian Schiffmann 11.2007
502! **************************************************************************************************
503 SUBROUTINE write_vibrations_molden(input, particles, freq, eigen_vec, intensities, calc_intens, &
504 dump_only_positive, logger, list)
505
506 TYPE(section_vals_type), POINTER :: input
507 TYPE(particle_type), DIMENSION(:), POINTER :: particles
508 REAL(kind=dp), DIMENSION(:) :: freq
509 REAL(kind=dp), DIMENSION(:, :) :: eigen_vec
510 REAL(kind=dp), DIMENSION(:), POINTER :: intensities
511 LOGICAL, INTENT(in) :: calc_intens, dump_only_positive
512 TYPE(cp_logger_type), POINTER :: logger
513 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: list
514
515 CHARACTER(len=*), PARAMETER :: routinen = 'write_vibrations_molden'
516
517 CHARACTER(LEN=2) :: element_symbol
518 INTEGER :: handle, i, iw, j, k, l, z
519 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
520 REAL(kind=dp) :: fint
521
522 CALL timeset(routinen, handle)
523
524 iw = cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
525 extension=".mol", file_status='REPLACE')
526
527 IF (iw > 0) THEN
528 cpassert(mod(SIZE(eigen_vec, 1), 3) == 0)
529 cpassert(SIZE(freq, 1) == SIZE(eigen_vec, 2))
530 ALLOCATE (my_list(SIZE(particles)))
531 ! Either we have a list of the subset of mobile atoms,
532 ! Or the eigenvectors must span the full space (all atoms)
533 IF (PRESENT(list)) THEN
534 my_list(:) = 0
535 DO i = 1, SIZE(list)
536 my_list(list(i)) = i
537 END DO
538 ELSE
539 cpassert(SIZE(particles) == SIZE(eigen_vec, 1)/3)
540 DO i = 1, SIZE(my_list)
541 my_list(i) = i
542 END DO
543 END IF
544 WRITE (iw, '(T2,A)') "[Molden Format]"
545 WRITE (iw, '(T2,A)') "[Atoms] AU"
546 DO i = 1, SIZE(particles)
547 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
548 element_symbol=element_symbol)
549 CALL get_ptable_info(element_symbol, number=z)
550
551 WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
552 element_symbol, i, z, particles(i)%r(:)
553
554 END DO
555 WRITE (iw, '(T2,A)') "[FREQ]"
556 DO i = 1, SIZE(freq, 1)
557 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(T5,F12.6)') freq(i)
558 END DO
559 WRITE (iw, '(T2,A)') "[FR-COORD]"
560 DO i = 1, SIZE(particles)
561 CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
562 element_symbol=element_symbol)
563 WRITE (iw, '(T2,A2,3X,3(F12.6,3X))') &
564 element_symbol, particles(i)%r(:)
565 END DO
566 WRITE (iw, '(T2,A)') "[FR-NORM-COORD]"
567 l = 0
568 DO i = 1, SIZE(eigen_vec, 2)
569 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) THEN
570 l = l + 1
571 WRITE (iw, '(T2,A,1X,I6)') "vibration", l
572 DO j = 1, SIZE(particles)
573 IF (my_list(j) /= 0) THEN
574 k = (my_list(j) - 1)*3
575 WRITE (iw, '(T2,3(F12.6,3X))') eigen_vec(k + 1, i), eigen_vec(k + 2, i), eigen_vec(k + 3, i)
576 ELSE
577 WRITE (iw, '(T2,3(F12.6,3X))') 0.0_dp, 0.0_dp, 0.0_dp
578 END IF
579 END DO
580 END IF
581 END DO
582 IF (calc_intens) THEN
583 fint = massunit
584 ! intensity units are a.u./amu
585 WRITE (iw, '(T2,A)') "[INT]"
586 DO i = 1, SIZE(intensities)
587 IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(3X,F18.6)') fint*intensities(i)**2
588 END DO
589 END IF
590 DEALLOCATE (my_list)
591 END IF
592 CALL cp_print_key_finished_output(iw, logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
593
594 CALL timestop(handle)
595
596 END SUBROUTINE write_vibrations_molden
597
598END 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:214
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
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_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:1295
input_constants::gto_spherical
integer, parameter, public gto_spherical
Definition input_constants.F:1295
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, unoccupied_orbs, unoccupied_evals)
Write out the MOs in molden format for visualisation.
Definition molden_utils.F:70
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:505
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::massunit
real(kind=dp), parameter, public massunit
Definition physcon.F:141
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_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
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:128
particle_types::particle_type
Definition particle_types.F:35
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