Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Write wfx file, works as interface to chargemol and multiwfn
10 : !> \note Works only for HF and KS-type wavefunctions, execution is aborted otherwise
11 : !> \par History
12 : !> 06.2020 created [Hossam Elgabarty]
13 : !> 01.2021 modified [Hossam]
14 : !> \author Hossam Elgabarty
15 : ! **************************************************************************************************
16 : MODULE qs_chargemol
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind,&
19 : get_atomic_kind_set
20 : USE basis_set_types, ONLY: get_gto_basis_set,&
21 : gto_basis_set_type
22 : USE cell_types, ONLY: cell_type,&
23 : get_cell
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
26 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
27 : USE cp_files, ONLY: open_file
28 : USE cp_fm_types, ONLY: cp_fm_get_info,&
29 : cp_fm_get_submatrix
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_get_default_io_unit,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_print_key_generate_filename
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type,&
36 : section_vals_val_get
37 : USE kinds, ONLY: default_path_length,&
38 : dp
39 : USE kpoint_types, ONLY: kpoint_type
40 : USE machine, ONLY: m_timestamp,&
41 : timestamp_length
42 : USE orbital_pointers, ONLY: nco,&
43 : nso
44 : USE orbital_symbols, ONLY: cgf_symbol,&
45 : sgf_symbol
46 : USE orbital_transformation_matrices, ONLY: orbtramat
47 : USE particle_types, ONLY: particle_type
48 : USE periodic_table, ONLY: get_ptable_info
49 : USE qs_energy_types, ONLY: qs_energy_type
50 : USE qs_environment_types, ONLY: get_qs_env,&
51 : qs_environment_type
52 : USE qs_force_types, ONLY: qs_force_type
53 : USE qs_kind_types, ONLY: get_qs_kind,&
54 : get_qs_kind_set,&
55 : qs_kind_type
56 : USE qs_mo_types, ONLY: mo_set_type
57 : USE qs_scf_types, ONLY: qs_scf_env_type
58 : USE qs_subsys_types, ONLY: qs_subsys_type
59 : USE scf_control_types, ONLY: scf_control_type
60 : #include "./base/base_uses.f90"
61 :
62 : IMPLICIT NONE
63 : PRIVATE
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_chargemol'
66 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
67 :
68 : INTEGER, PARAMETER :: chargemol_lmax = 5
69 :
70 : PUBLIC :: write_wfx
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param qs_env ...
77 : !> \param dft_section ...
78 : ! **************************************************************************************************
79 2 : SUBROUTINE write_wfx(qs_env, dft_section)
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(section_vals_type), POINTER :: dft_section
82 :
83 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_wfx'
84 :
85 2 : CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol
86 : CHARACTER(len=2) :: asym
87 2 : CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols
88 2 : CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol
89 : CHARACTER(LEN=default_path_length) :: filename
90 : CHARACTER(LEN=timestamp_length) :: timestamp
91 : INTEGER :: handle, i, iatom, icgf, ico, ikind, iloop, imo, ipgf, irow, iset, isgf, ishell, &
92 : ispin, iwfx, lshell, max_l, ncgf, ncol_global, nelec_alpha, nelec_beta, nkpoint, noccup, &
93 : nrow_global, nset, nsgf, nspins, num_atoms, output_unit, roks_high, roks_low, tot_npgf
94 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomic_numbers, kind_of
95 2 : INTEGER, DIMENSION(:), POINTER :: lmax, npgf, nshell
96 2 : INTEGER, DIMENSION(:, :), POINTER :: l
97 : LOGICAL :: do_kpoints, newl, periodic, unrestricted
98 : REAL(KIND=dp) :: zeta
99 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atoms_cart, cmatrix, smatrix
100 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: core_charges
101 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
102 2 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
103 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
104 : TYPE(cell_type), POINTER :: cell
105 : TYPE(cp_logger_type), POINTER :: logger
106 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
107 : TYPE(dft_control_type), POINTER :: dft_control
108 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
109 : TYPE(kpoint_type), POINTER :: kpoint_env
110 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
111 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
112 : TYPE(qs_energy_type), POINTER :: energy
113 2 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
114 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
115 : TYPE(qs_scf_env_type), POINTER :: scf_env
116 : TYPE(qs_subsys_type), POINTER :: subsys
117 : TYPE(scf_control_type), POINTER :: scf_control
118 : TYPE(section_vals_type), POINTER :: chargemol_section
119 :
120 : ! !--------------------------------------------------------------------------------------------!
121 :
122 2 : CALL timeset(routineN, handle)
123 :
124 2 : logger => cp_get_default_logger()
125 2 : output_unit = cp_logger_get_default_io_unit(logger)
126 :
127 2 : IF (output_unit > 0) THEN
128 : WRITE (output_unit, '(/,T2,A)') &
129 1 : '!-----------------------------------------------------------------------------!'
130 1 : WRITE (output_unit, '(T32,A)') "Writing Chargemol wfx file..."
131 : WRITE (output_unit, '(T2,A)') &
132 1 : '!-----------------------------------------------------------------------------!'
133 : END IF
134 :
135 : ! Collect environment info
136 2 : CPASSERT(ASSOCIATED(qs_env))
137 : CALL get_qs_env(qs_env, cell=cell, mos=mos, particle_set=particle_set, &
138 : qs_kind_set=qs_kind_set, scf_env=scf_env, scf_control=scf_control, &
139 : dft_control=dft_control, energy=energy, force=force, subsys=subsys, &
140 : atomic_kind_set=atomic_kind_set, do_kpoints=do_kpoints, &
141 2 : kpoints=kpoint_env, matrix_s=matrix_s)
142 :
143 2 : nkpoint = kpoint_env%nkp
144 :
145 2 : IF (scf_control%use_ot) THEN
146 0 : CPABORT("OT is incompatible with .wfx output. Switch off OT.")
147 : END IF
148 :
149 2 : IF (dft_control%roks) THEN
150 : CALL cp_abort(__LOCATION__, "The output of chargemol .wfx files is currently "// &
151 0 : "implemented only for (un)restricted HF and Kohn-Sham-like methods.")
152 : END IF
153 :
154 2 : IF (dft_control%lsd) THEN
155 : unrestricted = .TRUE.
156 : ELSE
157 2 : unrestricted = .FALSE.
158 : END IF
159 :
160 2 : chargemol_section => section_vals_get_subs_vals(dft_section, "PRINT%CHARGEMOL")
161 :
162 2 : CALL section_vals_val_get(chargemol_section, "PERIODIC", l_val=periodic)
163 :
164 : !!---------------------------------------------------------------------------------!
165 : ! output unit
166 : !!---------------------------------------------------------------------------------!
167 :
168 : filename = cp_print_key_generate_filename(logger, chargemol_section, &
169 2 : extension=".wfx", my_local=.FALSE.)
170 :
171 : CALL open_file(filename, file_status="UNKNOWN", &
172 : file_action="WRITE", &
173 2 : unit_number=iwfx)
174 :
175 : !!---------------------------------------------------------------------------------!
176 :
177 2 : IF (iwfx > 0) THEN
178 :
179 2 : nspins = SIZE(mos)
180 2 : IF (nspins == 1) THEN
181 2 : nelec_alpha = mos(1)%nelectron/2
182 2 : nelec_beta = nelec_alpha
183 : ELSE
184 0 : nelec_alpha = mos(1)%nelectron
185 0 : nelec_beta = mos(2)%nelectron
186 : END IF
187 :
188 2 : IF (dft_control%roks) THEN
189 0 : IF (mos(1)%homo > mos(2)%homo) THEN
190 : roks_high = 1
191 : roks_low = 2
192 : ELSE
193 0 : roks_high = 2
194 0 : roks_low = 1
195 : END IF
196 : END IF
197 :
198 : !!---------------------------------------------------------------------------------!
199 : ! Initial parsing of topology and basis set, check maximum l .le. chargemol_lmax
200 : !!---------------------------------------------------------------------------------!
201 2 : num_atoms = SIZE(particle_set)
202 6 : ALLOCATE (atoms_cart(3, num_atoms))
203 6 : ALLOCATE (atom_symbols(num_atoms))
204 6 : ALLOCATE (atomic_numbers(num_atoms))
205 6 : ALLOCATE (core_charges(num_atoms))
206 :
207 2 : max_l = 0
208 4 : DO iatom = 1, num_atoms
209 2 : NULLIFY (orb_basis_set)
210 8 : atoms_cart(1:3, iatom) = particle_set(iatom)%r(1:3)
211 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
212 2 : element_symbol=asym, kind_number=ikind)
213 2 : atom_symbols(iatom) = asym
214 2 : CALL get_ptable_info(asym, number=atomic_numbers(iatom))
215 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
216 2 : core_charge=core_charges(iatom))
217 2 : IF (ASSOCIATED(orb_basis_set)) THEN
218 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
219 2 : nset=nset, lmax=lmax)
220 :
221 6 : DO iset = 1, nset
222 6 : max_l = MAX(max_l, lmax(iset))
223 : END DO
224 : ELSE
225 0 : CPABORT("Unknown basis set type. ")
226 : END IF
227 6 : IF (max_l > chargemol_lmax) THEN
228 : CALL cp_abort(__LOCATION__, "Chargemol supports basis sets with l"// &
229 0 : " up to 5 only (h angular functions).")
230 : END IF
231 : END DO
232 :
233 : ! !===========================================================================!
234 : ! Header comment and Title
235 : ! !===========================================================================!
236 2 : WRITE (iwfx, "(A)") "# Chargemol input file generated by CP2K "
237 2 : CALL m_timestamp(timestamp)
238 2 : WRITE (iwfx, "(A,/)") "# Creation date "//timestamp(:19)
239 :
240 2 : WRITE (iwfx, "(A)") "<Title>"
241 2 : WRITE (iwfx, *) TRIM(logger%iter_info%project_name)
242 2 : WRITE (iwfx, "(A,/)") "</Title>"
243 :
244 : ! !===========================================================================!
245 : ! Keywords
246 : ! TODO: check: always GTO??
247 : ! !===========================================================================!
248 2 : WRITE (iwfx, "(A)") "<Keywords>"
249 2 : WRITE (iwfx, "(A)") "GTO"
250 2 : WRITE (iwfx, "(A,/)") "</Keywords>"
251 :
252 : ! !===========================================================================!
253 : ! Model
254 : ! !===========================================================================!
255 2 : WRITE (iwfx, "(A)") "#<Model>"
256 2 : IF (unrestricted) THEN
257 0 : WRITE (iwfx, "(A)") "#Unrestricted Kohn-Sham"
258 : ELSE
259 2 : WRITE (iwfx, "(A)") "#Restricted Kohn-Sham"
260 : END IF
261 2 : WRITE (iwfx, "(A,/)") "#</Model>"
262 :
263 : ! !===========================================================================!
264 : ! Number of nuclei
265 : ! !===========================================================================!
266 2 : WRITE (iwfx, "(A)") "<Number of Nuclei>"
267 2 : WRITE (iwfx, "(I6)") num_atoms
268 2 : WRITE (iwfx, "(A,/)") "</Number of Nuclei>"
269 :
270 : ! !===========================================================================!
271 : ! Number of Occupied MOs
272 : ! !===========================================================================!
273 2 : WRITE (iwfx, "(A)") "<Number of Occupied Molecular Orbitals>"
274 2 : noccup = 0
275 2 : IF (dft_control%roks .AND. nspins == 2) THEN
276 0 : noccup = MAX(mos(1)%homo, mos(2)%homo)
277 : ELSE
278 4 : DO ispin = 1, dft_control%nspins
279 4 : noccup = noccup + mos(ispin)%homo
280 : END DO
281 : END IF
282 2 : WRITE (iwfx, "(2I6)") noccup
283 2 : WRITE (iwfx, "(A,/)") "</Number of Occupied Molecular Orbitals>"
284 :
285 : ! !===========================================================================!
286 : ! Number of Perturbations
287 : ! !===========================================================================!
288 2 : WRITE (iwfx, "(A)") "<Number of Perturbations>"
289 2 : WRITE (iwfx, "(I6)") 0
290 2 : WRITE (iwfx, "(A,/)") "</Number of Perturbations>"
291 :
292 : ! !===========================================================================!
293 : ! Total charge
294 : ! !===========================================================================!
295 2 : WRITE (iwfx, "(A)") "<Net Charge>"
296 2 : WRITE (iwfx, "(F8.4)") REAL(dft_control%charge, KIND=dp)
297 2 : WRITE (iwfx, "(A,/)") "</Net Charge>"
298 :
299 : ! !===========================================================================!
300 : ! Number of electrons
301 : ! !===========================================================================!
302 2 : WRITE (iwfx, "(A)") "<Number of Electrons>"
303 2 : WRITE (iwfx, "(I6)") scf_env%nelectron
304 2 : WRITE (iwfx, "(A,/)") "</Number of Electrons>"
305 :
306 : !===========================================================================!
307 : ! Number of Alpha electrons
308 : !===========================================================================!
309 2 : WRITE (iwfx, "(A)") "<Number of Alpha Electrons>"
310 2 : WRITE (iwfx, "(I6)") nelec_alpha
311 2 : WRITE (iwfx, "(A,/)") "</Number of Alpha Electrons>"
312 :
313 : !===========================================================================!
314 : ! Number of Beta electrons
315 : !===========================================================================!
316 2 : WRITE (iwfx, "(A)") "<Number of Beta Electrons>"
317 2 : WRITE (iwfx, "(I6)") nelec_beta
318 2 : WRITE (iwfx, "(A,/)") "</Number of Beta Electrons>"
319 :
320 : !===========================================================================!
321 : ! Spin multiplicity
322 : !===========================================================================!
323 2 : WRITE (iwfx, "(A)") "<Electron Spin Multiplicity>"
324 2 : WRITE (iwfx, "(I4)") dft_control%multiplicity
325 2 : WRITE (iwfx, "(A,/)") "</Electron Spin Multiplicity>"
326 :
327 : ! !===========================================================================!
328 : ! Number of Core Electrons
329 : ! !===========================================================================!
330 2 : WRITE (iwfx, "(A)") "<Number of Core Electrons>"
331 6 : WRITE (iwfx, "(F16.10)") SUM(atomic_numbers) - SUM(core_charges)
332 2 : WRITE (iwfx, "(A,/)") "</Number of Core Electrons>"
333 :
334 : ! !===========================================================================!
335 : ! Nuclear Names
336 : ! !===========================================================================!
337 2 : WRITE (iwfx, "(A)") "<Nuclear Names>"
338 4 : DO iatom = 1, num_atoms
339 4 : WRITE (iwfx, "(A4)") atom_symbols(iatom)
340 : END DO
341 2 : WRITE (iwfx, "(A,/)") "</Nuclear Names>"
342 :
343 : ! !===========================================================================!
344 : ! Atomic Numbers
345 : ! !===========================================================================!
346 2 : WRITE (iwfx, "(A)") "<Atomic Numbers>"
347 4 : DO iatom = 1, num_atoms
348 4 : WRITE (iwfx, "(I4)") atomic_numbers(iatom)
349 : END DO
350 2 : WRITE (iwfx, "(A,/)") "</Atomic Numbers>"
351 :
352 : ! !===========================================================================!
353 : ! Nuclear charges
354 : ! !===========================================================================!
355 2 : WRITE (iwfx, "(A)") "<Nuclear Charges>"
356 4 : DO iatom = 1, num_atoms
357 4 : WRITE (iwfx, "(F12.8)") core_charges(iatom)
358 : END DO
359 2 : WRITE (iwfx, "(A,/)") "</Nuclear Charges>"
360 :
361 : ! !===========================================================================!
362 : ! Nuclear coordinates
363 : ! !===========================================================================!
364 2 : WRITE (iwfx, "(A)") "<Nuclear Cartesian Coordinates>"
365 4 : DO iatom = 1, num_atoms
366 4 : WRITE (iwfx, "(3ES26.16E3)") atoms_cart(1:3, iatom)
367 : END DO
368 2 : WRITE (iwfx, "(A,/)") "</Nuclear Cartesian Coordinates>"
369 :
370 : ! !===========================================================================!
371 : ! Periodic boundary conditions
372 : ! !===========================================================================!
373 2 : IF (periodic) THEN
374 2 : CALL get_cell(cell)
375 8 : IF (SUM(cell%perd) == 0) THEN
376 : CALL cp_warn(__LOCATION__, "Writing of periodic cell information"// &
377 0 : " requested in input, but system is not periodic, ignoring request.")
378 : ELSE
379 2 : WRITE (iwfx, "(A)") "<Number of Translation Vectors>"
380 8 : WRITE (iwfx, "(I3)") SUM(cell%perd)
381 2 : WRITE (iwfx, "(A,/)") "</Number of Translation Vectors>"
382 :
383 2 : WRITE (iwfx, "(A)") "<Translation Vectors>"
384 8 : DO iatom = 1, 3
385 8 : IF (cell%perd(iatom) == 1) THEN
386 6 : WRITE (iwfx, "(3F12.6)") cell%hmat(1:3, iatom)
387 : END IF
388 : END DO
389 2 : WRITE (iwfx, "(A,/)") "</Translation Vectors>"
390 : END IF
391 2 : WRITE (iwfx, "(A)") "<Number of Kpoints>"
392 2 : WRITE (iwfx, "(I3)") 1
393 2 : WRITE (iwfx, "(A,/)") "</Number of Kpoints>"
394 2 : WRITE (iwfx, "(A)") "<Kpoint Weights>"
395 2 : WRITE (iwfx, "(I3)") 1
396 2 : WRITE (iwfx, "(A,/)") "</Kpoint Weights>"
397 2 : WRITE (iwfx, "(A)") "<Kpoint Fractional Coordinates>"
398 2 : WRITE (iwfx, "(F8.6)") 0.0
399 2 : WRITE (iwfx, "(A,/)") "</Kpoint Fractional Coordinates>"
400 : END IF
401 :
402 : ! !===========================================================================!
403 : ! Primitive centers
404 : ! !===========================================================================!
405 2 : WRITE (iwfx, "(A)") "<Primitive Centers>"
406 2 : tot_npgf = 0
407 4 : DO iatom = 1, num_atoms
408 2 : iloop = 1
409 2 : NULLIFY (orb_basis_set)
410 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
411 2 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
412 6 : IF (ASSOCIATED(orb_basis_set)) THEN
413 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
414 : nset=nset, &
415 : npgf=npgf, &
416 : nshell=nshell, &
417 : l=l, &
418 2 : gcc=gcc)
419 :
420 6 : DO iset = 1, nset
421 16 : DO ishell = 1, nshell(iset)
422 10 : lshell = l(ishell, iset)
423 42 : DO ico = 1, nco(lshell)
424 114 : DO ipgf = 1, npgf(iset)
425 76 : tot_npgf = tot_npgf + 1
426 76 : IF (MOD(iloop + 10, 10) /= 0) THEN
427 70 : WRITE (iwfx, "(I4)", advance="no") iatom
428 70 : newl = .TRUE.
429 : ELSE
430 6 : WRITE (iwfx, "(I4)") iatom
431 6 : newl = .FALSE.
432 : END IF
433 104 : iloop = iloop + 1
434 : !END IF
435 : END DO
436 : END DO
437 : END DO
438 : END DO
439 2 : IF (newl) THEN
440 2 : WRITE (iwfx, *)
441 : END IF
442 : ELSE
443 0 : CPABORT("Unknown basis set type. ")
444 : END IF
445 : END DO
446 2 : WRITE (iwfx, "(A,/)") "</Primitive Centers>"
447 :
448 : ! !===========================================================================!
449 : ! Number of primitives
450 : ! !===========================================================================!
451 2 : WRITE (iwfx, "(A)") "<Number of Primitives>"
452 2 : WRITE (iwfx, "(I8)") tot_npgf
453 2 : WRITE (iwfx, "(A,/)") "</Number of Primitives>"
454 :
455 : ! !===========================================================================!
456 : ! Primitive Types
457 : ! !===========================================================================!
458 2 : WRITE (iwfx, "(A)") "<Primitive Types>"
459 4 : DO iatom = 1, num_atoms
460 2 : iloop = 1
461 2 : NULLIFY (orb_basis_set)
462 2 : NULLIFY (bcgf_symbol)
463 2 : NULLIFY (gcc)
464 2 : NULLIFY (zet)
465 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
466 2 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
467 2 : IF (ASSOCIATED(orb_basis_set)) THEN
468 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
469 : nset=nset, &
470 : nshell=nshell, &
471 : l=l, &
472 : npgf=npgf, &
473 : ncgf=ncgf, &
474 : zet=zet, &
475 : cgf_symbol=bcgf_symbol, &
476 : sgf_symbol=bsgf_symbol, &
477 2 : gcc=gcc)
478 :
479 2 : icgf = 1
480 6 : DO iset = 1, nset
481 16 : DO ishell = 1, nshell(iset)
482 10 : lshell = l(ishell, iset)
483 42 : DO ico = 1, nco(lshell)
484 104 : DO ipgf = 1, orb_basis_set%npgf(iset)
485 76 : IF (MOD(iloop + 10, 10) /= 0) THEN
486 : WRITE (iwfx, '(I6)', advance="no") &
487 70 : pgf_type_chargemol(bcgf_symbol(icgf) (3:))
488 70 : newl = .TRUE.
489 : ELSE
490 : WRITE (iwfx, '(I6)') &
491 6 : pgf_type_chargemol(bcgf_symbol(icgf) (3:))
492 6 : newl = .FALSE.
493 : END IF
494 104 : iloop = iloop + 1
495 : !END IF
496 : END DO
497 38 : icgf = icgf + 1
498 : END DO !ico
499 : END DO ! ishell
500 : END DO ! iset
501 :
502 : ELSE
503 0 : CPABORT("Unknown basis set type.")
504 : END IF
505 6 : IF (newl) THEN
506 2 : WRITE (iwfx, *)
507 : END IF
508 : END DO ! iatom
509 2 : WRITE (iwfx, "(A,/)") "</Primitive Types>"
510 :
511 : ! !===========================================================================!
512 : ! Primitive Exponents
513 : ! !===========================================================================!
514 : !! NOTE: CP2K orders the atomic orbitals as follows:
515 : !! Cartesian orbitals: x (slowest loop), y, z (fastest loop)
516 : !! Spherical orbitals: m = -l, ..., 0, ..., +l
517 :
518 2 : WRITE (iwfx, "(A)") "<Primitive Exponents>"
519 4 : DO iatom = 1, num_atoms
520 2 : NULLIFY (orb_basis_set)
521 2 : NULLIFY (gcc)
522 2 : NULLIFY (zet)
523 2 : iloop = 1
524 :
525 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
526 2 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
527 6 : IF (ASSOCIATED(orb_basis_set)) THEN
528 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
529 : nset=nset, &
530 : npgf=npgf, &
531 : nshell=nshell, &
532 : l=l, &
533 : zet=zet, &
534 2 : gcc=gcc)
535 :
536 6 : DO iset = 1, nset
537 16 : DO ishell = 1, nshell(iset)
538 10 : lshell = l(ishell, iset)
539 42 : DO ico = 1, nco(lshell)
540 114 : DO ipgf = 1, npgf(iset)
541 76 : IF (MOD(iloop + 5, 5) /= 0) THEN
542 62 : WRITE (iwfx, "(ES20.10)", advance="no") zet(ipgf, iset)
543 62 : newl = .TRUE.
544 : ELSE
545 14 : WRITE (iwfx, "(ES20.10)") zet(ipgf, iset)
546 14 : newl = .FALSE.
547 : END IF
548 104 : iloop = iloop + 1
549 : !END IF
550 : END DO
551 : END DO
552 : END DO
553 : END DO
554 2 : IF (newl) THEN
555 2 : WRITE (iwfx, *)
556 : END IF
557 : ELSE
558 0 : CPABORT("Unknown basis set type. ")
559 : END IF
560 : END DO
561 2 : WRITE (iwfx, "(A,/)") "</Primitive Exponents>"
562 :
563 : ! !===========================================================================!
564 : ! MO occupation numbers
565 : ! nonesense for roks at the moment
566 : ! !===========================================================================!
567 2 : WRITE (iwfx, "(A)") "<Molecular Orbital Occupation Numbers>"
568 2 : IF (.NOT. dft_control%roks) THEN
569 4 : DO ispin = 1, nspins
570 12 : DO imo = 1, mos(ispin)%homo
571 : WRITE (iwfx, "(f8.6)") &
572 10 : mos(ispin)%occupation_numbers(imo)
573 : END DO
574 : END DO
575 : ELSE
576 0 : DO imo = 1, mos(roks_low)%homo
577 : WRITE (iwfx, "(*(f8.6,1x))") &
578 : mos(roks_low)%occupation_numbers(imo) &
579 0 : + mos(roks_low)%occupation_numbers(imo)
580 : END DO
581 0 : DO imo = mos(roks_low)%homo + 1, mos(roks_high)%homo
582 : WRITE (iwfx, "(f8.6)") &
583 0 : mos(roks_high)%occupation_numbers(imo)
584 : END DO
585 : END IF
586 2 : WRITE (iwfx, "(A,/)") "</Molecular Orbital Occupation Numbers>"
587 :
588 : ! !===========================================================================!
589 : ! MO energies
590 : ! !===========================================================================!
591 2 : WRITE (iwfx, "(A)") "<Molecular Orbital Energies>"
592 4 : DO ispin = 1, nspins
593 12 : DO imo = 1, mos(ispin)%homo
594 10 : WRITE (iwfx, "(ES20.10)") mos(ispin)%eigenvalues(imo)
595 : END DO
596 : END DO
597 2 : WRITE (iwfx, "(A,/)") "</Molecular Orbital Energies>"
598 :
599 : ! !===========================================================================!
600 : ! MO Spin types
601 : ! nonesense for ROKS
602 : ! !===========================================================================!
603 2 : WRITE (iwfx, "(A)") "<Molecular Orbital Spin Types>"
604 10 : DO imo = 1, mos(1)%homo
605 10 : IF (nspins == 1) THEN
606 8 : WRITE (iwfx, '(A15)') "Alpha and Beta"
607 0 : ELSEIF (dft_control%uks) THEN
608 0 : WRITE (iwfx, '(A6)') "Alpha"
609 : ELSE
610 : CALL cp_abort(__LOCATION__, "This wavefunction type is currently"// &
611 0 : " not supported by the chargemol interface.")
612 : END IF
613 : END DO
614 2 : IF (nspins == 2) THEN
615 0 : DO imo = 1, mos(2)%homo
616 0 : WRITE (iwfx, '(A5)') "Beta"
617 : END DO
618 : END IF
619 2 : WRITE (iwfx, "(A,/)") "</Molecular Orbital Spin Types>"
620 :
621 : ! !===========================================================================!
622 : ! Kpoint index of orbitals
623 : ! !===========================================================================!
624 2 : IF (periodic) THEN
625 2 : WRITE (iwfx, "(A)") "<Kpoint Index for Orbitals>"
626 4 : DO ispin = 1, nspins
627 12 : DO imo = 1, mos(ispin)%homo
628 10 : WRITE (iwfx, "(I6)") 1
629 : END DO
630 : END DO
631 2 : WRITE (iwfx, "(A,/)") "</Kpoint Index for Orbitals>"
632 : END IF
633 :
634 : ! !===========================================================================!
635 : ! MO Primitive Coefficients
636 : ! !===========================================================================!
637 2 : WRITE (iwfx, "(A)") "<Molecular Orbital Primitive Coefficients>"
638 2 : NULLIFY (orb_basis_set)
639 2 : NULLIFY (gcc)
640 :
641 2 : IF (mos(1)%use_mo_coeff_b) THEN
642 0 : DO ispin = 1, SIZE(mos)
643 0 : IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
644 0 : CPASSERT(.FALSE.)
645 : END IF
646 : CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
647 0 : mos(ispin)%mo_coeff)
648 : END DO
649 : END IF
650 :
651 4 : DO ispin = 1, nspins
652 : CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
653 : nrow_global=nrow_global, &
654 2 : ncol_global=ncol_global)
655 :
656 8 : ALLOCATE (smatrix(nrow_global, ncol_global))
657 114 : smatrix = 0.0_dp
658 :
659 2 : CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
660 :
661 2 : CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
662 :
663 8 : ALLOCATE (cmatrix(ncgf, ncol_global))
664 :
665 122 : cmatrix = 0.0_dp
666 :
667 : ! get MO coefficients in terms of contracted cartesian functions
668 2 : icgf = 1
669 2 : isgf = 1
670 4 : DO iatom = 1, num_atoms
671 2 : NULLIFY (orb_basis_set)
672 2 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
673 : CALL get_qs_kind(qs_kind_set(ikind), &
674 2 : basis_set=orb_basis_set)
675 6 : IF (ASSOCIATED(orb_basis_set)) THEN
676 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
677 : nset=nset, &
678 : nshell=nshell, &
679 2 : l=l)
680 6 : DO iset = 1, nset
681 16 : DO ishell = 1, nshell(iset)
682 10 : lshell = l(ishell, iset)
683 : CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, &
684 : nso(lshell), 1.0_dp, &
685 : orbtramat(lshell)%c2s, nso(lshell), &
686 : smatrix(isgf, 1), nsgf, 0.0_dp, &
687 10 : cmatrix(icgf, 1), ncgf)
688 : !IF (iatom==1 .and. debug_this_module) THEN
689 : ! print *, lshell
690 : ! print *, "size ", size(orbtramat(lshell)%s2c,1),",",size(orbtramat(lshell)%s2c,2)
691 : ! DO itest = 1, size(orbtramat(lshell)%s2c,1)
692 : ! print *, orbtramat(lshell)%s2c(itest,:)
693 : ! END DO
694 : ! print *, " "
695 : ! DO itest = 1, 5
696 : ! print *, my_d_orbtramat(itest,:)
697 : ! END DO
698 : !END IF
699 10 : icgf = icgf + nco(lshell)
700 14 : isgf = isgf + nso(lshell)
701 : END DO
702 : END DO
703 : ELSE
704 0 : CPABORT("Unknown basis set type. ")
705 : END IF
706 : END DO ! iatom
707 :
708 : ! Now write MO coefficients in terms of cartesian primitives
709 10 : DO imo = 1, mos(ispin)%homo
710 :
711 8 : WRITE (iwfx, "(A)") "<MO Number>"
712 8 : IF (nspins == 1 .OR. ispin == 1) THEN
713 8 : WRITE (iwfx, "(I6)") imo
714 : ELSE
715 0 : WRITE (iwfx, "(I6)") imo + mos(1)%homo
716 : END IF
717 8 : WRITE (iwfx, "(A)") "</MO Number>"
718 :
719 8 : irow = 1 ! row of cmatrix, each row corresponds to
720 : ! a contracted Cartesian function
721 18 : DO iatom = 1, num_atoms
722 8 : iloop = 1
723 8 : NULLIFY (orb_basis_set)
724 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
725 8 : kind_number=ikind)
726 : CALL get_qs_kind(qs_kind_set(ikind), &
727 8 : basis_set=orb_basis_set)
728 8 : IF (ASSOCIATED(orb_basis_set)) THEN
729 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
730 : nset=nset, &
731 : nshell=nshell, &
732 : gcc=gcc, &
733 8 : l=l)
734 :
735 8 : icgf = 1
736 24 : DO iset = 1, nset
737 64 : DO ishell = 1, nshell(iset)
738 40 : lshell = orb_basis_set%l(ishell, iset)
739 :
740 152 : DO ico = orb_basis_set%first_cgf(ishell, iset), &
741 56 : orb_basis_set%last_cgf(ishell, iset)
742 :
743 : !loop over primitives
744 416 : DO ipgf = 1, orb_basis_set%npgf(iset)
745 :
746 304 : zeta = orb_basis_set%zet(ipgf, iset)
747 :
748 304 : IF (MOD(iloop + 5, 5) /= 0) THEN
749 : WRITE (iwfx, '(ES20.10)', advance="no") &
750 : cmatrix(irow, imo)* &
751 : orb_basis_set%norm_cgf(ico)* &
752 248 : orb_basis_set%gcc(ipgf, ishell, iset)
753 248 : newl = .TRUE.
754 : ELSE
755 : WRITE (iwfx, '(ES20.10)') &
756 : cmatrix(irow, imo)* &
757 : orb_basis_set%norm_cgf(ico)* &
758 56 : orb_basis_set%gcc(ipgf, ishell, iset)
759 56 : newl = .FALSE.
760 : END IF
761 416 : iloop = iloop + 1
762 : END DO ! end loop over primitives
763 152 : irow = irow + 1
764 : END DO !ico
765 : END DO ! ishell
766 : END DO ! iset
767 :
768 : ELSE
769 0 : CPABORT("Unknown basis set type.")
770 : END IF
771 24 : IF (newl) THEN
772 8 : WRITE (iwfx, *)
773 : END IF
774 : END DO ! iatom
775 : !Write (iwfx,*)
776 : END DO ! imo
777 :
778 2 : IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
779 8 : IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
780 :
781 : END DO ! ispin
782 2 : WRITE (iwfx, "(A,/)") "</Molecular Orbital Primitive Coefficients>"
783 :
784 : ! !===========================================================================!
785 : ! Energy
786 : ! !===========================================================================!
787 2 : WRITE (iwfx, "(A)") "<Energy>"
788 2 : WRITE (iwfx, "(ES26.16E3)") energy%total
789 2 : WRITE (iwfx, "(A,/)") "</Energy>"
790 :
791 : ! !===========================================================================!
792 : ! Virial ratio
793 : ! !===========================================================================!
794 2 : WRITE (iwfx, "(A)") "<Virial Ratio (-V/T)>"
795 2 : WRITE (iwfx, "(ES20.12)") - 1*(energy%total - energy%kinetic)/energy%kinetic
796 2 : WRITE (iwfx, "(A,/)") "</Virial Ratio (-V/T)>"
797 :
798 : ! !===========================================================================!
799 : ! Force-related quantities
800 : ! inactivated for now
801 : ! !===========================================================================!
802 : IF (ASSOCIATED(force) .AND. debug_this_module) THEN
803 : WRITE (iwfx, "(A)") "<Nuclear Cartesian Energy Gradients>"
804 :
805 : CALL get_qs_env(qs_env, force=force, atomic_kind_set=atomic_kind_set)
806 : ALLOCATE (atom_of_kind(num_atoms))
807 : ALLOCATE (kind_of(num_atoms))
808 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
809 : atom_of_kind=atom_of_kind, &
810 : kind_of=kind_of)
811 :
812 : DO iatom = 1, num_atoms
813 : ikind = kind_of(iatom)
814 : i = atom_of_kind(iatom)
815 : WRITE (iwfx, "(3ES20.12E3)") force(ikind)%total(1:3, i)
816 : END DO
817 :
818 : WRITE (iwfx, "(A,/)") "<Nuclear Cartesian Energy Gradients>"
819 :
820 : ! W
821 : WRITE (iwfx, "(A)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
822 : WRITE (iwfx, "(A)") ""
823 : WRITE (iwfx, "(A,/)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
824 :
825 : ! Full Virial ratio
826 : WRITE (iwfx, "(A)") "<Full Virial Ratio, -(V - W)/T)>"
827 : WRITE (iwfx, "(A)") ""
828 : WRITE (iwfx, "(A,/)") "</Full Virial Ratio, -(V - W)/T)>"
829 :
830 : DEALLOCATE (atom_of_kind)
831 : DEALLOCATE (kind_of)
832 :
833 : END IF
834 :
835 : ! !===========================================================================!
836 : ! Core-density
837 : ! !===========================================================================!
838 : ! Additional Electron Density Function (EDF)
839 : ! WRITE (iwfx, "(A)") "#<Additional Electron Density Function (EDF)>"
840 : ! WRITE (iwfx, "(A,/)") "#</Additional Electron Density Function (EDF)>"
841 :
842 : ! !-----------------------------------!
843 : ! Done
844 : ! !-----------------------------------!
845 2 : DEALLOCATE (atoms_cart)
846 2 : DEALLOCATE (atom_symbols)
847 2 : DEALLOCATE (atomic_numbers)
848 2 : DEALLOCATE (core_charges)
849 :
850 : ELSE ! no output unit
851 0 : iwfx = -1
852 : END IF
853 :
854 2 : IF (output_unit > 0) THEN
855 : WRITE (output_unit, '(/,T2,A)') &
856 1 : '!--------------------------- Chargemol wfx written ---------------------------!'
857 : END IF
858 :
859 2 : CALL timestop(handle)
860 :
861 4 : END SUBROUTINE write_wfx
862 :
863 : ! **************************************************************************************************
864 : !> \brief ...
865 : !> \param l_symb ...
866 : !> \return ...
867 : ! **************************************************************************************************
868 76 : INTEGER FUNCTION pgf_type_chargemol(l_symb) RESULT(label)
869 : !INTEGER, INTENT(OUT) :: label
870 : CHARACTER(len=*), INTENT(IN) :: l_symb
871 :
872 : SELECT CASE (l_symb)
873 : CASE ("s")
874 16 : label = 1
875 : CASE ("px")
876 16 : label = 2
877 : CASE ("py")
878 16 : label = 3
879 : CASE ("pz")
880 16 : label = 4
881 : CASE ("dx2")
882 2 : label = 5
883 : CASE ("dy2")
884 2 : label = 6
885 : CASE ("dz2")
886 2 : label = 7
887 : CASE ("dxy")
888 2 : label = 8
889 : CASE ("dxz")
890 2 : label = 9
891 : CASE ("dyz")
892 2 : label = 10
893 : CASE ("fx3")
894 0 : label = 11
895 : CASE ("fy3")
896 0 : label = 12
897 : CASE ("fz3")
898 0 : label = 13
899 : CASE ("fx2y")
900 0 : label = 14
901 : CASE ("fx2z")
902 0 : label = 15
903 : CASE ("fy2z")
904 0 : label = 16
905 : CASE ("fxy2")
906 0 : label = 17
907 : CASE ("fxz2")
908 0 : label = 18
909 : CASE ("fyz2")
910 0 : label = 19
911 : CASE ("fxyz")
912 0 : label = 20
913 : CASE ("gx4")
914 0 : label = 21
915 : CASE ("gy4")
916 0 : label = 22
917 : CASE ("gz4")
918 0 : label = 23
919 : CASE ("gx3y")
920 0 : label = 24
921 : CASE ("gx3z")
922 0 : label = 25
923 : CASE ("gxy3")
924 0 : label = 26
925 : CASE ("gy3z")
926 0 : label = 27
927 : CASE ("gxz3")
928 0 : label = 28
929 : CASE ("gyz3")
930 0 : label = 29
931 : CASE ("gx2y2")
932 0 : label = 30
933 : CASE ("gx2z2")
934 0 : label = 31
935 : CASE ("gy2z2")
936 0 : label = 32
937 : CASE ("gx2yz")
938 0 : label = 33
939 : CASE ("gxy2z")
940 0 : label = 34
941 : CASE ("gxyz2")
942 0 : label = 35
943 : CASE ("hz5")
944 0 : label = 36
945 : CASE ("hyz4")
946 0 : label = 37
947 : CASE ("hy2z3")
948 0 : label = 38
949 : CASE ("hy3z2")
950 0 : label = 39
951 : CASE ("hy4z")
952 0 : label = 40
953 : CASE ("hy5")
954 0 : label = 41
955 : CASE ("hxz4")
956 0 : label = 42
957 : CASE ("hxyz3")
958 0 : label = 43
959 : CASE ("hxy2z2")
960 0 : label = 44
961 : CASE ("hxy3z")
962 0 : label = 45
963 : CASE ("hxy4")
964 0 : label = 46
965 : CASE ("hx2z3")
966 0 : label = 47
967 : CASE ("hx2yz2")
968 0 : label = 48
969 : CASE ("hx2y2z")
970 0 : label = 49
971 : CASE ("hx2y3")
972 0 : label = 50
973 : CASE ("hx3z2")
974 0 : label = 51
975 : CASE ("hx3yz")
976 0 : label = 52
977 : CASE ("hx3y2")
978 0 : label = 53
979 : CASE ("hx4z")
980 0 : label = 54
981 : CASE ("hx4y")
982 0 : label = 55
983 : CASE ("hx5")
984 0 : label = 56
985 : CASE default
986 76 : CPABORT("The chargemol interface supports basis functions up to l=5 only")
987 : !label = 99
988 : END SELECT
989 :
990 76 : END FUNCTION pgf_type_chargemol
991 :
992 : END MODULE qs_chargemol
993 :
|
