(git:8ebf9ad)
Loading...
Searching...
No Matches
lri_environment_init.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief initializes the environment for lri
10!> lri : local resolution of the identity
11!> \par History
12!> created [06.2015]
13!> \author Dorothea Golze
14! **************************************************************************************************
15MODULE lri_environment_init
16 USE ao_util, ONLY: exp_radius
17 USE atomic_kind_types, ONLY: atomic_kind_type,&
18 get_atomic_kind_set
19 USE basis_set_types, ONLY: copy_gto_basis_set,&
20 gto_basis_set_type
21 USE bibliography, ONLY: golze2017a,&
22 golze2017b,&
23 cite_reference
24 USE cp_control_types, ONLY: dft_control_type
25 USE generic_shg_integrals, ONLY: int_overlap_aba_shg
26 USE generic_shg_integrals_init, ONLY: contraction_matrix_shg,&
27 contraction_matrix_shg_mix,&
28 get_clebsch_gordon_coefficients
29 USE input_section_types, ONLY: section_vals_type,&
30 section_vals_val_get
31 USE kinds, ONLY: dp
32 USE lri_environment_types, ONLY: deallocate_bas_properties,&
33 lri_env_create,&
34 lri_environment_type
35 USE mathconstants, ONLY: fac,&
36 pi,&
37 rootpi
38 USE mathlib, ONLY: invert_matrix
39 USE qs_environment_types, ONLY: get_qs_env,&
40 qs_environment_type
41 USE qs_kind_types, ONLY: get_qs_kind,&
42 qs_kind_type
43#include "./base/base_uses.f90"
44
45 IMPLICIT NONE
46
47 PRIVATE
48
49! **************************************************************************************************
50
51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_environment_init'
52
53 PUBLIC :: lri_env_init, lri_env_basis, lri_basis_init
54
55! **************************************************************************************************
56
57CONTAINS
58
59! **************************************************************************************************
60!> \brief initializes the lri env
61!> \param lri_env ...
62!> \param lri_section ...
63! **************************************************************************************************
64 SUBROUTINE lri_env_init(lri_env, lri_section)
65
66 TYPE(lri_environment_type), POINTER :: lri_env
67 TYPE(section_vals_type), POINTER :: lri_section
68
69 REAL(kind=dp), DIMENSION(:), POINTER :: radii
70
71 NULLIFY (lri_env)
72 ALLOCATE (lri_env)
73 CALL lri_env_create(lri_env)
74
75 ! init keywords
76 ! use RI for local pp terms
77 CALL section_vals_val_get(lri_section, "RI_STATISTIC", &
78 l_val=lri_env%statistics)
79 ! use exact one-center terms
80 CALL section_vals_val_get(lri_section, "EXACT_1C_TERMS", &
81 l_val=lri_env%exact_1c_terms)
82 ! use RI for local pp terms
83 CALL section_vals_val_get(lri_section, "PPL_RI", &
84 l_val=lri_env%ppl_ri)
85 ! check for debug (OS scheme)
86 CALL section_vals_val_get(lri_section, "DEBUG_LRI_INTEGRALS", &
87 l_val=lri_env%debug)
88 ! integrals based on solid harmonic Gaussians
89 CALL section_vals_val_get(lri_section, "SHG_LRI_INTEGRALS", &
90 l_val=lri_env%use_shg_integrals)
91 ! how to calculate inverse/pseuodinverse of overlap
92 CALL section_vals_val_get(lri_section, "LRI_OVERLAP_MATRIX", &
93 i_val=lri_env%lri_overlap_inv)
94 CALL section_vals_val_get(lri_section, "MAX_CONDITION_NUM", &
95 r_val=lri_env%cond_max)
96 ! integrals threshold (aba, abb)
97 CALL section_vals_val_get(lri_section, "EPS_O3_INT", &
98 r_val=lri_env%eps_o3_int)
99 ! RI SINV
100 CALL section_vals_val_get(lri_section, "RI_SINV", &
101 c_val=lri_env%ri_sinv_app)
102 ! Distant Pair Approximation
103 CALL section_vals_val_get(lri_section, "DISTANT_PAIR_APPROXIMATION", &
104 l_val=lri_env%distant_pair_approximation)
105 CALL section_vals_val_get(lri_section, "DISTANT_PAIR_METHOD", &
106 c_val=lri_env%distant_pair_method)
107 CALL section_vals_val_get(lri_section, "DISTANT_PAIR_RADII", r_vals=radii)
108 cpassert(SIZE(radii) == 2)
109 cpassert(radii(2) > radii(1))
110 cpassert(radii(1) > 0.0_dp)
111 lri_env%r_in = radii(1)
112 lri_env%r_out = radii(2)
113
114 CALL cite_reference(golze2017b)
115 IF (lri_env%use_shg_integrals) CALL cite_reference(golze2017a)
116
117 END SUBROUTINE lri_env_init
118! **************************************************************************************************
119!> \brief initializes the lri env
120!> \param ri_type ...
121!> \param qs_env ...
122!> \param lri_env ...
123!> \param qs_kind_set ...
124! **************************************************************************************************
125 SUBROUTINE lri_env_basis(ri_type, qs_env, lri_env, qs_kind_set)
126
127 CHARACTER(len=*), INTENT(IN) :: ri_type
128 TYPE(qs_environment_type), POINTER :: qs_env
129 TYPE(lri_environment_type), POINTER :: lri_env
130 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
131
132 INTEGER :: i, i1, i2, iat, ikind, ip, ipgf, iset, ishell, jp, l, lmax_ikind_orb, &
133 lmax_ikind_ri, maxl_orb, maxl_ri, n1, n2, natom, nbas, nkind, nribas, nspin
134 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
135 REAL(kind=dp) :: gcc, rad, rai, raj, xradius, zeta
136 REAL(kind=dp), DIMENSION(:), POINTER :: int_aux, norm
137 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
138 TYPE(dft_control_type), POINTER :: dft_control
139 TYPE(gto_basis_set_type), POINTER :: orb_basis_set, ri_basis_set
140
141 ! initialize the basic basis sets (orb and ri)
142 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
143 nkind = SIZE(atomic_kind_set)
144 ALLOCATE (lri_env%orb_basis(nkind), lri_env%ri_basis(nkind))
145 maxl_orb = 0
146 maxl_ri = 0
147 DO ikind = 1, nkind
148 NULLIFY (orb_basis_set, ri_basis_set)
149 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
150 IF (ri_type == "LRI") THEN
151 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="LRI_AUX")
152 ELSE IF (ri_type == "P_LRI") THEN
153 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="P_LRI_AUX")
154 ELSE IF (ri_type == "RI") THEN
155 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC")
156 ELSE
157 cpabort('ri_type')
158 END IF
159 NULLIFY (lri_env%orb_basis(ikind)%gto_basis_set)
160 NULLIFY (lri_env%ri_basis(ikind)%gto_basis_set)
161 IF (ASSOCIATED(orb_basis_set)) THEN
162 CALL copy_gto_basis_set(orb_basis_set, lri_env%orb_basis(ikind)%gto_basis_set)
163 CALL copy_gto_basis_set(ri_basis_set, lri_env%ri_basis(ikind)%gto_basis_set)
164 END IF
165 lmax_ikind_orb = maxval(lri_env%orb_basis(ikind)%gto_basis_set%lmax)
166 lmax_ikind_ri = maxval(lri_env%ri_basis(ikind)%gto_basis_set%lmax)
167 maxl_orb = max(maxl_orb, lmax_ikind_orb)
168 maxl_ri = max(maxl_ri, lmax_ikind_ri)
169 END DO
170
171 IF ((ri_type == "LRI") .OR. (ri_type == "P_LRI")) THEN
172 ! CG coefficients needed for lri integrals
173 IF (ASSOCIATED(lri_env%cg_shg)) THEN
174 CALL get_clebsch_gordon_coefficients(lri_env%cg_shg%cg_coeff, &
175 lri_env%cg_shg%cg_none0_list, &
176 lri_env%cg_shg%ncg_none0, &
177 maxl_orb, maxl_ri)
178 END IF
179 CALL lri_basis_init(lri_env)
180 ! distant pair approximation
181 IF (lri_env%distant_pair_approximation) THEN
182 !
183 SELECT CASE (lri_env%distant_pair_method)
184 CASE ("EW")
185 ! equal weight of 0.5
186 CASE ("AW")
187 ALLOCATE (lri_env%aradius(nkind))
188 DO i = 1, nkind
189 orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
190 lri_env%aradius(i) = orb_basis_set%kind_radius
191 END DO
192 CASE ("SW")
193 ALLOCATE (lri_env%wbas(nkind))
194 DO i = 1, nkind
195 orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
196 n1 = orb_basis_set%nsgf
197 ALLOCATE (lri_env%wbas(i)%vec(n1))
198 DO iset = 1, orb_basis_set%nset
199 i1 = orb_basis_set%first_sgf(1, iset)
200 n2 = orb_basis_set%nshell(iset)
201 i2 = orb_basis_set%last_sgf(n2, iset)
202 lri_env%wbas(i)%vec(i1:i2) = orb_basis_set%set_radius(iset)
203 END DO
204 END DO
205 CASE ("LW")
206 ALLOCATE (lri_env%wbas(nkind))
207 DO i = 1, nkind
208 orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
209 n1 = orb_basis_set%nsgf
210 ALLOCATE (lri_env%wbas(i)%vec(n1))
211 DO iset = 1, orb_basis_set%nset
212 DO ishell = 1, orb_basis_set%nshell(iset)
213 i1 = orb_basis_set%first_sgf(ishell, iset)
214 i2 = orb_basis_set%last_sgf(ishell, iset)
215 l = orb_basis_set%l(ishell, iset)
216 xradius = 0.0_dp
217 DO ipgf = 1, orb_basis_set%npgf(iset)
218 gcc = orb_basis_set%gcc(ipgf, ishell, iset)
219 zeta = orb_basis_set%zet(ipgf, iset)
220 rad = exp_radius(l, zeta, 1.e-5_dp, gcc, rlow=xradius)
221 xradius = max(xradius, rad)
222 END DO
223 lri_env%wbas(i)%vec(i1:i2) = xradius
224 END DO
225 END DO
226 END DO
227 CASE DEFAULT
228 cpabort("Unknown DISTANT_PAIR_METHOD in LRI")
229 END SELECT
230 !
231 ALLOCATE (lri_env%wmat(nkind, nkind))
232 SELECT CASE (lri_env%distant_pair_method)
233 CASE ("EW")
234 ! equal weight of 0.5
235 DO i1 = 1, nkind
236 n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
237 DO i2 = 1, nkind
238 n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
239 ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
240 lri_env%wmat(i1, i2)%mat(:, :) = 0.5_dp
241 END DO
242 END DO
243 CASE ("AW")
244 DO i1 = 1, nkind
245 n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
246 DO i2 = 1, nkind
247 n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
248 ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
249 rai = lri_env%aradius(i1)**2
250 raj = lri_env%aradius(i2)**2
251 IF (raj > rai) THEN
252 lri_env%wmat(i1, i2)%mat(:, :) = 1.0_dp
253 ELSE
254 lri_env%wmat(i1, i2)%mat(:, :) = 0.0_dp
255 END IF
256 END DO
257 END DO
258 CASE ("SW", "LW")
259 DO i1 = 1, nkind
260 n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
261 DO i2 = 1, nkind
262 n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
263 ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
264 DO ip = 1, SIZE(lri_env%wbas(i1)%vec)
265 rai = lri_env%wbas(i1)%vec(ip)**2
266 DO jp = 1, SIZE(lri_env%wbas(i2)%vec)
267 raj = lri_env%wbas(i2)%vec(jp)**2
268 IF (raj > rai) THEN
269 lri_env%wmat(i1, i2)%mat(ip, jp) = 1.0_dp
270 ELSE
271 lri_env%wmat(i1, i2)%mat(ip, jp) = 0.0_dp
272 END IF
273 END DO
274 END DO
275 END DO
276 END DO
277 END SELECT
278 END IF
279 ELSE IF (ri_type == "RI") THEN
280 ALLOCATE (lri_env%ri_fit)
281 NULLIFY (lri_env%ri_fit%nvec)
282 NULLIFY (lri_env%ri_fit%bas_ptr)
283 CALL get_qs_env(qs_env=qs_env, natom=natom)
284 ! initialize pointers to RI basis vector
285 ALLOCATE (lri_env%ri_fit%bas_ptr(2, natom))
286 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
287 nbas = 0
288 DO iat = 1, natom
289 ikind = kind_of(iat)
290 nribas = lri_env%ri_basis(ikind)%gto_basis_set%nsgf
291 lri_env%ri_fit%bas_ptr(1, iat) = nbas + 1
292 lri_env%ri_fit%bas_ptr(2, iat) = nbas + nribas
293 nbas = nbas + nribas
294 END DO
295 ! initialize vector t
296 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
297 nspin = dft_control%nspins
298 ALLOCATE (lri_env%ri_fit%tvec(nbas, nspin), lri_env%ri_fit%rm1t(nbas, nspin))
299 ! initialize vector a, expansion of density
300 ALLOCATE (lri_env%ri_fit%avec(nbas, nspin))
301 ! initialize vector fout, R^(-1)*(f-p*n)
302 ALLOCATE (lri_env%ri_fit%fout(nbas, nspin))
303 ! initialize vector with RI basis integrated
304 NULLIFY (norm, int_aux)
305 nbas = lri_env%ri_fit%bas_ptr(2, natom)
306 ALLOCATE (lri_env%ri_fit%nvec(nbas), lri_env%ri_fit%rm1n(nbas))
307 ikind = 0
308 DO iat = 1, natom
309 IF (ikind /= kind_of(iat)) THEN
310 ikind = kind_of(iat)
311 ri_basis_set => lri_env%ri_basis(ikind)%gto_basis_set
312 IF (ASSOCIATED(norm)) DEALLOCATE (norm)
313 IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
314 CALL basis_norm_s_func(ri_basis_set, norm)
315 CALL basis_int(ri_basis_set, int_aux, norm)
316 END IF
317 nbas = SIZE(int_aux)
318 i1 = lri_env%ri_fit%bas_ptr(1, iat)
319 i2 = lri_env%ri_fit%bas_ptr(2, iat)
320 lri_env%ri_fit%nvec(i1:i2) = int_aux(1:nbas)
321 END DO
322 IF (ASSOCIATED(norm)) DEALLOCATE (norm)
323 IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
324 DEALLOCATE (kind_of)
325 ELSE
326 cpabort('ri_type')
327 END IF
328
329 END SUBROUTINE lri_env_basis
330
331! **************************************************************************************************
332!> \brief initializes the lri basis: calculates the norm, self-overlap
333!> and integral of the ri basis
334!> \param lri_env ...
335! **************************************************************************************************
336 SUBROUTINE lri_basis_init(lri_env)
337 TYPE(lri_environment_type), POINTER :: lri_env
338
339 INTEGER :: ikind, nkind
340 INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index
341 REAL(kind=dp) :: delta
342 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dovlp3
343 REAL(kind=dp), DIMENSION(:), POINTER :: orb_norm_r, ri_int_fbas, ri_norm_r, &
344 ri_norm_s
345 REAL(kind=dp), DIMENSION(:, :), POINTER :: orb_ovlp, ri_ovlp, ri_ovlp_inv
346 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: scon_orb, scon_ri
347 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix
348 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
349
350 IF (ASSOCIATED(lri_env)) THEN
351 IF (ASSOCIATED(lri_env%orb_basis)) THEN
352 cpassert(ASSOCIATED(lri_env%ri_basis))
353 nkind = SIZE(lri_env%orb_basis)
354 CALL deallocate_bas_properties(lri_env)
355 ALLOCATE (lri_env%bas_prop(nkind))
356 DO ikind = 1, nkind
357 NULLIFY (orb_basis, ri_basis)
358 orb_basis => lri_env%orb_basis(ikind)%gto_basis_set
359 IF (ASSOCIATED(orb_basis)) THEN
360 ri_basis => lri_env%ri_basis(ikind)%gto_basis_set
361 cpassert(ASSOCIATED(ri_basis))
362 NULLIFY (ri_norm_r)
363 CALL basis_norm_radial(ri_basis, ri_norm_r)
364 NULLIFY (orb_norm_r)
365 CALL basis_norm_radial(orb_basis, orb_norm_r)
366 NULLIFY (ri_norm_s)
367 CALL basis_norm_s_func(ri_basis, ri_norm_s)
368 NULLIFY (ri_int_fbas)
369 CALL basis_int(ri_basis, ri_int_fbas, ri_norm_s)
370 lri_env%bas_prop(ikind)%int_fbas => ri_int_fbas
371 NULLIFY (ri_ovlp)
372 CALL basis_ovlp(ri_basis, ri_ovlp, ri_norm_r)
373 lri_env%bas_prop(ikind)%ri_ovlp => ri_ovlp
374 NULLIFY (orb_ovlp)
375 CALL basis_ovlp(orb_basis, orb_ovlp, orb_norm_r)
376 lri_env%bas_prop(ikind)%orb_ovlp => orb_ovlp
377 NULLIFY (scon_ri)
378 CALL contraction_matrix_shg(ri_basis, scon_ri)
379 lri_env%bas_prop(ikind)%scon_ri => scon_ri
380 NULLIFY (scon_orb)
381 CALL contraction_matrix_shg(orb_basis, scon_orb)
382 lri_env%bas_prop(ikind)%scon_orb => scon_orb
383 NULLIFY (scon_mix)
384 CALL contraction_matrix_shg_mix(orb_basis, ri_basis, &
385 orb_index, ri_index, scon_mix)
386 lri_env%bas_prop(ikind)%scon_mix => scon_mix
387 lri_env%bas_prop(ikind)%orb_index => orb_index
388 lri_env%bas_prop(ikind)%ri_index => ri_index
389 ALLOCATE (lri_env%bas_prop(ikind)%ovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf))
390 ALLOCATE (dovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf, 3))
391 CALL int_overlap_aba_shg(lri_env%bas_prop(ikind)%ovlp3, dovlp3, [0.0_dp, 0.0_dp, 0.0_dp], &
392 orb_basis, orb_basis, ri_basis, scon_orb, &
393 scon_mix, orb_index, ri_index, &
394 lri_env%cg_shg%cg_coeff, &
395 lri_env%cg_shg%cg_none0_list, &
396 lri_env%cg_shg%ncg_none0, &
397 calculate_forces=.false.)
398 DEALLOCATE (orb_norm_r, ri_norm_r, ri_norm_s)
399 DEALLOCATE (dovlp3)
400 ALLOCATE (ri_ovlp_inv(ri_basis%nsgf, ri_basis%nsgf))
401 CALL invert_matrix(ri_ovlp, ri_ovlp_inv, delta, improve=.true.)
402 lri_env%bas_prop(ikind)%ri_ovlp_inv => ri_ovlp_inv
403 END IF
404 END DO
405 END IF
406 END IF
407
408 END SUBROUTINE lri_basis_init
409
410! **************************************************************************************************
411!> \brief normalization for a contracted Gaussian s-function,
412!> spherical = cartesian Gaussian for s-functions
413!> \param basis ...
414!> \param norm ...
415! **************************************************************************************************
416 SUBROUTINE basis_norm_s_func(basis, norm)
417
418 TYPE(gto_basis_set_type), POINTER :: basis
419 REAL(dp), DIMENSION(:), POINTER :: norm
420
421 INTEGER :: ipgf, iset, isgf, ishell, jpgf, l, nbas
422 REAL(kind=dp) :: aai, aaj, cci, ccj, expa, ppl
423
424 NULLIFY (norm)
425
426 nbas = basis%nsgf
427 ALLOCATE (norm(nbas))
428 norm = 0._dp
429
430 DO iset = 1, basis%nset
431 DO ishell = 1, basis%nshell(iset)
432 l = basis%l(ishell, iset)
433 IF (l /= 0) cycle
434 expa = 0.5_dp*real(2*l + 3, dp)
435 ppl = pi**(3._dp/2._dp)
436 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
437 DO ipgf = 1, basis%npgf(iset)
438 cci = basis%gcc(ipgf, ishell, iset)
439 aai = basis%zet(ipgf, iset)
440 DO jpgf = 1, basis%npgf(iset)
441 ccj = basis%gcc(jpgf, ishell, iset)
442 aaj = basis%zet(jpgf, iset)
443 norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
444 END DO
445 END DO
446 norm(isgf) = 1.0_dp/sqrt(norm(isgf))
447 END DO
448 END DO
449 END DO
450
451 END SUBROUTINE basis_norm_s_func
452
453! **************************************************************************************************
454!> \brief normalization for radial part of contracted spherical Gaussian
455!> functions
456!> \param basis ...
457!> \param norm ...
458! **************************************************************************************************
459 SUBROUTINE basis_norm_radial(basis, norm)
460
461 TYPE(gto_basis_set_type), POINTER :: basis
462 REAL(dp), DIMENSION(:), POINTER :: norm
463
464 INTEGER :: ipgf, iset, isgf, ishell, jpgf, l, nbas
465 REAL(kind=dp) :: aai, aaj, cci, ccj, expa, ppl
466
467 NULLIFY (norm)
468
469 nbas = basis%nsgf
470 ALLOCATE (norm(nbas))
471 norm = 0._dp
472
473 DO iset = 1, basis%nset
474 DO ishell = 1, basis%nshell(iset)
475 l = basis%l(ishell, iset)
476 expa = 0.5_dp*real(2*l + 3, dp)
477 ppl = fac(2*l + 2)*rootpi/2._dp**real(2*l + 3, dp)/fac(l + 1)
478 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
479 DO ipgf = 1, basis%npgf(iset)
480 cci = basis%gcc(ipgf, ishell, iset)
481 aai = basis%zet(ipgf, iset)
482 DO jpgf = 1, basis%npgf(iset)
483 ccj = basis%gcc(jpgf, ishell, iset)
484 aaj = basis%zet(jpgf, iset)
485 norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
486 END DO
487 END DO
488 norm(isgf) = 1.0_dp/sqrt(norm(isgf))
489 END DO
490 END DO
491 END DO
492
493 END SUBROUTINE basis_norm_radial
494
495!*****************************************************************************
496!> \brief integral over a single (contracted) lri auxiliary basis function,
497!> integral is zero for all but s-functions
498!> \param basis ...
499!> \param int_aux ...
500!> \param norm ...
501! **************************************************************************************************
502 SUBROUTINE basis_int(basis, int_aux, norm)
503
504 TYPE(gto_basis_set_type), POINTER :: basis
505 REAL(dp), DIMENSION(:), POINTER :: int_aux, norm
506
507 INTEGER :: ipgf, iset, isgf, ishell, l, nbas
508 REAL(kind=dp) :: aa, cc, pp
509
510 nbas = basis%nsgf
511 ALLOCATE (int_aux(nbas))
512 int_aux = 0._dp
513
514 DO iset = 1, basis%nset
515 DO ishell = 1, basis%nshell(iset)
516 l = basis%l(ishell, iset)
517 IF (l /= 0) cycle
518 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
519 DO ipgf = 1, basis%npgf(iset)
520 cc = basis%gcc(ipgf, ishell, iset)
521 aa = basis%zet(ipgf, iset)
522 pp = (pi/aa)**(3._dp/2._dp)
523 int_aux(isgf) = int_aux(isgf) + norm(isgf)*cc*pp
524 END DO
525 END DO
526 END DO
527 END DO
528
529 END SUBROUTINE basis_int
530
531!*****************************************************************************
532!> \brief self-overlap of lri basis for contracted spherical Gaussians.
533!> Overlap of radial part. Norm contains only normalization of radial
534!> part. Norm and overlap of spherical harmonics not explicitly
535!> calculated since this cancels for the self-overlap anyway.
536!> \param basis ...
537!> \param ovlp ...
538!> \param norm ...
539! **************************************************************************************************
540 SUBROUTINE basis_ovlp(basis, ovlp, norm)
541
542 TYPE(gto_basis_set_type), POINTER :: basis
543 REAL(dp), DIMENSION(:, :), POINTER :: ovlp
544 REAL(dp), DIMENSION(:), POINTER :: norm
545
546 INTEGER :: ipgf, iset, isgf, ishell, jpgf, jset, &
547 jsgf, jshell, l, li, lj, m_i, m_j, nbas
548 REAL(kind=dp) :: aai, aaj, cci, ccj, expa, norm_i, &
549 norm_j, oo, ppl
550
551 nbas = basis%nsgf
552 ALLOCATE (ovlp(nbas, nbas))
553 ovlp = 0._dp
554
555 DO iset = 1, basis%nset
556 DO ishell = 1, basis%nshell(iset)
557 li = basis%l(ishell, iset)
558 DO jset = 1, basis%nset
559 DO jshell = 1, basis%nshell(jset)
560 lj = basis%l(jshell, jset)
561 IF (li == lj) THEN
562 l = li
563 expa = 0.5_dp*real(2*l + 3, dp)
564 ppl = fac(2*l + 2)*rootpi/2._dp**real(2*l + 3, dp)/fac(l + 1)
565 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
566 m_i = basis%m(isgf)
567 DO jsgf = basis%first_sgf(jshell, jset), basis%last_sgf(jshell, jset)
568 m_j = basis%m(jsgf)
569 IF (m_i == m_j) THEN
570 DO ipgf = 1, basis%npgf(iset)
571 cci = basis%gcc(ipgf, ishell, iset)
572 aai = basis%zet(ipgf, iset)
573 norm_i = norm(isgf)
574 DO jpgf = 1, basis%npgf(jset)
575 ccj = basis%gcc(jpgf, jshell, jset)
576 aaj = basis%zet(jpgf, jset)
577 oo = 1._dp/(aai + aaj)**expa
578 norm_j = norm(jsgf)
579 ovlp(isgf, jsgf) = ovlp(isgf, jsgf) + norm_i*norm_j*ppl*cci*ccj*oo
580 END DO
581 END DO
582 END IF
583 END DO
584 END DO
585 END IF
586 END DO
587 END DO
588 END DO
589 END DO
590
591 END SUBROUTINE basis_ovlp
592
593END MODULE lri_environment_init
mathlib::invert_matrix
Definition mathlib.F:74
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition ao_util.F:96
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
basis_set_types
Definition basis_set_types.F:15
basis_set_types::copy_gto_basis_set
subroutine, public copy_gto_basis_set(basis_set_in, basis_set_out)
...
Definition basis_set_types.F:214
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::golze2017b
integer, save, public golze2017b
Definition bibliography.F:36
bibliography::golze2017a
integer, save, public golze2017a
Definition bibliography.F:36
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
generic_shg_integrals_init
Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation of contracte...
Definition generic_shg_integrals_init.F:16
generic_shg_integrals_init::contraction_matrix_shg_mix
subroutine, public contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix)
mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis at the same atom
Definition generic_shg_integrals_init.F:135
generic_shg_integrals_init::contraction_matrix_shg
subroutine, public contraction_matrix_shg(basis, scon_shg)
contraction matrix for SHG integrals
Definition generic_shg_integrals_init.F:49
generic_shg_integrals_init::get_clebsch_gordon_coefficients
subroutine, public get_clebsch_gordon_coefficients(my_cg, cg_none0_list, ncg_none0, maxl1, maxl2)
calculate the Clebsch-Gordon (CG) coefficients for expansion of the product of two spherical harmonic...
Definition generic_shg_integrals_init.F:308
generic_shg_integrals
Calculation of contracted, spherical Gaussian integrals using the solid harmonic Gaussian (SHG) integ...
Definition generic_shg_integrals.F:21
generic_shg_integrals::int_overlap_aba_shg
subroutine, public int_overlap_aba_shg(saba, dsaba, rab, oba, obb, fba, scon_obb, scona_mix, oba_index, fba_index, cg_coeff, cg_none0_list, ncg_none0, calculate_forces)
calculate integrals (a,b,fa)
Definition generic_shg_integrals.F:251
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
lri_environment_init
initializes the environment for lri lri : local resolution of the identity
Definition lri_environment_init.F:15
lri_environment_init::lri_basis_init
subroutine, public lri_basis_init(lri_env)
initializes the lri basis: calculates the norm, self-overlap and integral of the ri basis
Definition lri_environment_init.F:337
lri_environment_init::lri_env_init
subroutine, public lri_env_init(lri_env, lri_section)
initializes the lri env
Definition lri_environment_init.F:65
lri_environment_init::lri_env_basis
subroutine, public lri_env_basis(ri_type, qs_env, lri_env, qs_kind_set)
initializes the lri env
Definition lri_environment_init.F:126
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
lri_environment_types::deallocate_bas_properties
subroutine, public deallocate_bas_properties(lri_env)
deallocates one-center overlap integrals, integral of ri basis and scon matrices
Definition lri_environment_types.F:1060
lri_environment_types::lri_env_create
subroutine, public lri_env_create(lri_env)
creates and initializes an lri_env
Definition lri_environment_types.F:392
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:528
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::dft_control_type
Definition cp_control_types.F:628
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:271
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177