(git:34ef472)
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-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  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 
594 END MODULE lri_environment_init
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:212
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::golze2017b
integer, save, public golze2017b
Definition: bibliography.F:43
bibliography::golze2017a
integer, save, public golze2017a
Definition: bibliography.F:43
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:1040
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:1058
lri_environment_types::lri_env_create
subroutine, public lri_env_create(lri_env)
creates and initializes an lri_env
Definition: lri_environment_types.F:390
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451