Line data Source code
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 2- and 3-center electron repulsion integral routines based on libint2
10 : !> Currently available operators: Coulomb, Truncated Coulomb, Short Range (erfc), Overlap
11 : !> \author A. Bussy (05.2019)
12 : ! **************************************************************************************************
13 :
14 : MODULE libint_2c_3c
15 : USE gamma, ONLY: fgamma => fgamma_0
16 : USE input_constants, ONLY: do_potential_coulomb,&
17 : do_potential_id,&
18 : do_potential_long,&
19 : do_potential_mix_cl_trunc,&
20 : do_potential_short,&
21 : do_potential_truncated
22 : USE kinds, ONLY: default_path_length,&
23 : dp
24 : USE libint_wrapper, ONLY: cp_libint_get_2eri_derivs,&
25 : cp_libint_get_2eris,&
26 : cp_libint_get_3eri_derivs,&
27 : cp_libint_get_3eris,&
28 : cp_libint_set_params_eri,&
29 : cp_libint_set_params_eri_deriv,&
30 : cp_libint_t,&
31 : prim_data_f_size
32 : USE mathconstants, ONLY: pi
33 : USE orbital_pointers, ONLY: nco,&
34 : ncoset
35 : USE t_c_g0, ONLY: get_lmax_init,&
36 : t_c_g0_n
37 : #include "./base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_2c_3c'
43 :
44 : PUBLIC :: eri_3center, eri_2center, cutoff_screen_factor, libint_potential_type, &
45 : eri_3center_derivs, eri_2center_derivs, compare_potential_types
46 :
47 : ! For screening of integrals with a truncated potential, it is important to use a slightly larger
48 : ! cutoff radius due to the discontinuity of the truncated Coulomb potential at the cutoff radius.
49 : REAL(KIND=dp), PARAMETER :: cutoff_screen_factor = 1.0001_dp
50 :
51 : TYPE :: params_2c
52 : INTEGER :: m_max = 0
53 : REAL(dp) :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
54 : REAL(dp), DIMENSION(3) :: W = 0.0_dp
55 : REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
56 : END TYPE
57 :
58 : TYPE :: params_3c
59 : INTEGER :: m_max = 0
60 : REAL(dp) :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
61 : REAL(dp), DIMENSION(3) :: Q = 0.0_dp, W = 0.0_dp
62 : REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
63 : END TYPE
64 :
65 : ! A potential type based on the hfx_potential_type concept, but only for implemented
66 : ! 2- and 3-center operators
67 : TYPE :: libint_potential_type
68 : INTEGER :: potential_type = do_potential_coulomb
69 : REAL(dp) :: omega = 0.0_dp !SR: erfc(w*r)/r
70 : REAL(dp) :: cutoff_radius = 0.0_dp !TC cutoff/effective SR range
71 : CHARACTER(default_path_length) :: filename = ""
72 : REAL(dp) :: scale_coulomb = 0.0_dp ! Only, for WFC methods
73 : REAL(dp) :: scale_longrange = 0.0_dp ! Only, for WFC methods
74 : END TYPE
75 :
76 : CONTAINS
77 :
78 : ! **************************************************************************************************
79 : !> \brief Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian
80 : !> gaussian orbitals
81 : !> \param int_abc the integrals as array of cartesian orbitals (allocated before hand)
82 : !> \param la_min ...
83 : !> \param la_max ...
84 : !> \param npgfa ...
85 : !> \param zeta ...
86 : !> \param rpgfa ...
87 : !> \param ra ...
88 : !> \param lb_min ...
89 : !> \param lb_max ...
90 : !> \param npgfb ...
91 : !> \param zetb ...
92 : !> \param rpgfb ...
93 : !> \param rb ...
94 : !> \param lc_min ...
95 : !> \param lc_max ...
96 : !> \param npgfc ...
97 : !> \param zetc ...
98 : !> \param rpgfc ...
99 : !> \param rc ...
100 : !> \param dab ...
101 : !> \param dac ...
102 : !> \param dbc ...
103 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
104 : !> \param potential_parameter the info about the potential
105 : !> \param int_abc_ext the extremal value of int_abc, i.e., MAXVAL(ABS(int_abc))
106 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
107 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
108 : !> the latter must be initialized too
109 : ! **************************************************************************************************
110 2163603 : SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
111 2163603 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
112 2163603 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
113 : dab, dac, dbc, lib, potential_parameter, &
114 : int_abc_ext)
115 :
116 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: int_abc
117 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
118 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
119 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
120 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
121 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
122 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
123 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
124 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
125 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
126 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
127 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
128 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
129 : REAL(dp), INTENT(INOUT), OPTIONAL :: int_abc_ext
130 :
131 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, ipgf, j, &
132 : jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
133 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
134 2163603 : REAL(dp), DIMENSION(:), POINTER :: p_work
135 : TYPE(params_3c), POINTER :: params
136 :
137 2163603 : NULLIFY (params, p_work)
138 64908090 : ALLOCATE (params)
139 :
140 2163603 : dr_ab = 0.0_dp
141 2163603 : dr_bc = 0.0_dp
142 2163603 : dr_ac = 0.0_dp
143 :
144 2163603 : op = potential_parameter%potential_type
145 :
146 : IF (op == do_potential_truncated .OR. op == do_potential_short &
147 2163603 : .OR. op == do_potential_mix_cl_trunc) THEN
148 1347753 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
149 1347753 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
150 815850 : ELSEIF (op == do_potential_coulomb) THEN
151 19291 : dr_bc = 1000000.0_dp
152 19291 : dr_ac = 1000000.0_dp
153 : END IF
154 :
155 2163603 : IF (PRESENT(int_abc_ext)) THEN
156 2127951 : int_abc_ext = 0.0_dp
157 : END IF
158 :
159 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
160 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
161 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
162 :
163 : !Looping over the pgfs
164 6980806 : DO ipgf = 1, npgfa
165 4817203 : zeti = zeta(ipgf)
166 4817203 : a_start = (ipgf - 1)*ncoset(la_max)
167 :
168 20249252 : DO jpgf = 1, npgfb
169 :
170 : ! screening
171 13268446 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
172 :
173 6292252 : zetj = zetb(jpgf)
174 6292252 : b_start = (jpgf - 1)*ncoset(lb_max)
175 :
176 30315371 : DO kpgf = 1, npgfc
177 :
178 : ! screening
179 19205916 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
180 11921910 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
181 :
182 9427201 : zetk = zetc(kpgf)
183 9427201 : c_start = (kpgf - 1)*ncoset(lc_max)
184 :
185 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
186 : CALL set_params_3c(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
187 9427201 : potential_parameter=potential_parameter, params_out=params)
188 :
189 21517482 : DO li = la_min, la_max
190 12090281 : a_offset = a_start + ncoset(li - 1)
191 12090281 : ncoa = nco(li)
192 33828431 : DO lj = MAX(li, lb_min), lb_max
193 12310949 : b_offset = b_start + ncoset(lj - 1)
194 12310949 : ncob = nco(lj)
195 46255097 : DO lk = lc_min, lc_max
196 21853867 : c_offset = c_start + ncoset(lk - 1)
197 21853867 : ncoc = nco(lk)
198 :
199 21853867 : a_mysize(1) = ncoa*ncob*ncoc
200 21853867 : CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
201 :
202 34164816 : IF (PRESENT(int_abc_ext)) THEN
203 86109397 : DO k = 1, ncoc
204 64989989 : p1 = (k - 1)*ncob
205 230224999 : DO j = 1, ncob
206 144115602 : p2 = (p1 + j - 1)*ncoa
207 444475001 : DO i = 1, ncoa
208 235369410 : p3 = p2 + i
209 235369410 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
210 379485012 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
211 : END DO
212 : END DO
213 : END DO
214 : ELSE
215 3219446 : DO k = 1, ncoc
216 2484987 : p1 = (k - 1)*ncob
217 8203101 : DO j = 1, ncob
218 4983655 : p2 = (p1 + j - 1)*ncoa
219 14964443 : DO i = 1, ncoa
220 7495801 : p3 = p2 + i
221 12479456 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
222 : END DO
223 : END DO
224 : END DO
225 : END IF
226 :
227 : END DO !lk
228 : END DO !lj
229 : END DO !li
230 :
231 : !swap centers 3 and 4 to compute (c|ab) with lb < la
232 9427201 : CALL set_params_3c(lib, rb, ra, rc, params_in=params)
233 :
234 34808502 : DO lj = lb_min, lb_max
235 12112855 : b_offset = b_start + ncoset(lj - 1)
236 12112855 : ncob = nco(lj)
237 34935133 : DO li = MAX(lj + 1, la_min), la_max
238 3616362 : a_offset = a_start + ncoset(li - 1)
239 3616362 : ncoa = nco(li)
240 22633106 : DO lk = lc_min, lc_max
241 6903889 : c_offset = c_start + ncoset(lk - 1)
242 6903889 : ncoc = nco(lk)
243 :
244 6903889 : a_mysize(1) = ncoa*ncob*ncoc
245 6903889 : CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
246 :
247 10520251 : IF (PRESENT(int_abc_ext)) THEN
248 27954368 : DO k = 1, ncoc
249 21261697 : p1 = (k - 1)*ncoa
250 108278160 : DO i = 1, ncoa
251 80323792 : p2 = (p1 + i - 1)*ncob
252 204374325 : DO j = 1, ncob
253 102788836 : p3 = p2 + j
254 102788836 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
255 183112628 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
256 : END DO
257 : END DO
258 : END DO
259 : ELSE
260 933928 : DO k = 1, ncoc
261 722710 : p1 = (k - 1)*ncoa
262 3520291 : DO i = 1, ncoa
263 2586363 : p2 = (p1 + i - 1)*ncob
264 6422524 : DO j = 1, ncob
265 3113451 : p3 = p2 + j
266 5699814 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
267 : END DO
268 : END DO
269 : END DO
270 : END IF
271 :
272 : END DO !lk
273 : END DO !li
274 : END DO !lj
275 :
276 : END DO !kpgf
277 : END DO !jpgf
278 : END DO !ipgf
279 :
280 2163603 : DEALLOCATE (params)
281 :
282 2163603 : END SUBROUTINE eri_3center
283 :
284 : ! **************************************************************************************************
285 : !> \brief Sets the internals of the cp_libint_t object for integrals of type (k|ji)
286 : !> \param lib ..
287 : !> \param ri ...
288 : !> \param rj ...
289 : !> \param rk ...
290 : !> \param zeti ...
291 : !> \param zetj ...
292 : !> \param zetk ...
293 : !> \param li_max ...
294 : !> \param lj_max ...
295 : !> \param lk_max ...
296 : !> \param potential_parameter ...
297 : !> \param params_in external parameters to use for libint
298 : !> \param params_out returns the libint parameters computed based on the other arguments
299 : !> \note The use of params_in and params_out comes from the fact that one might have to swap
300 : !> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
301 : !> remain the same upon such a change => might avoid recomputing things over and over again
302 : ! **************************************************************************************************
303 18854402 : SUBROUTINE set_params_3c(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
304 : potential_parameter, params_in, params_out)
305 :
306 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
307 : REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
308 : REAL(dp), INTENT(IN), OPTIONAL :: zeti, zetj, zetk
309 : INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
310 : TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
311 : TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
312 :
313 : INTEGER :: l
314 : LOGICAL :: use_gamma
315 : REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
316 : prefac, R, S1234, T, tmp
317 18854402 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
318 : TYPE(params_3c), POINTER :: params
319 :
320 : !Assume that one of params_in or params_out is present, and that in the latter case, all
321 : !other optional arguments are here
322 :
323 : !The internal structure of libint2 is based on 4-center integrals
324 : !For 3-center, one of those is a dummy center
325 : !The integral is assumed to be (k|ji) where the centers are ordered as:
326 : !k -> 1, j -> 3 and i -> 4 (the center #2 is the dummy center)
327 :
328 : !If external parameters are given, just use them
329 18854402 : IF (PRESENT(params_in)) THEN
330 9427201 : params => params_in
331 :
332 : !If no external parameters to use, compute them
333 : ELSE
334 9427201 : params => params_out
335 :
336 : !Note: some variable of 4-center integrals simplify with a dummy center:
337 : ! P -> rk, gammap -> zetk
338 9427201 : params%m_max = li_max + lj_max + lk_max
339 9427201 : gammaq = zeti + zetj
340 9427201 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
341 9427201 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
342 :
343 37708804 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
344 37708804 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
345 9427201 : params%Rho = zetk*gammaq/(zetk + gammaq)
346 :
347 207398422 : params%Fm = 0.0_dp
348 9427201 : SELECT CASE (potential_parameter%potential_type)
349 : CASE (do_potential_coulomb)
350 1632720 : T = params%Rho*SUM((params%Q - rk)**2)
351 1632720 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
352 408180 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
353 :
354 408180 : CALL fgamma(params%m_max, T, params%Fm)
355 8979960 : params%Fm = prefac*params%Fm
356 : CASE (do_potential_truncated)
357 4002001 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
358 16008004 : T = params%Rho*SUM((params%Q - rk)**2)
359 16008004 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
360 4002001 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
361 :
362 4002001 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
363 4002001 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
364 4002001 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
365 88044022 : params%Fm = prefac*params%Fm
366 : CASE (do_potential_short)
367 8439944 : T = params%Rho*SUM((params%Q - rk)**2)
368 8439944 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
369 2109986 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
370 :
371 2109986 : CALL fgamma(params%m_max, T, params%Fm)
372 :
373 2109986 : omega2 = potential_parameter%omega**2
374 2109986 : omega_corr2 = omega2/(omega2 + params%Rho)
375 2109986 : omega_corr = SQRT(omega_corr2)
376 2109986 : T = T*omega_corr2
377 2109986 : ALLOCATE (Fm(prim_data_f_size))
378 :
379 2109986 : CALL fgamma(params%m_max, T, Fm)
380 2109986 : tmp = -omega_corr
381 11707924 : DO l = 1, params%m_max + 1
382 9597938 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
383 11707924 : tmp = tmp*omega_corr2
384 : END DO
385 46419692 : params%Fm = prefac*params%Fm
386 : CASE (do_potential_mix_cl_trunc)
387 1399165 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
388 5596660 : T = params%Rho*SUM((params%Q - rk)**2)
389 5596660 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
390 1399165 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
391 :
392 1399165 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
393 1399165 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
394 1399165 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
395 :
396 1399165 : ALLOCATE (Fm(prim_data_f_size))
397 1399165 : CALL fgamma(params%m_max, T, Fm)
398 5377685 : DO l = 1, params%m_max + 1
399 : params%Fm(l) = params%Fm(l) &
400 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
401 5377685 : - Fm(l)*potential_parameter%scale_longrange
402 : END DO
403 1399165 : DEALLOCATE (Fm)
404 :
405 1399165 : omega2 = potential_parameter%omega**2
406 1399165 : omega_corr2 = omega2/(omega2 + params%Rho)
407 1399165 : omega_corr = SQRT(omega_corr2)
408 1399165 : T = T*omega_corr2
409 :
410 1399165 : ALLOCATE (Fm(prim_data_f_size))
411 1399165 : CALL fgamma(params%m_max, T, Fm)
412 1399165 : tmp = omega_corr
413 5377685 : DO l = 1, params%m_max + 1
414 3978520 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
415 5377685 : tmp = tmp*omega_corr2
416 : END DO
417 30781630 : params%Fm = prefac*params%Fm
418 : CASE (do_potential_id)
419 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
420 10555083 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
421 1507869 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
422 :
423 33173118 : params%Fm(:) = prefac
424 : CASE DEFAULT
425 9427201 : CPABORT("Requested operator NYI")
426 : END SELECT
427 :
428 : END IF
429 :
430 : CALL cp_libint_set_params_eri(lib, rk, rk, rj, ri, params%ZetaInv, params%EtaInv, &
431 : params%ZetapEtaInv, params%Rho, rk, params%Q, params%W, &
432 18854402 : params%m_max, params%Fm)
433 :
434 18854402 : END SUBROUTINE set_params_3c
435 :
436 : ! **************************************************************************************************
437 : !> \brief Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given
438 : !> set of cartesian gaussian orbitals. Returns x,y,z derivatives for 1st and 2nd center
439 : !> \param der_abc_1 the derivatives for the 1st center (allocated before hand)
440 : !> \param der_abc_2 the derivatives for the 2nd center (allocated before hand)
441 : !> \param la_min ...
442 : !> \param la_max ...
443 : !> \param npgfa ...
444 : !> \param zeta ...
445 : !> \param rpgfa ...
446 : !> \param ra ...
447 : !> \param lb_min ...
448 : !> \param lb_max ...
449 : !> \param npgfb ...
450 : !> \param zetb ...
451 : !> \param rpgfb ...
452 : !> \param rb ...
453 : !> \param lc_min ...
454 : !> \param lc_max ...
455 : !> \param npgfc ...
456 : !> \param zetc ...
457 : !> \param rpgfc ...
458 : !> \param rc ...
459 : !> \param dab ...
460 : !> \param dac ...
461 : !> \param dbc ...
462 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
463 : !> \param potential_parameter the info about the potential
464 : !> \param der_abc_1_ext the extremal value of der_abc_1, i.e., MAXVAL(ABS(der_abc_1))
465 : !> \param der_abc_2_ext ...
466 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
467 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
468 : !> the latter must be initialized too. Note that the derivative wrt to the third center
469 : !> can be obtained via translational invariance
470 : ! **************************************************************************************************
471 342310 : SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
472 342310 : la_min, la_max, npgfa, zeta, rpgfa, ra, &
473 342310 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
474 342310 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
475 : dab, dac, dbc, lib, potential_parameter, &
476 : der_abc_1_ext, der_abc_2_ext)
477 :
478 : REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT) :: der_abc_1, der_abc_2
479 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
480 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
481 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
482 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
483 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
484 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
485 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
486 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
487 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
488 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
489 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
490 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
491 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: der_abc_1_ext, der_abc_2_ext
492 :
493 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, i_deriv, &
494 : ipgf, j, jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
495 : INTEGER, DIMENSION(3) :: permute_1, permute_2
496 : LOGICAL :: do_ext
497 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
498 : REAL(dp), DIMENSION(3) :: der_abc_1_ext_prv, der_abc_2_ext_prv
499 342310 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
500 : TYPE(params_3c), POINTER :: params
501 :
502 342310 : NULLIFY (params, p_deriv)
503 10269300 : ALLOCATE (params)
504 :
505 342310 : permute_1 = [4, 5, 6]
506 342310 : permute_2 = [7, 8, 9]
507 :
508 342310 : dr_ab = 0.0_dp
509 342310 : dr_bc = 0.0_dp
510 342310 : dr_ac = 0.0_dp
511 :
512 342310 : op = potential_parameter%potential_type
513 :
514 : IF (op == do_potential_truncated .OR. op == do_potential_short &
515 342310 : .OR. op == do_potential_mix_cl_trunc) THEN
516 119807 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
517 119807 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
518 222503 : ELSEIF (op == do_potential_coulomb) THEN
519 9988 : dr_bc = 1000000.0_dp
520 9988 : dr_ac = 1000000.0_dp
521 : END IF
522 :
523 342310 : do_ext = .FALSE.
524 342310 : IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
525 342310 : der_abc_1_ext_prv = 0.0_dp
526 342310 : der_abc_2_ext_prv = 0.0_dp
527 :
528 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
529 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
530 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
531 :
532 : !Looping over the pgfs
533 1180747 : DO ipgf = 1, npgfa
534 838437 : zeti = zeta(ipgf)
535 838437 : a_start = (ipgf - 1)*ncoset(la_max)
536 :
537 4746758 : DO jpgf = 1, npgfb
538 :
539 : ! screening
540 3566011 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
541 :
542 1129693 : zetj = zetb(jpgf)
543 1129693 : b_start = (jpgf - 1)*ncoset(lb_max)
544 :
545 8285048 : DO kpgf = 1, npgfc
546 :
547 : ! screening
548 6316918 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
549 2649783 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
550 :
551 1857820 : zetk = zetc(kpgf)
552 1857820 : c_start = (kpgf - 1)*ncoset(lc_max)
553 :
554 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
555 : CALL set_params_3c_deriv(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
556 1857820 : potential_parameter=potential_parameter, params_out=params)
557 :
558 4281757 : DO li = la_min, la_max
559 2423937 : a_offset = a_start + ncoset(li - 1)
560 2423937 : ncoa = nco(li)
561 6872693 : DO lj = MAX(li, lb_min), lb_max
562 2590936 : b_offset = b_start + ncoset(lj - 1)
563 2590936 : ncob = nco(lj)
564 8765957 : DO lk = lc_min, lc_max
565 3751084 : c_offset = c_start + ncoset(lk - 1)
566 3751084 : ncoc = nco(lk)
567 :
568 3751084 : a_mysize(1) = ncoa*ncob*ncoc
569 :
570 3751084 : CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
571 :
572 3751084 : IF (do_ext) THEN
573 15004336 : DO i_deriv = 1, 3
574 42533851 : DO k = 1, ncoc
575 27529515 : p1 = (k - 1)*ncob
576 90660276 : DO j = 1, ncob
577 51877509 : p2 = (p1 + j - 1)*ncoa
578 157726497 : DO i = 1, ncoa
579 78319473 : p3 = p2 + i
580 :
581 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
582 78319473 : p_deriv(p3, permute_2(i_deriv))
583 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
584 78319473 : ABS(p_deriv(p3, permute_2(i_deriv))))
585 :
586 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
587 78319473 : p_deriv(p3, permute_1(i_deriv))
588 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
589 130196982 : ABS(p_deriv(p3, permute_1(i_deriv))))
590 :
591 : END DO
592 : END DO
593 : END DO
594 : END DO
595 : ELSE
596 0 : DO i_deriv = 1, 3
597 0 : DO k = 1, ncoc
598 0 : p1 = (k - 1)*ncob
599 0 : DO j = 1, ncob
600 0 : p2 = (p1 + j - 1)*ncoa
601 0 : DO i = 1, ncoa
602 0 : p3 = p2 + i
603 :
604 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
605 0 : p_deriv(p3, permute_2(i_deriv))
606 :
607 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
608 0 : p_deriv(p3, permute_1(i_deriv))
609 : END DO
610 : END DO
611 : END DO
612 : END DO
613 : END IF
614 :
615 6342020 : DEALLOCATE (p_deriv)
616 : END DO !lk
617 : END DO !lj
618 : END DO !li
619 :
620 : !swap centers 3 and 4 to compute (c|ab) with lb < la
621 1857820 : CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
622 :
623 7856965 : DO lj = lb_min, lb_max
624 2433134 : b_offset = b_start + ncoset(lj - 1)
625 2433134 : ncob = nco(lj)
626 9411651 : DO li = MAX(lj + 1, la_min), la_max
627 661599 : a_offset = a_start + ncoset(li - 1)
628 661599 : ncoa = nco(li)
629 4065070 : DO lk = lc_min, lc_max
630 970337 : c_offset = c_start + ncoset(lk - 1)
631 970337 : ncoc = nco(lk)
632 :
633 970337 : a_mysize(1) = ncoa*ncob*ncoc
634 970337 : CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
635 :
636 970337 : IF (do_ext) THEN
637 3881348 : DO i_deriv = 1, 3
638 11231387 : DO k = 1, ncoc
639 7350039 : p1 = (k - 1)*ncoa
640 34394001 : DO i = 1, ncoa
641 24132951 : p2 = (p1 + i - 1)*ncob
642 59103765 : DO j = 1, ncob
643 27620775 : p3 = p2 + j
644 :
645 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
646 27620775 : p_deriv(p3, permute_1(i_deriv))
647 :
648 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
649 27620775 : ABS(p_deriv(p3, permute_1(i_deriv))))
650 :
651 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
652 27620775 : p_deriv(p3, permute_2(i_deriv))
653 :
654 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
655 51753726 : ABS(p_deriv(p3, permute_2(i_deriv))))
656 : END DO
657 : END DO
658 : END DO
659 : END DO
660 : ELSE
661 0 : DO i_deriv = 1, 3
662 0 : DO k = 1, ncoc
663 0 : p1 = (k - 1)*ncoa
664 0 : DO i = 1, ncoa
665 0 : p2 = (p1 + i - 1)*ncob
666 0 : DO j = 1, ncob
667 0 : p3 = p2 + j
668 :
669 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
670 0 : p_deriv(p3, permute_1(i_deriv))
671 :
672 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
673 0 : p_deriv(p3, permute_2(i_deriv))
674 : END DO
675 : END DO
676 : END DO
677 : END DO
678 : END IF
679 :
680 1631936 : DEALLOCATE (p_deriv)
681 : END DO !lk
682 : END DO !li
683 : END DO !lj
684 :
685 : END DO !kpgf
686 : END DO !jpgf
687 : END DO !ipgf
688 :
689 342310 : IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
690 342310 : IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
691 :
692 342310 : DEALLOCATE (params)
693 :
694 342310 : END SUBROUTINE eri_3center_derivs
695 :
696 : ! **************************************************************************************************
697 : !> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|ji)
698 : !> \param lib ..
699 : !> \param ri ...
700 : !> \param rj ...
701 : !> \param rk ...
702 : !> \param zeti ...
703 : !> \param zetj ...
704 : !> \param zetk ...
705 : !> \param li_max ...
706 : !> \param lj_max ...
707 : !> \param lk_max ...
708 : !> \param potential_parameter ...
709 : !> \param params_in ...
710 : !> \param params_out ...
711 : !> \note The use of params_in and params_out comes from the fact that one might have to swap
712 : !> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
713 : !> remain the same upon such a change => might avoid recomputing things over and over again
714 : ! **************************************************************************************************
715 3715640 : SUBROUTINE set_params_3c_deriv(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
716 : potential_parameter, params_in, params_out)
717 :
718 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
719 : REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
720 : REAL(dp), INTENT(IN) :: zeti, zetj, zetk
721 : INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
722 : TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
723 : TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
724 :
725 : INTEGER :: l
726 : LOGICAL :: use_gamma
727 : REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
728 : prefac, R, S1234, T, tmp
729 3715640 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
730 : TYPE(params_3c), POINTER :: params
731 :
732 3715640 : IF (PRESENT(params_in)) THEN
733 1857820 : params => params_in
734 :
735 : ELSE
736 1857820 : params => params_out
737 :
738 1857820 : params%m_max = li_max + lj_max + lk_max + 1
739 1857820 : gammaq = zeti + zetj
740 1857820 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
741 1857820 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
742 :
743 7431280 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
744 7431280 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
745 1857820 : params%Rho = zetk*gammaq/(zetk + gammaq)
746 :
747 40872040 : params%Fm = 0.0_dp
748 1857820 : SELECT CASE (potential_parameter%potential_type)
749 : CASE (do_potential_coulomb)
750 416424 : T = params%Rho*SUM((params%Q - rk)**2)
751 416424 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
752 104106 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
753 :
754 104106 : CALL fgamma(params%m_max, T, params%Fm)
755 2290332 : params%Fm = prefac*params%Fm
756 : CASE (do_potential_truncated)
757 1244984 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
758 4979936 : T = params%Rho*SUM((params%Q - rk)**2)
759 4979936 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
760 1244984 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
761 :
762 1244984 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
763 1244984 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
764 1244984 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
765 27389648 : params%Fm = prefac*params%Fm
766 : CASE (do_potential_short)
767 0 : T = params%Rho*SUM((params%Q - rk)**2)
768 0 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
769 0 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
770 :
771 0 : CALL fgamma(params%m_max, T, params%Fm)
772 :
773 0 : omega2 = potential_parameter%omega**2
774 0 : omega_corr2 = omega2/(omega2 + params%Rho)
775 0 : omega_corr = SQRT(omega_corr2)
776 0 : T = T*omega_corr2
777 0 : ALLOCATE (Fm(prim_data_f_size))
778 :
779 0 : CALL fgamma(params%m_max, T, Fm)
780 0 : tmp = -omega_corr
781 0 : DO l = 1, params%m_max + 1
782 0 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
783 0 : tmp = tmp*omega_corr2
784 : END DO
785 0 : params%Fm = prefac*params%Fm
786 : CASE (do_potential_mix_cl_trunc)
787 0 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
788 0 : T = params%Rho*SUM((params%Q - rk)**2)
789 0 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
790 0 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
791 :
792 0 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
793 0 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
794 0 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
795 :
796 0 : ALLOCATE (Fm(prim_data_f_size))
797 0 : CALL fgamma(params%m_max, T, Fm)
798 0 : DO l = 1, params%m_max + 1
799 : params%Fm(l) = params%Fm(l) &
800 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
801 0 : - Fm(l)*potential_parameter%scale_longrange
802 : END DO
803 0 : DEALLOCATE (Fm)
804 :
805 0 : omega2 = potential_parameter%omega**2
806 0 : omega_corr2 = omega2/(omega2 + params%Rho)
807 0 : omega_corr = SQRT(omega_corr2)
808 0 : T = T*omega_corr2
809 :
810 0 : ALLOCATE (Fm(prim_data_f_size))
811 0 : CALL fgamma(params%m_max, T, Fm)
812 0 : tmp = omega_corr
813 0 : DO l = 1, params%m_max + 1
814 0 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
815 0 : tmp = tmp*omega_corr2
816 : END DO
817 0 : params%Fm = prefac*params%Fm
818 : CASE (do_potential_id)
819 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
820 3561110 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
821 508730 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
822 :
823 11192060 : params%Fm(:) = prefac
824 : CASE DEFAULT
825 1857820 : CPABORT("Requested operator NYI")
826 : END SELECT
827 :
828 : END IF
829 :
830 : CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, ri, rk, &
831 : params%Q, params%W, zetk, 0.0_dp, zetj, zeti, params%ZetaInv, &
832 3715640 : params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
833 :
834 3715640 : END SUBROUTINE set_params_3c_deriv
835 :
836 : ! **************************************************************************************************
837 : !> \brief Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian
838 : !> gaussian orbitals
839 : !> \param int_ab the integrals as array of cartesian orbitals (allocated before hand)
840 : !> \param la_min ...
841 : !> \param la_max ...
842 : !> \param npgfa ...
843 : !> \param zeta ...
844 : !> \param rpgfa ...
845 : !> \param ra ...
846 : !> \param lb_min ...
847 : !> \param lb_max ...
848 : !> \param npgfb ...
849 : !> \param zetb ...
850 : !> \param rpgfb ...
851 : !> \param rb ...
852 : !> \param dab ...
853 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
854 : !> \param potential_parameter the info about the potential
855 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
856 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
857 : !> the latter must be initialized too
858 : ! **************************************************************************************************
859 398419 : SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
860 398419 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
861 : dab, lib, potential_parameter)
862 :
863 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int_ab
864 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
865 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
866 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
867 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
868 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
869 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
870 : REAL(dp), INTENT(IN) :: dab
871 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
872 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
873 :
874 : INTEGER :: a_mysize(1), a_offset, a_start, &
875 : b_offset, b_start, i, ipgf, j, jpgf, &
876 : li, lj, ncoa, ncob, p1, p2
877 : REAL(dp) :: dr_ab, zeti, zetj
878 398419 : REAL(dp), DIMENSION(:), POINTER :: p_work
879 :
880 398419 : NULLIFY (p_work)
881 :
882 398419 : dr_ab = 0.0_dp
883 :
884 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
885 180884 : potential_parameter%potential_type == do_potential_short .OR. &
886 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
887 230171 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
888 168248 : ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
889 1616 : dr_ab = 1000000.0_dp
890 : END IF
891 :
892 : !Looping over the pgfs
893 1292204 : DO ipgf = 1, npgfa
894 893785 : zeti = zeta(ipgf)
895 893785 : a_start = (ipgf - 1)*ncoset(la_max)
896 :
897 5721224 : DO jpgf = 1, npgfb
898 4429020 : zetj = zetb(jpgf)
899 4429020 : b_start = (jpgf - 1)*ncoset(lb_max)
900 :
901 : !screening
902 4429020 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
903 :
904 2856168 : CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
905 :
906 9051913 : DO li = la_min, la_max
907 5301960 : a_offset = a_start + ncoset(li - 1)
908 5301960 : ncoa = nco(li)
909 20231143 : DO lj = lb_min, lb_max
910 10500163 : b_offset = b_start + ncoset(lj - 1)
911 10500163 : ncob = nco(lj)
912 :
913 10500163 : a_mysize(1) = ncoa*ncob
914 10500163 : CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
915 :
916 47859327 : DO j = 1, ncob
917 32057204 : p1 = (j - 1)*ncoa
918 142254741 : DO i = 1, ncoa
919 99697374 : p2 = p1 + i
920 131754578 : int_ab(a_offset + i, b_offset + j) = p_work(p2)
921 : END DO
922 : END DO
923 :
924 : END DO
925 : END DO
926 :
927 : END DO
928 : END DO
929 :
930 398419 : END SUBROUTINE eri_2center
931 :
932 : ! **************************************************************************************************
933 : !> \brief Sets the internals of the cp_libint_t object for integrals of type (k|j)
934 : !> \param lib ..
935 : !> \param rj ...
936 : !> \param rk ...
937 : !> \param zetj ...
938 : !> \param zetk ...
939 : !> \param lj_max ...
940 : !> \param lk_max ...
941 : !> \param potential_parameter ...
942 : ! **************************************************************************************************
943 2856168 : SUBROUTINE set_params_2c(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
944 :
945 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
946 : REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
947 : REAL(dp), INTENT(IN) :: zetj, zetk
948 : INTEGER, INTENT(IN) :: lj_max, lk_max
949 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
950 :
951 : INTEGER :: l, op
952 : LOGICAL :: use_gamma
953 : REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
954 : R, T, tmp
955 2856168 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
956 : TYPE(params_2c) :: params
957 :
958 : !The internal structure of libint2 is based on 4-center integrals
959 : !For 2-center, two of those are dummy centers
960 : !The integral is assumed to be (k|j) where the centers are ordered as:
961 : !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
962 :
963 : !Note: some variable of 4-center integrals simplify due to dummy centers:
964 : ! P -> rk, gammap -> zetk
965 : ! Q -> rj, gammaq -> zetj
966 :
967 2856168 : op = potential_parameter%potential_type
968 2856168 : params%m_max = lj_max + lk_max
969 2856168 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
970 2856168 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
971 :
972 11424672 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
973 2856168 : params%Rho = zetk*zetj/(zetk + zetj)
974 :
975 62835696 : params%Fm = 0.0_dp
976 : SELECT CASE (op)
977 : CASE (do_potential_coulomb)
978 106992 : T = params%Rho*SUM((rj - rk)**2)
979 26748 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
980 26748 : CALL fgamma(params%m_max, T, params%Fm)
981 588456 : params%Fm = prefac*params%Fm
982 : CASE (do_potential_truncated)
983 168755 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
984 675020 : T = params%Rho*SUM((rj - rk)**2)
985 168755 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
986 :
987 168755 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
988 168755 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
989 168755 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
990 3712610 : params%Fm = prefac*params%Fm
991 : CASE (do_potential_short)
992 10017192 : T = params%Rho*SUM((rj - rk)**2)
993 2504298 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
994 :
995 2504298 : CALL fgamma(params%m_max, T, params%Fm)
996 :
997 2504298 : omega2 = potential_parameter%omega**2
998 2504298 : omega_corr2 = omega2/(omega2 + params%Rho)
999 2504298 : omega_corr = SQRT(omega_corr2)
1000 2504298 : T = T*omega_corr2
1001 2504298 : ALLOCATE (Fm(prim_data_f_size))
1002 :
1003 2504298 : CALL fgamma(params%m_max, T, Fm)
1004 2504298 : tmp = -omega_corr
1005 12305810 : DO l = 1, params%m_max + 1
1006 9801512 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
1007 12305810 : tmp = tmp*omega_corr2
1008 : END DO
1009 55094556 : params%Fm = prefac*params%Fm
1010 : CASE (do_potential_mix_cl_trunc)
1011 61464 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1012 245856 : T = params%Rho*SUM((rj - rk)**2)
1013 61464 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1014 :
1015 61464 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
1016 61464 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1017 61464 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1018 :
1019 61464 : ALLOCATE (Fm(prim_data_f_size))
1020 61464 : CALL fgamma(params%m_max, T, Fm)
1021 223681 : DO l = 1, params%m_max + 1
1022 : params%Fm(l) = params%Fm(l) &
1023 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
1024 223681 : - Fm(l)*potential_parameter%scale_longrange
1025 : END DO
1026 61464 : DEALLOCATE (Fm)
1027 :
1028 61464 : omega2 = potential_parameter%omega**2
1029 61464 : omega_corr2 = omega2/(omega2 + params%Rho)
1030 61464 : omega_corr = SQRT(omega_corr2)
1031 61464 : T = T*omega_corr2
1032 :
1033 61464 : ALLOCATE (Fm(prim_data_f_size))
1034 61464 : CALL fgamma(params%m_max, T, Fm)
1035 61464 : tmp = omega_corr
1036 223681 : DO l = 1, params%m_max + 1
1037 162217 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
1038 223681 : tmp = tmp*omega_corr2
1039 : END DO
1040 1352208 : params%Fm = prefac*params%Fm
1041 : CASE (do_potential_id)
1042 :
1043 379612 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
1044 2087866 : params%Fm(:) = prefac
1045 : CASE DEFAULT
1046 2856168 : CPABORT("Requested operator NYI")
1047 : END SELECT
1048 :
1049 : CALL cp_libint_set_params_eri(lib, rk, rk, rj, rj, params%ZetaInv, params%EtaInv, &
1050 : params%ZetapEtaInv, params%Rho, rk, rj, params%W, &
1051 2856168 : params%m_max, params%Fm)
1052 :
1053 74260368 : END SUBROUTINE set_params_2c
1054 :
1055 : ! **************************************************************************************************
1056 : !> \brief Helper function to compare libint_potential_types
1057 : !> \param potential1 first potential
1058 : !> \param potential2 second potential
1059 : !> \return Boolean whether both potentials are equal
1060 : ! **************************************************************************************************
1061 14370 : PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
1062 : TYPE(libint_potential_type), INTENT(IN) :: potential1, potential2
1063 : LOGICAL :: equals
1064 :
1065 14370 : IF (potential1%potential_type /= potential2%potential_type) THEN
1066 : equals = .FALSE.
1067 : ELSE
1068 13552 : equals = .TRUE.
1069 736 : SELECT CASE (potential1%potential_type)
1070 : CASE (do_potential_short, do_potential_long)
1071 736 : IF (potential1%omega /= potential2%omega) equals = .FALSE.
1072 : CASE (do_potential_truncated)
1073 14 : IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
1074 : CASE (do_potential_mix_cl_trunc)
1075 2 : IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
1076 2 : IF (potential1%omega /= potential2%omega) equals = .FALSE.
1077 2 : IF (potential1%scale_coulomb /= potential2%scale_coulomb) equals = .FALSE.
1078 13554 : IF (potential1%scale_longrange /= potential2%scale_longrange) equals = .FALSE.
1079 : END SELECT
1080 : END IF
1081 :
1082 14370 : END FUNCTION compare_potential_types
1083 :
1084 : !> \brief Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given
1085 : !> set of cartesian gaussian orbitals. Returns the derivatives wrt to the first center
1086 : !> \param der_ab the derivatives as array of cartesian orbitals (allocated before hand)
1087 : !> \param la_min ...
1088 : !> \param la_max ...
1089 : !> \param npgfa ...
1090 : !> \param zeta ...
1091 : !> \param rpgfa ...
1092 : !> \param ra ...
1093 : !> \param lb_min ...
1094 : !> \param lb_max ...
1095 : !> \param npgfb ...
1096 : !> \param zetb ...
1097 : !> \param rpgfb ...
1098 : !> \param rb ...
1099 : !> \param dab ...
1100 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
1101 : !> \param potential_parameter the info about the potential
1102 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
1103 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
1104 : !> the latter must be initialized too
1105 : ! **************************************************************************************************
1106 143159 : SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
1107 143159 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
1108 : dab, lib, potential_parameter)
1109 :
1110 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: der_ab
1111 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
1112 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1113 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
1114 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
1115 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1116 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
1117 : REAL(dp), INTENT(IN) :: dab
1118 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
1119 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1120 :
1121 : INTEGER :: a_mysize(1), a_offset, a_start, &
1122 : b_offset, b_start, i, i_deriv, ipgf, &
1123 : j, jpgf, li, lj, ncoa, ncob, p1, p2
1124 : INTEGER, DIMENSION(3) :: permute
1125 : REAL(dp) :: dr_ab, zeti, zetj
1126 143159 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
1127 :
1128 143159 : NULLIFY (p_deriv)
1129 :
1130 143159 : permute = [4, 5, 6]
1131 :
1132 143159 : dr_ab = 0.0_dp
1133 :
1134 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
1135 79391 : potential_parameter%potential_type == do_potential_short .OR. &
1136 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
1137 64460 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
1138 78699 : ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
1139 9831 : dr_ab = 1000000.0_dp
1140 : END IF
1141 :
1142 : !Looping over the pgfs
1143 635114 : DO ipgf = 1, npgfa
1144 491955 : zeti = zeta(ipgf)
1145 491955 : a_start = (ipgf - 1)*ncoset(la_max)
1146 :
1147 3845505 : DO jpgf = 1, npgfb
1148 3210391 : zetj = zetb(jpgf)
1149 3210391 : b_start = (jpgf - 1)*ncoset(lb_max)
1150 :
1151 : !screening
1152 3210391 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
1153 :
1154 2171616 : CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
1155 :
1156 6307276 : DO li = la_min, la_max
1157 3643705 : a_offset = a_start + ncoset(li - 1)
1158 3643705 : ncoa = nco(li)
1159 13017159 : DO lj = lb_min, lb_max
1160 6163063 : b_offset = b_start + ncoset(lj - 1)
1161 6163063 : ncob = nco(lj)
1162 :
1163 6163063 : a_mysize(1) = ncoa*ncob
1164 6163063 : CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
1165 :
1166 24652252 : DO i_deriv = 1, 3
1167 75297271 : DO j = 1, ncob
1168 50645019 : p1 = (j - 1)*ncoa
1169 208323612 : DO i = 1, ncoa
1170 139189404 : p2 = p1 + i
1171 189834423 : der_ab(a_offset + i, b_offset + j, i_deriv) = p_deriv(p2, permute(i_deriv))
1172 : END DO
1173 : END DO
1174 : END DO
1175 :
1176 9806768 : DEALLOCATE (p_deriv)
1177 : END DO
1178 : END DO
1179 :
1180 : END DO
1181 : END DO
1182 :
1183 143159 : END SUBROUTINE eri_2center_derivs
1184 :
1185 : ! **************************************************************************************************
1186 : !> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|j)
1187 : !> \param lib ..
1188 : !> \param rj ...
1189 : !> \param rk ...
1190 : !> \param zetj ...
1191 : !> \param zetk ...
1192 : !> \param lj_max ...
1193 : !> \param lk_max ...
1194 : !> \param potential_parameter ...
1195 : ! **************************************************************************************************
1196 2171616 : SUBROUTINE set_params_2c_deriv(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
1197 :
1198 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
1199 : REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
1200 : REAL(dp), INTENT(IN) :: zetj, zetk
1201 : INTEGER, INTENT(IN) :: lj_max, lk_max
1202 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1203 :
1204 : INTEGER :: l, op
1205 : LOGICAL :: use_gamma
1206 : REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
1207 : R, T, tmp
1208 2171616 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
1209 : TYPE(params_2c) :: params
1210 :
1211 : !The internal structure of libint2 is based on 4-center integrals
1212 : !For 2-center, two of those are dummy centers
1213 : !The integral is assumed to be (k|j) where the centers are ordered as:
1214 : !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
1215 :
1216 : !Note: some variable of 4-center integrals simplify due to dummy centers:
1217 : ! P -> rk, gammap -> zetk
1218 : ! Q -> rj, gammaq -> zetj
1219 :
1220 2171616 : op = potential_parameter%potential_type
1221 2171616 : params%m_max = lj_max + lk_max + 1
1222 2171616 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
1223 2171616 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
1224 :
1225 8686464 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
1226 2171616 : params%Rho = zetk*zetj/(zetk + zetj)
1227 :
1228 47775552 : params%Fm = 0.0_dp
1229 : SELECT CASE (op)
1230 : CASE (do_potential_coulomb)
1231 102916 : T = params%Rho*SUM((rj - rk)**2)
1232 25729 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1233 25729 : CALL fgamma(params%m_max, T, params%Fm)
1234 566038 : params%Fm = prefac*params%Fm
1235 : CASE (do_potential_truncated)
1236 61941 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1237 247764 : T = params%Rho*SUM((rj - rk)**2)
1238 61941 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1239 :
1240 61941 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
1241 61941 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1242 61941 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1243 1362702 : params%Fm = prefac*params%Fm
1244 : CASE (do_potential_short)
1245 8096600 : T = params%Rho*SUM((rj - rk)**2)
1246 2024150 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1247 :
1248 2024150 : CALL fgamma(params%m_max, T, params%Fm)
1249 :
1250 2024150 : omega2 = potential_parameter%omega**2
1251 2024150 : omega_corr2 = omega2/(omega2 + params%Rho)
1252 2024150 : omega_corr = SQRT(omega_corr2)
1253 2024150 : T = T*omega_corr2
1254 2024150 : ALLOCATE (Fm(prim_data_f_size))
1255 :
1256 2024150 : CALL fgamma(params%m_max, T, Fm)
1257 2024150 : tmp = -omega_corr
1258 11369176 : DO l = 1, params%m_max + 1
1259 9345026 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
1260 11369176 : tmp = tmp*omega_corr2
1261 : END DO
1262 44531300 : params%Fm = prefac*params%Fm
1263 : CASE (do_potential_mix_cl_trunc)
1264 2578 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1265 10312 : T = params%Rho*SUM((rj - rk)**2)
1266 2578 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1267 :
1268 2578 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
1269 2578 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1270 2578 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1271 :
1272 2578 : ALLOCATE (Fm(prim_data_f_size))
1273 2578 : CALL fgamma(params%m_max, T, Fm)
1274 8234 : DO l = 1, params%m_max + 1
1275 : params%Fm(l) = params%Fm(l) &
1276 : *(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) &
1277 8234 : - Fm(l)*potential_parameter%scale_longrange
1278 : END DO
1279 2578 : DEALLOCATE (Fm)
1280 :
1281 2578 : omega2 = potential_parameter%omega**2
1282 2578 : omega_corr2 = omega2/(omega2 + params%Rho)
1283 2578 : omega_corr = SQRT(omega_corr2)
1284 2578 : T = T*omega_corr2
1285 :
1286 2578 : ALLOCATE (Fm(prim_data_f_size))
1287 2578 : CALL fgamma(params%m_max, T, Fm)
1288 2578 : tmp = omega_corr
1289 8234 : DO l = 1, params%m_max + 1
1290 5656 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp*potential_parameter%scale_longrange
1291 8234 : tmp = tmp*omega_corr2
1292 : END DO
1293 56716 : params%Fm = prefac*params%Fm
1294 : CASE (do_potential_id)
1295 :
1296 228872 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
1297 1258796 : params%Fm(:) = prefac
1298 : CASE DEFAULT
1299 2171616 : CPABORT("Requested operator NYI")
1300 : END SELECT
1301 :
1302 : CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, rj, rk, rj, params%W, zetk, 0.0_dp, &
1303 : zetj, 0.0_dp, params%ZetaInv, params%EtaInv, &
1304 : params%ZetapEtaInv, params%Rho, &
1305 2171616 : params%m_max, params%Fm)
1306 :
1307 56462016 : END SUBROUTINE set_params_2c_deriv
1308 :
1309 0 : END MODULE libint_2c_3c
1310 :
|