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