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