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