Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
87 25068 : SUBROUTINE calculate_atomic_orbitals(atomic_kind, qs_kind, agrid, iunit, pmat, fmat, &
88 8356 : density, wavefunction, wfninfo, confine, xc_section, nocc)
89 : TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
90 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
91 : TYPE(grid_atom_type), OPTIONAL :: agrid
92 : INTEGER, INTENT(IN), OPTIONAL :: iunit
93 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
94 : POINTER :: pmat, fmat
95 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: density
96 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: wavefunction, wfninfo
97 : LOGICAL, INTENT(IN), OPTIONAL :: confine
98 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
99 : INTEGER, DIMENSION(:), OPTIONAL :: nocc
100 :
101 : INTEGER :: i, ii, j, k, k1, k2, l, ll, m, mb, mo, &
102 : nr, nset, nsgf, z
103 : INTEGER, DIMENSION(0:lmat) :: nbb
104 : INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
105 : INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
106 8356 : INTEGER, DIMENSION(:), POINTER :: nshell
107 8356 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, ls
108 : LOGICAL :: ecp_semi_local, ghost, has_pp, uks
109 : REAL(KIND=dp) :: ok, scal, zeff
110 : REAL(KIND=dp), DIMENSION(0:lmat, 10, 2) :: edelta
111 : TYPE(all_potential_type), POINTER :: all_potential
112 : TYPE(atom_basis_type), POINTER :: basis
113 : TYPE(atom_integrals), POINTER :: integrals
114 : TYPE(atom_orbitals), POINTER :: orbitals
115 : TYPE(atom_potential_type), POINTER :: potential
116 : TYPE(atom_type), POINTER :: atom
117 : TYPE(gth_potential_type), POINTER :: gth_potential
118 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
119 : TYPE(sgp_potential_type), POINTER :: sgp_potential
120 :
121 8356 : NULLIFY (atom)
122 8356 : CALL create_atom_type(atom)
123 :
124 8356 : IF (PRESENT(xc_section)) THEN
125 0 : atom%xc_section => xc_section
126 : ELSE
127 8356 : NULLIFY (atom%xc_section)
128 : END IF
129 :
130 8356 : CALL get_atomic_kind(atomic_kind, z=z)
131 8356 : NULLIFY (all_potential, gth_potential, sgp_potential, orb_basis_set)
132 : CALL get_qs_kind(qs_kind, zeff=zeff, &
133 : basis_set=orb_basis_set, &
134 : ghost=ghost, &
135 : all_potential=all_potential, &
136 : gth_potential=gth_potential, &
137 8356 : sgp_potential=sgp_potential)
138 :
139 8356 : has_pp = ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
140 :
141 8356 : atom%z = z
142 : CALL set_atom(atom, &
143 : pp_calc=has_pp, &
144 : do_zmp=.FALSE., &
145 : doread=.FALSE., &
146 : read_vxc=.FALSE., &
147 : relativistic=do_nonrel_atom, &
148 : coulomb_integral_type=do_numeric, &
149 8356 : exchange_integral_type=do_numeric)
150 :
151 45665540 : ALLOCATE (potential, integrals)
152 :
153 8356 : IF (PRESENT(confine)) THEN
154 0 : potential%confinement = confine
155 : ELSE
156 8356 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
157 7270 : potential%confinement = .TRUE.
158 : ELSE
159 1086 : potential%confinement = .FALSE.
160 : END IF
161 : END IF
162 8356 : potential%conf_type = poly_conf
163 8356 : potential%acon = 0.1_dp
164 8356 : potential%rcon = 2.0_dp*ptable(z)%vdw_radius*bohr
165 8356 : potential%scon = 2.0_dp
166 :
167 8356 : IF (ASSOCIATED(gth_potential)) THEN
168 7246 : potential%ppot_type = gth_pseudo
169 7246 : CALL get_potential(gth_potential, zeff=zeff)
170 7246 : CALL gth_potential_conversion(gth_potential, potential%gth_pot)
171 7246 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
172 1110 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
173 24 : CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
174 24 : IF (ecp_semi_local) THEN
175 12 : potential%ppot_type = ecp_pseudo
176 12 : CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
177 12 : potential%ecp_pot%symbol = ptable(z)%symbol
178 : ELSE
179 12 : potential%ppot_type = sgp_pseudo
180 12 : CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
181 12 : potential%sgp_pot%symbol = ptable(z)%symbol
182 : END IF
183 24 : CALL get_potential(sgp_potential, zeff=zeff)
184 24 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
185 : ELSE
186 1086 : potential%ppot_type = no_pseudo
187 1086 : CALL set_atom(atom, zcore=z, potential=potential)
188 : END IF
189 :
190 : NULLIFY (basis)
191 158764 : ALLOCATE (basis)
192 8356 : CALL set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid)
193 :
194 8356 : CALL set_atom(atom, basis=basis)
195 :
196 : ! optimization defaults
197 8356 : atom%optimization%damping = 0.2_dp
198 8356 : atom%optimization%eps_scf = 1.e-6_dp
199 8356 : atom%optimization%eps_diis = 100._dp
200 8356 : atom%optimization%max_iter = 50
201 8356 : atom%optimization%n_diis = 5
202 :
203 : ! set up the electronic state
204 : CALL init_atom_electronic_state(atomic_kind=atomic_kind, &
205 : qs_kind=qs_kind, &
206 : ncalc=ncalc, &
207 : ncore=ncore, &
208 : nelem=nelem, &
209 8356 : edelta=edelta)
210 :
211 : ! restricted or unrestricted?
212 1194908 : IF (SUM(ABS(edelta)) > 0.0_dp) THEN
213 56 : uks = .TRUE.
214 56 : CALL set_atom(atom, method_type=do_uks_atom)
215 : ELSE
216 8300 : uks = .FALSE.
217 8300 : CALL set_atom(atom, method_type=do_rks_atom)
218 : END IF
219 :
220 3033228 : ALLOCATE (atom%state)
221 :
222 593276 : atom%state%core = 0._dp
223 417800 : atom%state%core(0:lmat, 1:7) = REAL(ncore(0:lmat, 1:7), dp)
224 593276 : atom%state%occ = 0._dp
225 8356 : IF (uks) THEN
226 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp) + &
227 2800 : edelta(0:lmat, 1:7, 1) + edelta(0:lmat, 1:7, 2)
228 : ELSE
229 415000 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp)
230 : END IF
231 593276 : atom%state%occupation = 0._dp
232 58492 : DO l = 0, lmat
233 : k = 0
234 401088 : DO i = 1, 7
235 401088 : IF (ncalc(l, i) > 0) THEN
236 13255 : k = k + 1
237 13255 : IF (uks) THEN
238 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp) + &
239 104 : edelta(l, i, 1) + edelta(l, i, 2)
240 104 : atom%state%occa(l, k) = 0.5_dp*REAL(ncalc(l, i), dp) + edelta(l, i, 1)
241 104 : atom%state%occb(l, k) = 0.5_dp*REAL(ncalc(l, i), dp) + edelta(l, i, 2)
242 : ELSE
243 13151 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp)
244 : END IF
245 : END IF
246 : END DO
247 50136 : ok = REAL(2*l + 1, KIND=dp)
248 58492 : IF (uks) THEN
249 2688 : DO i = 1, 7
250 2352 : atom%state%occ(l, i) = MIN(atom%state%occ(l, i), 2.0_dp*ok)
251 2352 : atom%state%occa(l, i) = MIN(atom%state%occa(l, i), ok)
252 2352 : atom%state%occb(l, i) = MIN(atom%state%occb(l, i), ok)
253 2688 : atom%state%occupation(l, i) = atom%state%occa(l, i) + atom%state%occb(l, i)
254 : END DO
255 : ELSE
256 398400 : DO i = 1, 7
257 348600 : atom%state%occ(l, i) = MIN(atom%state%occ(l, i), 2.0_dp*ok)
258 398400 : atom%state%occupation(l, i) = MIN(atom%state%occupation(l, i), 2.0_dp*ok)
259 : END DO
260 : END IF
261 : END DO
262 8356 : IF (uks) THEN
263 3976 : atom%state%multiplicity = NINT(ABS(SUM(atom%state%occa - atom%state%occb)) + 1)
264 : ELSE
265 8300 : atom%state%multiplicity = -1
266 : END IF
267 :
268 8356 : atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
269 58492 : atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
270 8356 : atom%state%maxl_calc = atom%state%maxl_occ
271 58492 : atom%state%maxn_calc = atom%state%maxn_occ
272 :
273 : ! total number of occupied orbitals
274 8356 : IF (PRESENT(nocc) .AND. ghost) THEN
275 420 : nocc = 0
276 : ELSEIF (PRESENT(nocc)) THEN
277 24252 : nocc = 0
278 56588 : DO l = 0, lmat
279 396116 : DO k = 1, 7
280 388032 : IF (uks) THEN
281 2352 : IF (atom%state%occa(l, k) > 0.0_dp) THEN
282 74 : nocc(1) = nocc(1) + 2*l + 1
283 : END IF
284 2352 : IF (atom%state%occb(l, k) > 0.0_dp) THEN
285 74 : nocc(2) = nocc(2) + 2*l + 1
286 : END IF
287 : ELSE
288 337176 : IF (atom%state%occupation(l, k) > 0.0_dp) THEN
289 12981 : nocc(1) = nocc(1) + 2*l + 1
290 12981 : nocc(2) = nocc(2) + 2*l + 1
291 : END IF
292 : END IF
293 : END DO
294 : END DO
295 : END IF
296 :
297 : ! calculate integrals
298 : ! general integrals
299 : CALL atom_int_setup(integrals, basis, potential=atom%potential, &
300 : eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
301 8356 : eri_exchange=(atom%exchange_integral_type == do_analytic))
302 : ! potential
303 8356 : CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
304 : ! relativistic correction terms
305 8356 : NULLIFY (integrals%tzora, integrals%hdkh)
306 8356 : CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=REAL(atom%zcore, dp))
307 8356 : CALL set_atom(atom, integrals=integrals)
308 :
309 8356 : NULLIFY (orbitals)
310 58492 : mo = MAXVAL(atom%state%maxn_calc)
311 58492 : mb = MAXVAL(atom%basis%nbas)
312 8356 : CALL create_atom_orbs(orbitals, mb, mo)
313 8356 : CALL set_atom(atom, orbitals=orbitals)
314 :
315 8356 : IF (.NOT. ghost) THEN
316 8216 : IF (PRESENT(iunit)) THEN
317 8214 : CALL calculate_atom(atom, iunit)
318 : ELSE
319 2 : CALL calculate_atom(atom, -1)
320 : END IF
321 : END IF
322 :
323 8356 : IF (PRESENT(pmat)) THEN
324 : ! recover density matrix in CP2K/GPW order and normalization
325 : CALL get_gto_basis_set(orb_basis_set, &
326 8224 : nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
327 8224 : set_index = 0
328 8224 : shell_index = 0
329 8224 : nbb = 0
330 24988 : DO i = 1, nset
331 57574 : DO j = 1, nshell(i)
332 32586 : l = ls(j, i)
333 49350 : IF (l <= lmat) THEN
334 32586 : nbb(l) = nbb(l) + 1
335 32586 : k = nbb(l)
336 32586 : CPASSERT(k <= 100)
337 32586 : set_index(l, k) = i
338 32586 : shell_index(l, k) = j
339 : END IF
340 : END DO
341 : END DO
342 :
343 8224 : IF (ASSOCIATED(pmat)) THEN
344 0 : DEALLOCATE (pmat)
345 : END IF
346 41108 : ALLOCATE (pmat(nsgf, nsgf, 2))
347 2416312 : pmat = 0._dp
348 16448 : IF (.NOT. ghost) THEN
349 56588 : DO l = 0, lmat
350 48504 : ll = 2*l
351 88648 : DO k1 = 1, atom%basis%nbas(l)
352 148226 : DO k2 = 1, atom%basis%nbas(l)
353 67662 : scal = SQRT(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))/REAL(2*l + 1, KIND=dp)
354 67662 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
355 67662 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
356 99722 : IF (uks) THEN
357 1360 : DO m = 0, ll
358 964 : pmat(i + m, j + m, 1) = atom%orbitals%pmata(k1, k2, l)*scal
359 1360 : pmat(i + m, j + m, 2) = atom%orbitals%pmatb(k1, k2, l)*scal
360 : END DO
361 : ELSE
362 205754 : DO m = 0, ll
363 205754 : pmat(i + m, j + m, 1) = atom%orbitals%pmat(k1, k2, l)*scal
364 : END DO
365 : END IF
366 : END DO
367 : END DO
368 : END DO
369 8084 : IF (uks) THEN
370 10184 : pmat(:, :, 1) = pmat(:, :, 1) + pmat(:, :, 2)
371 10184 : pmat(:, :, 2) = pmat(:, :, 1) - 2.0_dp*pmat(:, :, 2)
372 : END IF
373 : END IF
374 : END IF
375 :
376 8356 : IF (PRESENT(fmat)) THEN
377 : ! recover fock matrix in CP2K/GPW order.
378 : ! Caution: Normalization is not take care of, so it's probably weird.
379 : CALL get_gto_basis_set(orb_basis_set, &
380 130 : nset=nset, nshell=nshell, l=ls, nsgf=nsgf, first_sgf=first_sgf)
381 130 : set_index = 0
382 130 : shell_index = 0
383 130 : nbb = 0
384 262 : DO i = 1, nset
385 716 : DO j = 1, nshell(i)
386 454 : l = ls(j, i)
387 586 : IF (l <= lmat) THEN
388 454 : nbb(l) = nbb(l) + 1
389 454 : k = nbb(l)
390 454 : CPASSERT(k <= 100)
391 454 : set_index(l, k) = i
392 454 : shell_index(l, k) = j
393 : END IF
394 : END DO
395 : END DO
396 130 : IF (uks) CPABORT("calculate_atomic_orbitals: only RKS is implemented")
397 130 : IF (ASSOCIATED(fmat)) CPABORT("fmat already associated")
398 130 : IF (.NOT. ASSOCIATED(atom%fmat)) CPABORT("atom%fmat not associated")
399 520 : ALLOCATE (fmat(nsgf, nsgf, 1))
400 9276 : fmat = 0.0_dp
401 260 : IF (.NOT. ghost) THEN
402 910 : DO l = 0, lmat
403 780 : ll = 2*l
404 1364 : DO k1 = 1, atom%basis%nbas(l)
405 2016 : DO k2 = 1, atom%basis%nbas(l)
406 782 : scal = SQRT(atom%integrals%ovlp(k1, k1, l)*atom%integrals%ovlp(k2, k2, l))
407 782 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
408 782 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
409 2614 : DO m = 0, ll
410 2160 : fmat(i + m, j + m, 1) = atom%fmat%op(k1, k2, l)/scal
411 : END DO
412 : END DO
413 : END DO
414 : END DO
415 : END IF
416 : END IF
417 :
418 8356 : nr = basis%grid%nr
419 :
420 8356 : IF (PRESENT(density)) THEN
421 2 : IF (ASSOCIATED(density)) DEALLOCATE (density)
422 6 : ALLOCATE (density(nr))
423 2 : IF (ghost) THEN
424 0 : density = 0.0_dp
425 : ELSE
426 2 : CALL atom_density(density, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ)
427 : END IF
428 : END IF
429 :
430 8356 : IF (PRESENT(wavefunction)) THEN
431 2 : CPASSERT(PRESENT(wfninfo))
432 2 : IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
433 2 : IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
434 14 : mo = SUM(atom%state%maxn_occ)
435 12 : ALLOCATE (wavefunction(nr, mo), wfninfo(2, mo))
436 3722 : wavefunction = 0.0_dp
437 2 : IF (.NOT. ghost) THEN
438 : ii = 0
439 14 : DO l = 0, lmat
440 18 : DO i = 1, atom%state%maxn_occ(l)
441 16 : IF (atom%state%occupation(l, i) > 0.0_dp) THEN
442 4 : ii = ii + 1
443 4 : wfninfo(1, ii) = atom%state%occupation(l, i)
444 4 : wfninfo(2, ii) = REAL(l, dp)
445 84 : DO j = 1, atom%basis%nbas(l)
446 : wavefunction(:, ii) = wavefunction(:, ii) + &
447 148724 : atom%orbitals%wfn(j, i, l)*basis%bf(:, j, l)
448 : END DO
449 : END IF
450 : END DO
451 : END DO
452 2 : CPASSERT(mo == ii)
453 : END IF
454 : END IF
455 :
456 : ! clean up
457 8356 : CALL atom_int_release(integrals)
458 8356 : CALL atom_ppint_release(integrals)
459 8356 : CALL atom_relint_release(integrals)
460 8356 : CALL release_atom_basis(basis)
461 8356 : CALL release_atom_potential(potential)
462 8356 : CALL release_atom_type(atom)
463 :
464 8356 : DEALLOCATE (potential, basis, integrals)
465 :
466 8356 : END SUBROUTINE calculate_atomic_orbitals
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param density ...
471 : !> \param atomic_kind ...
472 : !> \param qs_kind ...
473 : !> \param ngto ...
474 : !> \param iunit ...
475 : !> \param optbasis ... Default=T, if basis should be optimized, if not basis is given in input (density)
476 : !> \param allelectron ...
477 : !> \param confine ...
478 : ! **************************************************************************************************
479 90 : SUBROUTINE calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, &
480 : optbasis, allelectron, confine)
481 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: density
482 : TYPE(atomic_kind_type), POINTER :: atomic_kind
483 : TYPE(qs_kind_type), POINTER :: qs_kind
484 : INTEGER, INTENT(IN) :: ngto
485 : INTEGER, INTENT(IN), OPTIONAL :: iunit
486 : LOGICAL, INTENT(IN), OPTIONAL :: optbasis, allelectron, confine
487 :
488 : INTEGER, PARAMETER :: num_gto = 40
489 :
490 : INTEGER :: i, ii, iw, k, l, ll, m, mb, mo, ngp, nn, &
491 : nr, quadtype, relativistic, z
492 : INTEGER, DIMENSION(0:lmat) :: starti
493 : INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
494 90 : INTEGER, DIMENSION(:), POINTER :: econf
495 : LOGICAL :: do_basopt, ecp_semi_local
496 : REAL(KIND=dp) :: al, aval, cc, cval, ear, rk, xx, zeff
497 : REAL(KIND=dp), DIMENSION(num_gto+2) :: results
498 : TYPE(all_potential_type), POINTER :: all_potential
499 : TYPE(atom_basis_type), POINTER :: basis
500 : TYPE(atom_integrals), POINTER :: integrals
501 : TYPE(atom_orbitals), POINTER :: orbitals
502 : TYPE(atom_potential_type), POINTER :: potential
503 : TYPE(atom_type), POINTER :: atom
504 : TYPE(grid_atom_type), POINTER :: grid
505 : TYPE(gth_potential_type), POINTER :: gth_potential
506 : TYPE(sgp_potential_type), POINTER :: sgp_potential
507 :
508 90 : NULLIFY (atom)
509 90 : CALL create_atom_type(atom)
510 :
511 90 : CALL get_atomic_kind(atomic_kind, z=z)
512 90 : NULLIFY (all_potential, gth_potential)
513 : CALL get_qs_kind(qs_kind, zeff=zeff, &
514 : all_potential=all_potential, &
515 : gth_potential=gth_potential, &
516 90 : sgp_potential=sgp_potential)
517 :
518 90 : IF (PRESENT(iunit)) THEN
519 8 : iw = iunit
520 : ELSE
521 82 : iw = -1
522 : END IF
523 :
524 90 : IF (PRESENT(allelectron)) THEN
525 4 : IF (allelectron) THEN
526 4 : NULLIFY (gth_potential)
527 4 : zeff = z
528 : END IF
529 : END IF
530 :
531 90 : do_basopt = .TRUE.
532 90 : IF (PRESENT(optbasis)) THEN
533 46 : do_basopt = optbasis
534 : END IF
535 :
536 90 : CPASSERT(ngto <= num_gto)
537 :
538 90 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
539 : ! PP calculation are non-relativistic
540 84 : relativistic = do_nonrel_atom
541 : ELSE
542 : ! AE calculations use DKH2
543 6 : relativistic = do_dkh2_atom
544 : END IF
545 :
546 90 : atom%z = z
547 : CALL set_atom(atom, &
548 : pp_calc=(ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)), &
549 : method_type=do_rks_atom, &
550 : relativistic=relativistic, &
551 : coulomb_integral_type=do_numeric, &
552 96 : exchange_integral_type=do_numeric)
553 :
554 493470 : ALLOCATE (potential, basis, integrals)
555 :
556 90 : IF (PRESENT(confine)) THEN
557 90 : potential%confinement = confine
558 : ELSE
559 0 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
560 0 : potential%confinement = .TRUE.
561 : ELSE
562 0 : potential%confinement = .FALSE.
563 : END IF
564 : END IF
565 90 : potential%conf_type = barrier_conf
566 90 : potential%acon = 200._dp
567 90 : potential%rcon = 4.0_dp
568 90 : potential%scon = 8.0_dp
569 :
570 90 : IF (ASSOCIATED(gth_potential)) THEN
571 84 : potential%ppot_type = gth_pseudo
572 84 : CALL get_potential(gth_potential, zeff=zeff)
573 84 : CALL gth_potential_conversion(gth_potential, potential%gth_pot)
574 84 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
575 6 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
576 0 : CALL get_potential(sgp_potential, ecp_semi_local=ecp_semi_local)
577 0 : IF (ecp_semi_local) THEN
578 0 : potential%ppot_type = ecp_pseudo
579 0 : CALL ecp_potential_conversion(sgp_potential, potential%ecp_pot)
580 0 : potential%ecp_pot%symbol = ptable(z)%symbol
581 : ELSE
582 0 : potential%ppot_type = sgp_pseudo
583 0 : CALL sgp_potential_conversion(sgp_potential, potential%sgp_pot)
584 0 : potential%sgp_pot%symbol = ptable(z)%symbol
585 : END IF
586 0 : CALL get_potential(sgp_potential, zeff=zeff)
587 0 : CALL set_atom(atom, zcore=NINT(zeff), potential=potential)
588 : ELSE
589 6 : potential%ppot_type = no_pseudo
590 6 : CALL set_atom(atom, zcore=z, potential=potential)
591 : END IF
592 :
593 : ! atomic grid
594 90 : NULLIFY (grid)
595 90 : ngp = 400
596 90 : quadtype = do_gapw_log
597 90 : CALL allocate_grid_atom(grid)
598 90 : CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
599 90 : grid%nr = ngp
600 90 : basis%grid => grid
601 :
602 90 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
603 :
604 : ! fill in the basis data structures
605 90 : basis%eps_eig = 1.e-12_dp
606 90 : basis%basis_type = GTO_BASIS
607 90 : CALL Clementi_geobas(z, cval, aval, basis%nbas, starti)
608 630 : basis%nprim = basis%nbas
609 630 : m = MAXVAL(basis%nbas)
610 270 : ALLOCATE (basis%am(m, 0:lmat))
611 12222 : basis%am = 0._dp
612 630 : DO l = 0, lmat
613 3266 : DO i = 1, basis%nbas(l)
614 2636 : ll = i - 1 + starti(l)
615 3176 : basis%am(i, l) = aval*cval**(ll)
616 : END DO
617 : END DO
618 :
619 90 : basis%geometrical = .TRUE.
620 90 : basis%aval = aval
621 90 : basis%cval = cval
622 630 : basis%start = starti
623 :
624 : ! initialize basis function on a radial grid
625 90 : nr = basis%grid%nr
626 630 : m = MAXVAL(basis%nbas)
627 450 : ALLOCATE (basis%bf(nr, m, 0:lmat))
628 270 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
629 270 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
630 4649022 : basis%bf = 0._dp
631 4649022 : basis%dbf = 0._dp
632 4649022 : basis%ddbf = 0._dp
633 630 : DO l = 0, lmat
634 3266 : DO i = 1, basis%nbas(l)
635 2636 : al = basis%am(i, l)
636 1057576 : DO k = 1, nr
637 1054400 : rk = basis%grid%rad(k)
638 1054400 : ear = EXP(-al*basis%grid%rad(k)**2)
639 1054400 : basis%bf(k, i, l) = rk**l*ear
640 1054400 : basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
641 : basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
642 1057036 : 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
643 : END DO
644 : END DO
645 : END DO
646 :
647 90 : CALL set_atom(atom, basis=basis)
648 :
649 : ! optimization defaults
650 90 : atom%optimization%damping = 0.2_dp
651 90 : atom%optimization%eps_scf = 1.e-6_dp
652 90 : atom%optimization%eps_diis = 100._dp
653 90 : atom%optimization%max_iter = 50
654 90 : atom%optimization%n_diis = 5
655 :
656 90 : nelem = 0
657 90 : ncore = 0
658 90 : ncalc = 0
659 90 : IF (ASSOCIATED(gth_potential)) THEN
660 84 : CALL get_potential(gth_potential, elec_conf=econf)
661 84 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
662 6 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
663 0 : CALL get_potential(sgp_potential, elec_conf=econf)
664 0 : CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
665 : ELSE
666 30 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
667 24 : ll = 2*(2*l + 1)
668 24 : nn = ptable(z)%e_conv(l)
669 24 : ii = 0
670 6 : DO
671 24 : ii = ii + 1
672 24 : IF (nn <= ll) THEN
673 24 : nelem(l, ii) = nn
674 : EXIT
675 : ELSE
676 0 : nelem(l, ii) = ll
677 0 : nn = nn - ll
678 : END IF
679 : END DO
680 : END DO
681 426 : ncalc = nelem - ncore
682 : END IF
683 :
684 90 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
685 0 : nelem = 0
686 0 : ncore = 0
687 0 : ncalc = 0
688 : END IF
689 :
690 32670 : ALLOCATE (atom%state)
691 :
692 6390 : atom%state%core = 0._dp
693 4500 : atom%state%core(0:lmat, 1:7) = REAL(ncore(0:lmat, 1:7), dp)
694 6390 : atom%state%occ = 0._dp
695 4500 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp)
696 6390 : atom%state%occupation = 0._dp
697 90 : atom%state%multiplicity = -1
698 630 : DO l = 0, lmat
699 : k = 0
700 4410 : DO i = 1, 7
701 4320 : IF (ncalc(l, i) > 0) THEN
702 134 : k = k + 1
703 134 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp)
704 : END IF
705 : END DO
706 : END DO
707 :
708 90 : atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
709 630 : atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
710 90 : atom%state%maxl_calc = atom%state%maxl_occ
711 630 : atom%state%maxn_calc = atom%state%maxn_occ
712 :
713 : ! calculate integrals
714 : ! general integrals
715 : CALL atom_int_setup(integrals, basis, potential=atom%potential, &
716 : eri_coulomb=(atom%coulomb_integral_type == do_analytic), &
717 90 : eri_exchange=(atom%exchange_integral_type == do_analytic))
718 : ! potential
719 90 : CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
720 : ! relativistic correction terms
721 90 : NULLIFY (integrals%tzora, integrals%hdkh)
722 90 : CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=REAL(atom%zcore, dp))
723 90 : CALL set_atom(atom, integrals=integrals)
724 :
725 90 : NULLIFY (orbitals)
726 630 : mo = MAXVAL(atom%state%maxn_calc)
727 630 : mb = MAXVAL(atom%basis%nbas)
728 90 : CALL create_atom_orbs(orbitals, mb, mo)
729 90 : CALL set_atom(atom, orbitals=orbitals)
730 :
731 90 : CALL calculate_atom(atom, iw)
732 :
733 90 : IF (do_basopt) THEN
734 44 : CALL atom_fit_density(atom, ngto, 0, iw, results=results)
735 44 : xx = results(1)
736 44 : cc = results(2)
737 428 : DO i = 1, ngto
738 384 : density(i, 1) = xx*cc**i
739 428 : density(i, 2) = results(2 + i)
740 : END DO
741 : ELSE
742 46 : CALL atom_fit_density(atom, ngto, 0, iw, agto=density(:, 1), results=results)
743 352 : density(1:ngto, 2) = results(1:ngto)
744 : END IF
745 :
746 : ! clean up
747 90 : CALL atom_int_release(integrals)
748 90 : CALL atom_ppint_release(integrals)
749 90 : CALL atom_relint_release(integrals)
750 90 : CALL release_atom_basis(basis)
751 90 : CALL release_atom_potential(potential)
752 90 : CALL release_atom_type(atom)
753 :
754 90 : DEALLOCATE (potential, basis, integrals)
755 :
756 90 : END SUBROUTINE calculate_atomic_density
757 :
758 : ! **************************************************************************************************
759 : !> \brief ...
760 : !> \param atomic_kind ...
761 : !> \param qs_kind ...
762 : !> \param rel_control ...
763 : !> \param rtmat ...
764 : ! **************************************************************************************************
765 26 : SUBROUTINE calculate_atomic_relkin(atomic_kind, qs_kind, rel_control, rtmat)
766 : TYPE(atomic_kind_type), INTENT(IN) :: atomic_kind
767 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
768 : TYPE(rel_control_type), POINTER :: rel_control
769 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rtmat
770 :
771 : INTEGER :: i, ii, ipgf, j, k, k1, k2, l, ll, m, n, &
772 : ngp, nj, nn, nr, ns, nset, nsgf, &
773 : quadtype, relativistic, z
774 : INTEGER, DIMENSION(0:lmat, 10) :: ncalc, ncore, nelem
775 : INTEGER, DIMENSION(0:lmat, 100) :: set_index, shell_index
776 26 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
777 26 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, ls
778 : REAL(KIND=dp) :: al, alpha, ear, prefac, rk, zeff
779 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: omat
780 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
781 26 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
782 : TYPE(all_potential_type), POINTER :: all_potential
783 : TYPE(atom_basis_type), POINTER :: basis
784 : TYPE(atom_integrals), POINTER :: integrals
785 : TYPE(atom_potential_type), POINTER :: potential
786 : TYPE(atom_type), POINTER :: atom
787 : TYPE(grid_atom_type), POINTER :: grid
788 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
789 :
790 26 : IF (rel_control%rel_method == rel_none) RETURN
791 :
792 26 : NULLIFY (all_potential, orb_basis_set)
793 26 : CALL get_qs_kind(qs_kind, basis_set=orb_basis_set, all_potential=all_potential)
794 :
795 26 : CPASSERT(ASSOCIATED(orb_basis_set))
796 :
797 26 : IF (ASSOCIATED(all_potential)) THEN
798 : ! only all electron atoms will get the relativistic correction
799 :
800 26 : CALL get_atomic_kind(atomic_kind, z=z)
801 26 : CALL get_qs_kind(qs_kind, zeff=zeff)
802 26 : NULLIFY (atom)
803 26 : CALL create_atom_type(atom)
804 26 : NULLIFY (atom%xc_section)
805 26 : NULLIFY (atom%orbitals)
806 26 : atom%z = z
807 26 : alpha = SQRT(all_potential%alpha_core_charge)
808 :
809 : ! set the method flag
810 26 : SELECT CASE (rel_control%rel_method)
811 : CASE DEFAULT
812 0 : CPABORT("")
813 : CASE (rel_dkh)
814 26 : SELECT CASE (rel_control%rel_DKH_order)
815 : CASE DEFAULT
816 0 : CPABORT("")
817 : CASE (0)
818 0 : relativistic = do_dkh0_atom
819 : CASE (1)
820 0 : relativistic = do_dkh1_atom
821 : CASE (2)
822 8 : relativistic = do_dkh2_atom
823 : CASE (3)
824 14 : relativistic = do_dkh3_atom
825 : END SELECT
826 : CASE (rel_zora)
827 26 : SELECT CASE (rel_control%rel_zora_type)
828 : CASE DEFAULT
829 0 : CPABORT("")
830 : CASE (rel_zora_full)
831 0 : CPABORT("")
832 : CASE (rel_zora_mp)
833 0 : relativistic = do_zoramp_atom
834 : CASE (rel_sczora_mp)
835 12 : relativistic = do_sczoramp_atom
836 : END SELECT
837 : END SELECT
838 :
839 : CALL set_atom(atom, &
840 : pp_calc=.FALSE., &
841 : method_type=do_rks_atom, &
842 : relativistic=relativistic, &
843 : coulomb_integral_type=do_numeric, &
844 26 : exchange_integral_type=do_numeric)
845 :
846 142558 : ALLOCATE (potential, basis, integrals)
847 :
848 26 : potential%ppot_type = no_pseudo
849 26 : CALL set_atom(atom, zcore=z, potential=potential)
850 :
851 : CALL get_gto_basis_set(orb_basis_set, &
852 : nset=nset, nshell=nshell, npgf=npgf, lmin=lmin, lmax=lmax, l=ls, nsgf=nsgf, zet=zet, gcc=gcc, &
853 26 : first_sgf=first_sgf, last_sgf=last_sgf)
854 :
855 26 : NULLIFY (grid)
856 26 : ngp = 400
857 26 : quadtype = do_gapw_log
858 26 : CALL allocate_grid_atom(grid)
859 26 : CALL create_grid_atom(grid, ngp, 1, 1, 0, quadtype)
860 26 : grid%nr = ngp
861 26 : basis%grid => grid
862 :
863 26 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
864 26 : basis%basis_type = CGTO_BASIS
865 26 : basis%eps_eig = 1.e-12_dp
866 :
867 : ! fill in the basis data structures
868 26 : set_index = 0
869 26 : shell_index = 0
870 182 : basis%nprim = 0
871 182 : basis%nbas = 0
872 120 : DO i = 1, nset
873 188 : DO j = lmin(i), MIN(lmax(i), lmat)
874 188 : basis%nprim(j) = basis%nprim(j) + npgf(i)
875 : END DO
876 458 : DO j = 1, nshell(i)
877 338 : l = ls(j, i)
878 432 : IF (l <= lmat) THEN
879 338 : basis%nbas(l) = basis%nbas(l) + 1
880 338 : k = basis%nbas(l)
881 338 : CPASSERT(k <= 100)
882 338 : set_index(l, k) = i
883 338 : shell_index(l, k) = j
884 : END IF
885 : END DO
886 : END DO
887 :
888 182 : nj = MAXVAL(basis%nprim)
889 182 : ns = MAXVAL(basis%nbas)
890 78 : ALLOCATE (basis%am(nj, 0:lmat))
891 2150 : basis%am = 0._dp
892 130 : ALLOCATE (basis%cm(nj, ns, 0:lmat))
893 17810 : basis%cm = 0._dp
894 182 : DO j = 0, lmat
895 : nj = 0
896 : ns = 0
897 746 : DO i = 1, nset
898 720 : IF (j >= lmin(i) .AND. j <= lmax(i)) THEN
899 734 : DO ipgf = 1, npgf(i)
900 734 : basis%am(nj + ipgf, j) = zet(ipgf, i)
901 : END DO
902 432 : DO ii = 1, nshell(i)
903 432 : IF (ls(ii, i) == j) THEN
904 338 : ns = ns + 1
905 4146 : DO ipgf = 1, npgf(i)
906 4146 : basis%cm(nj + ipgf, ns, j) = gcc(ipgf, ii, i)
907 : END DO
908 : END IF
909 : END DO
910 94 : nj = nj + npgf(i)
911 : END IF
912 : END DO
913 : END DO
914 :
915 : ! Normalization as used in the atomic code
916 : ! We have to undo the Quickstep normalization
917 182 : DO j = 0, lmat
918 156 : prefac = 2.0_dp*SQRT(pi/dfac(2*j + 1))
919 822 : DO ipgf = 1, basis%nprim(j)
920 6092 : DO ii = 1, basis%nbas(j)
921 5936 : basis%cm(ipgf, ii, j) = prefac*basis%cm(ipgf, ii, j)
922 : END DO
923 : END DO
924 : END DO
925 :
926 : ! initialize basis function on a radial grid
927 26 : nr = basis%grid%nr
928 182 : m = MAXVAL(basis%nbas)
929 130 : ALLOCATE (basis%bf(nr, m, 0:lmat))
930 78 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
931 78 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
932 :
933 351458 : basis%bf = 0._dp
934 351458 : basis%dbf = 0._dp
935 351458 : basis%ddbf = 0._dp
936 182 : DO l = 0, lmat
937 822 : DO i = 1, basis%nprim(l)
938 640 : al = basis%am(i, l)
939 256796 : DO k = 1, nr
940 256000 : rk = basis%grid%rad(k)
941 256000 : ear = EXP(-al*basis%grid%rad(k)**2)
942 2375040 : DO j = 1, basis%nbas(l)
943 2118400 : basis%bf(k, j, l) = basis%bf(k, j, l) + rk**l*ear*basis%cm(i, j, l)
944 : basis%dbf(k, j, l) = basis%dbf(k, j, l) &
945 2118400 : + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis%cm(i, j, l)
946 : basis%ddbf(k, j, l) = basis%ddbf(k, j, l) + &
947 : (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)* &
948 2374400 : rk**(l) + 4._dp*al*rk**(l + 2))*ear*basis%cm(i, j, l)
949 : END DO
950 : END DO
951 : END DO
952 : END DO
953 :
954 26 : CALL set_atom(atom, basis=basis)
955 :
956 : ! optimization defaults
957 26 : atom%optimization%damping = 0.2_dp
958 26 : atom%optimization%eps_scf = 1.e-6_dp
959 26 : atom%optimization%eps_diis = 100._dp
960 26 : atom%optimization%max_iter = 50
961 26 : atom%optimization%n_diis = 5
962 :
963 : ! electronic state
964 26 : nelem = 0
965 26 : ncore = 0
966 26 : ncalc = 0
967 130 : DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
968 104 : ll = 2*(2*l + 1)
969 104 : nn = ptable(z)%e_conv(l)
970 104 : ii = 0
971 26 : DO
972 146 : ii = ii + 1
973 146 : IF (nn <= ll) THEN
974 104 : nelem(l, ii) = nn
975 : EXIT
976 : ELSE
977 42 : nelem(l, ii) = ll
978 42 : nn = nn - ll
979 : END IF
980 : END DO
981 : END DO
982 1846 : ncalc = nelem - ncore
983 :
984 26 : IF (qs_kind%ghost .OR. qs_kind%floating) THEN
985 : nelem = 0
986 0 : ncore = 0
987 0 : ncalc = 0
988 : END IF
989 :
990 9438 : ALLOCATE (atom%state)
991 :
992 1846 : atom%state%core = 0._dp
993 1300 : atom%state%core(0:lmat, 1:7) = REAL(ncore(0:lmat, 1:7), dp)
994 1846 : atom%state%occ = 0._dp
995 1300 : atom%state%occ(0:lmat, 1:7) = REAL(ncalc(0:lmat, 1:7), dp)
996 1846 : atom%state%occupation = 0._dp
997 26 : atom%state%multiplicity = -1
998 182 : DO l = 0, lmat
999 : k = 0
1000 1274 : DO i = 1, 7
1001 1248 : IF (ncalc(l, i) > 0) THEN
1002 88 : k = k + 1
1003 88 : atom%state%occupation(l, k) = REAL(ncalc(l, i), dp)
1004 : END IF
1005 : END DO
1006 : END DO
1007 :
1008 26 : atom%state%maxl_occ = get_maxl_occ(atom%state%occupation)
1009 182 : atom%state%maxn_occ = get_maxn_occ(atom%state%occupation)
1010 26 : atom%state%maxl_calc = atom%state%maxl_occ
1011 182 : atom%state%maxn_calc = atom%state%maxn_occ
1012 :
1013 : ! calculate integrals
1014 : ! general integrals
1015 26 : CALL atom_int_setup(integrals, basis)
1016 : ! potential
1017 26 : CALL atom_ppint_setup(integrals, basis, potential=atom%potential)
1018 : ! relativistic correction terms
1019 26 : NULLIFY (integrals%tzora, integrals%hdkh)
1020 : CALL atom_relint_setup(integrals, basis, atom%relativistic, zcore=REAL(atom%zcore, dp), &
1021 26 : alpha=alpha)
1022 26 : CALL set_atom(atom, integrals=integrals)
1023 :
1024 : ! for DKH we need erfc integrals to correct non-relativistic
1025 13742 : integrals%core = 0.0_dp
1026 182 : DO l = 0, lmat
1027 156 : n = integrals%n(l)
1028 156 : m = basis%nprim(l)
1029 452 : ALLOCATE (omat(m, m))
1030 :
1031 156 : CALL sg_erfc(omat(1:m, 1:m), l, alpha, basis%am(1:m, l), basis%am(1:m, l))
1032 156 : integrals%core(1:n, 1:n, l) = MATMUL(TRANSPOSE(basis%cm(1:m, 1:n, l)), &
1033 635662 : MATMUL(omat(1:m, 1:m), basis%cm(1:m, 1:n, l)))
1034 :
1035 182 : DEALLOCATE (omat)
1036 : END DO
1037 :
1038 : ! recover relativistic kinetic matrix in CP2K/GPW order and normalization
1039 26 : IF (ASSOCIATED(rtmat)) THEN
1040 0 : DEALLOCATE (rtmat)
1041 : END IF
1042 104 : ALLOCATE (rtmat(nsgf, nsgf))
1043 159614 : rtmat = 0._dp
1044 182 : DO l = 0, lmat
1045 156 : ll = 2*l
1046 520 : DO k1 = 1, basis%nbas(l)
1047 4792 : DO k2 = 1, basis%nbas(l)
1048 4298 : i = first_sgf(shell_index(l, k1), set_index(l, k1))
1049 4298 : j = first_sgf(shell_index(l, k2), set_index(l, k2))
1050 338 : SELECT CASE (atom%relativistic)
1051 : CASE DEFAULT
1052 0 : CPABORT("")
1053 : CASE (do_zoramp_atom, do_sczoramp_atom)
1054 14656 : DO m = 0, ll
1055 14656 : rtmat(i + m, j + m) = integrals%tzora(k1, k2, l)
1056 : END DO
1057 : CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
1058 4898 : DO m = 0, ll
1059 : rtmat(i + m, j + m) = integrals%hdkh(k1, k2, l) - integrals%kin(k1, k2, l) + &
1060 904 : atom%zcore*integrals%core(k1, k2, l)
1061 : END DO
1062 : END SELECT
1063 : END DO
1064 : END DO
1065 : END DO
1066 968 : DO k1 = 1, nsgf
1067 80762 : DO k2 = k1, nsgf
1068 79794 : rtmat(k1, k2) = 0.5_dp*(rtmat(k1, k2) + rtmat(k2, k1))
1069 80736 : rtmat(k2, k1) = rtmat(k1, k2)
1070 : END DO
1071 : END DO
1072 :
1073 : ! clean up
1074 26 : CALL atom_int_release(integrals)
1075 26 : CALL atom_ppint_release(integrals)
1076 26 : CALL atom_relint_release(integrals)
1077 26 : CALL release_atom_basis(basis)
1078 26 : CALL release_atom_potential(potential)
1079 26 : CALL release_atom_type(atom)
1080 :
1081 78 : DEALLOCATE (potential, basis, integrals)
1082 :
1083 : ELSE
1084 :
1085 0 : IF (ASSOCIATED(rtmat)) THEN
1086 0 : DEALLOCATE (rtmat)
1087 : END IF
1088 0 : NULLIFY (rtmat)
1089 :
1090 : END IF
1091 :
1092 52 : END SUBROUTINE calculate_atomic_relkin
1093 :
1094 : ! **************************************************************************************************
1095 : !> \brief ...
1096 : !> \param gth_potential ...
1097 : !> \param gth_atompot ...
1098 : ! **************************************************************************************************
1099 14664 : SUBROUTINE gth_potential_conversion(gth_potential, gth_atompot)
1100 : TYPE(gth_potential_type), POINTER :: gth_potential
1101 : TYPE(atom_gthpot_type) :: gth_atompot
1102 :
1103 : INTEGER :: i, j, l, lm, n, ne, nexp_lpot, nexp_lsd, &
1104 : nexp_nlcc
1105 7332 : INTEGER, DIMENSION(:), POINTER :: nct_lpot, nct_lsd, nct_nlcc, nppnl, &
1106 7332 : ppeconf
1107 : LOGICAL :: lpot_present, lsd_present, nlcc_present
1108 : REAL(KIND=dp) :: ac, zeff
1109 7332 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_lpot, alpha_lsd, alpha_nlcc, ap, ce
1110 7332 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_lpot, cval_lsd, cval_nlcc
1111 7332 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: hp
1112 :
1113 : CALL get_potential(gth_potential, &
1114 : zeff=zeff, &
1115 : elec_conf=ppeconf, &
1116 : alpha_core_charge=ac, &
1117 : nexp_ppl=ne, &
1118 : cexp_ppl=ce, &
1119 : lppnl=lm, &
1120 : nprj_ppnl=nppnl, &
1121 : alpha_ppnl=ap, &
1122 7332 : hprj_ppnl=hp)
1123 :
1124 7332 : gth_atompot%zion = zeff
1125 7332 : gth_atompot%rc = SQRT(0.5_dp/ac)
1126 7332 : gth_atompot%ncl = ne
1127 43992 : gth_atompot%cl(:) = 0._dp
1128 7332 : IF (ac > 0._dp) THEN
1129 21608 : DO i = 1, ne
1130 21608 : gth_atompot%cl(i) = ce(i)/(2._dp*ac)**(i - 1)
1131 : END DO
1132 : END IF
1133 : !extended type
1134 7332 : gth_atompot%lpotextended = .FALSE.
1135 7332 : gth_atompot%lsdpot = .FALSE.
1136 7332 : gth_atompot%nlcc = .FALSE.
1137 7332 : gth_atompot%nexp_lpot = 0
1138 7332 : gth_atompot%nexp_lsd = 0
1139 7332 : gth_atompot%nexp_nlcc = 0
1140 : CALL get_potential(gth_potential, &
1141 : lpot_present=lpot_present, &
1142 : lsd_present=lsd_present, &
1143 7332 : nlcc_present=nlcc_present)
1144 7332 : IF (lpot_present) THEN
1145 : CALL get_potential(gth_potential, &
1146 : nexp_lpot=nexp_lpot, &
1147 : alpha_lpot=alpha_lpot, &
1148 : nct_lpot=nct_lpot, &
1149 8 : cval_lpot=cval_lpot)
1150 8 : gth_atompot%lpotextended = .TRUE.
1151 8 : gth_atompot%nexp_lpot = nexp_lpot
1152 20 : gth_atompot%alpha_lpot(1:nexp_lpot) = SQRT(0.5_dp/alpha_lpot(1:nexp_lpot))
1153 20 : gth_atompot%nct_lpot(1:nexp_lpot) = nct_lpot(1:nexp_lpot)
1154 20 : DO j = 1, nexp_lpot
1155 12 : ac = alpha_lpot(j)
1156 68 : DO i = 1, 4
1157 60 : gth_atompot%cval_lpot(i, j) = cval_lpot(i, j)/(2._dp*ac)**(i - 1)
1158 : END DO
1159 : END DO
1160 : END IF
1161 7332 : IF (lsd_present) THEN
1162 : CALL get_potential(gth_potential, &
1163 : nexp_lsd=nexp_lsd, &
1164 : alpha_lsd=alpha_lsd, &
1165 : nct_lsd=nct_lsd, &
1166 0 : cval_lsd=cval_lsd)
1167 0 : gth_atompot%lsdpot = .TRUE.
1168 0 : gth_atompot%nexp_lsd = nexp_lsd
1169 0 : gth_atompot%alpha_lsd(1:nexp_lsd) = SQRT(0.5_dp/alpha_lsd(1:nexp_lsd))
1170 0 : gth_atompot%nct_lsd(1:nexp_lsd) = nct_lsd(1:nexp_lsd)
1171 0 : DO j = 1, nexp_lpot
1172 0 : ac = alpha_lsd(j)
1173 0 : DO i = 1, 4
1174 0 : gth_atompot%cval_lsd(i, j) = cval_lsd(i, j)/(2._dp*ac)**(i - 1)
1175 : END DO
1176 : END DO
1177 : END IF
1178 :
1179 : ! nonlocal part
1180 51324 : gth_atompot%nl(:) = 0
1181 51324 : gth_atompot%rcnl(:) = 0._dp
1182 931164 : gth_atompot%hnl(:, :, :) = 0._dp
1183 14912 : DO l = 0, lm
1184 7580 : n = nppnl(l)
1185 7580 : gth_atompot%nl(l) = n
1186 7580 : gth_atompot%rcnl(l) = SQRT(0.5_dp/ap(l))
1187 27058 : gth_atompot%hnl(1:n, 1:n, l) = hp(1:n, 1:n, l)
1188 : END DO
1189 :
1190 7332 : IF (nlcc_present) THEN
1191 : CALL get_potential(gth_potential, &
1192 : nexp_nlcc=nexp_nlcc, &
1193 : alpha_nlcc=alpha_nlcc, &
1194 : nct_nlcc=nct_nlcc, &
1195 18 : cval_nlcc=cval_nlcc)
1196 18 : gth_atompot%nlcc = .TRUE.
1197 18 : gth_atompot%nexp_nlcc = nexp_nlcc
1198 36 : gth_atompot%alpha_nlcc(1:nexp_nlcc) = alpha_nlcc(1:nexp_nlcc)
1199 36 : gth_atompot%nct_nlcc(1:nexp_nlcc) = nct_nlcc(1:nexp_nlcc)
1200 108 : gth_atompot%cval_nlcc(1:4, 1:nexp_nlcc) = cval_nlcc(1:4, 1:nexp_nlcc)
1201 : END IF
1202 :
1203 7332 : END SUBROUTINE gth_potential_conversion
1204 :
1205 : ! **************************************************************************************************
1206 : !> \brief ...
1207 : !> \param sgp_potential ...
1208 : !> \param sgp_atompot ...
1209 : ! **************************************************************************************************
1210 36 : SUBROUTINE sgp_potential_conversion(sgp_potential, sgp_atompot)
1211 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1212 : TYPE(atom_sgppot_type) :: sgp_atompot
1213 :
1214 : INTEGER :: lm, n
1215 12 : INTEGER, DIMENSION(:), POINTER :: ppeconf
1216 : LOGICAL :: nlcc_present
1217 : REAL(KIND=dp) :: ac, zeff
1218 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: ap, ce
1219 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hhp
1220 12 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ccp
1221 :
1222 : CALL get_potential(sgp_potential, &
1223 : name=sgp_atompot%pname, &
1224 : zeff=zeff, &
1225 : elec_conf=ppeconf, &
1226 12 : alpha_core_charge=ac)
1227 12 : sgp_atompot%zion = zeff
1228 12 : sgp_atompot%ac_local = ac
1229 60 : sgp_atompot%econf(0:3) = ppeconf(0:3)
1230 : CALL get_potential(sgp_potential, lmax=lm, &
1231 : is_nonlocal=sgp_atompot%is_nonlocal, &
1232 12 : n_nonlocal=n, a_nonlocal=ap, h_nonlocal=hhp, c_nonlocal=ccp)
1233 : ! nonlocal
1234 48 : sgp_atompot%has_nonlocal = ANY(sgp_atompot%is_nonlocal)
1235 12 : sgp_atompot%lmax = lm
1236 12 : IF (sgp_atompot%has_nonlocal) THEN
1237 6 : CPASSERT(n <= SIZE(sgp_atompot%a_nonlocal))
1238 6 : sgp_atompot%n_nonlocal = n
1239 54 : sgp_atompot%a_nonlocal(1:n) = ap(1:n)
1240 60 : sgp_atompot%h_nonlocal(1:n, 0:lm) = hhp(1:n, 0:lm)
1241 444 : sgp_atompot%c_nonlocal(1:n, 1:n, 0:lm) = ccp(1:n, 1:n, 0:lm)
1242 : END IF
1243 : ! local
1244 12 : CALL get_potential(sgp_potential, n_local=n, a_local=ap, c_local=ce)
1245 12 : CPASSERT(n <= SIZE(sgp_atompot%a_local))
1246 12 : sgp_atompot%n_local = n
1247 156 : sgp_atompot%a_local(1:n) = ap(1:n)
1248 156 : sgp_atompot%c_local(1:n) = ce(1:n)
1249 : ! NLCC
1250 : CALL get_potential(sgp_potential, has_nlcc=nlcc_present, &
1251 12 : n_nlcc=n, a_nlcc=ap, c_nlcc=ce)
1252 12 : IF (nlcc_present) THEN
1253 0 : sgp_atompot%has_nlcc = .TRUE.
1254 0 : CPASSERT(n <= SIZE(sgp_atompot%a_nlcc))
1255 0 : sgp_atompot%n_nlcc = n
1256 0 : sgp_atompot%a_nlcc(1:n) = ap(1:n)
1257 0 : sgp_atompot%c_nlcc(1:n) = ce(1:n)
1258 : ELSE
1259 12 : sgp_atompot%has_nlcc = .FALSE.
1260 : END IF
1261 :
1262 12 : END SUBROUTINE sgp_potential_conversion
1263 :
1264 : ! **************************************************************************************************
1265 : !> \brief ...
1266 : !> \param sgp_potential ...
1267 : !> \param ecp_atompot ...
1268 : ! **************************************************************************************************
1269 24 : SUBROUTINE ecp_potential_conversion(sgp_potential, ecp_atompot)
1270 : TYPE(sgp_potential_type), POINTER :: sgp_potential
1271 : TYPE(atom_ecppot_type) :: ecp_atompot
1272 :
1273 12 : INTEGER, DIMENSION(:), POINTER :: ppeconf
1274 : LOGICAL :: ecp_local, ecp_semi_local
1275 : REAL(KIND=dp) :: zeff
1276 :
1277 12 : CALL get_potential(sgp_potential, ecp_local=ecp_local, ecp_semi_local=ecp_semi_local)
1278 12 : CPASSERT(ecp_semi_local .AND. ecp_local)
1279 : CALL get_potential(sgp_potential, &
1280 : name=ecp_atompot%pname, &
1281 : zeff=zeff, &
1282 12 : elec_conf=ppeconf)
1283 12 : ecp_atompot%zion = zeff
1284 60 : ecp_atompot%econf(0:3) = ppeconf(0:3)
1285 12 : CALL get_potential(sgp_potential, lmax=ecp_atompot%lmax)
1286 : ! local
1287 : CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, &
1288 12 : aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc)
1289 : ! nonlocal
1290 : CALL get_potential(sgp_potential, npot=ecp_atompot%npot, nrpot=ecp_atompot%nrpot, &
1291 12 : apot=ecp_atompot%apot, bpot=ecp_atompot%bpot)
1292 :
1293 12 : END SUBROUTINE ecp_potential_conversion
1294 : ! **************************************************************************************************
1295 :
1296 156 : END MODULE atom_kind_orbitals
|
