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