(git:d18deda)
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-2025 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 ALLOCATE (kind_of(natom))
287 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
288 nbas = 0
289 DO iat = 1, natom
290 ikind = kind_of(iat)
291 nribas = lri_env%ri_basis(ikind)%gto_basis_set%nsgf
292 lri_env%ri_fit%bas_ptr(1, iat) = nbas + 1
293 lri_env%ri_fit%bas_ptr(2, iat) = nbas + nribas
294 nbas = nbas + nribas
295 END DO
296 ! initialize vector t
297 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
298 nspin = dft_control%nspins
299 ALLOCATE (lri_env%ri_fit%tvec(nbas, nspin), lri_env%ri_fit%rm1t(nbas, nspin))
300 ! initialize vector a, expansion of density
301 ALLOCATE (lri_env%ri_fit%avec(nbas, nspin))
302 ! initialize vector fout, R^(-1)*(f-p*n)
303 ALLOCATE (lri_env%ri_fit%fout(nbas, nspin))
304 ! initialize vector with RI basis integrated
305 NULLIFY (norm, int_aux)
306 nbas = lri_env%ri_fit%bas_ptr(2, natom)
307 ALLOCATE (lri_env%ri_fit%nvec(nbas), lri_env%ri_fit%rm1n(nbas))
308 ikind = 0
309 DO iat = 1, natom
310 IF (ikind /= kind_of(iat)) THEN
311 ikind = kind_of(iat)
312 ri_basis_set => lri_env%ri_basis(ikind)%gto_basis_set
313 IF (ASSOCIATED(norm)) DEALLOCATE (norm)
314 IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
315 CALL basis_norm_s_func(ri_basis_set, norm)
316 CALL basis_int(ri_basis_set, int_aux, norm)
317 END IF
318 nbas = SIZE(int_aux)
319 i1 = lri_env%ri_fit%bas_ptr(1, iat)
320 i2 = lri_env%ri_fit%bas_ptr(2, iat)
321 lri_env%ri_fit%nvec(i1:i2) = int_aux(1:nbas)
322 END DO
323 IF (ASSOCIATED(norm)) DEALLOCATE (norm)
324 IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
325 DEALLOCATE (kind_of)
326 ELSE
327 cpabort('ri_type')
328 END IF
329
330 END SUBROUTINE lri_env_basis
331
332! **************************************************************************************************
333!> \brief initializes the lri basis: calculates the norm, self-overlap
334!> and integral of the ri basis
335!> \param lri_env ...
336! **************************************************************************************************
337 SUBROUTINE lri_basis_init(lri_env)
338 TYPE(lri_environment_type), POINTER :: lri_env
339
340 INTEGER :: ikind, nkind
341 INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index
342 REAL(kind=dp) :: delta
343 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dovlp3
344 REAL(kind=dp), DIMENSION(:), POINTER :: orb_norm_r, ri_int_fbas, ri_norm_r, &
345 ri_norm_s
346 REAL(kind=dp), DIMENSION(:, :), POINTER :: orb_ovlp, ri_ovlp, ri_ovlp_inv
347 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: scon_orb, scon_ri
348 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix
349 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
350
351 IF (ASSOCIATED(lri_env)) THEN
352 IF (ASSOCIATED(lri_env%orb_basis)) THEN
353 cpassert(ASSOCIATED(lri_env%ri_basis))
354 nkind = SIZE(lri_env%orb_basis)
355 CALL deallocate_bas_properties(lri_env)
356 ALLOCATE (lri_env%bas_prop(nkind))
357 DO ikind = 1, nkind
358 NULLIFY (orb_basis, ri_basis)
359 orb_basis => lri_env%orb_basis(ikind)%gto_basis_set
360 IF (ASSOCIATED(orb_basis)) THEN
361 ri_basis => lri_env%ri_basis(ikind)%gto_basis_set
362 cpassert(ASSOCIATED(ri_basis))
363 NULLIFY (ri_norm_r)
364 CALL basis_norm_radial(ri_basis, ri_norm_r)
365 NULLIFY (orb_norm_r)
366 CALL basis_norm_radial(orb_basis, orb_norm_r)
367 NULLIFY (ri_norm_s)
368 CALL basis_norm_s_func(ri_basis, ri_norm_s)
369 NULLIFY (ri_int_fbas)
370 CALL basis_int(ri_basis, ri_int_fbas, ri_norm_s)
371 lri_env%bas_prop(ikind)%int_fbas => ri_int_fbas
372 NULLIFY (ri_ovlp)
373 CALL basis_ovlp(ri_basis, ri_ovlp, ri_norm_r)
374 lri_env%bas_prop(ikind)%ri_ovlp => ri_ovlp
375 NULLIFY (orb_ovlp)
376 CALL basis_ovlp(orb_basis, orb_ovlp, orb_norm_r)
377 lri_env%bas_prop(ikind)%orb_ovlp => orb_ovlp
378 NULLIFY (scon_ri)
379 CALL contraction_matrix_shg(ri_basis, scon_ri)
380 lri_env%bas_prop(ikind)%scon_ri => scon_ri
381 NULLIFY (scon_orb)
382 CALL contraction_matrix_shg(orb_basis, scon_orb)
383 lri_env%bas_prop(ikind)%scon_orb => scon_orb
384 NULLIFY (scon_mix)
385 CALL contraction_matrix_shg_mix(orb_basis, ri_basis, &
386 orb_index, ri_index, scon_mix)
387 lri_env%bas_prop(ikind)%scon_mix => scon_mix
388 lri_env%bas_prop(ikind)%orb_index => orb_index
389 lri_env%bas_prop(ikind)%ri_index => ri_index
390 ALLOCATE (lri_env%bas_prop(ikind)%ovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf))
391 ALLOCATE (dovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf, 3))
392 CALL int_overlap_aba_shg(lri_env%bas_prop(ikind)%ovlp3, dovlp3, (/0.0_dp, 0.0_dp, 0.0_dp/), &
393 orb_basis, orb_basis, ri_basis, scon_orb, &
394 scon_mix, orb_index, ri_index, &
395 lri_env%cg_shg%cg_coeff, &
396 lri_env%cg_shg%cg_none0_list, &
397 lri_env%cg_shg%ncg_none0, &
398 calculate_forces=.false.)
399 DEALLOCATE (orb_norm_r, ri_norm_r, ri_norm_s)
400 DEALLOCATE (dovlp3)
401 ALLOCATE (ri_ovlp_inv(ri_basis%nsgf, ri_basis%nsgf))
402 CALL invert_matrix(ri_ovlp, ri_ovlp_inv, delta, improve=.true.)
403 lri_env%bas_prop(ikind)%ri_ovlp_inv => ri_ovlp_inv
404 END IF
405 END DO
406 END IF
407 END IF
408
409 END SUBROUTINE lri_basis_init
410
411! **************************************************************************************************
412!> \brief normalization for a contracted Gaussian s-function,
413!> spherical = cartesian Gaussian for s-functions
414!> \param basis ...
415!> \param norm ...
416! **************************************************************************************************
417 SUBROUTINE basis_norm_s_func(basis, norm)
418
419 TYPE(gto_basis_set_type), POINTER :: basis
420 REAL(dp), DIMENSION(:), POINTER :: norm
421
422 INTEGER :: ipgf, iset, isgf, ishell, jpgf, l, nbas
423 REAL(kind=dp) :: aai, aaj, cci, ccj, expa, ppl
424
425 NULLIFY (norm)
426
427 nbas = basis%nsgf
428 ALLOCATE (norm(nbas))
429 norm = 0._dp
430
431 DO iset = 1, basis%nset
432 DO ishell = 1, basis%nshell(iset)
433 l = basis%l(ishell, iset)
434 IF (l /= 0) cycle
435 expa = 0.5_dp*real(2*l + 3, dp)
436 ppl = pi**(3._dp/2._dp)
437 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
438 DO ipgf = 1, basis%npgf(iset)
439 cci = basis%gcc(ipgf, ishell, iset)
440 aai = basis%zet(ipgf, iset)
441 DO jpgf = 1, basis%npgf(iset)
442 ccj = basis%gcc(jpgf, ishell, iset)
443 aaj = basis%zet(jpgf, iset)
444 norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
445 END DO
446 END DO
447 norm(isgf) = 1.0_dp/sqrt(norm(isgf))
448 END DO
449 END DO
450 END DO
451
452 END SUBROUTINE basis_norm_s_func
453
454! **************************************************************************************************
455!> \brief normalization for radial part of contracted spherical Gaussian
456!> functions
457!> \param basis ...
458!> \param norm ...
459! **************************************************************************************************
460 SUBROUTINE basis_norm_radial(basis, norm)
461
462 TYPE(gto_basis_set_type), POINTER :: basis
463 REAL(dp), DIMENSION(:), POINTER :: norm
464
465 INTEGER :: ipgf, iset, isgf, ishell, jpgf, l, nbas
466 REAL(kind=dp) :: aai, aaj, cci, ccj, expa, ppl
467
468 NULLIFY (norm)
469
470 nbas = basis%nsgf
471 ALLOCATE (norm(nbas))
472 norm = 0._dp
473
474 DO iset = 1, basis%nset
475 DO ishell = 1, basis%nshell(iset)
476 l = basis%l(ishell, iset)
477 expa = 0.5_dp*real(2*l + 3, dp)
478 ppl = fac(2*l + 2)*rootpi/2._dp**real(2*l + 3, dp)/fac(l + 1)
479 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
480 DO ipgf = 1, basis%npgf(iset)
481 cci = basis%gcc(ipgf, ishell, iset)
482 aai = basis%zet(ipgf, iset)
483 DO jpgf = 1, basis%npgf(iset)
484 ccj = basis%gcc(jpgf, ishell, iset)
485 aaj = basis%zet(jpgf, iset)
486 norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
487 END DO
488 END DO
489 norm(isgf) = 1.0_dp/sqrt(norm(isgf))
490 END DO
491 END DO
492 END DO
493
494 END SUBROUTINE basis_norm_radial
495
496!*****************************************************************************
497!> \brief integral over a single (contracted) lri auxiliary basis function,
498!> integral is zero for all but s-functions
499!> \param basis ...
500!> \param int_aux ...
501!> \param norm ...
502! **************************************************************************************************
503 SUBROUTINE basis_int(basis, int_aux, norm)
504
505 TYPE(gto_basis_set_type), POINTER :: basis
506 REAL(dp), DIMENSION(:), POINTER :: int_aux, norm
507
508 INTEGER :: ipgf, iset, isgf, ishell, l, nbas
509 REAL(kind=dp) :: aa, cc, pp
510
511 nbas = basis%nsgf
512 ALLOCATE (int_aux(nbas))
513 int_aux = 0._dp
514
515 DO iset = 1, basis%nset
516 DO ishell = 1, basis%nshell(iset)
517 l = basis%l(ishell, iset)
518 IF (l /= 0) cycle
519 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
520 DO ipgf = 1, basis%npgf(iset)
521 cc = basis%gcc(ipgf, ishell, iset)
522 aa = basis%zet(ipgf, iset)
523 pp = (pi/aa)**(3._dp/2._dp)
524 int_aux(isgf) = int_aux(isgf) + norm(isgf)*cc*pp
525 END DO
526 END DO
527 END DO
528 END DO
529
530 END SUBROUTINE basis_int
531
532!*****************************************************************************
533!> \brief self-overlap of lri basis for contracted spherical Gaussians.
534!> Overlap of radial part. Norm contains only normalization of radial
535!> part. Norm and overlap of spherical harmonics not explicitly
536!> calculated since this cancels for the self-overlap anyway.
537!> \param basis ...
538!> \param ovlp ...
539!> \param norm ...
540! **************************************************************************************************
541 SUBROUTINE basis_ovlp(basis, ovlp, norm)
542
543 TYPE(gto_basis_set_type), POINTER :: basis
544 REAL(dp), DIMENSION(:, :), POINTER :: ovlp
545 REAL(dp), DIMENSION(:), POINTER :: norm
546
547 INTEGER :: ipgf, iset, isgf, ishell, jpgf, jset, &
548 jsgf, jshell, l, li, lj, m_i, m_j, nbas
549 REAL(kind=dp) :: aai, aaj, cci, ccj, expa, norm_i, &
550 norm_j, oo, ppl
551
552 nbas = basis%nsgf
553 ALLOCATE (ovlp(nbas, nbas))
554 ovlp = 0._dp
555
556 DO iset = 1, basis%nset
557 DO ishell = 1, basis%nshell(iset)
558 li = basis%l(ishell, iset)
559 DO jset = 1, basis%nset
560 DO jshell = 1, basis%nshell(jset)
561 lj = basis%l(jshell, jset)
562 IF (li == lj) THEN
563 l = li
564 expa = 0.5_dp*real(2*l + 3, dp)
565 ppl = fac(2*l + 2)*rootpi/2._dp**real(2*l + 3, dp)/fac(l + 1)
566 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
567 m_i = basis%m(isgf)
568 DO jsgf = basis%first_sgf(jshell, jset), basis%last_sgf(jshell, jset)
569 m_j = basis%m(jsgf)
570 IF (m_i == m_j) THEN
571 DO ipgf = 1, basis%npgf(iset)
572 cci = basis%gcc(ipgf, ishell, iset)
573 aai = basis%zet(ipgf, iset)
574 norm_i = norm(isgf)
575 DO jpgf = 1, basis%npgf(jset)
576 ccj = basis%gcc(jpgf, jshell, jset)
577 aaj = basis%zet(jpgf, jset)
578 oo = 1._dp/(aai + aaj)**expa
579 norm_j = norm(jsgf)
580 ovlp(isgf, jsgf) = ovlp(isgf, jsgf) + norm_i*norm_j*ppl*cci*ccj*oo
581 END DO
582 END DO
583 END IF
584 END DO
585 END DO
586 END IF
587 END DO
588 END DO
589 END DO
590 END DO
591
592 END SUBROUTINE basis_ovlp
593
594END MODULE lri_environment_init
mathlib::invert_matrix
Definition mathlib.F:70
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:213
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:338
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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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:454
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:570
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:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171