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