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