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 admm_types, ONLY: admm_type
14 : USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
15 : admm_uncorrect_for_eigenvalues
16 : USE atomic_kind_types, ONLY: get_atomic_kind
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE cell_types, ONLY: cell_type
20 : USE cp_array_utils, ONLY: cp_1d_r_p_type
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
23 : dbcsr_type
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
25 : USE cp_fm_types, ONLY: cp_fm_get_info,&
26 : cp_fm_get_submatrix,&
27 : cp_fm_type
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_type
30 : USE cp_output_handling, ONLY: cp_p_file,&
31 : cp_print_key_finished_output,&
32 : cp_print_key_should_output,&
33 : cp_print_key_unit_nr
34 : USE input_constants, ONLY: gto_cartesian,&
35 : gto_spherical
36 : USE input_section_types, ONLY: section_vals_type,&
37 : section_vals_val_get
38 : USE kinds, ONLY: dp
39 : USE mathconstants, ONLY: pi
40 : USE orbital_pointers, ONLY: nco,&
41 : nso
42 : USE orbital_transformation_matrices, ONLY: orbtramat
43 : USE particle_types, ONLY: particle_type
44 : USE periodic_table, ONLY: get_ptable_info
45 : USE physcon, ONLY: angstrom,&
46 : massunit
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : USE qs_kind_types, ONLY: get_qs_kind,&
50 : get_qs_kind_set,&
51 : qs_kind_type
52 : USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
53 : USE qs_mo_types, ONLY: get_mo_set,&
54 : mo_set_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molden_utils'
61 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
62 :
63 : INTEGER, PARAMETER :: molden_lmax = 4
64 : INTEGER, PARAMETER :: molden_ncomax = (molden_lmax + 1)*(molden_lmax + 2)/2 ! 15
65 :
66 : PUBLIC :: write_vibrations_molden, write_mos_molden
67 :
68 : CONTAINS
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 0 : 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 0 : IF (unit_choice == 2) THEN
84 0 : scale_factor = angstrom
85 0 : WRITE (iw, '(T2,A)') "[Cell] Angs"
86 : ELSE
87 0 : scale_factor = 1.0_dp
88 0 : WRITE (iw, '(T2,A)') "[Cell] AU"
89 : END IF
90 : WRITE (iw, '(T2,3(F12.6,3X))') &
91 0 : 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 0 : 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 0 : cell%hmat(1, 3)*scale_factor, cell%hmat(2, 3)*scale_factor, cell%hmat(3, 3)*scale_factor
96 0 : 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 12021 : SUBROUTINE write_mos_molden(mos, qs_kind_set, particle_set, print_section, cell, &
112 12021 : 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 12021 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
134 12021 : 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 12021 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmatrix, smatrix
141 12021 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
142 12021 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
143 12021 : 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 12021 : 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 12021 : CALL timeset(routineN, handle)
153 :
154 12021 : logger => cp_get_default_logger()
155 12021 : 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 22 : extension=".molden", file_status='REPLACE')
159 :
160 22 : print_warn = .TRUE.
161 :
162 22 : CALL section_vals_val_get(print_section, "UNIT", i_val=unit_choice)
163 22 : IF (unit_choice == 2) THEN
164 : scale_factor = angstrom
165 : ELSE
166 22 : scale_factor = 1.0_dp
167 : END IF
168 :
169 22 : CALL section_vals_val_get(print_section, "NDIGITS", i_val=ndigits)
170 22 : ndigits = MIN(MAX(3, ndigits), 30)
171 22 : WRITE (UNIT=fmtstr1, FMT='("(I6,1X,ES",I0,".",I0,")")') ndigits + 7, ndigits
172 22 : WRITE (UNIT=fmtstr2, FMT='("((T51,2F",I0,".",I0,"))")') ndigits + 10, ndigits
173 :
174 22 : CALL section_vals_val_get(print_section, "GTO_KIND", i_val=gto_kind)
175 22 : CALL section_vals_val_get(print_section, "WRITE_CELL", l_val=write_cell)
176 22 : CALL section_vals_val_get(print_section, "WRITE_PSEUDO", l_val=write_pseudo)
177 22 : CALL section_vals_val_get(print_section, "MARK_GHOST", l_val=mark_ghost)
178 :
179 22 : 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 0 : DO ispin = 1, SIZE(mos)
183 0 : CPASSERT(ASSOCIATED(mos(ispin)%mo_coeff_b))
184 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
185 0 : mos(ispin)%mo_coeff) !fm->dbcsr
186 : END DO
187 : END IF
188 :
189 22 : IF (iw > 0) THEN
190 11 : WRITE (iw, '(T2,A)') "[Molden Format]"
191 11 : IF (write_cell) THEN
192 0 : CPASSERT(PRESENT(cell))
193 0 : CPASSERT(ASSOCIATED(cell))
194 0 : CALL write_cell_molden(iw, cell, unit_choice)
195 : END IF
196 11 : IF (unit_choice == 2) THEN
197 0 : WRITE (iw, '(T2,A)') "[Atoms] Angs"
198 : ELSE
199 11 : WRITE (iw, '(T2,A)') "[Atoms] AU"
200 : END IF
201 144 : DO i = 1, SIZE(particle_set)
202 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
203 133 : element_symbol=element_symbol)
204 133 : CALL get_ptable_info(element_symbol, number=z)
205 133 : IF (mark_ghost) THEN
206 0 : CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_atom)
207 0 : IF (ghost_atom) z = 0
208 : END IF
209 :
210 : WRITE (iw, '(T2,A2,I6,I6,3X,3(F12.6,3X))') &
211 676 : element_symbol, i, z, particle_set(i)%r(:)*scale_factor
212 : END DO
213 11 : IF (write_pseudo) THEN
214 0 : WRITE (iw, '(T2,A)') "[Pseudo]"
215 0 : DO i = 1, SIZE(particle_set)
216 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
217 0 : element_symbol=element_symbol)
218 0 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
219 : WRITE (iw, '(T2,A2,I6,I6)') &
220 0 : element_symbol, i, NINT(zeff)
221 : END DO
222 : END IF
223 :
224 11 : WRITE (iw, '(T2,A)') "[GTO]"
225 :
226 144 : DO i = 1, SIZE(particle_set)
227 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=ikind, &
228 133 : element_symbol=element_symbol)
229 133 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
230 277 : IF (ASSOCIATED(orb_basis_set)) THEN
231 133 : 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 133 : gcc=gcc)
239 :
240 432 : DO iset = 1, nset
241 765 : DO ishell = 1, nshell(iset)
242 333 : lshell = l(ishell, iset)
243 632 : IF (lshell <= molden_lmax) THEN
244 : WRITE (UNIT=iw, FMT='(T25,A2,4X,I4,4X,F4.2)') &
245 333 : 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 333 : prefac = 2_dp**lshell*(2/pi)**0.75_dp
250 333 : expzet = 0.25_dp*(2*lshell + 3.0_dp)
251 : WRITE (UNIT=iw, FMT=fmtstr2) &
252 2156 : (zet(ipgf, iset), gcc(ipgf, ishell, iset)/(prefac*zet(ipgf, iset)**expzet), &
253 2489 : ipgf=1, npgf(iset))
254 : ELSE
255 0 : IF (print_warn) THEN
256 : CALL cp_warn(__LOCATION__, &
257 0 : "MOLDEN format does not support Gaussian orbitals with l > 4.")
258 0 : print_warn = .FALSE.
259 : END IF
260 : END IF
261 : END DO
262 : END DO
263 :
264 133 : WRITE (iw, '(A4)') " "
265 :
266 : END IF
267 :
268 : END DO
269 :
270 11 : IF (gto_kind == gto_spherical) THEN
271 11 : WRITE (iw, '(T2,A)') "[5D7F]"
272 11 : WRITE (iw, '(T2,A)') "[9G]"
273 : END IF
274 :
275 11 : 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 22 : IF (iw > 0) THEN
299 11 : 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 0 : [molden_ncomax, molden_lmax + 1])
310 11 : 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 11 : [molden_ncomax, molden_lmax + 1])
320 : END IF
321 : END IF
322 :
323 52 : DO ispin = 1, SIZE(mos)
324 30 : do_calc_energies = .FALSE.
325 30 : IF (PRESENT(calc_energies)) do_calc_energies = calc_energies
326 :
327 30 : IF (PRESENT(qs_env) .AND. do_calc_energies) THEN
328 4 : CALL get_qs_env(qs_env, matrix_ks=ks, dft_control=dft_control)
329 :
330 4 : matrix_ks => ks(ispin)%matrix
331 :
332 : ! With ADMM, we have to modify the Kohn-Sham matrix
333 4 : IF (dft_control%do_admm) THEN
334 0 : CALL get_qs_env(qs_env, admm_env=admm_env)
335 0 : CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks)
336 : END IF
337 :
338 4 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues)
339 :
340 4 : IF (ASSOCIATED(qs_env%mo_derivs)) THEN
341 0 : mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
342 : ELSE
343 4 : 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 4 : co_rotate_dbcsr=mo_coeff_deriv)
351 :
352 : ! With ADMM, we have to undo the modification of the Kohn-Sham matrix
353 4 : IF (dft_control%do_admm) THEN
354 0 : 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 30 : ncol_global=ncol_global)
361 120 : ALLOCATE (smatrix(nrow_global, ncol_global))
362 30 : CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
363 :
364 30 : IF (iw > 0) THEN
365 15 : IF (gto_kind == gto_cartesian) THEN
366 0 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
367 :
368 0 : ALLOCATE (cmatrix(ncgf, ncgf))
369 :
370 0 : cmatrix = 0.0_dp
371 :
372 : ! Transform spherical MOs to Cartesian MOs
373 :
374 0 : icgf = 1
375 0 : isgf = 1
376 0 : DO iatom = 1, SIZE(particle_set)
377 0 : NULLIFY (orb_basis_set)
378 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
379 : CALL get_qs_kind(qs_kind_set(ikind), &
380 0 : basis_set=orb_basis_set)
381 0 : 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 0 : l=l)
386 0 : DO iset = 1, nset
387 0 : DO ishell = 1, nshell(iset)
388 0 : 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 0 : cmatrix(icgf, 1), ncgf)
393 0 : icgf = icgf + nco(lshell)
394 0 : isgf = isgf + nso(lshell)
395 : END DO
396 : END DO
397 : END IF
398 : END DO ! iatom
399 : END IF
400 :
401 81 : DO icol = 1, mos(ispin)%nmo
402 : ! index of the first basis function for the given atom, set, and shell
403 66 : 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 66 : irow_in = 1
409 :
410 66 : WRITE (iw, '(A,ES20.10)') 'Ene=', mos(ispin)%eigenvalues(icol)
411 66 : IF (ispin < 2) THEN
412 61 : WRITE (iw, '(A)') 'Spin= Alpha'
413 : ELSE
414 5 : WRITE (iw, '(A)') 'Spin= Beta'
415 : END IF
416 66 : WRITE (iw, '(A,F12.7)') 'Occup=', mos(ispin)%occupation_numbers(icol)
417 :
418 539 : DO iatom = 1, SIZE(particle_set)
419 458 : NULLIFY (orb_basis_set)
420 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
421 458 : element_symbol=element_symbol, kind_number=ikind)
422 : CALL get_qs_kind(qs_kind_set(ikind), &
423 458 : basis_set=orb_basis_set)
424 982 : 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 458 : l=l)
429 :
430 458 : IF (gto_kind == gto_cartesian) THEN
431 : ! ----------------------------------------------
432 : ! Use cartesian (6D, 10F, 15G) representation.
433 : ! ----------------------------------------------
434 0 : icgf = 1
435 0 : DO iset = 1, nset
436 0 : DO ishell = 1, nshell(iset)
437 0 : lshell = l(ishell, iset)
438 :
439 0 : IF (lshell <= molden_lmax) THEN
440 : CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
441 0 : cmatrix(irow:irow + nco(lshell) - 1, icol))
442 0 : irow_in = irow_in + nco(lshell)
443 : END IF
444 :
445 0 : irow = irow + nco(lshell)
446 : END DO ! ishell
447 : END DO
448 :
449 458 : ELSE IF (gto_kind == gto_spherical) THEN
450 : ! ----------------------------------------------
451 : ! Use spherical (5D, 7F, 9G) representation.
452 : ! ----------------------------------------------
453 1704 : DO iset = 1, nset
454 3130 : DO ishell = 1, nshell(iset)
455 1426 : lshell = l(ishell, iset)
456 :
457 1426 : IF (lshell <= molden_lmax) THEN
458 : CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
459 1426 : smatrix(irow:irow + nso(lshell) - 1, icol))
460 1426 : irow_in = irow_in + nso(lshell)
461 : END IF
462 :
463 2672 : 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 30 : IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
474 82 : IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
475 : END DO
476 :
477 : ! Write unoccupied (virtual) orbitals if provided; only used with OT
478 22 : IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN
479 0 : DO ispin = 1, SIZE(unoccupied_orbs)
480 : CALL cp_fm_get_info(unoccupied_orbs(ispin), &
481 : nrow_global=nrow_global, &
482 0 : ncol_global=numos)
483 0 : ALLOCATE (smatrix(nrow_global, numos))
484 0 : CALL cp_fm_get_submatrix(unoccupied_orbs(ispin), smatrix)
485 :
486 0 : IF (iw > 0) THEN
487 0 : IF (gto_kind == gto_cartesian) THEN
488 0 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
489 0 : ALLOCATE (cmatrix(ncgf, numos))
490 0 : cmatrix = 0.0_dp
491 :
492 0 : icgf = 1
493 0 : isgf = 1
494 0 : DO iatom = 1, SIZE(particle_set)
495 0 : NULLIFY (orb_basis_set)
496 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
497 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
498 0 : IF (ASSOCIATED(orb_basis_set)) THEN
499 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
500 0 : nset=nset, nshell=nshell, l=l)
501 0 : DO iset = 1, nset
502 0 : DO ishell = 1, nshell(iset)
503 0 : 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 0 : cmatrix(icgf, 1), ncgf)
508 0 : icgf = icgf + nco(lshell)
509 0 : isgf = isgf + nso(lshell)
510 : END DO
511 : END DO
512 : END IF
513 : END DO
514 : END IF
515 :
516 0 : DO icol = 1, numos
517 0 : irow = 1
518 0 : irow_in = 1
519 :
520 0 : WRITE (iw, '(A,ES20.10)') 'Ene=', unoccupied_evals(ispin)%array(icol)
521 0 : IF (ispin < 2) THEN
522 0 : WRITE (iw, '(A)') 'Spin= Alpha'
523 : ELSE
524 0 : WRITE (iw, '(A)') 'Spin= Beta'
525 : END IF
526 0 : WRITE (iw, '(A,F12.7)') 'Occup=', 0.0_dp
527 :
528 0 : DO iatom = 1, SIZE(particle_set)
529 0 : NULLIFY (orb_basis_set)
530 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
531 0 : element_symbol=element_symbol, kind_number=ikind)
532 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
533 0 : IF (ASSOCIATED(orb_basis_set)) THEN
534 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
535 0 : nset=nset, nshell=nshell, l=l)
536 :
537 0 : IF (gto_kind == gto_cartesian) THEN
538 0 : icgf = 1
539 0 : DO iset = 1, nset
540 0 : DO ishell = 1, nshell(iset)
541 0 : lshell = l(ishell, iset)
542 0 : IF (lshell <= molden_lmax) THEN
543 : CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
544 0 : cmatrix(irow:irow + nco(lshell) - 1, icol))
545 0 : irow_in = irow_in + nco(lshell)
546 : END IF
547 0 : irow = irow + nco(lshell)
548 : END DO
549 : END DO
550 0 : ELSE IF (gto_kind == gto_spherical) THEN
551 0 : DO iset = 1, nset
552 0 : DO ishell = 1, nshell(iset)
553 0 : lshell = l(ishell, iset)
554 0 : IF (lshell <= molden_lmax) THEN
555 : CALL print_coeffs(iw, fmtstr1, ndigits, irow_in, orbmap(:, lshell), &
556 0 : smatrix(irow:irow + nso(lshell) - 1, icol))
557 0 : irow_in = irow_in + nso(lshell)
558 : END IF
559 0 : 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 0 : IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
570 0 : IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
571 : END DO ! ispin
572 : END IF
573 :
574 22 : CALL cp_print_key_finished_output(iw, logger, print_section, "")
575 :
576 : END IF
577 :
578 12021 : CALL timestop(handle)
579 :
580 24042 : 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 1426 : 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 22816 : DO orbital = 1, molden_ncomax
601 22816 : IF (orbmap(orbital) /= 0) THEN
602 2686 : IF (ABS(mo_coeff(orbmap(orbital))) >= 10.0_dp**(-ndigits)) THEN
603 1692 : WRITE (iw, fmtstr1) irow_in + orbital - 1, mo_coeff(orbmap(orbital))
604 : END IF
605 : END IF
606 : END DO
607 :
608 1426 : 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 60 : 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 60 : INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
642 : LOGICAL :: write_cell
643 : REAL(KIND=dp) :: fint
644 :
645 60 : CALL timeset(routineN, handle)
646 :
647 : iw = cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
648 60 : extension=".mol", file_status='REPLACE')
649 :
650 60 : IF (iw > 0) THEN
651 30 : CPASSERT(MOD(SIZE(eigen_vec, 1), 3) == 0)
652 30 : CPASSERT(SIZE(freq, 1) == SIZE(eigen_vec, 2))
653 90 : 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 30 : IF (PRESENT(list)) THEN
657 16 : my_list(:) = 0
658 60 : DO i = 1, SIZE(list)
659 60 : my_list(list(i)) = i
660 : END DO
661 : ELSE
662 14 : CPASSERT(SIZE(particles) == SIZE(eigen_vec, 1)/3)
663 534 : DO i = 1, SIZE(my_list)
664 534 : my_list(i) = i
665 : END DO
666 : END IF
667 30 : WRITE (iw, '(T2,A)') "[Molden Format]"
668 : CALL section_vals_val_get(input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB%WRITE_CELL", &
669 30 : l_val=write_cell)
670 30 : IF (write_cell) THEN
671 0 : CPASSERT(PRESENT(cell))
672 0 : CPASSERT(ASSOCIATED(cell))
673 0 : CALL write_cell_molden(iw, cell, 1)
674 : END IF
675 30 : WRITE (iw, '(T2,A)') "[Atoms] AU"
676 597 : DO i = 1, SIZE(particles)
677 : CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
678 567 : element_symbol=element_symbol)
679 567 : CALL get_ptable_info(element_symbol, number=z)
680 :
681 : WRITE (iw, '(T2,A2,I8,I8,3X,3(F12.6,3X))') &
682 2298 : element_symbol, i, z, particles(i)%r(:)
683 :
684 : END DO
685 30 : WRITE (iw, '(T2,A)') "[FREQ]"
686 183 : DO i = 1, SIZE(freq, 1)
687 183 : IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) WRITE (iw, '(T5,F12.6)') freq(i)
688 : END DO
689 30 : WRITE (iw, '(T2,A)') "[FR-COORD]"
690 597 : DO i = 1, SIZE(particles)
691 : CALL get_atomic_kind(atomic_kind=particles(i)%atomic_kind, &
692 567 : element_symbol=element_symbol)
693 : WRITE (iw, '(T2,A2,3X,3(F12.6,3X))') &
694 2298 : element_symbol, particles(i)%r(:)
695 : END DO
696 30 : WRITE (iw, '(T2,A)') "[FR-NORM-COORD]"
697 30 : l = 0
698 183 : DO i = 1, SIZE(eigen_vec, 2)
699 183 : IF ((.NOT. dump_only_positive) .OR. (freq(i) >= 0._dp)) THEN
700 153 : l = l + 1
701 153 : WRITE (iw, '(T2,A,1X,I6)') "vibration", l
702 4758 : DO j = 1, SIZE(particles)
703 4758 : IF (my_list(j) /= 0) THEN
704 4593 : k = (my_list(j) - 1)*3
705 4593 : 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 12 : 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 30 : IF (calc_intens) THEN
713 19 : fint = massunit
714 : ! intensity units are a.u./amu
715 19 : WRITE (iw, '(T2,A)') "[INT]"
716 128 : DO i = 1, SIZE(intensities)
717 128 : 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 60 : DEALLOCATE (my_list)
721 : END IF
722 60 : CALL cp_print_key_finished_output(iw, logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
723 :
724 60 : CALL timestop(handle)
725 :
726 60 : END SUBROUTINE write_vibrations_molden
727 :
728 : END MODULE molden_utils
|
