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