Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_short,&
20 : do_potential_truncated
21 : USE kinds, ONLY: default_path_length,&
22 : dp
23 : USE libint_wrapper, ONLY: cp_libint_get_2eri_derivs,&
24 : cp_libint_get_2eris,&
25 : cp_libint_get_3eri_derivs,&
26 : cp_libint_get_3eris,&
27 : cp_libint_set_params_eri,&
28 : cp_libint_set_params_eri_deriv,&
29 : cp_libint_t,&
30 : prim_data_f_size
31 : USE mathconstants, ONLY: pi
32 : USE orbital_pointers, ONLY: nco,&
33 : ncoset
34 : USE t_c_g0, ONLY: get_lmax_init,&
35 : t_c_g0_n
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libint_2c_3c'
42 :
43 : PUBLIC :: eri_3center, eri_2center, cutoff_screen_factor, libint_potential_type, &
44 : eri_3center_derivs, eri_2center_derivs, compare_potential_types
45 :
46 : ! For screening of integrals with a truncated potential, it is important to use a slightly larger
47 : ! cutoff radius due to the discontinuity of the truncated Coulomb potential at the cutoff radius.
48 : REAL(KIND=dp), PARAMETER :: cutoff_screen_factor = 1.0001_dp
49 :
50 : TYPE :: params_2c
51 : INTEGER :: m_max = 0
52 : REAL(dp) :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
53 : REAL(dp), DIMENSION(3) :: W = 0.0_dp
54 : REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
55 : END TYPE
56 :
57 : TYPE :: params_3c
58 : INTEGER :: m_max = 0
59 : REAL(dp) :: ZetaInv = 0.0_dp, EtaInv = 0.0_dp, ZetapEtaInv = 0.0_dp, Rho = 0.0_dp
60 : REAL(dp), DIMENSION(3) :: Q = 0.0_dp, W = 0.0_dp
61 : REAL(dp), DIMENSION(prim_data_f_size) :: Fm = 0.0_dp
62 : END TYPE
63 :
64 : ! A potential type based on the hfx_potential_type concept, but only for implemented
65 : ! 2- and 3-center operators
66 : TYPE :: libint_potential_type
67 : INTEGER :: potential_type = do_potential_coulomb
68 : REAL(dp) :: omega = 0.0_dp !SR: erfc(w*r)/r
69 : REAL(dp) :: cutoff_radius = 0.0_dp !TC cutoff/effective SR range
70 : CHARACTER(default_path_length) :: filename = ""
71 : REAL(dp) :: scale_coulomb = 0.0_dp ! Only, for WFC methods
72 : REAL(dp) :: scale_longrange = 0.0_dp ! Only, for WFC methods
73 : END TYPE
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian
79 : !> gaussian orbitals
80 : !> \param int_abc the integrals as array of cartesian orbitals (allocated before hand)
81 : !> \param la_min ...
82 : !> \param la_max ...
83 : !> \param npgfa ...
84 : !> \param zeta ...
85 : !> \param rpgfa ...
86 : !> \param ra ...
87 : !> \param lb_min ...
88 : !> \param lb_max ...
89 : !> \param npgfb ...
90 : !> \param zetb ...
91 : !> \param rpgfb ...
92 : !> \param rb ...
93 : !> \param lc_min ...
94 : !> \param lc_max ...
95 : !> \param npgfc ...
96 : !> \param zetc ...
97 : !> \param rpgfc ...
98 : !> \param rc ...
99 : !> \param dab ...
100 : !> \param dac ...
101 : !> \param dbc ...
102 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
103 : !> \param potential_parameter the info about the potential
104 : !> \param int_abc_ext the extremal value of int_abc, i.e., MAXVAL(ABS(int_abc))
105 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
106 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
107 : !> the latter must be initialized too
108 : ! **************************************************************************************************
109 1951407 : SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
110 1951407 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
111 1951407 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
112 : dab, dac, dbc, lib, potential_parameter, &
113 : int_abc_ext)
114 :
115 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: int_abc
116 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
117 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
118 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
119 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
120 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
121 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
122 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
123 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
124 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
125 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
126 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
127 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
128 : REAL(dp), INTENT(INOUT), OPTIONAL :: int_abc_ext
129 :
130 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, ipgf, j, &
131 : jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
132 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
133 1951407 : REAL(dp), DIMENSION(:), POINTER :: p_work
134 : TYPE(params_3c), POINTER :: params
135 :
136 1951407 : NULLIFY (params, p_work)
137 58542210 : ALLOCATE (params)
138 :
139 1951407 : dr_ab = 0.0_dp
140 1951407 : dr_bc = 0.0_dp
141 1951407 : dr_ac = 0.0_dp
142 :
143 1951407 : op = potential_parameter%potential_type
144 :
145 1951407 : IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
146 1153869 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
147 1153869 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
148 797538 : ELSEIF (op == do_potential_coulomb) THEN
149 16371 : dr_bc = 1000000.0_dp
150 16371 : dr_ac = 1000000.0_dp
151 : END IF
152 :
153 1951407 : IF (PRESENT(int_abc_ext)) THEN
154 1921035 : int_abc_ext = 0.0_dp
155 : END IF
156 :
157 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
158 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
159 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
160 :
161 : !Looping over the pgfs
162 6039027 : DO ipgf = 1, npgfa
163 4087620 : zeti = zeta(ipgf)
164 4087620 : a_start = (ipgf - 1)*ncoset(la_max)
165 :
166 16792346 : DO jpgf = 1, npgfb
167 :
168 : ! screening
169 10753319 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
170 :
171 5615963 : zetj = zetb(jpgf)
172 5615963 : b_start = (jpgf - 1)*ncoset(lb_max)
173 :
174 24740818 : DO kpgf = 1, npgfc
175 :
176 : ! screening
177 15037235 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
178 9795764 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
179 :
180 8028652 : zetk = zetc(kpgf)
181 8028652 : c_start = (kpgf - 1)*ncoset(lc_max)
182 :
183 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
184 : CALL set_params_3c(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
185 8028652 : potential_parameter=potential_parameter, params_out=params)
186 :
187 18504084 : DO li = la_min, la_max
188 10475432 : a_offset = a_start + ncoset(li - 1)
189 10475432 : ncoa = nco(li)
190 29354814 : DO lj = MAX(li, lb_min), lb_max
191 10850730 : b_offset = b_start + ncoset(lj - 1)
192 10850730 : ncob = nco(lj)
193 40482590 : DO lk = lc_min, lc_max
194 19156428 : c_offset = c_start + ncoset(lk - 1)
195 19156428 : ncoc = nco(lk)
196 :
197 19156428 : a_mysize(1) = ncoa*ncob*ncoc
198 19156428 : CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
199 :
200 30007158 : IF (PRESENT(int_abc_ext)) THEN
201 75605016 : DO k = 1, ncoc
202 57081646 : p1 = (k - 1)*ncob
203 204101895 : DO j = 1, ncob
204 128496879 : p2 = (p1 + j - 1)*ncoa
205 397775278 : DO i = 1, ncoa
206 212196753 : p3 = p2 + i
207 212196753 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
208 340693632 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
209 : END DO
210 : END DO
211 : END DO
212 : ELSE
213 2750604 : DO k = 1, ncoc
214 2117546 : p1 = (k - 1)*ncob
215 6941472 : DO j = 1, ncob
216 4190868 : p2 = (p1 + j - 1)*ncoa
217 12548432 : DO i = 1, ncoa
218 6240018 : p3 = p2 + i
219 10430886 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
220 : END DO
221 : END DO
222 : END DO
223 : END IF
224 :
225 : END DO !lk
226 : END DO !lj
227 : END DO !li
228 :
229 : !swap centers 3 and 4 to compute (c|ab) with lb < la
230 8028652 : CALL set_params_3c(lib, rb, ra, rc, params_in=params)
231 :
232 29281319 : DO lj = lb_min, lb_max
233 10499348 : b_offset = b_start + ncoset(lj - 1)
234 10499348 : ncob = nco(lj)
235 28767251 : DO li = MAX(lj + 1, la_min), la_max
236 3230668 : a_offset = a_start + ncoset(li - 1)
237 3230668 : ncoa = nco(li)
238 19863387 : DO lk = lc_min, lc_max
239 6133371 : c_offset = c_start + ncoset(lk - 1)
240 6133371 : ncoc = nco(lk)
241 :
242 6133371 : a_mysize(1) = ncoa*ncob*ncoc
243 6133371 : CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
244 :
245 9364039 : IF (PRESENT(int_abc_ext)) THEN
246 24893413 : DO k = 1, ncoc
247 18939768 : p1 = (k - 1)*ncoa
248 96962009 : DO i = 1, ncoa
249 72068596 : p2 = (p1 + i - 1)*ncob
250 183760388 : DO j = 1, ncob
251 92752024 : p3 = p2 + j
252 92752024 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
253 164820620 : int_abc_ext = MAX(int_abc_ext, ABS(p_work(p3)))
254 : END DO
255 : END DO
256 : END DO
257 : ELSE
258 785845 : DO k = 1, ncoc
259 606119 : p1 = (k - 1)*ncoa
260 2922451 : DO i = 1, ncoa
261 2136606 : p2 = (p1 + i - 1)*ncob
262 5275235 : DO j = 1, ncob
263 2532510 : p3 = p2 + j
264 4669116 : int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
265 : END DO
266 : END DO
267 : END DO
268 : END IF
269 :
270 : END DO !lk
271 : END DO !li
272 : END DO !lj
273 :
274 : END DO !kpgf
275 : END DO !jpgf
276 : END DO !ipgf
277 :
278 1951407 : DEALLOCATE (params)
279 :
280 1951407 : END SUBROUTINE eri_3center
281 :
282 : ! **************************************************************************************************
283 : !> \brief Sets the internals of the cp_libint_t object for integrals of type (k|ji)
284 : !> \param lib ..
285 : !> \param ri ...
286 : !> \param rj ...
287 : !> \param rk ...
288 : !> \param zeti ...
289 : !> \param zetj ...
290 : !> \param zetk ...
291 : !> \param li_max ...
292 : !> \param lj_max ...
293 : !> \param lk_max ...
294 : !> \param potential_parameter ...
295 : !> \param params_in external parameters to use for libint
296 : !> \param params_out returns the libint parameters computed based on the other arguments
297 : !> \note The use of params_in and params_out comes from the fact that one might have to swap
298 : !> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
299 : !> remain the same upon such a change => might avoid recomputing things over and over again
300 : ! **************************************************************************************************
301 16057304 : SUBROUTINE set_params_3c(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
302 : potential_parameter, params_in, params_out)
303 :
304 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
305 : REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
306 : REAL(dp), INTENT(IN), OPTIONAL :: zeti, zetj, zetk
307 : INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
308 : TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
309 : TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
310 :
311 : INTEGER :: l
312 : LOGICAL :: use_gamma
313 : REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
314 : prefac, R, S1234, T, tmp
315 16057304 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
316 : TYPE(params_3c), POINTER :: params
317 :
318 : !Assume that one of params_in or params_out is present, and that in the latter case, all
319 : !other optional arguments are here
320 :
321 : !The internal structure of libint2 is based on 4-center integrals
322 : !For 3-center, one of those is a dummy center
323 : !The integral is assumed to be (k|ji) where the centers are ordered as:
324 : !k -> 1, j -> 3 and i -> 4 (the center #2 is the dummy center)
325 :
326 : !If external parameters are given, just use them
327 16057304 : IF (PRESENT(params_in)) THEN
328 8028652 : params => params_in
329 :
330 : !If no external parameters to use, compute them
331 : ELSE
332 8028652 : params => params_out
333 :
334 : !Note: some variable of 4-center integrals simplify with a dummy center:
335 : ! P -> rk, gammap -> zetk
336 8028652 : params%m_max = li_max + lj_max + lk_max
337 8028652 : gammaq = zeti + zetj
338 8028652 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
339 8028652 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
340 :
341 32114608 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
342 32114608 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
343 8028652 : params%Rho = zetk*gammaq/(zetk + gammaq)
344 :
345 176630344 : params%Fm = 0.0_dp
346 8028652 : SELECT CASE (potential_parameter%potential_type)
347 : CASE (do_potential_coulomb)
348 1455148 : T = params%Rho*SUM((params%Q - rk)**2)
349 1455148 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
350 363787 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
351 :
352 363787 : CALL fgamma(params%m_max, T, params%Fm)
353 8003314 : params%Fm = prefac*params%Fm
354 : CASE (do_potential_truncated)
355 4118195 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
356 16472780 : T = params%Rho*SUM((params%Q - rk)**2)
357 16472780 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
358 4118195 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
359 :
360 4118195 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
361 4118195 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
362 4118195 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
363 90600290 : params%Fm = prefac*params%Fm
364 : CASE (do_potential_short)
365 8439944 : T = params%Rho*SUM((params%Q - rk)**2)
366 8439944 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
367 2109986 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
368 :
369 2109986 : CALL fgamma(params%m_max, T, params%Fm)
370 :
371 2109986 : omega2 = potential_parameter%omega**2
372 2109986 : omega_corr2 = omega2/(omega2 + params%Rho)
373 2109986 : omega_corr = SQRT(omega_corr2)
374 2109986 : T = T*omega_corr2
375 2109986 : ALLOCATE (Fm(prim_data_f_size))
376 :
377 2109986 : CALL fgamma(params%m_max, T, Fm)
378 2109986 : tmp = -omega_corr
379 11707924 : DO l = 1, params%m_max + 1
380 9597938 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
381 11707924 : tmp = tmp*omega_corr2
382 : END DO
383 46419692 : params%Fm = prefac*params%Fm
384 : CASE (do_potential_id)
385 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
386 10056788 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
387 1436684 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
388 :
389 31607048 : params%Fm(:) = prefac
390 : CASE DEFAULT
391 8028652 : CPABORT("Requested operator NYI")
392 : END SELECT
393 :
394 : END IF
395 :
396 : CALL cp_libint_set_params_eri(lib, rk, rk, rj, ri, params%ZetaInv, params%EtaInv, &
397 : params%ZetapEtaInv, params%Rho, rk, params%Q, params%W, &
398 16057304 : params%m_max, params%Fm)
399 :
400 16057304 : END SUBROUTINE set_params_3c
401 :
402 : ! **************************************************************************************************
403 : !> \brief Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given
404 : !> set of cartesian gaussian orbitals. Returns x,y,z derivatives for 1st and 2nd center
405 : !> \param der_abc_1 the derivatives for the 1st center (allocated before hand)
406 : !> \param der_abc_2 the derivatives for the 2nd center (allocated before hand)
407 : !> \param la_min ...
408 : !> \param la_max ...
409 : !> \param npgfa ...
410 : !> \param zeta ...
411 : !> \param rpgfa ...
412 : !> \param ra ...
413 : !> \param lb_min ...
414 : !> \param lb_max ...
415 : !> \param npgfb ...
416 : !> \param zetb ...
417 : !> \param rpgfb ...
418 : !> \param rb ...
419 : !> \param lc_min ...
420 : !> \param lc_max ...
421 : !> \param npgfc ...
422 : !> \param zetc ...
423 : !> \param rpgfc ...
424 : !> \param rc ...
425 : !> \param dab ...
426 : !> \param dac ...
427 : !> \param dbc ...
428 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
429 : !> \param potential_parameter the info about the potential
430 : !> \param der_abc_1_ext the extremal value of der_abc_1, i.e., MAXVAL(ABS(der_abc_1))
431 : !> \param der_abc_2_ext ...
432 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
433 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
434 : !> the latter must be initialized too. Note that the derivative wrt to the third center
435 : !> can be obtained via translational invariance
436 : ! **************************************************************************************************
437 337392 : SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
438 337392 : la_min, la_max, npgfa, zeta, rpgfa, ra, &
439 337392 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
440 337392 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
441 : dab, dac, dbc, lib, potential_parameter, &
442 : der_abc_1_ext, der_abc_2_ext)
443 :
444 : REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT) :: der_abc_1, der_abc_2
445 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
446 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
447 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
448 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
449 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
450 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
451 : INTEGER, INTENT(IN) :: lc_min, lc_max, npgfc
452 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc
453 : REAL(dp), DIMENSION(3), INTENT(IN) :: rc
454 : REAL(KIND=dp), INTENT(IN) :: dab, dac, dbc
455 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
456 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
457 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: der_abc_1_ext, der_abc_2_ext
458 :
459 : INTEGER :: a_mysize(1), a_offset, a_start, b_offset, b_start, c_offset, c_start, i, i_deriv, &
460 : ipgf, j, jpgf, k, kpgf, li, lj, lk, ncoa, ncob, ncoc, op, p1, p2, p3
461 : INTEGER, DIMENSION(3) :: permute_1, permute_2
462 : LOGICAL :: do_ext
463 : REAL(dp) :: dr_ab, dr_ac, dr_bc, zeti, zetj, zetk
464 : REAL(dp), DIMENSION(3) :: der_abc_1_ext_prv, der_abc_2_ext_prv
465 337392 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
466 : TYPE(params_3c), POINTER :: params
467 :
468 337392 : NULLIFY (params, p_deriv)
469 10121760 : ALLOCATE (params)
470 :
471 337392 : permute_1 = [4, 5, 6]
472 337392 : permute_2 = [7, 8, 9]
473 :
474 337392 : dr_ab = 0.0_dp
475 337392 : dr_bc = 0.0_dp
476 337392 : dr_ac = 0.0_dp
477 :
478 337392 : op = potential_parameter%potential_type
479 :
480 337392 : IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
481 119807 : dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
482 119807 : dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
483 217585 : ELSEIF (op == do_potential_coulomb) THEN
484 9988 : dr_bc = 1000000.0_dp
485 9988 : dr_ac = 1000000.0_dp
486 : END IF
487 :
488 337392 : do_ext = .FALSE.
489 337392 : IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .TRUE.
490 337392 : der_abc_1_ext_prv = 0.0_dp
491 337392 : der_abc_2_ext_prv = 0.0_dp
492 :
493 : !Note: we want to compute all possible integrals based on the 3-centers (ab|c) before
494 : ! having to switch to (ba|c) (or the other way around) due to angular momenta in libint
495 : ! For a triplet of centers (k|ji), we can only compute integrals for which lj >= li
496 :
497 : !Looping over the pgfs
498 1168658 : DO ipgf = 1, npgfa
499 831266 : zeti = zeta(ipgf)
500 831266 : a_start = (ipgf - 1)*ncoset(la_max)
501 :
502 4724795 : DO jpgf = 1, npgfb
503 :
504 : ! screening
505 3556137 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
506 :
507 1121923 : zetj = zetb(jpgf)
508 1121923 : b_start = (jpgf - 1)*ncoset(lb_max)
509 :
510 8250173 : DO kpgf = 1, npgfc
511 :
512 : ! screening
513 6296984 : IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) CYCLE
514 2636987 : IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) CYCLE
515 :
516 1847476 : zetk = zetc(kpgf)
517 1847476 : c_start = (kpgf - 1)*ncoset(lc_max)
518 :
519 : !start with all the (c|ba) integrals (standard order) and keep to lb >= la
520 : CALL set_params_3c_deriv(lib, ra, rb, rc, zeti, zetj, zetk, la_max, lb_max, lc_max, &
521 1847476 : potential_parameter=potential_parameter, params_out=params)
522 :
523 4257889 : DO li = la_min, la_max
524 2410413 : a_offset = a_start + ncoset(li - 1)
525 2410413 : ncoa = nco(li)
526 6833201 : DO lj = MAX(li, lb_min), lb_max
527 2575312 : b_offset = b_start + ncoset(lj - 1)
528 2575312 : ncob = nco(lj)
529 8721185 : DO lk = lc_min, lc_max
530 3735460 : c_offset = c_start + ncoset(lk - 1)
531 3735460 : ncoc = nco(lk)
532 :
533 3735460 : a_mysize(1) = ncoa*ncob*ncoc
534 :
535 3735460 : CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
536 :
537 3735460 : IF (do_ext) THEN
538 14941840 : DO i_deriv = 1, 3
539 42343153 : DO k = 1, ncoc
540 27401313 : p1 = (k - 1)*ncob
541 90165768 : DO j = 1, ncob
542 51558075 : p2 = (p1 + j - 1)*ncoa
543 156643779 : DO i = 1, ncoa
544 77684391 : p3 = p2 + i
545 :
546 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
547 77684391 : p_deriv(p3, permute_2(i_deriv))
548 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
549 77684391 : ABS(p_deriv(p3, permute_2(i_deriv))))
550 :
551 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
552 77684391 : p_deriv(p3, permute_1(i_deriv))
553 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
554 129242466 : ABS(p_deriv(p3, permute_1(i_deriv))))
555 :
556 : END DO
557 : END DO
558 : END DO
559 : END DO
560 : ELSE
561 0 : DO i_deriv = 1, 3
562 0 : DO k = 1, ncoc
563 0 : p1 = (k - 1)*ncob
564 0 : DO j = 1, ncob
565 0 : p2 = (p1 + j - 1)*ncoa
566 0 : DO i = 1, ncoa
567 0 : p3 = p2 + i
568 :
569 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
570 0 : p_deriv(p3, permute_2(i_deriv))
571 :
572 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
573 0 : p_deriv(p3, permute_1(i_deriv))
574 : END DO
575 : END DO
576 : END DO
577 : END DO
578 : END IF
579 :
580 6310772 : DEALLOCATE (p_deriv)
581 : END DO !lk
582 : END DO !lj
583 : END DO !li
584 :
585 : !swap centers 3 and 4 to compute (c|ab) with lb < la
586 1847476 : CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
587 :
588 7823151 : DO lj = lb_min, lb_max
589 2419538 : b_offset = b_start + ncoset(lj - 1)
590 2419538 : ncob = nco(lj)
591 9374293 : DO li = MAX(lj + 1, la_min), la_max
592 657771 : a_offset = a_start + ncoset(li - 1)
593 657771 : ncoa = nco(li)
594 4043818 : DO lk = lc_min, lc_max
595 966509 : c_offset = c_start + ncoset(lk - 1)
596 966509 : ncoc = nco(lk)
597 :
598 966509 : a_mysize(1) = ncoa*ncob*ncoc
599 966509 : CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
600 :
601 966509 : IF (do_ext) THEN
602 3866036 : DO i_deriv = 1, 3
603 11172347 : DO k = 1, ncoc
604 7306311 : p1 = (k - 1)*ncoa
605 34168005 : DO i = 1, ncoa
606 23962167 : p2 = (p1 + i - 1)*ncob
607 58639269 : DO j = 1, ncob
608 27370791 : p3 = p2 + j
609 :
610 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
611 27370791 : p_deriv(p3, permute_1(i_deriv))
612 :
613 : der_abc_1_ext_prv(i_deriv) = MAX(der_abc_1_ext_prv(i_deriv), &
614 27370791 : ABS(p_deriv(p3, permute_1(i_deriv))))
615 :
616 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
617 27370791 : p_deriv(p3, permute_2(i_deriv))
618 :
619 : der_abc_2_ext_prv(i_deriv) = MAX(der_abc_2_ext_prv(i_deriv), &
620 51332958 : ABS(p_deriv(p3, permute_2(i_deriv))))
621 : END DO
622 : END DO
623 : END DO
624 : END DO
625 : ELSE
626 0 : DO i_deriv = 1, 3
627 0 : DO k = 1, ncoc
628 0 : p1 = (k - 1)*ncoa
629 0 : DO i = 1, ncoa
630 0 : p2 = (p1 + i - 1)*ncob
631 0 : DO j = 1, ncob
632 0 : p3 = p2 + j
633 :
634 : der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
635 0 : p_deriv(p3, permute_1(i_deriv))
636 :
637 : der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
638 0 : p_deriv(p3, permute_2(i_deriv))
639 : END DO
640 : END DO
641 : END DO
642 : END DO
643 : END IF
644 :
645 1624280 : DEALLOCATE (p_deriv)
646 : END DO !lk
647 : END DO !li
648 : END DO !lj
649 :
650 : END DO !kpgf
651 : END DO !jpgf
652 : END DO !ipgf
653 :
654 337392 : IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
655 337392 : IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
656 :
657 337392 : DEALLOCATE (params)
658 :
659 337392 : END SUBROUTINE eri_3center_derivs
660 :
661 : ! **************************************************************************************************
662 : !> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|ji)
663 : !> \param lib ..
664 : !> \param ri ...
665 : !> \param rj ...
666 : !> \param rk ...
667 : !> \param zeti ...
668 : !> \param zetj ...
669 : !> \param zetk ...
670 : !> \param li_max ...
671 : !> \param lj_max ...
672 : !> \param lk_max ...
673 : !> \param potential_parameter ...
674 : !> \param params_in ...
675 : !> \param params_out ...
676 : !> \note The use of params_in and params_out comes from the fact that one might have to swap
677 : !> centers 3 and 4 because of angular momenta and pretty much all the parameters of libint
678 : !> remain the same upon such a change => might avoid recomputing things over and over again
679 : ! **************************************************************************************************
680 3694952 : SUBROUTINE set_params_3c_deriv(lib, ri, rj, rk, zeti, zetj, zetk, li_max, lj_max, lk_max, &
681 : potential_parameter, params_in, params_out)
682 :
683 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
684 : REAL(dp), DIMENSION(3), INTENT(IN) :: ri, rj, rk
685 : REAL(dp), INTENT(IN) :: zeti, zetj, zetk
686 : INTEGER, INTENT(IN), OPTIONAL :: li_max, lj_max, lk_max
687 : TYPE(libint_potential_type), INTENT(IN), OPTIONAL :: potential_parameter
688 : TYPE(params_3c), OPTIONAL, POINTER :: params_in, params_out
689 :
690 : INTEGER :: l
691 : LOGICAL :: use_gamma
692 : REAL(dp) :: gammaq, omega2, omega_corr, omega_corr2, &
693 : prefac, R, S1234, T, tmp
694 3694952 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
695 : TYPE(params_3c), POINTER :: params
696 :
697 3694952 : IF (PRESENT(params_in)) THEN
698 1847476 : params => params_in
699 :
700 : ELSE
701 1847476 : params => params_out
702 :
703 1847476 : params%m_max = li_max + lj_max + lk_max + 1
704 1847476 : gammaq = zeti + zetj
705 1847476 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
706 1847476 : params%ZetapEtaInv = 1._dp/(zetk + gammaq)
707 :
708 7389904 : params%Q = (zeti*ri + zetj*rj)*params%EtaInv
709 7389904 : params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
710 1847476 : params%Rho = zetk*gammaq/(zetk + gammaq)
711 :
712 40644472 : params%Fm = 0.0_dp
713 1847476 : SELECT CASE (potential_parameter%potential_type)
714 : CASE (do_potential_coulomb)
715 416424 : T = params%Rho*SUM((params%Q - rk)**2)
716 416424 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
717 104106 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
718 :
719 104106 : CALL fgamma(params%m_max, T, params%Fm)
720 2290332 : params%Fm = prefac*params%Fm
721 : CASE (do_potential_truncated)
722 1244984 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
723 4979936 : T = params%Rho*SUM((params%Q - rk)**2)
724 4979936 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
725 1244984 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
726 :
727 1244984 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
728 1244984 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
729 1244984 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
730 27389648 : params%Fm = prefac*params%Fm
731 : CASE (do_potential_short)
732 0 : T = params%Rho*SUM((params%Q - rk)**2)
733 0 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2))
734 0 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)*S1234
735 :
736 0 : CALL fgamma(params%m_max, T, params%Fm)
737 :
738 0 : omega2 = potential_parameter%omega**2
739 0 : omega_corr2 = omega2/(omega2 + params%Rho)
740 0 : omega_corr = SQRT(omega_corr2)
741 0 : T = T*omega_corr2
742 0 : ALLOCATE (Fm(prim_data_f_size))
743 :
744 0 : CALL fgamma(params%m_max, T, Fm)
745 0 : tmp = -omega_corr
746 0 : DO l = 1, params%m_max + 1
747 0 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
748 0 : tmp = tmp*omega_corr2
749 : END DO
750 0 : params%Fm = prefac*params%Fm
751 : CASE (do_potential_id)
752 : S1234 = EXP(-zeti*zetj*params%EtaInv*SUM((rj - ri)**2) &
753 3488702 : - gammaq*zetk*params%ZetapEtaInv*SUM((params%Q - rk)**2))
754 498386 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*S1234
755 :
756 10964492 : params%Fm(:) = prefac
757 : CASE DEFAULT
758 1847476 : CPABORT("Requested operator NYI")
759 : END SELECT
760 :
761 : END IF
762 :
763 : CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, ri, rk, &
764 : params%Q, params%W, zetk, 0.0_dp, zetj, zeti, params%ZetaInv, &
765 3694952 : params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
766 :
767 3694952 : END SUBROUTINE set_params_3c_deriv
768 :
769 : ! **************************************************************************************************
770 : !> \brief Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian
771 : !> gaussian orbitals
772 : !> \param int_ab the integrals as array of cartesian orbitals (allocated before hand)
773 : !> \param la_min ...
774 : !> \param la_max ...
775 : !> \param npgfa ...
776 : !> \param zeta ...
777 : !> \param rpgfa ...
778 : !> \param ra ...
779 : !> \param lb_min ...
780 : !> \param lb_max ...
781 : !> \param npgfb ...
782 : !> \param zetb ...
783 : !> \param rpgfb ...
784 : !> \param rb ...
785 : !> \param dab ...
786 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
787 : !> \param potential_parameter the info about the potential
788 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
789 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
790 : !> the latter must be initialized too
791 : ! **************************************************************************************************
792 384757 : SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
793 384757 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
794 : dab, lib, potential_parameter)
795 :
796 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int_ab
797 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
798 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
799 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
800 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
801 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
802 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
803 : REAL(dp), INTENT(IN) :: dab
804 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
805 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
806 :
807 : INTEGER :: a_mysize(1), a_offset, a_start, &
808 : b_offset, b_start, i, ipgf, j, jpgf, &
809 : li, lj, ncoa, ncob, p1, p2
810 : REAL(dp) :: dr_ab, zeti, zetj
811 384757 : REAL(dp), DIMENSION(:), POINTER :: p_work
812 :
813 384757 : NULLIFY (p_work)
814 :
815 384757 : dr_ab = 0.0_dp
816 :
817 384757 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
818 : potential_parameter%potential_type == do_potential_short) THEN
819 220087 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
820 164670 : ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
821 1584 : dr_ab = 1000000.0_dp
822 : END IF
823 :
824 : !Looping over the pgfs
825 1223368 : DO ipgf = 1, npgfa
826 838611 : zeti = zeta(ipgf)
827 838611 : a_start = (ipgf - 1)*ncoset(la_max)
828 :
829 5430244 : DO jpgf = 1, npgfb
830 4206876 : zetj = zetb(jpgf)
831 4206876 : b_start = (jpgf - 1)*ncoset(lb_max)
832 :
833 : !screening
834 4206876 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
835 :
836 2788540 : CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
837 :
838 8820531 : DO li = la_min, la_max
839 5193380 : a_offset = a_start + ncoset(li - 1)
840 5193380 : ncoa = nco(li)
841 19716184 : DO lj = lb_min, lb_max
842 10315928 : b_offset = b_start + ncoset(lj - 1)
843 10315928 : ncob = nco(lj)
844 :
845 10315928 : a_mysize(1) = ncoa*ncob
846 10315928 : CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
847 :
848 47036899 : DO j = 1, ncob
849 31527591 : p1 = (j - 1)*ncoa
850 139963128 : DO i = 1, ncoa
851 98119609 : p2 = p1 + i
852 129647200 : int_ab(a_offset + i, b_offset + j) = p_work(p2)
853 : END DO
854 : END DO
855 :
856 : END DO
857 : END DO
858 :
859 : END DO
860 : END DO
861 :
862 384757 : END SUBROUTINE eri_2center
863 :
864 : ! **************************************************************************************************
865 : !> \brief Sets the internals of the cp_libint_t object for integrals of type (k|j)
866 : !> \param lib ..
867 : !> \param rj ...
868 : !> \param rk ...
869 : !> \param zetj ...
870 : !> \param zetk ...
871 : !> \param lj_max ...
872 : !> \param lk_max ...
873 : !> \param potential_parameter ...
874 : ! **************************************************************************************************
875 2788540 : SUBROUTINE set_params_2c(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
876 :
877 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
878 : REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
879 : REAL(dp), INTENT(IN) :: zetj, zetk
880 : INTEGER, INTENT(IN) :: lj_max, lk_max
881 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
882 :
883 : INTEGER :: l, op
884 : LOGICAL :: use_gamma
885 : REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
886 : R, T, tmp
887 2788540 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
888 : TYPE(params_2c) :: params
889 :
890 : !The internal structure of libint2 is based on 4-center integrals
891 : !For 2-center, two of those are dummy centers
892 : !The integral is assumed to be (k|j) where the centers are ordered as:
893 : !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
894 :
895 : !Note: some variable of 4-center integrals simplify due to dummy centers:
896 : ! P -> rk, gammap -> zetk
897 : ! Q -> rj, gammaq -> zetj
898 :
899 2788540 : op = potential_parameter%potential_type
900 2788540 : params%m_max = lj_max + lk_max
901 2788540 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
902 2788540 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
903 :
904 11154160 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
905 2788540 : params%Rho = zetk*zetj/(zetk + zetj)
906 :
907 61347880 : params%Fm = 0.0_dp
908 : SELECT CASE (op)
909 : CASE (do_potential_coulomb)
910 97168 : T = params%Rho*SUM((rj - rk)**2)
911 24292 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
912 24292 : CALL fgamma(params%m_max, T, params%Fm)
913 534424 : params%Fm = prefac*params%Fm
914 : CASE (do_potential_truncated)
915 169427 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
916 677708 : T = params%Rho*SUM((rj - rk)**2)
917 169427 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
918 :
919 169427 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
920 169427 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
921 169427 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
922 3727394 : params%Fm = prefac*params%Fm
923 : CASE (do_potential_short)
924 10017192 : T = params%Rho*SUM((rj - rk)**2)
925 2504298 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
926 :
927 2504298 : CALL fgamma(params%m_max, T, params%Fm)
928 :
929 2504298 : omega2 = potential_parameter%omega**2
930 2504298 : omega_corr2 = omega2/(omega2 + params%Rho)
931 2504298 : omega_corr = SQRT(omega_corr2)
932 2504298 : T = T*omega_corr2
933 2504298 : ALLOCATE (Fm(prim_data_f_size))
934 :
935 2504298 : CALL fgamma(params%m_max, T, Fm)
936 2504298 : tmp = -omega_corr
937 12305810 : DO l = 1, params%m_max + 1
938 9801512 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
939 12305810 : tmp = tmp*omega_corr2
940 : END DO
941 55094556 : params%Fm = prefac*params%Fm
942 : CASE (do_potential_id)
943 :
944 362092 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
945 1991506 : params%Fm(:) = prefac
946 : CASE DEFAULT
947 2788540 : CPABORT("Requested operator NYI")
948 : END SELECT
949 :
950 : CALL cp_libint_set_params_eri(lib, rk, rk, rj, rj, params%ZetaInv, params%EtaInv, &
951 : params%ZetapEtaInv, params%Rho, rk, rj, params%W, &
952 2788540 : params%m_max, params%Fm)
953 :
954 72502040 : END SUBROUTINE set_params_2c
955 :
956 : ! **************************************************************************************************
957 : !> \brief Helper function to compare libint_potential_types
958 : !> \param potential1 first potential
959 : !> \param potential2 second potential
960 : !> \return Boolean whether both potentials are equal
961 : ! **************************************************************************************************
962 13828 : PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
963 : TYPE(libint_potential_type), INTENT(IN) :: potential1, potential2
964 : LOGICAL :: equals
965 :
966 13828 : IF (potential1%potential_type /= potential2%potential_type) THEN
967 : equals = .FALSE.
968 : ELSE
969 13052 : equals = .TRUE.
970 734 : SELECT CASE (potential1%potential_type)
971 : CASE (do_potential_short, do_potential_long)
972 734 : IF (potential1%omega /= potential2%omega) equals = .FALSE.
973 : CASE (do_potential_truncated)
974 13052 : IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .FALSE.
975 : END SELECT
976 : END IF
977 :
978 13828 : END FUNCTION compare_potential_types
979 :
980 : !> \brief Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given
981 : !> set of cartesian gaussian orbitals. Returns the derivatives wrt to the first center
982 : !> \param der_ab the derivatives as array of cartesian orbitals (allocated before hand)
983 : !> \param la_min ...
984 : !> \param la_max ...
985 : !> \param npgfa ...
986 : !> \param zeta ...
987 : !> \param rpgfa ...
988 : !> \param ra ...
989 : !> \param lb_min ...
990 : !> \param lb_max ...
991 : !> \param npgfb ...
992 : !> \param zetb ...
993 : !> \param rpgfb ...
994 : !> \param rb ...
995 : !> \param dab ...
996 : !> \param lib the libint_t object for evaluation (assume that it is initialized outside)
997 : !> \param potential_parameter the info about the potential
998 : !> \note Prior to calling this routine, the cp_libint_t type passed as argument must be initialized,
999 : !> the libint library must be static initialized, and in case of truncated Coulomb operator,
1000 : !> the latter must be initialized too
1001 : ! **************************************************************************************************
1002 137567 : SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
1003 137567 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
1004 : dab, lib, potential_parameter)
1005 :
1006 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: der_ab
1007 : INTEGER, INTENT(IN) :: la_min, la_max, npgfa
1008 : REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa
1009 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra
1010 : INTEGER, INTENT(IN) :: lb_min, lb_max, npgfb
1011 : REAL(dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb
1012 : REAL(dp), DIMENSION(3), INTENT(IN) :: rb
1013 : REAL(dp), INTENT(IN) :: dab
1014 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
1015 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1016 :
1017 : INTEGER :: a_mysize(1), a_offset, a_start, &
1018 : b_offset, b_start, i, i_deriv, ipgf, &
1019 : j, jpgf, li, lj, ncoa, ncob, p1, p2
1020 : INTEGER, DIMENSION(3) :: permute
1021 : REAL(dp) :: dr_ab, zeti, zetj
1022 137567 : REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
1023 :
1024 137567 : NULLIFY (p_deriv)
1025 :
1026 137567 : permute = [4, 5, 6]
1027 :
1028 137567 : dr_ab = 0.0_dp
1029 :
1030 137567 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
1031 : potential_parameter%potential_type == do_potential_short) THEN
1032 61518 : dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
1033 76049 : ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
1034 9831 : dr_ab = 1000000.0_dp
1035 : END IF
1036 :
1037 : !Looping over the pgfs
1038 620990 : DO ipgf = 1, npgfa
1039 483423 : zeti = zeta(ipgf)
1040 483423 : a_start = (ipgf - 1)*ncoset(la_max)
1041 :
1042 3811929 : DO jpgf = 1, npgfb
1043 3190939 : zetj = zetb(jpgf)
1044 3190939 : b_start = (jpgf - 1)*ncoset(lb_max)
1045 :
1046 : !screening
1047 3190939 : IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) CYCLE
1048 :
1049 2166116 : CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
1050 :
1051 6287744 : DO li = la_min, la_max
1052 3638205 : a_offset = a_start + ncoset(li - 1)
1053 3638205 : ncoa = nco(li)
1054 12986707 : DO lj = lb_min, lb_max
1055 6157563 : b_offset = b_start + ncoset(lj - 1)
1056 6157563 : ncob = nco(lj)
1057 :
1058 6157563 : a_mysize(1) = ncoa*ncob
1059 6157563 : CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
1060 :
1061 24630252 : DO i_deriv = 1, 3
1062 75240423 : DO j = 1, ncob
1063 50610171 : p1 = (j - 1)*ncoa
1064 208174620 : DO i = 1, ncoa
1065 139091760 : p2 = p1 + i
1066 189701931 : der_ab(a_offset + i, b_offset + j, i_deriv) = p_deriv(p2, permute(i_deriv))
1067 : END DO
1068 : END DO
1069 : END DO
1070 :
1071 9795768 : DEALLOCATE (p_deriv)
1072 : END DO
1073 : END DO
1074 :
1075 : END DO
1076 : END DO
1077 :
1078 137567 : END SUBROUTINE eri_2center_derivs
1079 :
1080 : ! **************************************************************************************************
1081 : !> \brief Sets the internals of the cp_libint_t object for derivatives of integrals of type (k|j)
1082 : !> \param lib ..
1083 : !> \param rj ...
1084 : !> \param rk ...
1085 : !> \param zetj ...
1086 : !> \param zetk ...
1087 : !> \param lj_max ...
1088 : !> \param lk_max ...
1089 : !> \param potential_parameter ...
1090 : ! **************************************************************************************************
1091 2166116 : SUBROUTINE set_params_2c_deriv(lib, rj, rk, zetj, zetk, lj_max, lk_max, potential_parameter)
1092 :
1093 : TYPE(cp_libint_t), INTENT(INOUT) :: lib
1094 : REAL(dp), DIMENSION(3), INTENT(IN) :: rj, rk
1095 : REAL(dp), INTENT(IN) :: zetj, zetk
1096 : INTEGER, INTENT(IN) :: lj_max, lk_max
1097 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1098 :
1099 : INTEGER :: l, op
1100 : LOGICAL :: use_gamma
1101 : REAL(dp) :: omega2, omega_corr, omega_corr2, prefac, &
1102 : R, T, tmp
1103 2166116 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Fm
1104 : TYPE(params_2c) :: params
1105 :
1106 : !The internal structure of libint2 is based on 4-center integrals
1107 : !For 2-center, two of those are dummy centers
1108 : !The integral is assumed to be (k|j) where the centers are ordered as:
1109 : !k -> 1, j -> 3 and (the centers #2 & #4 are dummy centers)
1110 :
1111 : !Note: some variable of 4-center integrals simplify due to dummy centers:
1112 : ! P -> rk, gammap -> zetk
1113 : ! Q -> rj, gammaq -> zetj
1114 :
1115 2166116 : op = potential_parameter%potential_type
1116 2166116 : params%m_max = lj_max + lk_max + 1
1117 2166116 : params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
1118 2166116 : params%ZetapEtaInv = 1._dp/(zetk + zetj)
1119 :
1120 8664464 : params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
1121 2166116 : params%Rho = zetk*zetj/(zetk + zetj)
1122 :
1123 47654552 : params%Fm = 0.0_dp
1124 : SELECT CASE (op)
1125 : CASE (do_potential_coulomb)
1126 102916 : T = params%Rho*SUM((rj - rk)**2)
1127 25729 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1128 25729 : CALL fgamma(params%m_max, T, params%Fm)
1129 566038 : params%Fm = prefac*params%Fm
1130 : CASE (do_potential_truncated)
1131 60843 : R = potential_parameter%cutoff_radius*SQRT(params%Rho)
1132 243372 : T = params%Rho*SUM((rj - rk)**2)
1133 60843 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1134 :
1135 60843 : CPASSERT(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
1136 60843 : CALL t_c_g0_n(params%Fm, use_gamma, R, T, params%m_max)
1137 60843 : IF (use_gamma) CALL fgamma(params%m_max, T, params%Fm)
1138 1338546 : params%Fm = prefac*params%Fm
1139 : CASE (do_potential_short)
1140 8096600 : T = params%Rho*SUM((rj - rk)**2)
1141 2024150 : prefac = 2._dp*pi/params%Rho*SQRT((pi*params%ZetapEtaInv)**3)
1142 :
1143 2024150 : CALL fgamma(params%m_max, T, params%Fm)
1144 :
1145 2024150 : omega2 = potential_parameter%omega**2
1146 2024150 : omega_corr2 = omega2/(omega2 + params%Rho)
1147 2024150 : omega_corr = SQRT(omega_corr2)
1148 2024150 : T = T*omega_corr2
1149 2024150 : ALLOCATE (Fm(prim_data_f_size))
1150 :
1151 2024150 : CALL fgamma(params%m_max, T, Fm)
1152 2024150 : tmp = -omega_corr
1153 11369176 : DO l = 1, params%m_max + 1
1154 9345026 : params%Fm(l) = params%Fm(l) + Fm(l)*tmp
1155 11369176 : tmp = tmp*omega_corr2
1156 : END DO
1157 44531300 : params%Fm = prefac*params%Fm
1158 : CASE (do_potential_id)
1159 :
1160 221576 : prefac = SQRT((pi*params%ZetapEtaInv)**3)*EXP(-zetj*zetk*params%ZetapEtaInv*SUM((rk - rj)**2))
1161 1218668 : params%Fm(:) = prefac
1162 : CASE DEFAULT
1163 2166116 : CPABORT("Requested operator NYI")
1164 : END SELECT
1165 :
1166 : CALL cp_libint_set_params_eri_deriv(lib, rk, rk, rj, rj, rk, rj, params%W, zetk, 0.0_dp, &
1167 : zetj, 0.0_dp, params%ZetaInv, params%EtaInv, &
1168 : params%ZetapEtaInv, params%Rho, &
1169 2166116 : params%m_max, params%Fm)
1170 :
1171 56319016 : END SUBROUTINE set_params_2c_deriv
1172 :
1173 0 : END MODULE libint_2c_3c
1174 :
|
