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