(git:374b731)
Loading...
Searching...
No Matches
libint_2c_3c.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief 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
14MODULE 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
75CONTAINS
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 SUBROUTINE eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
110 lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
111 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 REAL(dp), DIMENSION(:), POINTER :: p_work
134 TYPE(params_3c), POINTER :: params
135
136 NULLIFY (params, p_work)
137 ALLOCATE (params)
138
139 dr_ab = 0.0_dp
140 dr_bc = 0.0_dp
141 dr_ac = 0.0_dp
142
143 op = potential_parameter%potential_type
144
145 IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
146 dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
147 dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
148 ELSEIF (op == do_potential_coulomb) THEN
149 dr_bc = 1000000.0_dp
150 dr_ac = 1000000.0_dp
151 END IF
152
153 IF (PRESENT(int_abc_ext)) THEN
154 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 DO ipgf = 1, npgfa
163 zeti = zeta(ipgf)
164 a_start = (ipgf - 1)*ncoset(la_max)
165
166 DO jpgf = 1, npgfb
167
168 ! screening
169 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
170
171 zetj = zetb(jpgf)
172 b_start = (jpgf - 1)*ncoset(lb_max)
173
174 DO kpgf = 1, npgfc
175
176 ! screening
177 IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) cycle
178 IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) cycle
179
180 zetk = zetc(kpgf)
181 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 potential_parameter=potential_parameter, params_out=params)
186
187 DO li = la_min, la_max
188 a_offset = a_start + ncoset(li - 1)
189 ncoa = nco(li)
190 DO lj = max(li, lb_min), lb_max
191 b_offset = b_start + ncoset(lj - 1)
192 ncob = nco(lj)
193 DO lk = lc_min, lc_max
194 c_offset = c_start + ncoset(lk - 1)
195 ncoc = nco(lk)
196
197 a_mysize(1) = ncoa*ncob*ncoc
198 CALL cp_libint_get_3eris(li, lj, lk, lib, p_work, a_mysize)
199
200 IF (PRESENT(int_abc_ext)) THEN
201 DO k = 1, ncoc
202 p1 = (k - 1)*ncob
203 DO j = 1, ncob
204 p2 = (p1 + j - 1)*ncoa
205 DO i = 1, ncoa
206 p3 = p2 + i
207 int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
208 int_abc_ext = max(int_abc_ext, abs(p_work(p3)))
209 END DO
210 END DO
211 END DO
212 ELSE
213 DO k = 1, ncoc
214 p1 = (k - 1)*ncob
215 DO j = 1, ncob
216 p2 = (p1 + j - 1)*ncoa
217 DO i = 1, ncoa
218 p3 = p2 + i
219 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 CALL set_params_3c(lib, rb, ra, rc, params_in=params)
231
232 DO lj = lb_min, lb_max
233 b_offset = b_start + ncoset(lj - 1)
234 ncob = nco(lj)
235 DO li = max(lj + 1, la_min), la_max
236 a_offset = a_start + ncoset(li - 1)
237 ncoa = nco(li)
238 DO lk = lc_min, lc_max
239 c_offset = c_start + ncoset(lk - 1)
240 ncoc = nco(lk)
241
242 a_mysize(1) = ncoa*ncob*ncoc
243 CALL cp_libint_get_3eris(lj, li, lk, lib, p_work, a_mysize)
244
245 IF (PRESENT(int_abc_ext)) THEN
246 DO k = 1, ncoc
247 p1 = (k - 1)*ncoa
248 DO i = 1, ncoa
249 p2 = (p1 + i - 1)*ncob
250 DO j = 1, ncob
251 p3 = p2 + j
252 int_abc(a_offset + i, b_offset + j, c_offset + k) = p_work(p3)
253 int_abc_ext = max(int_abc_ext, abs(p_work(p3)))
254 END DO
255 END DO
256 END DO
257 ELSE
258 DO k = 1, ncoc
259 p1 = (k - 1)*ncoa
260 DO i = 1, ncoa
261 p2 = (p1 + i - 1)*ncob
262 DO j = 1, ncob
263 p3 = p2 + j
264 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 DEALLOCATE (params)
279
280 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 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 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 IF (PRESENT(params_in)) THEN
328 params => params_in
329
330 !If no external parameters to use, compute them
331 ELSE
332 params => params_out
333
334 !Note: some variable of 4-center integrals simplify with a dummy center:
335 ! P -> rk, gammap -> zetk
336 params%m_max = li_max + lj_max + lk_max
337 gammaq = zeti + zetj
338 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
339 params%ZetapEtaInv = 1._dp/(zetk + gammaq)
340
341 params%Q = (zeti*ri + zetj*rj)*params%EtaInv
342 params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
343 params%Rho = zetk*gammaq/(zetk + gammaq)
344
345 params%Fm = 0.0_dp
346 SELECT CASE (potential_parameter%potential_type)
347 CASE (do_potential_coulomb)
348 t = params%Rho*sum((params%Q - rk)**2)
349 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
350 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
351
352 CALL fgamma(params%m_max, t, params%Fm)
353 params%Fm = prefac*params%Fm
354 CASE (do_potential_truncated)
355 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
356 t = params%Rho*sum((params%Q - rk)**2)
357 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
358 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
359
360 cpassert(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
361 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
362 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
363 params%Fm = prefac*params%Fm
364 CASE (do_potential_short)
365 t = params%Rho*sum((params%Q - rk)**2)
366 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
367 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
368
369 CALL fgamma(params%m_max, t, params%Fm)
370
371 omega2 = potential_parameter%omega**2
372 omega_corr2 = omega2/(omega2 + params%Rho)
373 omega_corr = sqrt(omega_corr2)
374 t = t*omega_corr2
375 ALLOCATE (fm(prim_data_f_size))
376
377 CALL fgamma(params%m_max, t, fm)
378 tmp = -omega_corr
379 DO l = 1, params%m_max + 1
380 params%Fm(l) = params%Fm(l) + fm(l)*tmp
381 tmp = tmp*omega_corr2
382 END DO
383 params%Fm = prefac*params%Fm
384 CASE (do_potential_id)
385 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2) &
386 - gammaq*zetk*params%ZetapEtaInv*sum((params%Q - rk)**2))
387 prefac = sqrt((pi*params%ZetapEtaInv)**3)*s1234
388
389 params%Fm(:) = prefac
390 CASE DEFAULT
391 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 params%m_max, params%Fm)
399
400 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 SUBROUTINE eri_3center_derivs(der_abc_1, der_abc_2, &
438 la_min, la_max, npgfa, zeta, rpgfa, ra, &
439 lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
440 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 REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
466 TYPE(params_3c), POINTER :: params
467
468 NULLIFY (params, p_deriv)
469 ALLOCATE (params)
470
471 permute_1 = [4, 5, 6]
472 permute_2 = [7, 8, 9]
473
474 dr_ab = 0.0_dp
475 dr_bc = 0.0_dp
476 dr_ac = 0.0_dp
477
478 op = potential_parameter%potential_type
479
480 IF (op == do_potential_truncated .OR. op == do_potential_short) THEN
481 dr_bc = potential_parameter%cutoff_radius*cutoff_screen_factor
482 dr_ac = potential_parameter%cutoff_radius*cutoff_screen_factor
483 ELSEIF (op == do_potential_coulomb) THEN
484 dr_bc = 1000000.0_dp
485 dr_ac = 1000000.0_dp
486 END IF
487
488 do_ext = .false.
489 IF (PRESENT(der_abc_1_ext) .OR. PRESENT(der_abc_2_ext)) do_ext = .true.
490 der_abc_1_ext_prv = 0.0_dp
491 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 DO ipgf = 1, npgfa
499 zeti = zeta(ipgf)
500 a_start = (ipgf - 1)*ncoset(la_max)
501
502 DO jpgf = 1, npgfb
503
504 ! screening
505 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
506
507 zetj = zetb(jpgf)
508 b_start = (jpgf - 1)*ncoset(lb_max)
509
510 DO kpgf = 1, npgfc
511
512 ! screening
513 IF (rpgfb(jpgf) + rpgfc(kpgf) + dr_bc < dbc) cycle
514 IF (rpgfa(ipgf) + rpgfc(kpgf) + dr_ac < dac) cycle
515
516 zetk = zetc(kpgf)
517 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 potential_parameter=potential_parameter, params_out=params)
522
523 DO li = la_min, la_max
524 a_offset = a_start + ncoset(li - 1)
525 ncoa = nco(li)
526 DO lj = max(li, lb_min), lb_max
527 b_offset = b_start + ncoset(lj - 1)
528 ncob = nco(lj)
529 DO lk = lc_min, lc_max
530 c_offset = c_start + ncoset(lk - 1)
531 ncoc = nco(lk)
532
533 a_mysize(1) = ncoa*ncob*ncoc
534
535 CALL cp_libint_get_3eri_derivs(li, lj, lk, lib, p_deriv, a_mysize)
536
537 IF (do_ext) THEN
538 DO i_deriv = 1, 3
539 DO k = 1, ncoc
540 p1 = (k - 1)*ncob
541 DO j = 1, ncob
542 p2 = (p1 + j - 1)*ncoa
543 DO i = 1, ncoa
544 p3 = p2 + i
545
546 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
547 p_deriv(p3, permute_2(i_deriv))
548 der_abc_1_ext_prv(i_deriv) = max(der_abc_1_ext_prv(i_deriv), &
549 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 p_deriv(p3, permute_1(i_deriv))
553 der_abc_2_ext_prv(i_deriv) = max(der_abc_2_ext_prv(i_deriv), &
554 abs(p_deriv(p3, permute_1(i_deriv))))
555
556 END DO
557 END DO
558 END DO
559 END DO
560 ELSE
561 DO i_deriv = 1, 3
562 DO k = 1, ncoc
563 p1 = (k - 1)*ncob
564 DO j = 1, ncob
565 p2 = (p1 + j - 1)*ncoa
566 DO i = 1, ncoa
567 p3 = p2 + i
568
569 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
570 p_deriv(p3, permute_2(i_deriv))
571
572 der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
573 p_deriv(p3, permute_1(i_deriv))
574 END DO
575 END DO
576 END DO
577 END DO
578 END IF
579
580 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 CALL set_params_3c_deriv(lib, rb, ra, rc, zetj, zeti, zetk, params_in=params)
587
588 DO lj = lb_min, lb_max
589 b_offset = b_start + ncoset(lj - 1)
590 ncob = nco(lj)
591 DO li = max(lj + 1, la_min), la_max
592 a_offset = a_start + ncoset(li - 1)
593 ncoa = nco(li)
594 DO lk = lc_min, lc_max
595 c_offset = c_start + ncoset(lk - 1)
596 ncoc = nco(lk)
597
598 a_mysize(1) = ncoa*ncob*ncoc
599 CALL cp_libint_get_3eri_derivs(lj, li, lk, lib, p_deriv, a_mysize)
600
601 IF (do_ext) THEN
602 DO i_deriv = 1, 3
603 DO k = 1, ncoc
604 p1 = (k - 1)*ncoa
605 DO i = 1, ncoa
606 p2 = (p1 + i - 1)*ncob
607 DO j = 1, ncob
608 p3 = p2 + j
609
610 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
611 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 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 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 abs(p_deriv(p3, permute_2(i_deriv))))
621 END DO
622 END DO
623 END DO
624 END DO
625 ELSE
626 DO i_deriv = 1, 3
627 DO k = 1, ncoc
628 p1 = (k - 1)*ncoa
629 DO i = 1, ncoa
630 p2 = (p1 + i - 1)*ncob
631 DO j = 1, ncob
632 p3 = p2 + j
633
634 der_abc_1(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
635 p_deriv(p3, permute_1(i_deriv))
636
637 der_abc_2(a_offset + i, b_offset + j, c_offset + k, i_deriv) = &
638 p_deriv(p3, permute_2(i_deriv))
639 END DO
640 END DO
641 END DO
642 END DO
643 END IF
644
645 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 IF (PRESENT(der_abc_1_ext)) der_abc_1_ext = der_abc_1_ext_prv
655 IF (PRESENT(der_abc_2_ext)) der_abc_2_ext = der_abc_2_ext_prv
656
657 DEALLOCATE (params)
658
659 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: fm
695 TYPE(params_3c), POINTER :: params
696
697 IF (PRESENT(params_in)) THEN
698 params => params_in
699
700 ELSE
701 params => params_out
702
703 params%m_max = li_max + lj_max + lk_max + 1
704 gammaq = zeti + zetj
705 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/gammaq
706 params%ZetapEtaInv = 1._dp/(zetk + gammaq)
707
708 params%Q = (zeti*ri + zetj*rj)*params%EtaInv
709 params%W = (zetk*rk + gammaq*params%Q)*params%ZetapEtaInv
710 params%Rho = zetk*gammaq/(zetk + gammaq)
711
712 params%Fm = 0.0_dp
713 SELECT CASE (potential_parameter%potential_type)
714 CASE (do_potential_coulomb)
715 t = params%Rho*sum((params%Q - rk)**2)
716 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
717 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
718
719 CALL fgamma(params%m_max, t, params%Fm)
720 params%Fm = prefac*params%Fm
721 CASE (do_potential_truncated)
722 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
723 t = params%Rho*sum((params%Q - rk)**2)
724 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
725 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
726
727 cpassert(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
728 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
729 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
730 params%Fm = prefac*params%Fm
731 CASE (do_potential_short)
732 t = params%Rho*sum((params%Q - rk)**2)
733 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2))
734 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)*s1234
735
736 CALL fgamma(params%m_max, t, params%Fm)
737
738 omega2 = potential_parameter%omega**2
739 omega_corr2 = omega2/(omega2 + params%Rho)
740 omega_corr = sqrt(omega_corr2)
741 t = t*omega_corr2
742 ALLOCATE (fm(prim_data_f_size))
743
744 CALL fgamma(params%m_max, t, fm)
745 tmp = -omega_corr
746 DO l = 1, params%m_max + 1
747 params%Fm(l) = params%Fm(l) + fm(l)*tmp
748 tmp = tmp*omega_corr2
749 END DO
750 params%Fm = prefac*params%Fm
751 CASE (do_potential_id)
752 s1234 = exp(-zeti*zetj*params%EtaInv*sum((rj - ri)**2) &
753 - gammaq*zetk*params%ZetapEtaInv*sum((params%Q - rk)**2))
754 prefac = sqrt((pi*params%ZetapEtaInv)**3)*s1234
755
756 params%Fm(:) = prefac
757 CASE DEFAULT
758 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 params%EtaInv, params%ZetapEtaInv, params%Rho, params%m_max, params%Fm)
766
767 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 SUBROUTINE eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
793 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 REAL(dp), DIMENSION(:), POINTER :: p_work
812
813 NULLIFY (p_work)
814
815 dr_ab = 0.0_dp
816
817 IF (potential_parameter%potential_type == do_potential_truncated .OR. &
818 potential_parameter%potential_type == do_potential_short) THEN
819 dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
820 ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
821 dr_ab = 1000000.0_dp
822 END IF
823
824 !Looping over the pgfs
825 DO ipgf = 1, npgfa
826 zeti = zeta(ipgf)
827 a_start = (ipgf - 1)*ncoset(la_max)
828
829 DO jpgf = 1, npgfb
830 zetj = zetb(jpgf)
831 b_start = (jpgf - 1)*ncoset(lb_max)
832
833 !screening
834 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
835
836 CALL set_params_2c(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
837
838 DO li = la_min, la_max
839 a_offset = a_start + ncoset(li - 1)
840 ncoa = nco(li)
841 DO lj = lb_min, lb_max
842 b_offset = b_start + ncoset(lj - 1)
843 ncob = nco(lj)
844
845 a_mysize(1) = ncoa*ncob
846 CALL cp_libint_get_2eris(li, lj, lib, p_work, a_mysize)
847
848 DO j = 1, ncob
849 p1 = (j - 1)*ncoa
850 DO i = 1, ncoa
851 p2 = p1 + i
852 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 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 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 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 op = potential_parameter%potential_type
900 params%m_max = lj_max + lk_max
901 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
902 params%ZetapEtaInv = 1._dp/(zetk + zetj)
903
904 params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
905 params%Rho = zetk*zetj/(zetk + zetj)
906
907 params%Fm = 0.0_dp
908 SELECT CASE (op)
909 CASE (do_potential_coulomb)
910 t = params%Rho*sum((rj - rk)**2)
911 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
912 CALL fgamma(params%m_max, t, params%Fm)
913 params%Fm = prefac*params%Fm
914 CASE (do_potential_truncated)
915 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
916 t = params%Rho*sum((rj - rk)**2)
917 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
918
919 cpassert(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
920 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
921 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
922 params%Fm = prefac*params%Fm
923 CASE (do_potential_short)
924 t = params%Rho*sum((rj - rk)**2)
925 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
926
927 CALL fgamma(params%m_max, t, params%Fm)
928
929 omega2 = potential_parameter%omega**2
930 omega_corr2 = omega2/(omega2 + params%Rho)
931 omega_corr = sqrt(omega_corr2)
932 t = t*omega_corr2
933 ALLOCATE (fm(prim_data_f_size))
934
935 CALL fgamma(params%m_max, t, fm)
936 tmp = -omega_corr
937 DO l = 1, params%m_max + 1
938 params%Fm(l) = params%Fm(l) + fm(l)*tmp
939 tmp = tmp*omega_corr2
940 END DO
941 params%Fm = prefac*params%Fm
942 CASE (do_potential_id)
943
944 prefac = sqrt((pi*params%ZetapEtaInv)**3)*exp(-zetj*zetk*params%ZetapEtaInv*sum((rk - rj)**2))
945 params%Fm(:) = prefac
946 CASE DEFAULT
947 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 params%m_max, params%Fm)
953
954 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 PURE FUNCTION compare_potential_types(potential1, potential2) RESULT(equals)
963 TYPE(libint_potential_type), INTENT(IN) :: potential1, potential2
964 LOGICAL :: equals
965
966 IF (potential1%potential_type /= potential2%potential_type) THEN
967 equals = .false.
968 ELSE
969 equals = .true.
970 SELECT CASE (potential1%potential_type)
971 CASE (do_potential_short, do_potential_long)
972 IF (potential1%omega /= potential2%omega) equals = .false.
973 CASE (do_potential_truncated)
974 IF (potential1%cutoff_radius /= potential2%cutoff_radius) equals = .false.
975 END SELECT
976 END IF
977
978 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 SUBROUTINE eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
1003 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 REAL(dp), DIMENSION(:, :), POINTER :: p_deriv
1023
1024 NULLIFY (p_deriv)
1025
1026 permute = [4, 5, 6]
1027
1028 dr_ab = 0.0_dp
1029
1030 IF (potential_parameter%potential_type == do_potential_truncated .OR. &
1031 potential_parameter%potential_type == do_potential_short) THEN
1032 dr_ab = potential_parameter%cutoff_radius*cutoff_screen_factor
1033 ELSEIF (potential_parameter%potential_type == do_potential_coulomb) THEN
1034 dr_ab = 1000000.0_dp
1035 END IF
1036
1037 !Looping over the pgfs
1038 DO ipgf = 1, npgfa
1039 zeti = zeta(ipgf)
1040 a_start = (ipgf - 1)*ncoset(la_max)
1041
1042 DO jpgf = 1, npgfb
1043 zetj = zetb(jpgf)
1044 b_start = (jpgf - 1)*ncoset(lb_max)
1045
1046 !screening
1047 IF (rpgfa(ipgf) + rpgfb(jpgf) + dr_ab < dab) cycle
1048
1049 CALL set_params_2c_deriv(lib, ra, rb, zeti, zetj, la_max, lb_max, potential_parameter)
1050
1051 DO li = la_min, la_max
1052 a_offset = a_start + ncoset(li - 1)
1053 ncoa = nco(li)
1054 DO lj = lb_min, lb_max
1055 b_offset = b_start + ncoset(lj - 1)
1056 ncob = nco(lj)
1057
1058 a_mysize(1) = ncoa*ncob
1059 CALL cp_libint_get_2eri_derivs(li, lj, lib, p_deriv, a_mysize)
1060
1061 DO i_deriv = 1, 3
1062 DO j = 1, ncob
1063 p1 = (j - 1)*ncoa
1064 DO i = 1, ncoa
1065 p2 = p1 + i
1066 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 DEALLOCATE (p_deriv)
1072 END DO
1073 END DO
1074
1075 END DO
1076 END DO
1077
1078 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 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 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 op = potential_parameter%potential_type
1116 params%m_max = lj_max + lk_max + 1
1117 params%ZetaInv = 1._dp/zetk; params%EtaInv = 1._dp/zetj
1118 params%ZetapEtaInv = 1._dp/(zetk + zetj)
1119
1120 params%W = (zetk*rk + zetj*rj)*params%ZetapEtaInv
1121 params%Rho = zetk*zetj/(zetk + zetj)
1122
1123 params%Fm = 0.0_dp
1124 SELECT CASE (op)
1125 CASE (do_potential_coulomb)
1126 t = params%Rho*sum((rj - rk)**2)
1127 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1128 CALL fgamma(params%m_max, t, params%Fm)
1129 params%Fm = prefac*params%Fm
1130 CASE (do_potential_truncated)
1131 r = potential_parameter%cutoff_radius*sqrt(params%Rho)
1132 t = params%Rho*sum((rj - rk)**2)
1133 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1134
1135 cpassert(get_lmax_init() .GE. params%m_max) !check if truncated coulomb init correctly
1136 CALL t_c_g0_n(params%Fm, use_gamma, r, t, params%m_max)
1137 IF (use_gamma) CALL fgamma(params%m_max, t, params%Fm)
1138 params%Fm = prefac*params%Fm
1139 CASE (do_potential_short)
1140 t = params%Rho*sum((rj - rk)**2)
1141 prefac = 2._dp*pi/params%Rho*sqrt((pi*params%ZetapEtaInv)**3)
1142
1143 CALL fgamma(params%m_max, t, params%Fm)
1144
1145 omega2 = potential_parameter%omega**2
1146 omega_corr2 = omega2/(omega2 + params%Rho)
1147 omega_corr = sqrt(omega_corr2)
1148 t = t*omega_corr2
1149 ALLOCATE (fm(prim_data_f_size))
1150
1151 CALL fgamma(params%m_max, t, fm)
1152 tmp = -omega_corr
1153 DO l = 1, params%m_max + 1
1154 params%Fm(l) = params%Fm(l) + fm(l)*tmp
1155 tmp = tmp*omega_corr2
1156 END DO
1157 params%Fm = prefac*params%Fm
1158 CASE (do_potential_id)
1159
1160 prefac = sqrt((pi*params%ZetapEtaInv)**3)*exp(-zetj*zetk*params%ZetapEtaInv*sum((rk - rj)**2))
1161 params%Fm(:) = prefac
1162 CASE DEFAULT
1163 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 params%m_max, params%Fm)
1170
1171 END SUBROUTINE set_params_2c_deriv
1172
1173END MODULE libint_2c_3c
1174
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::fgamma_0
subroutine, public fgamma_0(nmax, t, f)
Calculation of the incomplete Gamma function F(t) for multicenter integrals over Gaussian functions....
Definition gamma.F:154
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:801
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:801
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:801
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:801
input_constants::do_potential_long
integer, parameter, public do_potential_long
Definition input_constants.F:801
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::eri_2center
subroutine, public eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian gaussian orbita...
Definition libint_2c_3c.F:795
libint_2c_3c::compare_potential_types
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare libint_potential_types.
Definition libint_2c_3c.F:963
libint_2c_3c::eri_3center
subroutine, public eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, int_abc_ext)
Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian gaussian orbit...
Definition libint_2c_3c.F:114
libint_2c_3c::cutoff_screen_factor
real(kind=dp), parameter, public cutoff_screen_factor
Definition libint_2c_3c.F:48
libint_2c_3c::eri_2center_derivs
subroutine, public eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given set of cartes...
Definition libint_2c_3c.F:1005
libint_2c_3c::eri_3center_derivs
subroutine, public eri_3center_derivs(der_abc_1, der_abc_2, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, der_abc_1_ext, der_abc_2_ext)
Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given set of carte...
Definition libint_2c_3c.F:443
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::cp_libint_set_params_eri_deriv
subroutine, public cp_libint_set_params_eri_deriv(libint, a, b, c, d, p, q, w, zeta_a, zeta_b, zeta_c, zeta_d, zetainv, etainv, zetapetainv, rho, m_max, f)
Definition libint_wrapper.F:334
libint_wrapper::cp_libint_get_2eri_derivs
subroutine, public cp_libint_get_2eri_derivs(n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:956
libint_wrapper::cp_libint_set_params_eri
subroutine, public cp_libint_set_params_eri(libint, a, b, c, d, zetainv, etainv, zetapetainv, rho, p, q, w, m_max, f)
Definition libint_wrapper.F:583
libint_wrapper::cp_libint_get_3eris
subroutine, public cp_libint_get_3eris(n_c, n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:879
libint_wrapper::cp_libint_get_2eris
subroutine, public cp_libint_get_2eris(n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:998
libint_wrapper::prim_data_f_size
integer, parameter, public prim_data_f_size
Definition libint_wrapper.F:63
libint_wrapper::cp_libint_get_3eri_derivs
subroutine, public cp_libint_get_3eri_derivs(n_c, n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:913
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:57
t_c_g0::t_c_g0_n
subroutine, public t_c_g0_n(res, use_gamma, r, t, nderiv)
...
Definition t_c_g0.F:84
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition t_c_g0.F:1464
libint_2c_3c::libint_potential_type
Definition libint_2c_3c.F:66
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:69