Line data Source code
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 Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation
10 : !> of contracted, spherical Gaussian integrals using the solid harmonics. Initialization of
11 : !> the contraction matrices
12 : !> \par History
13 : !> created [08.2016]
14 : !> \author Dorothea Golze
15 : ! **************************************************************************************************
16 : MODULE generic_shg_integrals_init
17 : USE basis_set_types, ONLY: gto_basis_set_type
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: fac,&
20 : ifac,&
21 : pi
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals_init'
29 :
30 : PUBLIC :: contraction_matrix_shg, contraction_matrix_shg_mix, contraction_matrix_shg_rx2m
31 :
32 : ! **************************************************************************************************
33 :
34 : CONTAINS
35 :
36 : ! **************************************************************************************************
37 : !> \brief contraction matrix for SHG integrals
38 : !> \param basis ...
39 : !> \param scon_shg contraction matrix
40 : ! **************************************************************************************************
41 15393246 : SUBROUTINE contraction_matrix_shg(basis, scon_shg)
42 :
43 : TYPE(gto_basis_set_type), POINTER :: basis
44 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_shg
45 :
46 : INTEGER :: ipgf, iset, ishell, l, maxpgf, maxshell, &
47 : nset
48 15393246 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
49 : REAL(KIND=dp) :: aif, gcc
50 15393246 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: norm
51 15393246 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
52 :
53 15393246 : nset = basis%nset
54 15393246 : npgf => basis%npgf
55 15393246 : nshell => basis%nshell
56 15393246 : zet => basis%zet
57 :
58 15393246 : maxpgf = SIZE(basis%gcc, 1)
59 15393246 : maxshell = SIZE(basis%gcc, 2)
60 61572984 : ALLOCATE (norm(basis%nset, maxshell))
61 76966230 : ALLOCATE (scon_shg(maxpgf, maxshell, nset))
62 86076268 : scon_shg = 0.0_dp
63 :
64 15393246 : CALL basis_norm_shg(basis, norm)
65 :
66 38949850 : DO iset = 1, nset
67 62509444 : DO ishell = 1, nshell(iset)
68 23559594 : l = basis%l(ishell, iset)
69 70679416 : DO ipgf = 1, npgf(iset)
70 23563218 : aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l)
71 23563218 : gcc = basis%gcc(ipgf, ishell, iset)
72 47122812 : scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif
73 : END DO
74 : END DO
75 : END DO
76 :
77 15393246 : DEALLOCATE (norm)
78 :
79 15393246 : END SUBROUTINE contraction_matrix_shg
80 :
81 : !***************************************************************************************************
82 : !> \brief normalization solid harmonic Gaussians (SHG)
83 : !> \param basis ...
84 : !> \param norm ...
85 : ! **************************************************************************************************
86 15393546 : SUBROUTINE basis_norm_shg(basis, norm)
87 :
88 : TYPE(gto_basis_set_type), POINTER :: basis
89 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: norm
90 :
91 : INTEGER :: ipgf, iset, ishell, jpgf, l
92 : REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl
93 :
94 54355594 : norm = 0._dp
95 :
96 38952086 : DO iset = 1, basis%nset
97 62516644 : DO ishell = 1, basis%nshell(iset)
98 23564558 : l = basis%l(ishell, iset)
99 23564558 : expa = 0.5_dp*REAL(2*l + 3, dp)
100 23564558 : ppl = fac(2*l + 2)*pi**(1.5_dp)/fac(l + 1)
101 23564558 : ppl = ppl/(2._dp**REAL(2*l + 1, dp))
102 23564558 : ppl = ppl/REAL(2*l + 1, dp)
103 47136364 : DO ipgf = 1, basis%npgf(iset)
104 23571806 : cci = basis%gcc(ipgf, ishell, iset)
105 23571806 : aai = basis%zet(ipgf, iset)
106 70756426 : DO jpgf = 1, basis%npgf(iset)
107 23620062 : ccj = basis%gcc(jpgf, ishell, iset)
108 23620062 : aaj = basis%zet(jpgf, iset)
109 47191868 : norm(iset, ishell) = norm(iset, ishell) + cci*ccj*ppl/(aai + aaj)**expa
110 : END DO
111 : END DO
112 47123098 : norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell))
113 : END DO
114 : END DO
115 :
116 15393546 : END SUBROUTINE basis_norm_shg
117 :
118 : ! **************************************************************************************************
119 : !> \brief mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis
120 : !> at the same atom
121 : !> \param orb_basis orbital basis
122 : !> \param ri_basis ...
123 : !> \param orb_index index for orbital basis
124 : !> \param ri_index index for ri basis
125 : !> \param scon_mix mixed contraction matrix
126 : ! **************************************************************************************************
127 150 : SUBROUTINE contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix)
128 :
129 : TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
130 : INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index
131 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix
132 :
133 : INTEGER :: forb, fri, iil, il, ipgf, iset, ishell, jpgf, jset, jshell, l, l1, l2, lmax_orb, &
134 : lmax_ri, maxpgf_orb, maxpgf_ri, maxshell_orb, maxshell_ri, nf_orb, nf_ri, nl, nl_max, &
135 : nset_orb, nset_ri
136 150 : INTEGER, DIMENSION(:), POINTER :: npgf_orb, npgf_ri, nshell_orb, nshell_ri
137 : REAL(KIND=dp) :: cjf, const, const1, const2, gcc_orb, &
138 : gcc_ri, prefac, scon1, scon2, zet
139 150 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac
140 150 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_orb, norm_ri, zet_orb, zet_ri
141 :
142 150 : nset_orb = orb_basis%nset
143 150 : npgf_orb => orb_basis%npgf
144 150 : nshell_orb => orb_basis%nshell
145 150 : zet_orb => orb_basis%zet
146 150 : nset_ri = ri_basis%nset
147 150 : npgf_ri => ri_basis%npgf
148 150 : nshell_ri => ri_basis%nshell
149 150 : zet_ri => ri_basis%zet
150 :
151 150 : maxpgf_orb = SIZE(orb_basis%gcc, 1)
152 150 : maxshell_orb = SIZE(orb_basis%gcc, 2)
153 600 : ALLOCATE (norm_orb(nset_orb, maxshell_orb))
154 150 : maxpgf_ri = SIZE(ri_basis%gcc, 1)
155 150 : maxshell_ri = SIZE(ri_basis%gcc, 2)
156 600 : ALLOCATE (norm_ri(nset_ri, maxshell_ri))
157 :
158 150 : CALL basis_norm_shg(orb_basis, norm_orb)
159 150 : CALL basis_norm_shg(ri_basis, norm_ri)
160 :
161 750 : ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb))
162 750 : ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri))
163 :
164 : !** index orbital basis set
165 150 : nf_orb = 0
166 344 : DO iset = 1, nset_orb
167 864 : DO ishell = 1, nshell_orb(iset)
168 4094 : DO ipgf = 1, npgf_orb(iset)
169 3380 : nf_orb = nf_orb + 1
170 3900 : orb_index(ipgf, ishell, iset) = nf_orb
171 : END DO
172 : END DO
173 : END DO
174 :
175 : !** index ri basis set
176 : nf_ri = 0
177 1892 : DO iset = 1, nset_ri
178 6336 : DO ishell = 1, nshell_ri(iset)
179 11394 : DO ipgf = 1, npgf_ri(iset)
180 5208 : nf_ri = nf_ri + 1
181 9652 : ri_index(ipgf, ishell, iset) = nf_ri
182 : END DO
183 : END DO
184 : END DO
185 :
186 344 : lmax_orb = MAXVAL(orb_basis%lmax)
187 1892 : lmax_ri = MAXVAL(ri_basis%lmax)
188 150 : nl_max = INT((lmax_orb + lmax_ri)/2) + 1
189 900 : ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max))
190 2043224 : scon_mix = 0.0_dp
191 :
192 450 : ALLOCATE (shg_fac(0:nl_max - 1))
193 150 : shg_fac(0) = 1.0_dp
194 :
195 344 : DO iset = 1, nset_orb
196 864 : DO ishell = 1, nshell_orb(iset)
197 520 : l1 = orb_basis%l(ishell, iset)
198 520 : const1 = SQRT(1.0_dp/REAL(2*l1 + 1, dp))
199 7054 : DO jset = 1, nset_ri
200 23662 : DO jshell = 1, nshell_ri(jset)
201 16802 : l2 = ri_basis%l(jshell, jset)
202 16802 : const2 = SQRT(1.0_dp/REAL(2*l2 + 1, dp))
203 16802 : nl = INT((l1 + l2)/2)
204 16802 : IF (l1 == 0 .OR. l2 == 0) nl = 0
205 48748 : DO il = 0, nl
206 25606 : l = l1 + l2 - 2*il
207 25606 : const = const1*const2*2.0_dp*SQRT(pi*REAL(2*l + 1, dp))
208 37330 : DO iil = 1, il
209 : shg_fac(iil) = fac(l + iil - 1)*ifac(l)*REAL(l, dp) &
210 37330 : *fac(il)/fac(il - iil)/fac(iil)
211 : END DO
212 216144 : DO ipgf = 1, npgf_orb(iset)
213 173736 : forb = orb_index(ipgf, ishell, iset)
214 173736 : gcc_orb = orb_basis%gcc(ipgf, ishell, iset)
215 173736 : scon1 = norm_orb(iset, ishell)*gcc_orb
216 397072 : DO jpgf = 1, npgf_ri(jset)
217 197730 : fri = ri_index(jpgf, jshell, jset)
218 197730 : gcc_ri = ri_basis%gcc(jpgf, jshell, jset)
219 197730 : scon2 = norm_ri(jset, jshell)*gcc_ri
220 197730 : zet = zet_orb(ipgf, iset) + zet_ri(jpgf, jset)
221 197730 : cjf = 1.0_dp/((2._dp*zet)**l)
222 197730 : prefac = const*cjf*scon1*scon2
223 655656 : DO iil = 0, il
224 481920 : scon_mix(iil + 1, fri, forb, il + 1) = prefac*shg_fac(iil)/zet**iil
225 : END DO
226 : END DO
227 : END DO
228 : END DO
229 : END DO
230 : END DO
231 : END DO
232 : END DO
233 :
234 150 : DEALLOCATE (norm_orb, norm_ri, shg_fac)
235 :
236 150 : END SUBROUTINE contraction_matrix_shg_mix
237 :
238 : ! **************************************************************************************************
239 : !> \brief ...
240 : !> \param basis ...
241 : !> \param m ...
242 : !> \param scon_shg ...
243 : !> \param scon_rx2m ...
244 : ! **************************************************************************************************
245 2 : SUBROUTINE contraction_matrix_shg_rx2m(basis, m, scon_shg, scon_rx2m)
246 :
247 : TYPE(gto_basis_set_type), POINTER :: basis
248 : INTEGER, INTENT(IN) :: m
249 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_shg
250 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: scon_rx2m
251 :
252 : INTEGER :: ipgf, iset, ishell, j, l, maxpgf, &
253 : maxshell, nset
254 2 : INTEGER, DIMENSION(:), POINTER :: npgf, nshell
255 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac
256 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
257 :
258 2 : npgf => basis%npgf
259 2 : nshell => basis%nshell
260 2 : zet => basis%zet
261 2 : nset = basis%nset
262 :
263 2 : maxpgf = SIZE(basis%gcc, 1)
264 2 : maxshell = SIZE(basis%gcc, 2)
265 12 : ALLOCATE (scon_rx2m(maxpgf, m + 1, maxshell, nset))
266 2 : scon_rx2m = 0.0_dp
267 6 : ALLOCATE (shg_fac(0:m))
268 2 : shg_fac(0) = 1.0_dp
269 :
270 22 : DO iset = 1, nset
271 88 : DO ishell = 1, nshell(iset)
272 66 : l = basis%l(ishell, iset)
273 264 : DO j = 1, m
274 : shg_fac(j) = fac(l + j - 1)*ifac(l)*REAL(l, dp) &
275 264 : *fac(m)/fac(m - j)/fac(j)
276 : END DO
277 152 : DO ipgf = 1, npgf(iset)
278 396 : DO j = 0, m
279 : scon_rx2m(ipgf, j + 1, ishell, iset) = scon_shg(ipgf, ishell, iset) &
280 330 : *shg_fac(j)/zet(ipgf, iset)**j
281 : END DO
282 : END DO
283 : END DO
284 : END DO
285 :
286 2 : DEALLOCATE (shg_fac)
287 :
288 2 : END SUBROUTINE contraction_matrix_shg_rx2m
289 :
290 : END MODULE generic_shg_integrals_init
|
