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 calculate the orbitals for a given atomic kind type
10 : ! **************************************************************************************************
11 : MODULE atom_kind_orbitals
12 : USE ai_onecenter, ONLY: sg_erfc
13 : USE atom_electronic_structure, ONLY: calculate_atom
14 : USE atom_fit, ONLY: atom_fit_density
15 : USE atom_operators, ONLY: atom_int_release,&
16 : atom_int_setup,&
17 : atom_ppint_release,&
18 : atom_ppint_setup,&
19 : atom_relint_release,&
20 : atom_relint_setup
21 : USE atom_set_basis, ONLY: set_kind_basis_atomic
22 : USE atom_types, ONLY: &
23 : CGTO_BASIS, Clementi_geobas, GTO_BASIS, atom_basis_type, atom_ecppot_type, &
24 : atom_gthpot_type, atom_integrals, atom_orbitals, atom_potential_type, atom_sgppot_type, &
25 : atom_type, create_atom_orbs, create_atom_type, lmat, release_atom_basis, &
26 : release_atom_potential, release_atom_type, set_atom
27 : USE atom_utils, ONLY: atom_density,&
28 : get_maxl_occ,&
29 : get_maxn_occ
30 : USE atomic_kind_types, ONLY: atomic_kind_type,&
31 : get_atomic_kind
32 : USE basis_set_types, ONLY: get_gto_basis_set,&
33 : gto_basis_set_type
34 : USE external_potential_types, ONLY: all_potential_type,&
35 : get_potential,&
36 : gth_potential_type,&
37 : sgp_potential_type
38 : USE input_constants, ONLY: &
39 : barrier_conf, do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, &
40 : do_gapw_log, do_nonrel_atom, do_numeric, do_rks_atom, do_sczoramp_atom, do_uks_atom, &
41 : do_zoramp_atom, ecp_pseudo, gth_pseudo, no_pseudo, poly_conf, rel_dkh, rel_none, &
42 : rel_sczora_mp, rel_zora, rel_zora_full, rel_zora_mp, sgp_pseudo
43 : USE input_section_types, ONLY: section_vals_type
44 : USE kinds, ONLY: dp
45 : USE mathconstants, ONLY: dfac,&
46 : pi
47 : USE periodic_table, ONLY: ptable
48 : USE physcon, ONLY: bohr
49 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
50 : create_grid_atom,&
51 : grid_atom_type
52 : USE qs_kind_types, ONLY: get_qs_kind,&
53 : init_atom_electronic_state,&
54 : qs_kind_type,&
55 : set_pseudo_state
56 : USE rel_control_types, ONLY: rel_control_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_kind_orbitals'
64 :
65 : PUBLIC :: calculate_atomic_orbitals, calculate_atomic_density, &
66 : calculate_atomic_relkin, gth_potential_conversion
67 :
68 : ! **************************************************************************************************
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief ...
74 : !> \param atomic_kind ...
75 : !> \param qs_kind ...
76 : !> \param agrid ...
77 : !> \param iunit ...
78 : !> \param pmat ...
79 : !> \param fmat ...
80 : !> \param density ...
81 : !> \param wavefunction ...
82 : !> \param wfninfo ...
83 : !> \param confine ...
84 : !> \param xc_section ...
85 : !> \param nocc ...
86 : !> \param which_l ...
87 : !> \param which_n ...
88 : !> \param proj_shell_charge ...
89 : !> \param ao_coef ...
90 : ! **************************************************************************************************
91 30708 : SUBROUTINE calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, &
92 10236 : density, wavefunction, wfninfo, confine, xc_section, nocc, &
93 : which_l, which_n, proj_shell_charge, ao_coef)
94 : TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
95 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
96 : TYPE(grid_atom_type), OPTIONAL :: agrid
97 : INTEGER, INTENT(IN), OPTIONAL :: iunit
98 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
99 : POINTER :: pmat, fmat
100 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: density
101 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: wavefunction, wfninfo
102 : LOGICAL, INTENT(IN), OPTIONAL :: confine
103 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
104 : INTEGER, DIMENSION(:), OPTIONAL :: nocc
105 : INTEGER, INTENT(IN), OPTIONAL :: which_l, which_n
106 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: proj_shell_charge(4)
107 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
108 : OPTIONAL, POINTER :: ao_coef
109 :
110 : INTEGER :: i, ii, j, k, k1, k2, l, ll, m, mb, mo, &
111 : nr, nset, nsgf, projector_index, &
112 : valence_n, z
113 : INTEGER, DIMENSION(0:lmat) :: nbb
114 : INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
115 : INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
116 10236 : INTEGER, DIMENSION(:), POINTER :: nshell
117 10236 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, ls
118 : LOGICAL :: ecp_semi_local, ghost, has_pp, uks
119 : REAL(KIND=dp) :: ok, qleft, scal, take, zeff
120 : REAL(KIND=dp), DIMENSION(0:lmat, 10) :: occupation_neutral
121 : REAL(KIND=dp), DIMENSION(0:lmat, 10, 2) :: edelta
122 : TYPE(all_potential_type), POINTER :: all_potential
123 : TYPE(atom_basis_type), POINTER :: basis
124 : TYPE(atom_integrals), POINTER :: integrals
125 : TYPE(atom_orbitals), POINTER :: orbitals
126 : TYPE(atom_potential_type), POINTER :: potential
127 : TYPE(atom_type), POINTER :: atom
128 : TYPE(gth_potential_type), POINTER :: gth_potential
129 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
130 : TYPE(sgp_potential_type), POINTER :: sgp_potential
131 :
132 : IF (PRESENT(ao_coef) .OR. PRESENT(which_l) .OR. &
133 10206 : PRESENT(which_n) .OR. PRESENT(proj_shell_charge)) THEN
134 : IF (.NOT. PRESENT(ao_coef) .OR. .NOT. PRESENT(which_l) .OR. &
135 30 : .NOT. PRESENT(which_n) .OR. .NOT. PRESENT(proj_shell_charge)) THEN
136 0 : CPABORT("ao_coef, which_l, which_n, and proj_shell_charge must be provided together")
137 : END IF
138 : END IF
139 :
140 10236 : NULLIFY (atom)
141 10236 : CALL create_atom_type(atom)
142 :
143 10236 : IF (PRESENT(xc_section)) THEN
144 0 : atom%xc_section => xc_section
145 : ELSE
146 10236 : NULLIFY (atom%xc_section)
147 : END IF
148 :
149 10236 : CALL get_atomic_kind(atomic_kind, z=z)
150 10236 : NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set)
151 : CALL get_qs_kind(qs_kind, zeff=zeff, &
152 : basis_set=orb_basis_set, &
153 : ghost=ghost, &
154 : all_potential=all_potential, &
155 : gth_potential=gth_potential, &
156 10236 : sgp_potential=sgp_potential)
157 :
158 10236 : has_pp = ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
159 :
160 10236 : atom%z = z
161 : CALL set_atom(atom, &
162 : pp_calc=has_pp, &
163 : do_zmp=.FALSE., &
164 : doread=.FALSE., &
165 : read_vxc=.FALSE., &
166 : relativistic=do_nonrel_atom, &
167 : coulomb_integral_type=do_numeric, &
168 10236 : exchange_integral_type=do_numeric)
169 :
170 57229476 : ALLOCATE (potential, integrals)
171 :
172 10236 : IF (PRESENT(confine)) THEN
173 0 : potential%confinement = confine
174 : ELSE
175 10236 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
176 8800 : potential%confinement = .TRUE.
177 : ELSE
178 1436 : potential%confinement = .FALSE.
179 : END IF
180 : END IF
181 10236 : potential%conf_type = poly_conf
182 10236 : potential%acon = 0.1_dp
183 10236 : potential%rcon = 2.0_dp*ptable(z)%vdw_radius*bohr
184 10236 : potential%scon = 2.0_dp
185 :
186 10236 : IF (ASSOCIATED(gth_potential)) THEN
187 8708 : potential%ppot_type = gth_pseudo
188 8708 : CALL get_potential(gth_potential, zeff=zeff)
189 8708 : CALL gth_potential_conversion(gth_potential, potential%gth_pot)
190 8708 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
191 1528 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
192 92 : CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
193 92 : IF (ecp_semi_local) THEN
194 80 : potential%ppot_type = ecp_pseudo
195 80 : CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
196 80 : potential%ecp_pot%symbol = ptable(z)%symbol
197 : ELSE
198 12 : potential%ppot_type = sgp_pseudo
199 12 : CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
200 12 : potential%sgp_pot%symbol = ptable(z)%symbol
201 : END IF
202 92 : CALL get_potential(sgp_potential, zeff=zeff)
203 92 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
204 : ELSE
205 1436 : potential%ppot_type = no_pseudo
206 1436 : CALL set_atom(atom, zcore=z, potential=potential)
207 : END IF
208 :
209 : NULLIFY (basis)
210 194484 : ALLOCATE (basis)
211 :
212 10236 : CALL set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid)
213 :
214 10236 : CALL set_atom(atom, basis=basis)
215 :
216 : ! optimization defaults
217 10236 : atom%optimization%damping = 0.2_dp
218 10236 : atom%optimization%eps_scf = 1.e-6_dp
219 10236 : atom%optimization%eps_diis = 100._dp
220 10236 : atom%optimization%max_iter = 50
221 10236 : atom%optimization%n_diis = 5
222 :
223 : ! set up the electronic state
224 : CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
225 : qs_kind=qs_kind, &
226 : ncalc=ncalc, &
227 : ncore=ncore, &
228 : nelem=nelem, &
229 : edelta=edelta, &
230 10236 : rks=PRESENT(ao_coef))
231 :
232 : ! restricted or unrestricted?
233 1463748 : IF (SUM(ABS(edelta)) > 0.0_dp) THEN
234 64 : uks = .TRUE.
235 64 : CALL set_atom(atom, method_type=do_uks_atom)
236 : ELSE
237 10172 : uks = .FALSE.
238 10172 : CALL set_atom(atom, method_type=do_rks_atom)
239 : END IF
240 :
241 3715668 : ALLOCATE (atom%state)
242 :
243 726756 : atom%state%core = 0._dp
244 511800 : atom%state%core(0:lmat, 1:7) = REAL(ncore(0:lmat, 1:7), dp)
245 726756 : atom%state%occ = 0._dp
246 10236 : IF (uks) THEN
247 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp) + &
248 3200 : edelta(0:lmat, 1:7, 1) + edelta(0:lmat, 1:7, 2)
249 : ELSE
250 508600 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp)
251 : END IF
252 726756 : atom%state%occupation = 0._dp
253 71652 : DO l = 0, lmat
254 : k = 0
255 491328 : DO i = 1, 7
256 491328 : IF (ncalc(l, i) > 0) THEN
257 16203 : k = k + 1
258 16203 : IF (uks) THEN
259 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp) + &
260 118 : edelta(l, i, 1) + edelta(l, i, 2)
261 118 : atom%state%occa(l, k) = 0.5_dp*REAL(ncalc(l, i), dp) + edelta(l, i, 1)
262 118 : atom%state%occb(l, k) = 0.5_dp*REAL(ncalc(l, i), dp) + edelta(l, i, 2)
263 : ELSE
264 16085 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp)
265 : END IF
266 : END IF
267 : END DO
268 61416 : ok = REAL(2*l + 1, KIND=dp)
269 71652 : IF (uks) THEN
270 3072 : DO i = 1, 7
271 2688 : atom%state%occ(l, i) = MIN(atom%state%occ(l, i), 2.0_dp*ok)
272 2688 : atom%state%occa(l, i) = MIN(atom%state%occa(l, i), ok)
273 2688 : atom%state%occb(l, i) = MIN(atom%state%occb(l, i), ok)
274 3072 : atom%state%occupation(l, i) = atom%state%occa(l, i) + atom%state%occb(l, i)
275 : END DO
276 : ELSE
277 488256 : DO i = 1, 7
278 427224 : atom%state%occ(l, i) = MIN(atom%state%occ(l, i), 2.0_dp*ok)
279 488256 : atom%state%occupation(l, i) = MIN(atom%state%occupation(l, i), 2.0_dp*ok)
280 : END DO
281 : END IF
282 : END DO
283 :
284 10236 : IF (PRESENT(proj_shell_charge)) THEN
285 : ! Save the neutral/original calculated shell list.
286 2130 : occupation_neutral(0:lmat, 1:10) = atom%state%occupation(0:lmat, 1:10)
287 30 : IF (uks) THEN
288 : CALL cp_abort(__LOCATION__, "proj_shell_charge in"// &
289 0 : " calculate_atomic_orbitals is implemented only for RKS")
290 : END IF
291 150 : DO l = 0, MIN(3, lmat)
292 120 : qleft = proj_shell_charge(l + 1)
293 150 : IF (qleft > 0.0_dp) THEN
294 : ! Remove electrons from the highest shell in this angular momentum channel.
295 210 : DO k = 7, 1, -1
296 240 : IF (atom%state%occupation(l, k) > 0.0_dp) THEN
297 30 : take = MIN(qleft, atom%state%occupation(l, k))
298 30 : atom%state%occupation(l, k) = atom%state%occupation(l, k) - take
299 30 : qleft = qleft - take
300 30 : IF (qleft <= 1.0e-12_dp) EXIT
301 : END IF
302 : END DO
303 30 : IF (qleft > 1.0e-10_dp) THEN
304 : CALL cp_abort(__LOCATION__, "proj_shell_charge removes"// &
305 0 : " more electrons than available in this angular channel")
306 : END IF
307 : ! Keep atom%state%occ consistent with atom%state%occupation.
308 : qleft = proj_shell_charge(l + 1)
309 210 : DO k = 7, 1, -1
310 240 : IF (atom%state%occ(l, k) > 0.0_dp) THEN
311 30 : take = MIN(qleft, atom%state%occ(l, k))
312 30 : atom%state%occ(l, k) = atom%state%occ(l, k) - take
313 30 : qleft = qleft - take
314 30 : IF (qleft <= 1.0e-12_dp) EXIT
315 : END IF
316 : END DO
317 90 : ELSE IF (qleft < 0.0_dp) THEN
318 : CALL cp_abort(__LOCATION__, "Negative proj_shell_charge"// &
319 0 : " is not supported in this implementation")
320 : END IF
321 : END DO
322 : END IF
323 :
324 10236 : IF (uks) THEN
325 4544 : atom%state%multiplicity = NINT(ABS(SUM(atom%state%occa - atom%state%occb)) + 1)
326 : ELSE
327 10172 : atom%state%multiplicity = -1
328 : END IF
329 :
330 10236 : atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
331 71652 : atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
332 10236 : IF (PRESENT(proj_shell_charge)) THEN
333 30 : atom%state%maxl_calc = get_maxl_occ(occupation_neutral)
334 210 : atom%state%maxn_calc = get_maxn_occ(occupation_neutral)
335 : ELSE
336 10206 : atom%state%maxl_calc = atom%state%maxl_occ
337 71442 : atom%state%maxn_calc = atom%state%maxn_occ
338 : END IF
339 :
340 : ! total number of occupied orbitals
341 10236 : IF (PRESENT(nocc) .AND. ghost) THEN
342 444 : nocc = 0
343 : ELSE IF (PRESENT(nocc)) THEN
344 29742 : nocc = 0
345 69398 : DO l = 0, lmat
346 485786 : DO k = 1, 7
347 475872 : IF (uks) THEN
348 2688 : IF (atom%state%occa(l, k) > 0.0_dp) THEN
349 88 : nocc(1) = nocc(1) + 2*l + 1
350 : END IF
351 2688 : IF (atom%state%occb(l, k) > 0.0_dp) THEN
352 82 : nocc(2) = nocc(2) + 2*l + 1
353 : END IF
354 : ELSE
355 413700 : IF (atom%state%occupation(l, k) > 0.0_dp) THEN
356 15857 : nocc(1) = nocc(1) + 2*l + 1
357 15857 : nocc(2) = nocc(2) + 2*l + 1
358 : END IF
359 : END IF
360 : END DO
361 : END DO
362 : END IF
363 :
364 : ! calculate integrals
365 : ! general integrals
366 : CALL atom_int_setup(integrals, basis, potential=atom%potential, &
367 : eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
368 10236 : eri_exchange=(atom%exchange_integral_type == do_analytic))
369 : ! potential
370 10236 : CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
371 : ! relativistic correction terms
372 10236 : NULLIFY (integrals%tzora, integrals%hdkh)
373 10236 : CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=REAL(atom%zcore, dp))
374 10236 : CALL set_atom(atom, integrals=integrals)
375 :
376 10236 : NULLIFY (orbitals)
377 71652 : mo = MAXVAL(atom%state%maxn_calc)
378 71652 : mb = MAXVAL(atom%basis%nbas)
379 10236 : CALL create_atom_orbs(orbitals, mb, mo)
380 10236 : CALL set_atom(atom, orbitals=orbitals)
381 :
382 10236 : IF (.NOT. ghost) THEN
383 10088 : IF (PRESENT(iunit)) THEN
384 10052 : CALL calculate_atom(atom, iunit)
385 : ELSE
386 36 : CALL calculate_atom(atom, -1)
387 : END IF
388 : END IF
389 10236 : IF (PRESENT(pmat)) THEN
390 : ! recover density matrix in CP2K/GPW order and normalization
391 : CALL get_gto_basis_set(orb_basis_set, &
392 10062 : nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
393 10062 : set_index = 0
394 10062 : shell_index = 0
395 10062 : nbb = 0
396 30368 : DO i = 1, nset
397 69460 : DO j = 1, nshell(i)
398 39092 : l = ls(j, i)
399 59398 : IF (l <= lmat) THEN
400 39092 : nbb(l) = nbb(l) + 1
401 39092 : k = nbb(l)
402 39092 : CPASSERT(k <= 100)
403 39092 : set_index(l, k) = i
404 39092 : shell_index(l, k) = j
405 : END IF
406 : END DO
407 : END DO
408 :
409 10062 : IF (ASSOCIATED(pmat)) THEN
410 0 : DEALLOCATE (pmat)
411 : END IF
412 50298 : ALLOCATE (pmat(nsgf, nsgf, 2))
413 2766610 : pmat = 0._dp
414 20124 : IF (.NOT. ghost) THEN
415 69398 : DO l = 0, lmat
416 59484 : ll = 2*l
417 107936 : DO k1 = 1, atom%basis%nbas(l)
418 177970 : DO k2 = 1, atom%basis%nbas(l)
419 79948 : scal = SQRT(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))/REAL(2*l + 1, KIND=dp)
420 79948 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
421 79948 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
422 118486 : IF (uks) THEN
423 1564 : DO m = 0, ll
424 1102 : pmat(i + m, j + m, 1) = atom%orbitals%pmata(k1, k2, l)*scal
425 1564 : pmat(i + m, j + m, 2) = atom%orbitals%pmatb(k1, k2, l)*scal
426 : END DO
427 : ELSE
428 240946 : DO m = 0, ll
429 240946 : pmat(i + m, j + m, 1) = atom%orbitals%pmat(k1, k2, l)*scal
430 : END DO
431 : END IF
432 : END DO
433 : END DO
434 : END DO
435 9914 : IF (uks) THEN
436 11776 : pmat(:, :, 1) = pmat(:, :, 1) + pmat(:, :, 2)
437 11776 : pmat(:, :, 2) = pmat(:, :, 1) - 2.0_dp*pmat(:, :, 2)
438 : END IF
439 : END IF
440 : END IF
441 :
442 10236 : IF (PRESENT(fmat)) THEN
443 : ! recover fock matrix in CP2K/GPW order.
444 : ! Caution: Normalization is not take care of, so it's probably weird.
445 : CALL get_gto_basis_set(orb_basis_set, &
446 138 : nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
447 138 : set_index = 0
448 138 : shell_index = 0
449 138 : nbb = 0
450 278 : DO i = 1, nset
451 764 : DO j = 1, nshell(i)
452 486 : l = ls(j, i)
453 626 : IF (l <= lmat) THEN
454 486 : nbb(l) = nbb(l) + 1
455 486 : k = nbb(l)
456 486 : CPASSERT(k <= 100)
457 486 : set_index(l, k) = i
458 486 : shell_index(l, k) = j
459 : END IF
460 : END DO
461 : END DO
462 138 : IF (uks) CPABORT("calculate_atomic_orbitals: only RKS is implemented")
463 138 : IF (ASSOCIATED(fmat)) CPABORT("fmat already associated")
464 138 : IF (.NOT. ASSOCIATED(atom%fmat)) CPABORT("atom%fmat not associated")
465 552 : ALLOCATE (fmat(nsgf, nsgf, 1))
466 10140 : fmat = 0.0_dp
467 276 : IF (.NOT. ghost) THEN
468 966 : DO l = 0, lmat
469 828 : ll = 2*l
470 1452 : DO k1 = 1, atom%basis%nbas(l)
471 2152 : DO k2 = 1, atom%basis%nbas(l)
472 838 : scal = SQRT(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))
473 838 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
474 838 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
475 2814 : DO m = 0, ll
476 2328 : fmat(i + m, j + m, 1) = atom%fmat%op(k1, k2, l)/scal
477 : END DO
478 : END DO
479 : END DO
480 : END DO
481 : END IF
482 : END IF
483 :
484 10236 : nr = basis%grid%nr
485 :
486 10236 : IF (PRESENT(density)) THEN
487 6 : IF (ASSOCIATED(density)) DEALLOCATE (density)
488 18 : ALLOCATE (density(nr))
489 6 : IF (ghost) THEN
490 0 : density = 0.0_dp
491 : ELSE
492 6 : CALL atom_density(density, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ)
493 : END IF
494 : END IF
495 :
496 10236 : IF (PRESENT(wavefunction)) THEN
497 6 : CPASSERT(PRESENT(wfninfo))
498 6 : IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
499 6 : IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
500 42 : mo = SUM(atom%state%maxn_occ)
501 36 : ALLOCATE (wavefunction(nr, mo), wfninfo(2, mo))
502 14886 : wavefunction = 0.0_dp
503 6 : IF (.NOT. ghost) THEN
504 : ii = 0
505 42 : DO l = 0, lmat
506 58 : DO i = 1, atom%state%maxn_occ(l)
507 52 : IF (atom%state%occupation(l, i) > 0.0_dp) THEN
508 16 : ii = ii + 1
509 16 : wfninfo(1, ii) = atom%state%occupation(l, i)
510 16 : wfninfo(2, ii) = REAL(l, dp)
511 336 : DO j = 1, atom%basis%nbas(l)
512 : wavefunction(:, ii) = wavefunction(:, ii) + &
513 594896 : atom%orbitals%wfn(j, i, l)*basis%bf(:, j, l)
514 : END DO
515 : END IF
516 : END DO
517 : END DO
518 6 : CPASSERT(mo == ii)
519 : END IF
520 : END IF
521 :
522 10236 : IF (PRESENT(ao_coef)) THEN
523 30 : CPASSERT(PRESENT(which_l))
524 30 : CPASSERT(PRESENT(which_n))
525 30 : CPASSERT(which_l >= 0 .AND. which_l <= lmat)
526 330 : valence_n = COUNT(ncore(which_l, :) > 0) + which_l + 1
527 30 : projector_index = which_n - valence_n + 1
528 30 : IF (projector_index < 1 .OR. projector_index > SIZE(atom%orbitals%wfn, 2)) THEN
529 : CALL cp_abort(__LOCATION__, &
530 : "The requested tensorial atomic projector N is outside "// &
531 0 : "the available orbital range for the selected L.")
532 : END IF
533 30 : IF (ASSOCIATED(ao_coef)) DEALLOCATE (ao_coef)
534 90 : ALLOCATE (ao_coef(SIZE(atom%orbitals%wfn(:, 1, which_l))))
535 150 : ao_coef(:) = atom%orbitals%wfn(:, projector_index, which_l)
536 90 : IF (MAXVAL(ABS(ao_coef)) <= 100.0_dp*EPSILON(1.0_dp)) THEN
537 : CALL cp_abort(__LOCATION__, &
538 : "The requested atomic projector is zero. Check the "// &
539 0 : "specified N and L quantum numbers.")
540 : END IF
541 90 : DO j = 1, SIZE(ao_coef(:))
542 90 : ao_coef(j) = ao_coef(j)*SQRT(atom%integrals%ovlp(j, j, which_l))
543 : END DO
544 : END IF
545 :
546 : ! clean up
547 10236 : CALL atom_int_release(integrals)
548 10236 : CALL atom_ppint_release(integrals)
549 10236 : CALL atom_relint_release(integrals)
550 10236 : CALL release_atom_basis(basis)
551 10236 : CALL release_atom_potential(potential)
552 10236 : CALL release_atom_type(atom)
553 :
554 10236 : DEALLOCATE (potential, basis, integrals)
555 :
556 10236 : END SUBROUTINE calculate_atomic_orbitals
557 :
558 : ! **************************************************************************************************
559 : !> \brief ...
560 : !> \param density ...
561 : !> \param atomic_kind ...
562 : !> \param qs_kind ...
563 : !> \param ngto ...
564 : !> \param iunit ...
565 : !> \param optbasis ... Default=T, if basis should be optimized, if not basis is given in input (density)
566 : !> \param allelectron ...
567 : !> \param confine ...
568 : ! **************************************************************************************************
569 66 : SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, &
570 : optbasis, allelectron, confine)
571 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: density
572 : TYPE(atomic_kind_type), POINTER :: atomic_kind
573 : TYPE(qs_kind_type), POINTER :: qs_kind
574 : INTEGER, INTENT(IN) :: ngto
575 : INTEGER, INTENT(IN), OPTIONAL :: iunit
576 : LOGICAL, INTENT(IN), OPTIONAL :: optbasis, allelectron, confine
577 :
578 : INTEGER, PARAMETER :: num_gto = 40
579 :
580 : INTEGER :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, &
581 : nr, quadtype, relativistic, z
582 : INTEGER, DIMENSION(0:lmat) :: starti
583 : INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
584 66 : INTEGER, DIMENSION(:), POINTER :: econf
585 : LOGICAL :: do_basopt, ecp_semi_local, monovalent
586 : REAL(KIND=dp) :: al, aval, cc, cval, ear, rk, xx, zeff
587 : REAL(KIND=dp), DIMENSION(num_gto+2) :: results
588 : TYPE(all_potential_type), POINTER :: all_potential
589 : TYPE(atom_basis_type), POINTER :: basis
590 : TYPE(atom_integrals), POINTER :: integrals
591 : TYPE(atom_orbitals), POINTER :: orbitals
592 : TYPE(atom_potential_type), POINTER :: potential
593 : TYPE(atom_type), POINTER :: atom
594 : TYPE(grid_atom_type), POINTER :: grid
595 : TYPE(gth_potential_type), POINTER :: gth_potential
596 : TYPE(sgp_potential_type), POINTER :: sgp_potential
597 :
598 66 : NULLIFY (atom)
599 66 : CALL create_atom_type(atom)
600 :
601 66 : CALL get_atomic_kind(atomic_kind, z=z)
602 66 : NULLIFY (all_potential, gth_potential)
603 : CALL get_qs_kind(qs_kind, zeff=zeff, &
604 : all_potential=all_potential, &
605 : gth_potential=gth_potential, &
606 : sgp_potential=sgp_potential, &
607 66 : monovalent=monovalent)
608 :
609 66 : IF (PRESENT(iunit)) THEN
610 8 : iw = iunit
611 : ELSE
612 58 : iw = -1
613 : END IF
614 :
615 66 : IF (PRESENT(allelectron)) THEN
616 4 : IF (allelectron) THEN
617 4 : NULLIFY (gth_potential)
618 4 : zeff = z
619 : END IF
620 : END IF
621 :
622 66 : do_basopt = .TRUE.
623 66 : IF (PRESENT(optbasis)) THEN
624 22 : do_basopt = optbasis
625 : END IF
626 :
627 66 : CPASSERT(ngto <= num_gto)
628 :
629 66 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
630 : ! PP calculation are non-relativistic
631 60 : relativistic = do_nonrel_atom
632 : ELSE
633 : ! AE calculations use DKH2
634 6 : relativistic = do_dkh2_atom
635 : END IF
636 :
637 66 : atom%z = z
638 : CALL set_atom(atom, &
639 : pp_calc=(ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)), &
640 : method_type=do_rks_atom, &
641 : relativistic=relativistic, &
642 : coulomb_integral_type=do_numeric, &
643 72 : exchange_integral_type=do_numeric)
644 :
645 370194 : ALLOCATE (potential, basis, integrals)
646 :
647 66 : IF (PRESENT(confine)) THEN
648 66 : potential%confinement = confine
649 : ELSE
650 0 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
651 0 : potential%confinement = .TRUE.
652 : ELSE
653 0 : potential%confinement = .FALSE.
654 : END IF
655 : END IF
656 66 : potential%conf_type = barrier_conf
657 66 : potential%acon = 200._dp
658 66 : potential%rcon = 4.0_dp
659 66 : potential%scon = 8.0_dp
660 :
661 66 : IF (ASSOCIATED(gth_potential)) THEN
662 60 : potential%ppot_type = gth_pseudo
663 60 : CALL get_potential(gth_potential, zeff=zeff)
664 60 : CALL gth_potential_conversion(gth_potential, potential%gth_pot)
665 60 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
666 6 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
667 0 : CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
668 0 : IF (ecp_semi_local) THEN
669 0 : potential%ppot_type = ecp_pseudo
670 0 : CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
671 0 : potential%ecp_pot%symbol = ptable(z)%symbol
672 : ELSE
673 0 : potential%ppot_type = sgp_pseudo
674 0 : CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
675 0 : potential%sgp_pot%symbol = ptable(z)%symbol
676 : END IF
677 0 : CALL get_potential(sgp_potential, zeff=zeff)
678 0 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
679 : ELSE
680 6 : potential%ppot_type = no_pseudo
681 6 : CALL set_atom(atom, zcore=z, potential=potential)
682 : END IF
683 :
684 : ! atomic grid
685 66 : NULLIFY (grid)
686 66 : ngp = 400
687 66 : quadtype = do_gapw_log
688 66 : CALL allocate_grid_atom(grid)
689 66 : CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
690 66 : grid%nr = ngp
691 66 : basis%grid => grid
692 :
693 66 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
694 :
695 : ! fill in the basis data structures
696 66 : basis%eps_eig = 1.e-12_dp
697 66 : basis%basis_type = GTO_BASIS
698 66 : CALL Clementi_geobas(z, cval, aval, basis%nbas, starti)
699 462 : basis%nprim = basis%nbas
700 462 : m = MAXVAL(basis%nbas)
701 198 : ALLOCATE (basis%am(m, 0:lmat))
702 8886 : basis%am = 0._dp
703 462 : DO l = 0, lmat
704 2314 : DO i = 1, basis%nbas(l)
705 1852 : ll = i - 1 + starti(l)
706 2248 : basis%am(i, l) = aval*cval**(ll)
707 : END DO
708 : END DO
709 :
710 66 : basis%geometrical = .TRUE.
711 66 : basis%aval = aval
712 66 : basis%cval = cval
713 462 : basis%start = starti
714 :
715 : ! initialize basis function on a radial grid
716 66 : nr = basis%grid%nr
717 462 : m = MAXVAL(basis%nbas)
718 330 : ALLOCATE (basis%bf(nr, m, 0:lmat))
719 198 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
720 198 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
721 3378486 : basis%bf = 0._dp
722 3378486 : basis%dbf = 0._dp
723 3378486 : basis%ddbf = 0._dp
724 462 : DO l = 0, lmat
725 2314 : DO i = 1, basis%nbas(l)
726 1852 : al = basis%am(i, l)
727 743048 : DO k = 1, nr
728 740800 : rk = basis%grid%rad(k)
729 740800 : ear = EXP(-al*basis%grid%rad(k)**2)
730 740800 : basis%bf(k, i, l) = rk**l*ear
731 740800 : basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
732 : basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
733 742652 : 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
734 : END DO
735 : END DO
736 : END DO
737 :
738 66 : CALL set_atom(atom, basis=basis)
739 :
740 : ! optimization defaults
741 66 : atom%optimization%damping = 0.2_dp
742 66 : atom%optimization%eps_scf = 1.e-6_dp
743 66 : atom%optimization%eps_diis = 100._dp
744 66 : atom%optimization%max_iter = 50
745 66 : atom%optimization%n_diis = 5
746 :
747 66 : nelem = 0
748 66 : ncore = 0
749 66 : ncalc = 0
750 66 : IF (monovalent) THEN
751 0 : ncalc(0, 1) = 1
752 0 : nelem(0, 1) = 1
753 66 : ELSE IF (ASSOCIATED(gth_potential)) THEN
754 60 : CALL get_potential(gth_potential, elec_conf=econf)
755 60 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
756 6 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
757 0 : CALL get_potential(sgp_potential, elec_conf=econf)
758 0 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
759 : ELSE
760 30 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
761 24 : ll = 2*(2*l + 1)
762 24 : nn = ptable(z)%e_conv(l)
763 24 : ii = 0
764 6 : DO
765 24 : ii = ii + 1
766 24 : IF (nn <= ll) THEN
767 24 : nelem(l, ii) = nn
768 : EXIT
769 : ELSE
770 0 : nelem(l, ii) = ll
771 0 : nn = nn - ll
772 : END IF
773 : END DO
774 : END DO
775 426 : ncalc = nelem - ncore
776 : END IF
777 :
778 66 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
779 0 : nelem = 0
780 0 : ncore = 0
781 0 : ncalc = 0
782 : END IF
783 :
784 23958 : ALLOCATE (atom%state)
785 :
786 4686 : atom%state%core = 0._dp
787 3300 : atom%state%core(0:lmat, 1:7) = REAL(ncore(0:lmat, 1:7), dp)
788 4686 : atom%state%occ = 0._dp
789 3300 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp)
790 4686 : atom%state%occupation = 0._dp
791 66 : atom%state%multiplicity = -1
792 462 : DO l = 0, lmat
793 : k = 0
794 3234 : DO i = 1, 7
795 3168 : IF (ncalc(l, i) > 0) THEN
796 94 : k = k + 1
797 94 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp)
798 : END IF
799 : END DO
800 : END DO
801 :
802 66 : atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
803 462 : atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
804 66 : atom%state%maxl_calc = atom%state%maxl_occ
805 462 : atom%state%maxn_calc = atom%state%maxn_occ
806 :
807 : ! calculate integrals
808 : ! general integrals
809 : CALL atom_int_setup(integrals, basis, potential=atom%potential, &
810 : eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
811 66 : eri_exchange=(atom%exchange_integral_type == do_analytic))
812 : ! potential
813 66 : CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
814 : ! relativistic correction terms
815 66 : NULLIFY (integrals%tzora, integrals%hdkh)
816 66 : CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=REAL(atom%zcore, dp))
817 66 : CALL set_atom(atom, integrals=integrals)
818 :
819 66 : NULLIFY (orbitals)
820 462 : mo = MAXVAL(atom%state%maxn_calc)
821 462 : mb = MAXVAL(atom%basis%nbas)
822 66 : CALL create_atom_orbs(orbitals, mb, mo)
823 66 : CALL set_atom(atom, orbitals=orbitals)
824 :
825 66 : CALL calculate_atom(atom, iw)
826 :
827 66 : IF (do_basopt) THEN
828 44 : CALL atom_fit_density(atom, ngto, 0, iw, results=results)
829 44 : xx = results(1)
830 44 : cc = results(2)
831 428 : DO i = 1, ngto
832 384 : density(i, 1) = xx*cc**i
833 428 : density(i, 2) = results(2 + i)
834 : END DO
835 : ELSE
836 22 : CALL atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results)
837 168 : density(1:ngto, 2) = results(1:ngto)
838 : END IF
839 :
840 : ! clean up
841 66 : CALL atom_int_release(integrals)
842 66 : CALL atom_ppint_release(integrals)
843 66 : CALL atom_relint_release(integrals)
844 66 : CALL release_atom_basis(basis)
845 66 : CALL release_atom_potential(potential)
846 66 : CALL release_atom_type(atom)
847 :
848 66 : DEALLOCATE (potential, basis, integrals)
849 :
850 66 : END SUBROUTINE calculate_atomic_density
851 :
852 : ! **************************************************************************************************
853 : !> \brief ...
854 : !> \param atomic_kind ...
855 : !> \param qs_kind ...
856 : !> \param rel_control ...
857 : !> \param rtmat ...
858 : ! **************************************************************************************************
859 28 : SUBROUTINE calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
860 : TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
861 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
862 : TYPE(rel_control_type), POINTER :: rel_control
863 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rtmat
864 :
865 : INTEGER :: i, ii, ipgf, j, k, k1, k2, l, ll, m, n, &
866 : ngp, nj, nn, nr, ns, nset, nsgf, &
867 : quadtype, relativistic, z
868 : INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
869 : INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
870 28 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
871 28 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, ls
872 : REAL(KIND=dp) :: al, alpha, ear, prefac, rk, zeff
873 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
874 28 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
875 28 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
876 : TYPE(all_potential_type), POINTER :: all_potential
877 : TYPE(atom_basis_type), POINTER :: basis
878 : TYPE(atom_integrals), POINTER :: integrals
879 : TYPE(atom_potential_type), POINTER :: potential
880 : TYPE(atom_type), POINTER :: atom
881 : TYPE(grid_atom_type), POINTER :: grid
882 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
883 :
884 28 : IF (rel_control%rel_method == rel_none) RETURN
885 :
886 28 : NULLIFY (all_potential, orb_basis_set)
887 28 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, all_potential=all_potential)
888 :
889 28 : CPASSERT(ASSOCIATED(orb_basis_set))
890 :
891 28 : IF (ASSOCIATED(all_potential)) THEN
892 : ! only all electron atoms will get the relativistic correction
893 :
894 28 : CALL get_atomic_kind(atomic_kind, z=z)
895 28 : CALL get_qs_kind(qs_kind, zeff=zeff)
896 28 : NULLIFY (atom)
897 28 : CALL create_atom_type(atom)
898 28 : NULLIFY (atom%xc_section)
899 28 : NULLIFY (atom%orbitals)
900 28 : atom%z = z
901 28 : alpha = SQRT(all_potential%alpha_core_charge)
902 :
903 : ! set the method flag
904 28 : SELECT CASE (rel_control%rel_method)
905 : CASE DEFAULT
906 0 : CPABORT("Unknown relativistic method for calculate_atomic_relkin")
907 : CASE (rel_dkh)
908 28 : SELECT CASE (rel_control%rel_DKH_order)
909 : CASE DEFAULT
910 0 : CPABORT("Unknown DKH order for calculate_atomic_relkin")
911 : CASE (0)
912 0 : relativistic = do_dkh0_atom
913 : CASE (1)
914 0 : relativistic = do_dkh1_atom
915 : CASE (2)
916 8 : relativistic = do_dkh2_atom
917 : CASE (3)
918 16 : relativistic = do_dkh3_atom
919 : END SELECT
920 : CASE (rel_zora)
921 28 : SELECT CASE (rel_control%rel_zora_type)
922 : CASE DEFAULT
923 0 : CPABORT("Unknown ZORA type for calculate_atomic_relkin")
924 : CASE (rel_zora_full)
925 0 : CPABORT("ZORA full not yet implemented in calculate_atomic_relkin")
926 : CASE (rel_zora_mp)
927 0 : relativistic = do_zoramp_atom
928 : CASE (rel_sczora_mp)
929 12 : relativistic = do_sczoramp_atom
930 : END SELECT
931 : END SELECT
932 :
933 : CALL set_atom(atom, &
934 : pp_calc=.FALSE., &
935 : method_type=do_rks_atom, &
936 : relativistic=relativistic, &
937 : coulomb_integral_type=do_numeric, &
938 28 : exchange_integral_type=do_numeric)
939 :
940 157052 : ALLOCATE (potential, basis, integrals)
941 :
942 28 : potential%ppot_type = no_pseudo
943 28 : CALL set_atom(atom, zcore=z, potential=potential)
944 :
945 : CALL get_gto_basis_set(orb_basis_set, &
946 : nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
947 28 : first_sgf=first_sgf, last_sgf=last_sgf)
948 :
949 28 : NULLIFY (grid)
950 28 : ngp = 400
951 28 : quadtype = do_gapw_log
952 28 : CALL allocate_grid_atom(grid)
953 28 : CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
954 28 : grid%nr = ngp
955 28 : basis%grid => grid
956 :
957 28 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
958 28 : basis%basis_type = CGTO_BASIS
959 28 : basis%eps_eig = 1.e-12_dp
960 :
961 : ! fill in the basis data structures
962 28 : set_index = 0
963 28 : shell_index = 0
964 196 : basis%nprim = 0
965 196 : basis%nbas = 0
966 138 : DO i = 1, nset
967 220 : DO j = lmin(i), MIN(lmax(i), lmat)
968 220 : basis%nprim(j) = basis%nprim(j) + npgf(i)
969 : END DO
970 492 : DO j = 1, nshell(i)
971 354 : l = ls(j, i)
972 464 : IF (l <= lmat) THEN
973 354 : basis%nbas(l) = basis%nbas(l) + 1
974 354 : k = basis%nbas(l)
975 354 : CPASSERT(k <= 100)
976 354 : set_index(l, k) = i
977 354 : shell_index(l, k) = j
978 : END IF
979 : END DO
980 : END DO
981 :
982 196 : nj = MAXVAL(basis%nprim)
983 196 : ns = MAXVAL(basis%nbas)
984 84 : ALLOCATE (basis%am(nj, 0:lmat))
985 2308 : basis%am = 0._dp
986 140 : ALLOCATE (basis%cm(nj, ns, 0:lmat))
987 18448 : basis%cm = 0._dp
988 196 : DO j = 0, lmat
989 : nj = 0
990 : ns = 0
991 856 : DO i = 1, nset
992 828 : IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
993 792 : DO ipgf = 1, npgf(i)
994 792 : basis%am(nj + ipgf, j) = zet(ipgf, i)
995 : END DO
996 464 : DO ii = 1, nshell(i)
997 464 : IF (ls(ii, i) == j) THEN
998 354 : ns = ns + 1
999 4204 : DO ipgf = 1, npgf(i)
1000 4204 : basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
1001 : END DO
1002 : END IF
1003 : END DO
1004 110 : nj = nj + npgf(i)
1005 : END IF
1006 : END DO
1007 : END DO
1008 :
1009 : ! Normalization as used in the atomic code
1010 : ! We have to undo the Quickstep normalization
1011 196 : DO j = 0, lmat
1012 168 : prefac = 2.0_dp*SQRT(pi/dfac(2*j + 1))
1013 878 : DO ipgf = 1, basis%nprim(j)
1014 6292 : DO ii = 1, basis%nbas(j)
1015 6124 : basis%cm(ipgf, ii, j) = prefac*basis%cm(ipgf, ii, j)
1016 : END DO
1017 : END DO
1018 : END DO
1019 :
1020 : ! initialize basis function on a radial grid
1021 28 : nr = basis%grid%nr
1022 196 : m = MAXVAL(basis%nbas)
1023 140 : ALLOCATE (basis%bf(nr, m, 0:lmat))
1024 84 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
1025 84 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
1026 :
1027 370720 : basis%bf = 0._dp
1028 370720 : basis%dbf = 0._dp
1029 370720 : basis%ddbf = 0._dp
1030 196 : DO l = 0, lmat
1031 878 : DO i = 1, basis%nprim(l)
1032 682 : al = basis%am(i, l)
1033 273650 : DO k = 1, nr
1034 272800 : rk = basis%grid%rad(k)
1035 272800 : ear = EXP(-al*basis%grid%rad(k)**2)
1036 2450282 : DO j = 1, basis%nbas(l)
1037 2176800 : basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
1038 : basis%dbf(k, j, l) = basis%dbf(k, j, l) &
1039 2176800 : + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
1040 : basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
1041 : (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)* &
1042 2449600 : rk**(l) + 4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
1043 : END DO
1044 : END DO
1045 : END DO
1046 : END DO
1047 :
1048 28 : CALL set_atom(atom, basis=basis)
1049 :
1050 : ! optimization defaults
1051 28 : atom%optimization%damping = 0.2_dp
1052 28 : atom%optimization%eps_scf = 1.e-6_dp
1053 28 : atom%optimization%eps_diis = 100._dp
1054 28 : atom%optimization%max_iter = 50
1055 28 : atom%optimization%n_diis = 5
1056 :
1057 : ! electronic state
1058 28 : nelem = 0
1059 28 : ncore = 0
1060 28 : ncalc = 0
1061 140 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
1062 112 : ll = 2*(2*l + 1)
1063 112 : nn = ptable(z)%e_conv(l)
1064 112 : ii = 0
1065 28 : DO
1066 160 : ii = ii + 1
1067 160 : IF (nn <= ll) THEN
1068 112 : nelem(l, ii) = nn
1069 : EXIT
1070 : ELSE
1071 48 : nelem(l, ii) = ll
1072 48 : nn = nn - ll
1073 : END IF
1074 : END DO
1075 : END DO
1076 1988 : ncalc = nelem - ncore
1077 :
1078 28 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
1079 : nelem = 0
1080 0 : ncore = 0
1081 0 : ncalc = 0
1082 : END IF
1083 :
1084 10164 : ALLOCATE (atom%state)
1085 :
1086 1988 : atom%state%core = 0._dp
1087 1400 : atom%state%core(0:lmat, 1:7) = REAL(ncore(0:lmat, 1:7), dp)
1088 1988 : atom%state%occ = 0._dp
1089 1400 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp)
1090 1988 : atom%state%occupation = 0._dp
1091 28 : atom%state%multiplicity = -1
1092 196 : DO l = 0, lmat
1093 : k = 0
1094 1372 : DO i = 1, 7
1095 1344 : IF (ncalc(l, i) > 0) THEN
1096 98 : k = k + 1
1097 98 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp)
1098 : END IF
1099 : END DO
1100 : END DO
1101 :
1102 28 : atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
1103 196 : atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
1104 28 : atom%state%maxl_calc = atom%state%maxl_occ
1105 196 : atom%state%maxn_calc = atom%state%maxn_occ
1106 :
1107 : ! calculate integrals
1108 : ! general integrals
1109 28 : CALL atom_int_setup(integrals, basis)
1110 : ! potential
1111 28 : CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
1112 : ! relativistic correction terms
1113 28 : NULLIFY (integrals%tzora, integrals%hdkh)
1114 : CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=REAL(atom%zcore, dp), &
1115 28 : alpha=alpha)
1116 28 : CALL set_atom(atom, integrals=integrals)
1117 :
1118 : ! for DKH we need erfc integrals to correct non-relativistic
1119 13996 : integrals%core = 0.0_dp
1120 196 : DO l = 0, lmat
1121 168 : n = integrals%n(l)
1122 168 : m = basis%nprim(l)
1123 488 : ALLOCATE (omat(m, m))
1124 :
1125 168 : CALL sg_erfc(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
1126 168 : integrals%core(1:n, 1:n, l) = MATMUL(TRANSPOSE(basis%cm(1:m, 1:n, l)), &
1127 432156 : MATMUL(omat(1:m, 1:m), basis%cm(1:m, 1:n, l)))
1128 :
1129 196 : DEALLOCATE (omat)
1130 : END DO
1131 :
1132 : ! recover relativistic kinetic matrix in CP2K/GPW order and normalization
1133 28 : IF (ASSOCIATED(rtmat)) THEN
1134 0 : DEALLOCATE (rtmat)
1135 : END IF
1136 112 : ALLOCATE (rtmat(nsgf, nsgf))
1137 160300 : rtmat = 0._dp
1138 196 : DO l = 0, lmat
1139 168 : ll = 2*l
1140 550 : DO k1 = 1, basis%nbas(l)
1141 4872 : DO k2 = 1, basis%nbas(l)
1142 4350 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
1143 4350 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
1144 354 : SELECT CASE (atom%relativistic)
1145 : CASE DEFAULT
1146 0 : CPABORT("Unknown relativistic type for calculate_atomic_relkin")
1147 : CASE (do_zoramp_atom, do_sczoramp_atom)
1148 14656 : DO m = 0, ll
1149 14656 : rtmat(i + m, j + m) = integrals%tzora(k1, k2, l)
1150 : END DO
1151 : CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
1152 5046 : DO m = 0, ll
1153 : rtmat(i + m, j + m) = integrals%hdkh(k1, k2, l) - integrals%kin(k1, k2, l) + &
1154 1052 : atom%zcore*integrals%core(k1, k2, l)
1155 : END DO
1156 : END SELECT
1157 : END DO
1158 : END DO
1159 : END DO
1160 1006 : DO k1 = 1, nsgf
1161 81142 : DO k2 = k1, nsgf
1162 80136 : rtmat(k1, k2) = 0.5_dp*(rtmat(k1, k2) + rtmat(k2, k1))
1163 81114 : rtmat(k2, k1) = rtmat(k1, k2)
1164 : END DO
1165 : END DO
1166 :
1167 : ! clean up
1168 28 : CALL atom_int_release(integrals)
1169 28 : CALL atom_ppint_release(integrals)
1170 28 : CALL atom_relint_release(integrals)
1171 28 : CALL release_atom_basis(basis)
1172 28 : CALL release_atom_potential(potential)
1173 28 : CALL release_atom_type(atom)
1174 :
1175 84 : DEALLOCATE (potential, basis, integrals)
1176 :
1177 : ELSE
1178 :
1179 0 : IF (ASSOCIATED(rtmat)) THEN
1180 0 : DEALLOCATE (rtmat)
1181 : END IF
1182 0 : NULLIFY (rtmat)
1183 :
1184 : END IF
1185 :
1186 56 : END SUBROUTINE calculate_atomic_relkin
1187 :
1188 : ! **************************************************************************************************
1189 : !> \brief ...
1190 : !> \param gth_potential ...
1191 : !> \param gth_atompot ...
1192 : ! **************************************************************************************************
1193 26322 : SUBROUTINE gth_potential_conversion(gth_potential, gth_atompot)
1194 : TYPE(gth_potential_type), POINTER :: gth_potential
1195 : TYPE(atom_gthpot_type) :: gth_atompot
1196 :
1197 : INTEGER :: i, j, l, lm, n, ne, nexp_lpot, nexp_lsd, &
1198 : nexp_nlcc
1199 8774 : INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nct_nlcc, nppnl, &
1200 8774 : ppeconf
1201 : LOGICAL :: lpot_present, lsd_present, nlcc_present, &
1202 : soc_present
1203 : REAL(KIND=dp) :: ac, zeff
1204 8774 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_lpot, alpha_lsd, alpha_nlcc, ap, ce
1205 8774 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, cval_lsd, cval_nlcc
1206 8774 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hp, kp
1207 :
1208 : CALL get_potential(gth_potential, &
1209 : zeff=zeff, &
1210 : elec_conf=ppeconf, &
1211 : alpha_core_charge=ac, &
1212 : nexp_ppl=ne, &
1213 : cexp_ppl=ce, &
1214 : lppnl=lm, &
1215 : nprj_ppnl=nppnl, &
1216 : alpha_ppnl=ap, &
1217 : kprj_ppnl=kp, &
1218 8774 : hprj_ppnl=hp)
1219 :
1220 8774 : gth_atompot%zion = zeff
1221 8774 : gth_atompot%rc = SQRT(0.5_dp/ac)
1222 8774 : gth_atompot%ncl = ne
1223 52644 : gth_atompot%cl(:) = 0._dp
1224 8774 : IF (ac > 0._dp) THEN
1225 25882 : DO i = 1, ne
1226 25882 : gth_atompot%cl(i) = ce(i)/(2._dp*ac)**(i - 1)
1227 : END DO
1228 : END IF
1229 : !extended type
1230 8774 : gth_atompot%lpotextended = .FALSE.
1231 8774 : gth_atompot%lsdpot = .FALSE.
1232 8774 : gth_atompot%nlcc = .FALSE.
1233 8774 : gth_atompot%nexp_lpot = 0
1234 8774 : gth_atompot%nexp_lsd = 0
1235 8774 : gth_atompot%nexp_nlcc = 0
1236 : CALL get_potential(gth_potential, &
1237 : lpot_present=lpot_present, &
1238 : lsd_present=lsd_present, &
1239 8774 : nlcc_present=nlcc_present)
1240 8774 : IF (lpot_present) THEN
1241 : CALL get_potential(gth_potential, &
1242 : nexp_lpot=nexp_lpot, &
1243 : alpha_lpot=alpha_lpot, &
1244 : nct_lpot=nct_lpot, &
1245 8 : cval_lpot=cval_lpot)
1246 8 : gth_atompot%lpotextended = .TRUE.
1247 8 : gth_atompot%nexp_lpot = nexp_lpot
1248 20 : gth_atompot%alpha_lpot(1:nexp_lpot) = SQRT(0.5_dp/alpha_lpot(1:nexp_lpot))
1249 20 : gth_atompot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot)
1250 20 : DO j = 1, nexp_lpot
1251 12 : ac = alpha_lpot(j)
1252 68 : DO i = 1, 4
1253 60 : gth_atompot%cval_lpot(i, j) = cval_lpot(i, j)/(2._dp*ac)**(i - 1)
1254 : END DO
1255 : END DO
1256 : END IF
1257 8774 : IF (lsd_present) THEN
1258 : CALL get_potential(gth_potential, &
1259 : nexp_lsd=nexp_lsd, &
1260 : alpha_lsd=alpha_lsd, &
1261 : nct_lsd=nct_lsd, &
1262 0 : cval_lsd=cval_lsd)
1263 0 : gth_atompot%lsdpot = .TRUE.
1264 0 : gth_atompot%nexp_lsd = nexp_lsd
1265 0 : gth_atompot%alpha_lsd(1:nexp_lsd) = SQRT(0.5_dp/alpha_lsd(1:nexp_lsd))
1266 0 : gth_atompot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd)
1267 0 : DO j = 1, nexp_lpot
1268 0 : ac = alpha_lsd(j)
1269 0 : DO i = 1, 4
1270 0 : gth_atompot%cval_lsd(i, j) = cval_lsd(i, j)/(2._dp*ac)**(i - 1)
1271 : END DO
1272 : END DO
1273 : END IF
1274 :
1275 : ! nonlocal part
1276 61418 : gth_atompot%nl(:) = 0
1277 61418 : gth_atompot%rcnl(:) = 0._dp
1278 1114298 : gth_atompot%hnl(:, :, :) = 0._dp
1279 17606 : DO l = 0, lm
1280 8832 : n = nppnl(l)
1281 8832 : gth_atompot%nl(l) = n
1282 8832 : gth_atompot%rcnl(l) = SQRT(0.5_dp/ap(l))
1283 31976 : gth_atompot%hnl(1:n, 1:n, l) = hp(1:n, 1:n, l)
1284 : END DO
1285 :
1286 : ! SOC
1287 8774 : CALL get_potential(gth_potential, soc_present=soc_present)
1288 8774 : gth_atompot%soc = soc_present
1289 1114298 : gth_atompot%knl = 0.0_dp
1290 8774 : IF (soc_present) THEN
1291 118 : DO l = 1, lm
1292 60 : n = nppnl(l)
1293 390 : gth_atompot%knl(1:n, 1:n, l) = kp(1:n, 1:n, l)
1294 : END DO
1295 : END IF
1296 :
1297 8774 : IF (nlcc_present) THEN
1298 : CALL get_potential(gth_potential, &
1299 : nexp_nlcc=nexp_nlcc, &
1300 : alpha_nlcc=alpha_nlcc, &
1301 : nct_nlcc=nct_nlcc, &
1302 32 : cval_nlcc=cval_nlcc)
1303 32 : gth_atompot%nlcc = .TRUE.
1304 32 : gth_atompot%nexp_nlcc = nexp_nlcc
1305 64 : gth_atompot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc)
1306 64 : gth_atompot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc)
1307 192 : gth_atompot%cval_nlcc(1:4, 1:nexp_nlcc) = cval_nlcc(1:4, 1:nexp_nlcc)
1308 : END IF
1309 :
1310 8774 : END SUBROUTINE gth_potential_conversion
1311 :
1312 : ! **************************************************************************************************
1313 : !> \brief ...
1314 : !> \param sgp_potential ...
1315 : !> \param sgp_atompot ...
1316 : ! **************************************************************************************************
1317 36 : SUBROUTINE sgp_potential_conversion(sgp_potential, sgp_atompot)
1318 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1319 : TYPE(atom_sgppot_type) :: sgp_atompot
1320 :
1321 : INTEGER :: lm, n
1322 12 : INTEGER, DIMENSION(:), POINTER :: ppeconf
1323 : LOGICAL :: nlcc_present
1324 : REAL(KIND=dp) :: ac, zeff
1325 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: ap, ce
1326 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hhp
1327 12 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ccp
1328 :
1329 : CALL get_potential(sgp_potential, &
1330 : name=sgp_atompot%pname, &
1331 : zeff=zeff, &
1332 : elec_conf=ppeconf, &
1333 12 : alpha_core_charge=ac)
1334 12 : sgp_atompot%zion = zeff
1335 12 : sgp_atompot%ac_local = ac
1336 60 : sgp_atompot%econf(0:3) = ppeconf(0:3)
1337 : CALL get_potential(sgp_potential, lmax=lm, &
1338 : is_nonlocal=sgp_atompot%is_nonlocal, &
1339 12 : n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
1340 : ! nonlocal
1341 48 : sgp_atompot%has_nonlocal = ANY(sgp_atompot%is_nonlocal)
1342 12 : sgp_atompot%lmax = lm
1343 12 : IF (sgp_atompot%has_nonlocal) THEN
1344 6 : CPASSERT(n <= SIZE(sgp_atompot%a_nonlocal))
1345 6 : sgp_atompot%n_nonlocal = n
1346 54 : sgp_atompot%a_nonlocal(1:n) = ap(1:n)
1347 60 : sgp_atompot%h_nonlocal(1:n, 0:lm) = hhp(1:n, 0:lm)
1348 444 : sgp_atompot%c_nonlocal(1:n, 1:n, 0:lm) = ccp(1:n, 1:n, 0:lm)
1349 : END IF
1350 : ! local
1351 12 : CALL get_potential(sgp_potential, n_local=n, a_local=ap, c_local=ce)
1352 12 : CPASSERT(n <= SIZE(sgp_atompot%a_local))
1353 12 : sgp_atompot%n_local = n
1354 156 : sgp_atompot%a_local(1:n) = ap(1:n)
1355 156 : sgp_atompot%c_local(1:n) = ce(1:n)
1356 : ! NLCC
1357 : CALL get_potential(sgp_potential, has_nlcc=nlcc_present, &
1358 12 : n_nlcc=n, a_nlcc=ap, c_nlcc=ce)
1359 12 : IF (nlcc_present) THEN
1360 0 : sgp_atompot%has_nlcc = .TRUE.
1361 0 : CPASSERT(n <= SIZE(sgp_atompot%a_nlcc))
1362 0 : sgp_atompot%n_nlcc = n
1363 0 : sgp_atompot%a_nlcc(1:n) = ap(1:n)
1364 0 : sgp_atompot%c_nlcc(1:n) = ce(1:n)
1365 : ELSE
1366 12 : sgp_atompot%has_nlcc = .FALSE.
1367 : END IF
1368 :
1369 12 : END SUBROUTINE sgp_potential_conversion
1370 :
1371 : ! **************************************************************************************************
1372 : !> \brief ...
1373 : !> \param sgp_potential ...
1374 : !> \param ecp_atompot ...
1375 : ! **************************************************************************************************
1376 160 : SUBROUTINE ecp_potential_conversion(sgp_potential, ecp_atompot)
1377 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1378 : TYPE(atom_ecppot_type) :: ecp_atompot
1379 :
1380 80 : INTEGER, DIMENSION(:), POINTER :: ppeconf
1381 : LOGICAL :: ecp_local, ecp_semi_local
1382 : REAL(KIND=dp) :: zeff
1383 :
1384 80 : CALL get_potential(sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
1385 80 : CPASSERT(ecp_semi_local .AND. ecp_local)
1386 : CALL get_potential(sgp_potential, &
1387 : name=ecp_atompot%pname, &
1388 : zeff=zeff, &
1389 80 : elec_conf=ppeconf)
1390 80 : ecp_atompot%zion = zeff
1391 400 : ecp_atompot%econf(0:3) = ppeconf(0:3)
1392 80 : CALL get_potential(sgp_potential, sl_lmax=ecp_atompot%lmax)
1393 : ! local
1394 : CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, &
1395 80 : aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc)
1396 : ! nonlocal
1397 : CALL get_potential(sgp_potential, npot=ecp_atompot%npot, nrpot=ecp_atompot%nrpot, &
1398 80 : apot=ecp_atompot%apot, bpot=ecp_atompot%bpot)
1399 :
1400 80 : END SUBROUTINE ecp_potential_conversion
1401 : ! **************************************************************************************************
1402 :
1403 168 : END MODULE atom_kind_orbitals
|
