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