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 initializes the environment for lri
10 : !> lri : local resolution of the identity
11 : !> \par History
12 : !> created [06.2015]
13 : !> \author Dorothea Golze
14 : ! **************************************************************************************************
15 : MODULE 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 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 : !> \brief initializes the lri env
61 : !> \param lri_env ...
62 : !> \param lri_section ...
63 : ! **************************************************************************************************
64 60 : 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 60 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii
70 :
71 : NULLIFY (lri_env)
72 0 : ALLOCATE (lri_env)
73 60 : 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 60 : l_val=lri_env%statistics)
79 : ! use exact one-center terms
80 : CALL section_vals_val_get(lri_section, "EXACT_1C_TERMS", &
81 60 : 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 60 : l_val=lri_env%ppl_ri)
85 : ! check for debug (OS scheme)
86 : CALL section_vals_val_get(lri_section, "DEBUG_LRI_INTEGRALS", &
87 60 : l_val=lri_env%debug)
88 : ! integrals based on solid harmonic Gaussians
89 : CALL section_vals_val_get(lri_section, "SHG_LRI_INTEGRALS", &
90 60 : 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 60 : i_val=lri_env%lri_overlap_inv)
94 : CALL section_vals_val_get(lri_section, "MAX_CONDITION_NUM", &
95 60 : r_val=lri_env%cond_max)
96 : ! integrals threshold (aba, abb)
97 : CALL section_vals_val_get(lri_section, "EPS_O3_INT", &
98 60 : r_val=lri_env%eps_o3_int)
99 : ! RI SINV
100 : CALL section_vals_val_get(lri_section, "RI_SINV", &
101 60 : c_val=lri_env%ri_sinv_app)
102 : ! Distant Pair Approximation
103 : CALL section_vals_val_get(lri_section, "DISTANT_PAIR_APPROXIMATION", &
104 60 : l_val=lri_env%distant_pair_approximation)
105 : CALL section_vals_val_get(lri_section, "DISTANT_PAIR_METHOD", &
106 60 : c_val=lri_env%distant_pair_method)
107 60 : CALL section_vals_val_get(lri_section, "DISTANT_PAIR_RADII", r_vals=radii)
108 60 : CPASSERT(SIZE(radii) == 2)
109 60 : CPASSERT(radii(2) > radii(1))
110 60 : CPASSERT(radii(1) > 0.0_dp)
111 60 : lri_env%r_in = radii(1)
112 60 : lri_env%r_out = radii(2)
113 :
114 60 : CALL cite_reference(Golze2017b)
115 60 : IF (lri_env%use_shg_integrals) CALL cite_reference(Golze2017a)
116 :
117 60 : 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 60 : 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 60 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
135 : REAL(KIND=dp) :: gcc, rad, rai, raj, xradius, zeta
136 60 : REAL(KIND=dp), DIMENSION(:), POINTER :: int_aux, norm
137 60 : 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 60 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
143 60 : nkind = SIZE(atomic_kind_set)
144 520 : ALLOCATE (lri_env%orb_basis(nkind), lri_env%ri_basis(nkind))
145 60 : maxl_orb = 0
146 60 : maxl_ri = 0
147 170 : DO ikind = 1, nkind
148 110 : NULLIFY (orb_basis_set, ri_basis_set)
149 110 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
150 110 : IF (ri_type == "LRI") THEN
151 90 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="LRI_AUX")
152 20 : ELSE IF (ri_type == "P_LRI") THEN
153 20 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="P_LRI_AUX")
154 0 : ELSE IF (ri_type == "RI") THEN
155 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC")
156 : ELSE
157 0 : CPABORT('ri_type')
158 : END IF
159 110 : NULLIFY (lri_env%orb_basis(ikind)%gto_basis_set)
160 110 : NULLIFY (lri_env%ri_basis(ikind)%gto_basis_set)
161 110 : IF (ASSOCIATED(orb_basis_set)) THEN
162 110 : CALL copy_gto_basis_set(orb_basis_set, lri_env%orb_basis(ikind)%gto_basis_set)
163 110 : CALL copy_gto_basis_set(ri_basis_set, lri_env%ri_basis(ikind)%gto_basis_set)
164 : END IF
165 264 : lmax_ikind_orb = MAXVAL(lri_env%orb_basis(ikind)%gto_basis_set%lmax)
166 1368 : lmax_ikind_ri = MAXVAL(lri_env%ri_basis(ikind)%gto_basis_set%lmax)
167 110 : maxl_orb = MAX(maxl_orb, lmax_ikind_orb)
168 170 : maxl_ri = MAX(maxl_ri, lmax_ikind_ri)
169 : END DO
170 :
171 60 : IF ((ri_type == "LRI") .OR. (ri_type == "P_LRI")) THEN
172 : ! CG coefficients needed for lri integrals
173 60 : 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 60 : maxl_orb, maxl_ri)
178 : END IF
179 60 : CALL lri_basis_init(lri_env)
180 : ! distant pair approximation
181 60 : IF (lri_env%distant_pair_approximation) THEN
182 : !
183 0 : SELECT CASE (lri_env%distant_pair_method)
184 : CASE ("EW")
185 : ! equal weight of 0.5
186 : CASE ("AW")
187 0 : ALLOCATE (lri_env%aradius(nkind))
188 0 : DO i = 1, nkind
189 0 : orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
190 0 : lri_env%aradius(i) = orb_basis_set%kind_radius
191 : END DO
192 : CASE ("SW")
193 0 : ALLOCATE (lri_env%wbas(nkind))
194 0 : DO i = 1, nkind
195 0 : orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
196 0 : n1 = orb_basis_set%nsgf
197 0 : ALLOCATE (lri_env%wbas(i)%vec(n1))
198 0 : DO iset = 1, orb_basis_set%nset
199 0 : i1 = orb_basis_set%first_sgf(1, iset)
200 0 : n2 = orb_basis_set%nshell(iset)
201 0 : i2 = orb_basis_set%last_sgf(n2, iset)
202 0 : lri_env%wbas(i)%vec(i1:i2) = orb_basis_set%set_radius(iset)
203 : END DO
204 : END DO
205 : CASE ("LW")
206 48 : ALLOCATE (lri_env%wbas(nkind))
207 46 : DO i = 1, nkind
208 30 : orb_basis_set => lri_env%orb_basis(i)%gto_basis_set
209 30 : n1 = orb_basis_set%nsgf
210 90 : ALLOCATE (lri_env%wbas(i)%vec(n1))
211 76 : DO iset = 1, orb_basis_set%nset
212 182 : DO ishell = 1, orb_basis_set%nshell(iset)
213 122 : i1 = orb_basis_set%first_sgf(ishell, iset)
214 122 : i2 = orb_basis_set%last_sgf(ishell, iset)
215 122 : l = orb_basis_set%l(ishell, iset)
216 122 : xradius = 0.0_dp
217 956 : DO ipgf = 1, orb_basis_set%npgf(iset)
218 834 : gcc = orb_basis_set%gcc(ipgf, ishell, iset)
219 834 : zeta = orb_basis_set%zet(ipgf, iset)
220 834 : rad = exp_radius(l, zeta, 1.e-5_dp, gcc, rlow=xradius)
221 956 : xradius = MAX(xradius, rad)
222 : END DO
223 430 : lri_env%wbas(i)%vec(i1:i2) = xradius
224 : END DO
225 : END DO
226 : END DO
227 : CASE DEFAULT
228 18 : CPABORT("Unknown DISTANT_PAIR_METHOD in LRI")
229 : END SELECT
230 : !
231 72 : ALLOCATE (lri_env%wmat(nkind, nkind))
232 18 : SELECT CASE (lri_env%distant_pair_method)
233 : CASE ("EW")
234 : ! equal weight of 0.5
235 6 : DO i1 = 1, nkind
236 4 : n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
237 14 : DO i2 = 1, nkind
238 8 : n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
239 32 : ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
240 732 : lri_env%wmat(i1, i2)%mat(:, :) = 0.5_dp
241 : END DO
242 : END DO
243 : CASE ("AW")
244 0 : DO i1 = 1, nkind
245 0 : n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
246 0 : DO i2 = 1, nkind
247 0 : n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
248 0 : ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
249 0 : rai = lri_env%aradius(i1)**2
250 0 : raj = lri_env%aradius(i2)**2
251 0 : IF (raj > rai) THEN
252 0 : lri_env%wmat(i1, i2)%mat(:, :) = 1.0_dp
253 : ELSE
254 0 : lri_env%wmat(i1, i2)%mat(:, :) = 0.0_dp
255 : END IF
256 : END DO
257 : END DO
258 : CASE ("SW", "LW")
259 48 : DO i1 = 1, nkind
260 30 : n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf
261 104 : DO i2 = 1, nkind
262 58 : n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf
263 232 : ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2))
264 618 : DO ip = 1, SIZE(lri_env%wbas(i1)%vec)
265 530 : rai = lri_env%wbas(i1)%vec(ip)**2
266 5462 : DO jp = 1, SIZE(lri_env%wbas(i2)%vec)
267 4874 : raj = lri_env%wbas(i2)%vec(jp)**2
268 5404 : IF (raj > rai) THEN
269 2000 : lri_env%wmat(i1, i2)%mat(ip, jp) = 1.0_dp
270 : ELSE
271 2874 : 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 0 : ELSE IF (ri_type == "RI") THEN
280 0 : ALLOCATE (lri_env%ri_fit)
281 0 : NULLIFY (lri_env%ri_fit%nvec)
282 0 : NULLIFY (lri_env%ri_fit%bas_ptr)
283 0 : CALL get_qs_env(qs_env=qs_env, natom=natom)
284 : ! initialize pointers to RI basis vector
285 0 : ALLOCATE (lri_env%ri_fit%bas_ptr(2, natom))
286 0 : ALLOCATE (kind_of(natom))
287 0 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
288 0 : nbas = 0
289 0 : DO iat = 1, natom
290 0 : ikind = kind_of(iat)
291 0 : nribas = lri_env%ri_basis(ikind)%gto_basis_set%nsgf
292 0 : lri_env%ri_fit%bas_ptr(1, iat) = nbas + 1
293 0 : lri_env%ri_fit%bas_ptr(2, iat) = nbas + nribas
294 0 : nbas = nbas + nribas
295 : END DO
296 : ! initialize vector t
297 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
298 0 : nspin = dft_control%nspins
299 0 : ALLOCATE (lri_env%ri_fit%tvec(nbas, nspin), lri_env%ri_fit%rm1t(nbas, nspin))
300 : ! initialize vector a, expansion of density
301 0 : ALLOCATE (lri_env%ri_fit%avec(nbas, nspin))
302 : ! initialize vector fout, R^(-1)*(f-p*n)
303 0 : ALLOCATE (lri_env%ri_fit%fout(nbas, nspin))
304 : ! initialize vector with RI basis integrated
305 0 : NULLIFY (norm, int_aux)
306 0 : nbas = lri_env%ri_fit%bas_ptr(2, natom)
307 0 : ALLOCATE (lri_env%ri_fit%nvec(nbas), lri_env%ri_fit%rm1n(nbas))
308 0 : ikind = 0
309 0 : DO iat = 1, natom
310 0 : IF (ikind /= kind_of(iat)) THEN
311 0 : ikind = kind_of(iat)
312 0 : ri_basis_set => lri_env%ri_basis(ikind)%gto_basis_set
313 0 : IF (ASSOCIATED(norm)) DEALLOCATE (norm)
314 0 : IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
315 0 : CALL basis_norm_s_func(ri_basis_set, norm)
316 0 : CALL basis_int(ri_basis_set, int_aux, norm)
317 : END IF
318 0 : nbas = SIZE(int_aux)
319 0 : i1 = lri_env%ri_fit%bas_ptr(1, iat)
320 0 : i2 = lri_env%ri_fit%bas_ptr(2, iat)
321 0 : lri_env%ri_fit%nvec(i1:i2) = int_aux(1:nbas)
322 : END DO
323 0 : IF (ASSOCIATED(norm)) DEALLOCATE (norm)
324 0 : IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux)
325 0 : DEALLOCATE (kind_of)
326 : ELSE
327 0 : CPABORT('ri_type')
328 : END IF
329 :
330 120 : 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 78 : SUBROUTINE lri_basis_init(lri_env)
338 : TYPE(lri_environment_type), POINTER :: lri_env
339 :
340 : INTEGER :: ikind, nkind
341 78 : INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index
342 : REAL(KIND=dp) :: delta
343 78 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dovlp3
344 78 : REAL(KIND=dp), DIMENSION(:), POINTER :: orb_norm_r, ri_int_fbas, ri_norm_r, &
345 78 : ri_norm_s
346 78 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: orb_ovlp, ri_ovlp, ri_ovlp_inv
347 78 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_orb, scon_ri
348 78 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix
349 : TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
350 :
351 78 : IF (ASSOCIATED(lri_env)) THEN
352 78 : IF (ASSOCIATED(lri_env%orb_basis)) THEN
353 78 : CPASSERT(ASSOCIATED(lri_env%ri_basis))
354 78 : nkind = SIZE(lri_env%orb_basis)
355 78 : CALL deallocate_bas_properties(lri_env)
356 362 : ALLOCATE (lri_env%bas_prop(nkind))
357 206 : DO ikind = 1, nkind
358 128 : NULLIFY (orb_basis, ri_basis)
359 128 : orb_basis => lri_env%orb_basis(ikind)%gto_basis_set
360 206 : IF (ASSOCIATED(orb_basis)) THEN
361 128 : ri_basis => lri_env%ri_basis(ikind)%gto_basis_set
362 128 : CPASSERT(ASSOCIATED(ri_basis))
363 128 : NULLIFY (ri_norm_r)
364 128 : CALL basis_norm_radial(ri_basis, ri_norm_r)
365 128 : NULLIFY (orb_norm_r)
366 128 : CALL basis_norm_radial(orb_basis, orb_norm_r)
367 128 : NULLIFY (ri_norm_s)
368 128 : CALL basis_norm_s_func(ri_basis, ri_norm_s)
369 128 : NULLIFY (ri_int_fbas)
370 128 : CALL basis_int(ri_basis, ri_int_fbas, ri_norm_s)
371 128 : lri_env%bas_prop(ikind)%int_fbas => ri_int_fbas
372 128 : NULLIFY (ri_ovlp)
373 128 : CALL basis_ovlp(ri_basis, ri_ovlp, ri_norm_r)
374 128 : lri_env%bas_prop(ikind)%ri_ovlp => ri_ovlp
375 128 : NULLIFY (orb_ovlp)
376 128 : CALL basis_ovlp(orb_basis, orb_ovlp, orb_norm_r)
377 128 : lri_env%bas_prop(ikind)%orb_ovlp => orb_ovlp
378 128 : NULLIFY (scon_ri)
379 128 : CALL contraction_matrix_shg(ri_basis, scon_ri)
380 128 : lri_env%bas_prop(ikind)%scon_ri => scon_ri
381 128 : NULLIFY (scon_orb)
382 128 : CALL contraction_matrix_shg(orb_basis, scon_orb)
383 128 : lri_env%bas_prop(ikind)%scon_orb => scon_orb
384 128 : NULLIFY (scon_mix)
385 : CALL contraction_matrix_shg_mix(orb_basis, ri_basis, &
386 128 : orb_index, ri_index, scon_mix)
387 128 : lri_env%bas_prop(ikind)%scon_mix => scon_mix
388 128 : lri_env%bas_prop(ikind)%orb_index => orb_index
389 128 : lri_env%bas_prop(ikind)%ri_index => ri_index
390 640 : ALLOCATE (lri_env%bas_prop(ikind)%ovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf))
391 768 : 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 128 : calculate_forces=.FALSE.)
399 128 : DEALLOCATE (orb_norm_r, ri_norm_r, ri_norm_s)
400 128 : DEALLOCATE (dovlp3)
401 512 : ALLOCATE (ri_ovlp_inv(ri_basis%nsgf, ri_basis%nsgf))
402 128 : CALL invert_matrix(ri_ovlp, ri_ovlp_inv, delta, improve=.TRUE.)
403 128 : 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 156 : 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 128 : 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 128 : NULLIFY (norm)
426 :
427 128 : nbas = basis%nsgf
428 384 : ALLOCATE (norm(nbas))
429 13778 : norm = 0._dp
430 :
431 1504 : DO iset = 1, basis%nset
432 5138 : DO ishell = 1, basis%nshell(iset)
433 3634 : l = basis%l(ishell, iset)
434 3634 : IF (l /= 0) CYCLE
435 1106 : expa = 0.5_dp*REAL(2*l + 3, dp)
436 1106 : ppl = pi**(3._dp/2._dp)
437 3588 : DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
438 2580 : DO ipgf = 1, basis%npgf(iset)
439 1474 : cci = basis%gcc(ipgf, ishell, iset)
440 1474 : aai = basis%zet(ipgf, iset)
441 6442 : DO jpgf = 1, basis%npgf(iset)
442 3862 : ccj = basis%gcc(jpgf, ishell, iset)
443 3862 : aaj = basis%zet(jpgf, iset)
444 5336 : norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
445 : END DO
446 : END DO
447 4740 : norm(isgf) = 1.0_dp/SQRT(norm(isgf))
448 : END DO
449 : END DO
450 : END DO
451 :
452 128 : 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 256 : 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 256 : NULLIFY (norm)
469 :
470 256 : nbas = basis%nsgf
471 768 : ALLOCATE (norm(nbas))
472 14884 : norm = 0._dp
473 :
474 1804 : DO iset = 1, basis%nset
475 5868 : DO ishell = 1, basis%nshell(iset)
476 4064 : l = basis%l(ishell, iset)
477 4064 : expa = 0.5_dp*REAL(2*l + 3, dp)
478 4064 : ppl = fac(2*l + 2)*rootpi/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
479 20240 : DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
480 36428 : DO ipgf = 1, basis%npgf(iset)
481 21800 : cci = basis%gcc(ipgf, ishell, iset)
482 21800 : aai = basis%zet(ipgf, iset)
483 105212 : DO jpgf = 1, basis%npgf(iset)
484 68784 : ccj = basis%gcc(jpgf, ishell, iset)
485 68784 : aaj = basis%zet(jpgf, iset)
486 90584 : norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa
487 : END DO
488 : END DO
489 18692 : norm(isgf) = 1.0_dp/SQRT(norm(isgf))
490 : END DO
491 : END DO
492 : END DO
493 :
494 256 : 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 128 : 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 128 : nbas = basis%nsgf
512 384 : ALLOCATE (int_aux(nbas))
513 13778 : int_aux = 0._dp
514 :
515 1504 : DO iset = 1, basis%nset
516 5138 : DO ishell = 1, basis%nshell(iset)
517 3634 : l = basis%l(ishell, iset)
518 3634 : IF (l /= 0) CYCLE
519 3588 : DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
520 6214 : DO ipgf = 1, basis%npgf(iset)
521 1474 : cc = basis%gcc(ipgf, ishell, iset)
522 1474 : aa = basis%zet(ipgf, iset)
523 1474 : pp = (pi/aa)**(3._dp/2._dp)
524 2580 : int_aux(isgf) = int_aux(isgf) + norm(isgf)*cc*pp
525 : END DO
526 : END DO
527 : END DO
528 : END DO
529 :
530 128 : 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 256 : 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 256 : nbas = basis%nsgf
553 1024 : ALLOCATE (ovlp(nbas, nbas))
554 2343004 : ovlp = 0._dp
555 :
556 1804 : DO iset = 1, basis%nset
557 5868 : DO ishell = 1, basis%nshell(iset)
558 4064 : li = basis%l(ishell, iset)
559 52122 : DO jset = 1, basis%nset
560 192430 : DO jshell = 1, basis%nshell(jset)
561 141856 : lj = basis%l(jshell, jset)
562 188366 : IF (li == lj) THEN
563 34508 : l = li
564 34508 : expa = 0.5_dp*REAL(2*l + 3, dp)
565 34508 : ppl = fac(2*l + 2)*rootpi/2._dp**REAL(2*l + 3, dp)/fac(l + 1)
566 155708 : DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset)
567 121200 : m_i = basis%m(isgf)
568 775640 : DO jsgf = basis%first_sgf(jshell, jset), basis%last_sgf(jshell, jset)
569 619932 : m_j = basis%m(jsgf)
570 741132 : IF (m_i == m_j) THEN
571 254220 : DO ipgf = 1, basis%npgf(iset)
572 133020 : cci = basis%gcc(ipgf, ishell, iset)
573 133020 : aai = basis%zet(ipgf, iset)
574 133020 : norm_i = norm(isgf)
575 463872 : DO jpgf = 1, basis%npgf(jset)
576 209652 : ccj = basis%gcc(jpgf, jshell, jset)
577 209652 : aaj = basis%zet(jpgf, jset)
578 209652 : oo = 1._dp/(aai + aaj)**expa
579 209652 : norm_j = norm(jsgf)
580 342672 : 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 256 : END SUBROUTINE basis_ovlp
593 :
594 : END MODULE lri_environment_init
|
