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