(git:374b731)
Loading...
Searching...
No Matches
hfx_libint_interface.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 Interface to the Libint-Library
10!> \par History
11!> 11.2006 created [Manuel Guidon]
12!> 11.2019 Fixed potential_id initial values (A. Bussy)
13!> \author Manuel Guidon
14! **************************************************************************************************
15MODULE hfx_libint_interface
16
17 USE cell_types, ONLY: cell_type,&
18 real_to_scaled
19 USE gamma, ONLY: fgamma => fgamma_0
20 USE hfx_contraction_methods, ONLY: contract
21 USE hfx_types, ONLY: hfx_pgf_product_list,&
22 hfx_potential_type
23 USE input_constants, ONLY: &
24 do_potential_coulomb, do_potential_gaussian, do_potential_id, do_potential_long, &
25 do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, do_potential_short, &
26 do_potential_truncated
27 USE kinds, ONLY: dp,&
28 int_8
29 USE libint_wrapper, ONLY: cp_libint_get_derivs,&
30 cp_libint_get_eris,&
31 cp_libint_set_params_eri,&
32 cp_libint_set_params_eri_deriv,&
33 cp_libint_set_params_eri_screen,&
34 cp_libint_t,&
35 get_ssss_f_val,&
36 prim_data_f_size
37 USE mathconstants, ONLY: pi
38 USE orbital_pointers, ONLY: nco
39 USE t_c_g0, ONLY: t_c_g0_n
40#include "./base/base_uses.f90"
41
42 IMPLICIT NONE
43 PRIVATE
44 PUBLIC :: evaluate_eri, &
45 evaluate_deriv_eri, &
46 evaluate_eri_screen
47
48 INTEGER, DIMENSION(12), PARAMETER :: full_perm1 = (/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12/)
49 INTEGER, DIMENSION(12), PARAMETER :: full_perm2 = (/4, 5, 6, 1, 2, 3, 7, 8, 9, 10, 11, 12/)
50 INTEGER, DIMENSION(12), PARAMETER :: full_perm3 = (/1, 2, 3, 4, 5, 6, 10, 11, 12, 7, 8, 9/)
51 INTEGER, DIMENSION(12), PARAMETER :: full_perm4 = (/4, 5, 6, 1, 2, 3, 10, 11, 12, 7, 8, 9/)
52 INTEGER, DIMENSION(12), PARAMETER :: full_perm5 = (/7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6/)
53 INTEGER, DIMENSION(12), PARAMETER :: full_perm6 = (/7, 8, 9, 10, 11, 12, 4, 5, 6, 1, 2, 3/)
54 INTEGER, DIMENSION(12), PARAMETER :: full_perm7 = (/10, 11, 12, 7, 8, 9, 1, 2, 3, 4, 5, 6/)
55 INTEGER, DIMENSION(12), PARAMETER :: full_perm8 = (/10, 11, 12, 7, 8, 9, 4, 5, 6, 1, 2, 3/)
56
57!***
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief Fill data structure used in libint
63!> \param A ...
64!> \param B ...
65!> \param C ...
66!> \param D ...
67!> \param Zeta_A ...
68!> \param Zeta_B ...
69!> \param Zeta_C ...
70!> \param Zeta_D ...
71!> \param m_max ...
72!> \param potential_parameter ...
73!> \param libint ...
74!> \param R11 ...
75!> \param R22 ...
76!> \par History
77!> 03.2007 created [Manuel Guidon]
78!> \author Manuel Guidon
79! **************************************************************************************************
80 SUBROUTINE build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, m_max, &
81 potential_parameter, libint, R11, R22)
82
83 REAL(KIND=dp) :: a(3), b(3), c(3), d(3)
84 REAL(KIND=dp), INTENT(IN) :: zeta_a, zeta_b, zeta_c, zeta_d
85 INTEGER, INTENT(IN) :: m_max
86 TYPE(hfx_potential_type) :: potential_parameter
87 TYPE(cp_libint_t) :: libint
88 REAL(dp) :: R11, R22
89
90 INTEGER :: i
91 LOGICAL :: use_gamma
92 REAL(KIND=dp) :: ab(3), ab2, cd(3), cd2, den, eta, etainv, factor, g(3), num, omega2, &
93 omega_corr, omega_corr2, p(3), pq(3), pq2, q(3), r, r1, r2, rho, rhoinv, s1234, ssss, t, &
94 tmp, w(3), zeta, zetainv, zetapetainv
95 REAL(KIND=dp), DIMENSION(prim_data_f_size) :: f, fm
96
97 zeta = zeta_a + zeta_b
98 zetainv = 1.0_dp/zeta
99 eta = zeta_c + zeta_d
100 etainv = 1.0_dp/eta
101 zetapetainv = zeta + eta
102 zetapetainv = 1.0_dp/zetapetainv
103 rho = zeta*eta*zetapetainv
104 rhoinv = 1.0_dp/rho
105
106 DO i = 1, 3
107 p(i) = (zeta_a*a(i) + zeta_b*b(i))*zetainv
108 q(i) = (zeta_c*c(i) + zeta_d*d(i))*etainv
109 ab(i) = a(i) - b(i)
110 cd(i) = c(i) - d(i)
111 pq(i) = p(i) - q(i)
112 w(i) = (zeta*p(i) + eta*q(i))*zetapetainv
113 END DO
114
115 ab2 = dot_product(ab, ab)
116 cd2 = dot_product(cd, cd)
117 pq2 = dot_product(pq, pq)
118
119 s1234 = exp((-zeta_a*zeta_b*zetainv*ab2) + (-zeta_c*zeta_d*etainv*cd2))
120 t = rho*pq2
121
122 SELECT CASE (potential_parameter%potential_type)
123 CASE (do_potential_truncated)
124 r = potential_parameter%cutoff_radius*sqrt(rho)
125 r1 = r11
126 r2 = r22
127 IF (pq2 > (r1 + r2 + potential_parameter%cutoff_radius)**2) THEN
128 RETURN
129 END IF
130 CALL t_c_g0_n(f(1), use_gamma, r, t, m_max)
131 IF (use_gamma) CALL fgamma(m_max, t, f(1))
132 factor = 2.0_dp*pi*rhoinv
133 CASE (do_potential_coulomb)
134 CALL fgamma(m_max, t, f(1))
135 factor = 2.0_dp*pi*rhoinv
136 CASE (do_potential_short)
137 r = potential_parameter%cutoff_radius*sqrt(rho)
138 r1 = r11
139 r2 = r22
140 IF (pq2 > (r1 + r2 + potential_parameter%cutoff_radius)**2) THEN
141 RETURN
142 END IF
143 CALL fgamma(m_max, t, f(1))
144 omega2 = potential_parameter%omega**2
145 omega_corr2 = omega2/(omega2 + rho)
146 omega_corr = sqrt(omega_corr2)
147 t = t*omega_corr2
148 CALL fgamma(m_max, t, fm)
149 tmp = -omega_corr
150 DO i = 1, m_max + 1
151 f(i) = f(i) + fm(i)*tmp
152 tmp = tmp*omega_corr2
153 END DO
154 factor = 2.0_dp*pi*rhoinv
155 CASE (do_potential_long)
156 omega2 = potential_parameter%omega**2
157 omega_corr2 = omega2/(omega2 + rho)
158 omega_corr = sqrt(omega_corr2)
159 t = t*omega_corr2
160 CALL fgamma(m_max, t, f(1))
161 tmp = omega_corr
162 DO i = 1, m_max + 1
163 f(i) = f(i)*tmp
164 tmp = tmp*omega_corr2
165 END DO
166 factor = 2.0_dp*pi*rhoinv
167 CASE (do_potential_mix_cl)
168 CALL fgamma(m_max, t, f(1))
169 omega2 = potential_parameter%omega**2
170 omega_corr2 = omega2/(omega2 + rho)
171 omega_corr = sqrt(omega_corr2)
172 t = t*omega_corr2
173 CALL fgamma(m_max, t, fm)
174 tmp = omega_corr
175 DO i = 1, m_max + 1
176 f(i) = f(i)*potential_parameter%scale_coulomb + fm(i)*tmp*potential_parameter%scale_longrange
177 tmp = tmp*omega_corr2
178 END DO
179 factor = 2.0_dp*pi*rhoinv
180 CASE (do_potential_mix_cl_trunc)
181
182 ! truncated
183 r = potential_parameter%cutoff_radius*sqrt(rho)
184 r1 = r11
185 r2 = r22
186 IF (pq2 > (r1 + r2 + potential_parameter%cutoff_radius)**2) THEN
187 RETURN
188 END IF
189 CALL t_c_g0_n(f(1), use_gamma, r, t, m_max)
190 IF (use_gamma) CALL fgamma(m_max, t, f(1))
191
192 ! Coulomb
193 CALL fgamma(m_max, t, fm)
194
195 DO i = 1, m_max + 1
196 f(i) = f(i)*(potential_parameter%scale_coulomb + potential_parameter%scale_longrange) - &
197 fm(i)*potential_parameter%scale_longrange
198 END DO
199
200 ! longrange
201 omega2 = potential_parameter%omega**2
202 omega_corr2 = omega2/(omega2 + rho)
203 omega_corr = sqrt(omega_corr2)
204 t = t*omega_corr2
205 CALL fgamma(m_max, t, fm)
206 tmp = omega_corr
207 DO i = 1, m_max + 1
208 f(i) = f(i) + fm(i)*tmp*potential_parameter%scale_longrange
209 tmp = tmp*omega_corr2
210 END DO
211 factor = 2.0_dp*pi*rhoinv
212
213 CASE (do_potential_gaussian)
214 omega2 = potential_parameter%omega**2
215 t = -omega2*t/(rho + omega2)
216 tmp = 1.0_dp
217 DO i = 1, m_max + 1
218 f(i) = exp(t)*tmp
219 tmp = tmp*omega2/(rho + omega2)
220 END DO
221 factor = (pi/(rho + omega2))**(1.5_dp)
222 CASE (do_potential_mix_lg)
223 omega2 = potential_parameter%omega**2
224 omega_corr2 = omega2/(omega2 + rho)
225 omega_corr = sqrt(omega_corr2)
226 t = t*omega_corr2
227 CALL fgamma(m_max, t, fm)
228 tmp = omega_corr*2.0_dp*pi*rhoinv*potential_parameter%scale_longrange
229 DO i = 1, m_max + 1
230 fm(i) = fm(i)*tmp
231 tmp = tmp*omega_corr2
232 END DO
233 t = rho*pq2
234 t = -omega2*t/(rho + omega2)
235 tmp = (pi/(rho + omega2))**(1.5_dp)*potential_parameter%scale_gaussian
236 DO i = 1, m_max + 1
237 f(i) = exp(t)*tmp + fm(i)
238 tmp = tmp*omega2/(rho + omega2)
239 END DO
240 factor = 1.0_dp
241 CASE (do_potential_id)
242 ssss = -zeta_a*zeta_b*zetainv*ab2
243
244 num = (zeta_a + zeta_b)*zeta_c
245 den = zeta_a + zeta_b + zeta_c
246 ssss = ssss - num/den*sum((p - c)**2)
247
248 g(:) = (zeta_a*a(:) + zeta_b*b(:) + zeta_c*c(:))/den
249 num = den*zeta_d
250 den = den + zeta_d
251 ssss = ssss - num/den*sum((g - d)**2)
252
253 f(:) = exp(ssss)
254 factor = 1.0_dp
255 IF (s1234 > epsilon(0.0_dp)) factor = 1.0_dp/s1234
256 END SELECT
257
258 tmp = (pi*zetapetainv)**3
259 factor = factor*s1234*sqrt(tmp)
260
261 DO i = 1, m_max + 1
262 f(i) = f(i)*factor
263 END DO
264
265 CALL cp_libint_set_params_eri_screen(libint, a, b, c, d, p, q, w, zetainv, etainv, zetapetainv, rho, m_max, f)
266
267 END SUBROUTINE build_quartet_data_screen
268
269! **************************************************************************************************
270!> \brief Fill data structure used in libderiv
271!> \param libint ...
272!> \param A ...
273!> \param B ...
274!> \param C ...
275!> \param D ...
276!> \param Zeta_A ...
277!> \param Zeta_B ...
278!> \param Zeta_C ...
279!> \param Zeta_D ...
280!> \param ZetaInv ...
281!> \param EtaInv ...
282!> \param ZetapEtaInv ...
283!> \param Rho ...
284!> \param RhoInv ...
285!> \param P ...
286!> \param Q ...
287!> \param W ...
288!> \param m_max ...
289!> \param F ...
290!> \par History
291!> 03.2007 created [Manuel Guidon]
292!> \author Manuel Guidon
293! **************************************************************************************************
294 SUBROUTINE build_deriv_data(libint, A, B, C, D, &
295 Zeta_A, Zeta_B, Zeta_C, Zeta_D, &
296 ZetaInv, EtaInv, ZetapEtaInv, Rho, RhoInv, &
297 P, Q, W, m_max, F)
298
299 TYPE(cp_libint_t) :: libint
300 REAL(KIND=dp), INTENT(IN) :: a(3), b(3), c(3), d(3)
301 REAL(dp), INTENT(IN) :: Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, &
302 EtaInv, ZetapEtaInv, Rho, RhoInv, &
303 P(3), Q(3), W(3)
304 INTEGER, INTENT(IN) :: m_max
305 REAL(KIND=dp), DIMENSION(:) :: f
306
307 mark_used(d) ! todo: fix
308 mark_used(zeta_c)
309 mark_used(rhoinv)
310
311 CALL cp_libint_set_params_eri_deriv(libint, a, b, c, d, p, q, w, zeta_a, zeta_b, zeta_c, zeta_d, &
312 zetainv, etainv, zetapetainv, rho, m_max, f)
313
314 END SUBROUTINE build_deriv_data
315
316! **************************************************************************************************
317!> \brief Evaluates derivatives of electron repulsion integrals for a primitive quartet
318!> \param deriv ...
319!> \param nproducts ...
320!> \param pgf_product_list ...
321!> \param n_a ...
322!> \param n_b ...
323!> \param n_c ...
324!> \param n_d ...
325!> \param ncoa ...
326!> \param ncob ...
327!> \param ncoc ...
328!> \param ncod ...
329!> \param nsgfa ...
330!> \param nsgfb ...
331!> \param nsgfc ...
332!> \param nsgfd ...
333!> \param primitives ...
334!> \param max_contraction ...
335!> \param tmp_max_all ...
336!> \param eps_schwarz ...
337!> \param neris ...
338!> \param Zeta_A ...
339!> \param Zeta_B ...
340!> \param Zeta_C ...
341!> \param Zeta_D ...
342!> \param ZetaInv ...
343!> \param EtaInv ...
344!> \param s_offset_a ...
345!> \param s_offset_b ...
346!> \param s_offset_c ...
347!> \param s_offset_d ...
348!> \param nl_a ...
349!> \param nl_b ...
350!> \param nl_c ...
351!> \param nl_d ...
352!> \param nsoa ...
353!> \param nsob ...
354!> \param nsoc ...
355!> \param nsod ...
356!> \param sphi_a ...
357!> \param sphi_b ...
358!> \param sphi_c ...
359!> \param sphi_d ...
360!> \param work ...
361!> \param work2 ...
362!> \param work_forces ...
363!> \param buffer1 ...
364!> \param buffer2 ...
365!> \param primitives_tmp ...
366!> \param use_virial ...
367!> \param work_virial ...
368!> \param work2_virial ...
369!> \param primitives_tmp_virial ...
370!> \param primitives_virial ...
371!> \param cell ...
372!> \param tmp_max_all_virial ...
373!> \par History
374!> 03.2007 created [Manuel Guidon]
375!> 08.2007 refactured permutation part [Manuel Guidon]
376!> \author Manuel Guidon
377! **************************************************************************************************
378 SUBROUTINE evaluate_deriv_eri(deriv, nproducts, pgf_product_list, &
379 n_a, n_b, n_c, n_d, &
380 ncoa, ncob, ncoc, ncod, &
381 nsgfa, nsgfb, nsgfc, nsgfd, &
382 primitives, max_contraction, tmp_max_all, &
383 eps_schwarz, neris, &
384 Zeta_A, Zeta_B, Zeta_C, Zeta_D, ZetaInv, EtaInv, &
385 s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
386 nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, &
387 sphi_a, sphi_b, sphi_c, sphi_d, &
388 work, work2, work_forces, buffer1, buffer2, primitives_tmp, &
389 use_virial, work_virial, work2_virial, primitives_tmp_virial, &
390 primitives_virial, cell, tmp_max_all_virial)
391
392 TYPE(cp_libint_t) :: deriv
393 INTEGER, INTENT(IN) :: nproducts
394 TYPE(hfx_pgf_product_list), DIMENSION(nproducts) :: pgf_product_list
395 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, &
396 ncod, nsgfa, nsgfb, nsgfc, nsgfd
397 REAL(dp), &
398 DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12) :: primitives
399 REAL(dp) :: max_contraction, tmp_max_all, eps_schwarz
400 INTEGER(int_8) :: neris
401 REAL(dp), INTENT(IN) :: zeta_a, zeta_b, zeta_c, zeta_d, zetainv, &
402 etainv
403 INTEGER :: s_offset_a, s_offset_b, s_offset_c, &
404 s_offset_d, nl_a, nl_b, nl_c, nl_d, &
405 nsoa, nsob, nsoc, nsod
406 REAL(dp), DIMENSION(ncoa, nsoa*nl_a) :: sphi_a
407 REAL(dp), DIMENSION(ncob, nsob*nl_b) :: sphi_b
408 REAL(dp), DIMENSION(ncoc, nsoc*nl_c) :: sphi_c
409 REAL(dp), DIMENSION(ncod, nsod*nl_d) :: sphi_d
410 REAL(dp) :: work(ncoa*ncob*ncoc*ncod, 12), work2(ncoa, ncob, ncoc, ncod, 12), &
411 work_forces(ncoa*ncob*ncoc*ncod, 12)
412 REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2
413 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*& nl_c, nsod*nl_d) :: primitives_tmp
414 LOGICAL, INTENT(IN) :: use_virial
415 REAL(dp) :: work_virial(ncoa*ncob*ncoc*ncod, 12, 3), &
416 work2_virial(ncoa, ncob, ncoc, ncod, 12, 3)
417 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*& nl_c, nsod*nl_d) :: primitives_tmp_virial
418 REAL(dp), &
419 DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd, 12, 3) :: primitives_virial
420 TYPE(cell_type), POINTER :: cell
421 REAL(dp) :: tmp_max_all_virial
422
423 INTEGER :: a_mysize(1), i, j, k, l, m, m_max, &
424 mysize, n, p1, p2, p3, p4, perm_case
425 REAL(dp) :: a(3), ab(3), b(3), c(3), cd(3), d(3), &
426 p(3), q(3), rho, rhoinv, scoord(12), &
427 tmp_max, tmp_max_virial, w(3), &
428 zetapetainv
429
430 m_max = n_a + n_b + n_c + n_d
431 m_max = m_max + 1
432 mysize = ncoa*ncob*ncoc*ncod
433 a_mysize = mysize
434
435 work = 0.0_dp
436 work2 = 0.0_dp
437
438 IF (use_virial) THEN
439 work_virial = 0.0_dp
440 work2_virial = 0.0_dp
441 END IF
442
443 perm_case = 1
444 IF (n_a < n_b) THEN
445 perm_case = perm_case + 1
446 END IF
447 IF (n_c < n_d) THEN
448 perm_case = perm_case + 2
449 END IF
450 IF (n_a + n_b > n_c + n_d) THEN
451 perm_case = perm_case + 4
452 END IF
453
454 SELECT CASE (perm_case)
455 CASE (1)
456 DO i = 1, nproducts
457 a = pgf_product_list(i)%ra
458 b = pgf_product_list(i)%rb
459 c = pgf_product_list(i)%rc
460 d = pgf_product_list(i)%rd
461 rho = pgf_product_list(i)%Rho
462 rhoinv = pgf_product_list(i)%RhoInv
463 p = pgf_product_list(i)%P
464 q = pgf_product_list(i)%Q
465 w = pgf_product_list(i)%W
466 ab = pgf_product_list(i)%AB
467 cd = pgf_product_list(i)%CD
468 zetapetainv = pgf_product_list(i)%ZetapEtaInv
469
470 CALL build_deriv_data(deriv, a, b, c, d, &
471 zeta_a, zeta_b, zeta_c, zeta_d, &
472 zetainv, etainv, zetapetainv, rho, rhoinv, &
473 p, q, w, m_max, pgf_product_list(i)%Fm)
474 CALL cp_libint_get_derivs(n_d, n_c, n_b, n_a, deriv, work_forces, a_mysize)
475 DO k = 4, 6
476 DO j = 1, mysize
477 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
478 work_forces(j, k + 3) + &
479 work_forces(j, k + 6))
480 END DO
481 END DO
482 DO k = 1, 12
483 DO j = 1, mysize
484 work(j, k) = work(j, k) + work_forces(j, k)
485 END DO
486 END DO
487 neris = neris + 12*mysize
488 IF (use_virial) THEN
489 CALL real_to_scaled(scoord(1:3), a, cell)
490 CALL real_to_scaled(scoord(4:6), b, cell)
491 CALL real_to_scaled(scoord(7:9), c, cell)
492 CALL real_to_scaled(scoord(10:12), d, cell)
493 DO k = 1, 12
494 DO j = 1, mysize
495 DO m = 1, 3
496 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
497 END DO
498 END DO
499 END DO
500 END IF
501 END DO
502
503 DO n = 1, 12
504 tmp_max = 0.0_dp
505 DO i = 1, mysize
506 tmp_max = max(tmp_max, abs(work(i, n)))
507 END DO
508 tmp_max = tmp_max*max_contraction
509 tmp_max_all = max(tmp_max_all, tmp_max)
510
511 DO i = 1, ncoa
512 p1 = (i - 1)*ncob
513 DO j = 1, ncob
514 p2 = (p1 + j - 1)*ncoc
515 DO k = 1, ncoc
516 p3 = (p2 + k - 1)*ncod
517 DO l = 1, ncod
518 p4 = p3 + l
519 work2(i, j, k, l, full_perm1(n)) = work(p4, n)
520 END DO
521 END DO
522 END DO
523 END DO
524 END DO
525 IF (use_virial) THEN
526 DO n = 1, 12
527 tmp_max_virial = 0.0_dp
528 DO i = 1, mysize
529 tmp_max_virial = max(tmp_max_virial, &
530 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
531 END DO
532 tmp_max_virial = tmp_max_virial*max_contraction
533 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
534
535 DO i = 1, ncoa
536 p1 = (i - 1)*ncob
537 DO j = 1, ncob
538 p2 = (p1 + j - 1)*ncoc
539 DO k = 1, ncoc
540 p3 = (p2 + k - 1)*ncod
541 DO l = 1, ncod
542 p4 = p3 + l
543 work2_virial(i, j, k, l, full_perm1(n), 1:3) = work_virial(p4, n, 1:3)
544 END DO
545 END DO
546 END DO
547 END DO
548 END DO
549 END IF
550 CASE (2)
551 DO i = 1, nproducts
552 a = pgf_product_list(i)%ra
553 b = pgf_product_list(i)%rb
554 c = pgf_product_list(i)%rc
555 d = pgf_product_list(i)%rd
556 rho = pgf_product_list(i)%Rho
557 rhoinv = pgf_product_list(i)%RhoInv
558 p = pgf_product_list(i)%P
559 q = pgf_product_list(i)%Q
560 w = pgf_product_list(i)%W
561 ab = pgf_product_list(i)%AB
562 cd = pgf_product_list(i)%CD
563 zetapetainv = pgf_product_list(i)%ZetapEtaInv
564
565 CALL build_deriv_data(deriv, b, a, c, d, &
566 zeta_b, zeta_a, zeta_c, zeta_d, &
567 zetainv, etainv, zetapetainv, rho, rhoinv, &
568 p, q, w, m_max, pgf_product_list(i)%Fm)
569 CALL cp_libint_get_derivs(n_d, n_c, n_a, n_b, deriv, work_forces, a_mysize)
570 DO k = 4, 6
571 DO j = 1, mysize
572 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
573 work_forces(j, k + 3) + &
574 work_forces(j, k + 6))
575 END DO
576 END DO
577 DO k = 1, 12
578 DO j = 1, mysize
579 work(j, k) = work(j, k) + work_forces(j, k)
580 END DO
581 END DO
582 neris = neris + 12*mysize
583 IF (use_virial) THEN
584 CALL real_to_scaled(scoord(1:3), b, cell)
585 CALL real_to_scaled(scoord(4:6), a, cell)
586 CALL real_to_scaled(scoord(7:9), c, cell)
587 CALL real_to_scaled(scoord(10:12), d, cell)
588 DO k = 1, 12
589 DO j = 1, mysize
590 DO m = 1, 3
591 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
592 END DO
593 END DO
594 END DO
595 END IF
596
597 END DO
598 DO n = 1, 12
599 tmp_max = 0.0_dp
600 DO i = 1, mysize
601 tmp_max = max(tmp_max, abs(work(i, n)))
602 END DO
603 tmp_max = tmp_max*max_contraction
604 tmp_max_all = max(tmp_max_all, tmp_max)
605
606 DO j = 1, ncob
607 p1 = (j - 1)*ncoa
608 DO i = 1, ncoa
609 p2 = (p1 + i - 1)*ncoc
610 DO k = 1, ncoc
611 p3 = (p2 + k - 1)*ncod
612 DO l = 1, ncod
613 p4 = p3 + l
614 work2(i, j, k, l, full_perm2(n)) = work(p4, n)
615 END DO
616 END DO
617 END DO
618 END DO
619 END DO
620 IF (use_virial) THEN
621 DO n = 1, 12
622 tmp_max_virial = 0.0_dp
623 DO i = 1, mysize
624 tmp_max_virial = max(tmp_max_virial, &
625 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
626 END DO
627 tmp_max_virial = tmp_max_virial*max_contraction
628 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
629
630 DO j = 1, ncob
631 p1 = (j - 1)*ncoa
632 DO i = 1, ncoa
633 p2 = (p1 + i - 1)*ncoc
634 DO k = 1, ncoc
635 p3 = (p2 + k - 1)*ncod
636 DO l = 1, ncod
637 p4 = p3 + l
638 work2_virial(i, j, k, l, full_perm2(n), 1:3) = work_virial(p4, n, 1:3)
639 END DO
640 END DO
641 END DO
642 END DO
643 END DO
644 END IF
645 CASE (3)
646 DO i = 1, nproducts
647 a = pgf_product_list(i)%ra
648 b = pgf_product_list(i)%rb
649 c = pgf_product_list(i)%rc
650 d = pgf_product_list(i)%rd
651 rho = pgf_product_list(i)%Rho
652 rhoinv = pgf_product_list(i)%RhoInv
653 p = pgf_product_list(i)%P
654 q = pgf_product_list(i)%Q
655 w = pgf_product_list(i)%W
656 ab = pgf_product_list(i)%AB
657 cd = pgf_product_list(i)%CD
658 zetapetainv = pgf_product_list(i)%ZetapEtaInv
659
660 CALL build_deriv_data(deriv, a, b, d, c, &
661 zeta_a, zeta_b, zeta_d, zeta_c, &
662 zetainv, etainv, zetapetainv, rho, rhoinv, &
663 p, q, w, m_max, pgf_product_list(i)%Fm)
664 CALL cp_libint_get_derivs(n_c, n_d, n_b, n_a, deriv, work_forces, a_mysize)
665 DO k = 4, 6
666 DO j = 1, mysize
667 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
668 work_forces(j, k + 3) + &
669 work_forces(j, k + 6))
670 END DO
671 END DO
672 DO k = 1, 12
673 DO j = 1, mysize
674 work(j, k) = work(j, k) + work_forces(j, k)
675 END DO
676 END DO
677 neris = neris + 12*mysize
678 IF (use_virial) THEN
679 CALL real_to_scaled(scoord(1:3), a, cell)
680 CALL real_to_scaled(scoord(4:6), b, cell)
681 CALL real_to_scaled(scoord(7:9), d, cell)
682 CALL real_to_scaled(scoord(10:12), c, cell)
683 DO k = 1, 12
684 DO j = 1, mysize
685 DO m = 1, 3
686 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
687 END DO
688 END DO
689 END DO
690 END IF
691
692 END DO
693 DO n = 1, 12
694 tmp_max = 0.0_dp
695 DO i = 1, mysize
696 tmp_max = max(tmp_max, abs(work(i, n)))
697 END DO
698 tmp_max = tmp_max*max_contraction
699 tmp_max_all = max(tmp_max_all, tmp_max)
700
701 DO i = 1, ncoa
702 p1 = (i - 1)*ncob
703 DO j = 1, ncob
704 p2 = (p1 + j - 1)*ncod
705 DO l = 1, ncod
706 p3 = (p2 + l - 1)*ncoc
707 DO k = 1, ncoc
708 p4 = p3 + k
709 work2(i, j, k, l, full_perm3(n)) = work(p4, n)
710 END DO
711 END DO
712 END DO
713 END DO
714 END DO
715 IF (use_virial) THEN
716 DO n = 1, 12
717 tmp_max_virial = 0.0_dp
718 DO i = 1, mysize
719 tmp_max_virial = max(tmp_max_virial, &
720 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
721 END DO
722 tmp_max_virial = tmp_max_virial*max_contraction
723 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
724
725 DO i = 1, ncoa
726 p1 = (i - 1)*ncob
727 DO j = 1, ncob
728 p2 = (p1 + j - 1)*ncod
729 DO l = 1, ncod
730 p3 = (p2 + l - 1)*ncoc
731 DO k = 1, ncoc
732 p4 = p3 + k
733 work2_virial(i, j, k, l, full_perm3(n), 1:3) = work_virial(p4, n, 1:3)
734 END DO
735 END DO
736 END DO
737 END DO
738 END DO
739 END IF
740 CASE (4)
741 DO i = 1, nproducts
742 a = pgf_product_list(i)%ra
743 b = pgf_product_list(i)%rb
744 c = pgf_product_list(i)%rc
745 d = pgf_product_list(i)%rd
746 rho = pgf_product_list(i)%Rho
747 rhoinv = pgf_product_list(i)%RhoInv
748 p = pgf_product_list(i)%P
749 q = pgf_product_list(i)%Q
750 w = pgf_product_list(i)%W
751 ab = pgf_product_list(i)%AB
752 cd = pgf_product_list(i)%CD
753 zetapetainv = pgf_product_list(i)%ZetapEtaInv
754 CALL build_deriv_data(deriv, b, a, d, c, &
755 zeta_b, zeta_a, zeta_d, zeta_c, &
756 zetainv, etainv, zetapetainv, rho, rhoinv, &
757 p, q, w, m_max, pgf_product_list(i)%Fm)
758 CALL cp_libint_get_derivs(n_c, n_d, n_a, n_b, deriv, work_forces, a_mysize)
759 DO k = 4, 6
760 DO j = 1, mysize
761 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
762 work_forces(j, k + 3) + &
763 work_forces(j, k + 6))
764 END DO
765 END DO
766 DO k = 1, 12
767 DO j = 1, mysize
768 work(j, k) = work(j, k) + work_forces(j, k)
769 END DO
770 END DO
771 neris = neris + 12*mysize
772 IF (use_virial) THEN
773 CALL real_to_scaled(scoord(1:3), b, cell)
774 CALL real_to_scaled(scoord(4:6), a, cell)
775 CALL real_to_scaled(scoord(7:9), d, cell)
776 CALL real_to_scaled(scoord(10:12), c, cell)
777 DO k = 1, 12
778 DO j = 1, mysize
779 DO m = 1, 3
780 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
781 END DO
782 END DO
783 END DO
784 END IF
785
786 END DO
787 DO n = 1, 12
788 tmp_max = 0.0_dp
789 DO i = 1, mysize
790 tmp_max = max(tmp_max, abs(work(i, n)))
791 END DO
792 tmp_max = tmp_max*max_contraction
793 tmp_max_all = max(tmp_max_all, tmp_max)
794
795 DO j = 1, ncob
796 p1 = (j - 1)*ncoa
797 DO i = 1, ncoa
798 p2 = (p1 + i - 1)*ncod
799 DO l = 1, ncod
800 p3 = (p2 + l - 1)*ncoc
801 DO k = 1, ncoc
802 p4 = p3 + k
803 work2(i, j, k, l, full_perm4(n)) = work(p4, n)
804 END DO
805 END DO
806 END DO
807 END DO
808 END DO
809 IF (use_virial) THEN
810 DO n = 1, 12
811 tmp_max_virial = 0.0_dp
812 DO i = 1, mysize
813 tmp_max_virial = max(tmp_max_virial, &
814 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
815 END DO
816 tmp_max_virial = tmp_max_virial*max_contraction
817 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
818
819 DO j = 1, ncob
820 p1 = (j - 1)*ncoa
821 DO i = 1, ncoa
822 p2 = (p1 + i - 1)*ncod
823 DO l = 1, ncod
824 p3 = (p2 + l - 1)*ncoc
825 DO k = 1, ncoc
826 p4 = p3 + k
827 work2_virial(i, j, k, l, full_perm4(n), 1:3) = work_virial(p4, n, 1:3)
828 END DO
829 END DO
830 END DO
831 END DO
832 END DO
833 END IF
834 CASE (5)
835 DO i = 1, nproducts
836 a = pgf_product_list(i)%ra
837 b = pgf_product_list(i)%rb
838 c = pgf_product_list(i)%rc
839 d = pgf_product_list(i)%rd
840 rho = pgf_product_list(i)%Rho
841 rhoinv = pgf_product_list(i)%RhoInv
842 p = pgf_product_list(i)%P
843 q = pgf_product_list(i)%Q
844 w = pgf_product_list(i)%W
845 ab = pgf_product_list(i)%AB
846 cd = pgf_product_list(i)%CD
847 zetapetainv = pgf_product_list(i)%ZetapEtaInv
848 CALL build_deriv_data(deriv, c, d, a, b, &
849 zeta_c, zeta_d, zeta_a, zeta_b, &
850 etainv, zetainv, zetapetainv, rho, rhoinv, &
851 q, p, w, m_max, pgf_product_list(i)%Fm)
852 CALL cp_libint_get_derivs(n_b, n_a, n_d, n_c, deriv, work_forces, a_mysize)
853 DO k = 4, 6
854 DO j = 1, mysize
855 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
856 work_forces(j, k + 3) + &
857 work_forces(j, k + 6))
858 END DO
859 END DO
860 DO k = 1, 12
861 DO j = 1, mysize
862 work(j, k) = work(j, k) + work_forces(j, k)
863 END DO
864 END DO
865 neris = neris + 12*mysize
866 IF (use_virial) THEN
867 CALL real_to_scaled(scoord(1:3), c, cell)
868 CALL real_to_scaled(scoord(4:6), d, cell)
869 CALL real_to_scaled(scoord(7:9), a, cell)
870 CALL real_to_scaled(scoord(10:12), b, cell)
871 DO k = 1, 12
872 DO j = 1, mysize
873 DO m = 1, 3
874 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
875 END DO
876 END DO
877 END DO
878 END IF
879
880 END DO
881 DO n = 1, 12
882 tmp_max = 0.0_dp
883 DO i = 1, mysize
884 tmp_max = max(tmp_max, abs(work(i, n)))
885 END DO
886 tmp_max = tmp_max*max_contraction
887 tmp_max_all = max(tmp_max_all, tmp_max)
888
889 DO k = 1, ncoc
890 p1 = (k - 1)*ncod
891 DO l = 1, ncod
892 p2 = (p1 + l - 1)*ncoa
893 DO i = 1, ncoa
894 p3 = (p2 + i - 1)*ncob
895 DO j = 1, ncob
896 p4 = p3 + j
897 work2(i, j, k, l, full_perm5(n)) = work(p4, n)
898 END DO
899 END DO
900 END DO
901 END DO
902 END DO
903 IF (use_virial) THEN
904 DO n = 1, 12
905 tmp_max_virial = 0.0_dp
906 DO i = 1, mysize
907 tmp_max_virial = max(tmp_max_virial, &
908 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
909 END DO
910 tmp_max_virial = tmp_max_virial*max_contraction
911 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
912
913 DO k = 1, ncoc
914 p1 = (k - 1)*ncod
915 DO l = 1, ncod
916 p2 = (p1 + l - 1)*ncoa
917 DO i = 1, ncoa
918 p3 = (p2 + i - 1)*ncob
919 DO j = 1, ncob
920 p4 = p3 + j
921 work2_virial(i, j, k, l, full_perm5(n), 1:3) = work_virial(p4, n, 1:3)
922 END DO
923 END DO
924 END DO
925 END DO
926 END DO
927 END IF
928 CASE (6)
929 DO i = 1, nproducts
930 a = pgf_product_list(i)%ra
931 b = pgf_product_list(i)%rb
932 c = pgf_product_list(i)%rc
933 d = pgf_product_list(i)%rd
934 rho = pgf_product_list(i)%Rho
935 rhoinv = pgf_product_list(i)%RhoInv
936 p = pgf_product_list(i)%P
937 q = pgf_product_list(i)%Q
938 w = pgf_product_list(i)%W
939 ab = pgf_product_list(i)%AB
940 cd = pgf_product_list(i)%CD
941 zetapetainv = pgf_product_list(i)%ZetapEtaInv
942
943 CALL build_deriv_data(deriv, c, d, b, a, &
944 zeta_c, zeta_d, zeta_b, zeta_a, &
945 etainv, zetainv, zetapetainv, rho, rhoinv, &
946 q, p, w, m_max, pgf_product_list(i)%Fm)
947 CALL cp_libint_get_derivs(n_a, n_b, n_d, n_c, deriv, work_forces, a_mysize)
948 DO k = 4, 6
949 DO j = 1, mysize
950 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
951 work_forces(j, k + 3) + &
952 work_forces(j, k + 6))
953 END DO
954 END DO
955 DO k = 1, 12
956 DO j = 1, mysize
957 work(j, k) = work(j, k) + work_forces(j, k)
958 END DO
959 END DO
960 neris = neris + 12*mysize
961 IF (use_virial) THEN
962 CALL real_to_scaled(scoord(1:3), c, cell)
963 CALL real_to_scaled(scoord(4:6), d, cell)
964 CALL real_to_scaled(scoord(7:9), b, cell)
965 CALL real_to_scaled(scoord(10:12), a, cell)
966 DO k = 1, 12
967 DO j = 1, mysize
968 DO m = 1, 3
969 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
970 END DO
971 END DO
972 END DO
973 END IF
974
975 END DO
976 DO n = 1, 12
977 tmp_max = 0.0_dp
978 DO i = 1, mysize
979 tmp_max = max(tmp_max, abs(work(i, n)))
980 END DO
981 tmp_max = tmp_max*max_contraction
982 tmp_max_all = max(tmp_max_all, tmp_max)
983
984 DO k = 1, ncoc
985 p1 = (k - 1)*ncod
986 DO l = 1, ncod
987 p2 = (p1 + l - 1)*ncob
988 DO j = 1, ncob
989 p3 = (p2 + j - 1)*ncoa
990 DO i = 1, ncoa
991 p4 = p3 + i
992 work2(i, j, k, l, full_perm6(n)) = work(p4, n)
993 END DO
994 END DO
995 END DO
996 END DO
997 END DO
998 IF (use_virial) THEN
999 DO n = 1, 12
1000 tmp_max_virial = 0.0_dp
1001 DO i = 1, mysize
1002 tmp_max_virial = max(tmp_max_virial, &
1003 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
1004 END DO
1005 tmp_max_virial = tmp_max_virial*max_contraction
1006 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
1007
1008 DO k = 1, ncoc
1009 p1 = (k - 1)*ncod
1010 DO l = 1, ncod
1011 p2 = (p1 + l - 1)*ncob
1012 DO j = 1, ncob
1013 p3 = (p2 + j - 1)*ncoa
1014 DO i = 1, ncoa
1015 p4 = p3 + i
1016 work2_virial(i, j, k, l, full_perm6(n), 1:3) = work_virial(p4, n, 1:3)
1017 END DO
1018 END DO
1019 END DO
1020 END DO
1021 END DO
1022 END IF
1023 CASE (7)
1024 DO i = 1, nproducts
1025 a = pgf_product_list(i)%ra
1026 b = pgf_product_list(i)%rb
1027 c = pgf_product_list(i)%rc
1028 d = pgf_product_list(i)%rd
1029 rho = pgf_product_list(i)%Rho
1030 rhoinv = pgf_product_list(i)%RhoInv
1031 p = pgf_product_list(i)%P
1032 q = pgf_product_list(i)%Q
1033 w = pgf_product_list(i)%W
1034 ab = pgf_product_list(i)%AB
1035 cd = pgf_product_list(i)%CD
1036 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1037
1038 CALL build_deriv_data(deriv, d, c, a, b, &
1039 zeta_d, zeta_c, zeta_a, zeta_b, &
1040 etainv, zetainv, zetapetainv, rho, rhoinv, &
1041 q, p, w, m_max, pgf_product_list(i)%Fm)
1042 CALL cp_libint_get_derivs(n_b, n_a, n_c, n_d, deriv, work_forces, a_mysize)
1043 DO k = 4, 6
1044 DO j = 1, mysize
1045 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
1046 work_forces(j, k + 3) + &
1047 work_forces(j, k + 6))
1048 END DO
1049 END DO
1050 DO k = 1, 12
1051 DO j = 1, mysize
1052 work(j, k) = work(j, k) + work_forces(j, k)
1053 END DO
1054 END DO
1055 neris = neris + 12*mysize
1056 IF (use_virial) THEN
1057 CALL real_to_scaled(scoord(1:3), d, cell)
1058 CALL real_to_scaled(scoord(4:6), c, cell)
1059 CALL real_to_scaled(scoord(7:9), a, cell)
1060 CALL real_to_scaled(scoord(10:12), b, cell)
1061 DO k = 1, 12
1062 DO j = 1, mysize
1063 DO m = 1, 3
1064 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
1065 END DO
1066 END DO
1067 END DO
1068 END IF
1069
1070 END DO
1071 DO n = 1, 12
1072 tmp_max = 0.0_dp
1073 DO i = 1, mysize
1074 tmp_max = max(tmp_max, abs(work(i, n)))
1075 END DO
1076 tmp_max = tmp_max*max_contraction
1077 tmp_max_all = max(tmp_max_all, tmp_max)
1078
1079 DO l = 1, ncod
1080 p1 = (l - 1)*ncoc
1081 DO k = 1, ncoc
1082 p2 = (p1 + k - 1)*ncoa
1083 DO i = 1, ncoa
1084 p3 = (p2 + i - 1)*ncob
1085 DO j = 1, ncob
1086 p4 = p3 + j
1087 work2(i, j, k, l, full_perm7(n)) = work(p4, n)
1088 END DO
1089 END DO
1090 END DO
1091 END DO
1092 END DO
1093 IF (use_virial) THEN
1094 DO n = 1, 12
1095 tmp_max_virial = 0.0_dp
1096 DO i = 1, mysize
1097 tmp_max_virial = max(tmp_max_virial, &
1098 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
1099 END DO
1100 tmp_max_virial = tmp_max_virial*max_contraction
1101 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
1102
1103 DO l = 1, ncod
1104 p1 = (l - 1)*ncoc
1105 DO k = 1, ncoc
1106 p2 = (p1 + k - 1)*ncoa
1107 DO i = 1, ncoa
1108 p3 = (p2 + i - 1)*ncob
1109 DO j = 1, ncob
1110 p4 = p3 + j
1111 work2_virial(i, j, k, l, full_perm7(n), 1:3) = work_virial(p4, n, 1:3)
1112 END DO
1113 END DO
1114 END DO
1115 END DO
1116 END DO
1117 END IF
1118 CASE (8)
1119 DO i = 1, nproducts
1120 a = pgf_product_list(i)%ra
1121 b = pgf_product_list(i)%rb
1122 c = pgf_product_list(i)%rc
1123 d = pgf_product_list(i)%rd
1124 rho = pgf_product_list(i)%Rho
1125 rhoinv = pgf_product_list(i)%RhoInv
1126 p = pgf_product_list(i)%P
1127 q = pgf_product_list(i)%Q
1128 w = pgf_product_list(i)%W
1129 ab = pgf_product_list(i)%AB
1130 cd = pgf_product_list(i)%CD
1131 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1132
1133 CALL build_deriv_data(deriv, d, c, b, a, &
1134 zeta_d, zeta_c, zeta_b, zeta_a, &
1135 etainv, zetainv, zetapetainv, rho, rhoinv, &
1136 q, p, w, m_max, pgf_product_list(i)%Fm)
1137 CALL cp_libint_get_derivs(n_a, n_b, n_c, n_d, deriv, work_forces, a_mysize)
1138 DO k = 4, 6
1139 DO j = 1, mysize
1140 work_forces(j, k) = -1.0_dp*(work_forces(j, k - 3) + &
1141 work_forces(j, k + 3) + &
1142 work_forces(j, k + 6))
1143 END DO
1144 END DO
1145 DO k = 1, 12
1146 DO j = 1, mysize
1147 work(j, k) = work(j, k) + work_forces(j, k)
1148 END DO
1149 END DO
1150 neris = neris + 12*mysize
1151 IF (use_virial) THEN
1152 CALL real_to_scaled(scoord(1:3), d, cell)
1153 CALL real_to_scaled(scoord(4:6), c, cell)
1154 CALL real_to_scaled(scoord(7:9), b, cell)
1155 CALL real_to_scaled(scoord(10:12), a, cell)
1156 DO k = 1, 12
1157 DO j = 1, mysize
1158 DO m = 1, 3
1159 work_virial(j, k, m) = work_virial(j, k, m) + work_forces(j, k)*scoord(int((k - 1)/3)*3 + m)
1160 END DO
1161 END DO
1162 END DO
1163 END IF
1164
1165 END DO
1166 DO n = 1, 12
1167 tmp_max = 0.0_dp
1168 DO i = 1, mysize
1169 tmp_max = max(tmp_max, abs(work(i, n)))
1170 END DO
1171 tmp_max = tmp_max*max_contraction
1172 tmp_max_all = max(tmp_max_all, tmp_max)
1173
1174 DO l = 1, ncod
1175 p1 = (l - 1)*ncoc
1176 DO k = 1, ncoc
1177 p2 = (p1 + k - 1)*ncob
1178 DO j = 1, ncob
1179 p3 = (p2 + j - 1)*ncoa
1180 DO i = 1, ncoa
1181 p4 = p3 + i
1182 work2(i, j, k, l, full_perm8(n)) = work(p4, n)
1183 END DO
1184 END DO
1185 END DO
1186 END DO
1187 END DO
1188 IF (use_virial) THEN
1189 DO n = 1, 12
1190 tmp_max_virial = 0.0_dp
1191 DO i = 1, mysize
1192 tmp_max_virial = max(tmp_max_virial, &
1193 abs(work_virial(i, n, 1)), abs(work_virial(i, n, 2)), abs(work_virial(i, n, 3)))
1194 END DO
1195 tmp_max_virial = tmp_max_virial*max_contraction
1196 tmp_max_all_virial = max(tmp_max_all_virial, tmp_max_virial)
1197
1198 DO l = 1, ncod
1199 p1 = (l - 1)*ncoc
1200 DO k = 1, ncoc
1201 p2 = (p1 + k - 1)*ncob
1202 DO j = 1, ncob
1203 p3 = (p2 + j - 1)*ncoa
1204 DO i = 1, ncoa
1205 p4 = p3 + i
1206 work2_virial(i, j, k, l, full_perm8(n), 1:3) = work_virial(p4, n, 1:3)
1207 END DO
1208 END DO
1209 END DO
1210 END DO
1211 END DO
1212 END IF
1213 END SELECT
1214
1215 IF (.NOT. use_virial) THEN
1216 IF (tmp_max_all < eps_schwarz) RETURN
1217 END IF
1218
1219 IF (tmp_max_all >= eps_schwarz) THEN
1220 DO i = 1, 12
1221 primitives_tmp(:, :, :, :) = 0.0_dp
1222 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
1223 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2(1, 1, 1, 1, i), &
1224 sphi_a, &
1225 sphi_b, &
1226 sphi_c, &
1227 sphi_d, &
1228 primitives_tmp(1, 1, 1, 1), &
1229 buffer1, buffer2)
1230
1231 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
1232 s_offset_b + 1:s_offset_b + nsob*nl_b, &
1233 s_offset_c + 1:s_offset_c + nsoc*nl_c, &
1234 s_offset_d + 1:s_offset_d + nsod*nl_d, i) = &
1235 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
1236 s_offset_b + 1:s_offset_b + nsob*nl_b, &
1237 s_offset_c + 1:s_offset_c + nsoc*nl_c, &
1238 s_offset_d + 1:s_offset_d + nsod*nl_d, i) + primitives_tmp(:, :, :, :)
1239 END DO
1240 END IF
1241
1242 IF (use_virial .AND. tmp_max_all_virial >= eps_schwarz) THEN
1243 DO i = 1, 12
1244 DO m = 1, 3
1245 primitives_tmp_virial(:, :, :, :) = 0.0_dp
1246 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
1247 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2_virial(1, 1, 1, 1, i, m), &
1248 sphi_a, &
1249 sphi_b, &
1250 sphi_c, &
1251 sphi_d, &
1252 primitives_tmp_virial(1, 1, 1, 1), &
1253 buffer1, buffer2)
1254
1255 primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
1256 s_offset_b + 1:s_offset_b + nsob*nl_b, &
1257 s_offset_c + 1:s_offset_c + nsoc*nl_c, &
1258 s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) = &
1259 primitives_virial(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
1260 s_offset_b + 1:s_offset_b + nsob*nl_b, &
1261 s_offset_c + 1:s_offset_c + nsoc*nl_c, &
1262 s_offset_d + 1:s_offset_d + nsod*nl_d, i, m) + primitives_tmp_virial(:, :, :, :)
1263 END DO
1264 END DO
1265 END IF
1266
1267 END SUBROUTINE evaluate_deriv_eri
1268
1269! **************************************************************************************************
1270!> \brief Evaluates max-abs values of electron repulsion integrals for a primitive quartet
1271!> \param libint ...
1272!> \param A ...
1273!> \param B ...
1274!> \param C ...
1275!> \param D ...
1276!> \param Zeta_A ...
1277!> \param Zeta_B ...
1278!> \param Zeta_C ...
1279!> \param Zeta_D ...
1280!> \param n_a ...
1281!> \param n_b ...
1282!> \param n_c ...
1283!> \param n_d ...
1284!> \param max_val ...
1285!> \param potential_parameter ...
1286!> \param R1 ...
1287!> \param R2 ...
1288!> \param p_work ...
1289!> \par History
1290!> 03.2007 created [Manuel Guidon]
1291!> 08.2007 refactured permutation part [Manuel Guidon]
1292!> \author Manuel Guidon
1293! **************************************************************************************************
1294 SUBROUTINE evaluate_eri_screen(libint, A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D, &
1295 n_a, n_b, n_c, n_d, &
1296 max_val, potential_parameter, R1, R2, &
1297 p_work)
1298
1299 TYPE(cp_libint_t) :: libint
1300 REAL(dp), INTENT(IN) :: a(3), b(3), c(3), d(3), zeta_a, zeta_b, &
1301 zeta_c, zeta_d
1302 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d
1303 REAL(dp), INTENT(INOUT) :: max_val
1304 TYPE(hfx_potential_type) :: potential_parameter
1305 REAL(dp) :: r1, r2
1306 REAL(dp), DIMENSION(:), POINTER :: p_work
1307
1308 INTEGER :: a_mysize(1), i, m_max, mysize, perm_case
1309
1310 m_max = n_a + n_b + n_c + n_d
1311 mysize = nco(n_a)*nco(n_b)*nco(n_c)*nco(n_d)
1312 a_mysize = mysize
1313
1314 IF (m_max /= 0) THEN
1315 perm_case = 1
1316 IF (n_a < n_b) THEN
1317 perm_case = perm_case + 1
1318 END IF
1319 IF (n_c < n_d) THEN
1320 perm_case = perm_case + 2
1321 END IF
1322 IF (n_a + n_b > n_c + n_d) THEN
1323 perm_case = perm_case + 4
1324 END IF
1325
1326 SELECT CASE (perm_case)
1327 CASE (1)
1328 CALL build_quartet_data_screen(a, b, c, d, zeta_a, zeta_b, zeta_c, zeta_d, m_max, &
1329 potential_parameter, libint, r1, r2)
1330
1331 CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
1332 DO i = 1, mysize
1333 max_val = max(max_val, abs(p_work(i)))
1334 END DO
1335 CASE (2)
1336 CALL build_quartet_data_screen(b, a, c, d, zeta_b, zeta_a, zeta_c, zeta_d, m_max, &
1337 potential_parameter, libint, r1, r2)
1338
1339 CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
1340 DO i = 1, mysize
1341 max_val = max(max_val, abs(p_work(i)))
1342 END DO
1343 CASE (3)
1344 CALL build_quartet_data_screen(a, b, d, c, zeta_a, zeta_b, zeta_d, zeta_c, m_max, &
1345 potential_parameter, libint, r1, r2)
1346
1347 CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
1348 DO i = 1, mysize
1349 max_val = max(max_val, abs(p_work(i)))
1350 END DO
1351 CASE (4)
1352 CALL build_quartet_data_screen(b, a, d, c, zeta_b, zeta_a, zeta_d, zeta_c, m_max, &
1353 potential_parameter, libint, r1, r2)
1354
1355 CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
1356 DO i = 1, mysize
1357 max_val = max(max_val, abs(p_work(i)))
1358 END DO
1359 CASE (5)
1360 CALL build_quartet_data_screen(c, d, a, b, zeta_c, zeta_d, zeta_a, zeta_b, m_max, &
1361 potential_parameter, libint, r1, r2)
1362
1363 CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
1364 DO i = 1, mysize
1365 max_val = max(max_val, abs(p_work(i)))
1366 END DO
1367 CASE (6)
1368 CALL build_quartet_data_screen(c, d, b, a, zeta_c, zeta_d, zeta_b, zeta_a, m_max, &
1369 potential_parameter, libint, r1, r2)
1370
1371 CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
1372 DO i = 1, mysize
1373 max_val = max(max_val, abs(p_work(i)))
1374 END DO
1375 CASE (7)
1376 CALL build_quartet_data_screen(d, c, a, b, zeta_d, zeta_c, zeta_a, zeta_b, m_max, &
1377 potential_parameter, libint, r1, r2)
1378
1379 CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
1380 DO i = 1, mysize
1381 max_val = max(max_val, abs(p_work(i)))
1382 END DO
1383 CASE (8)
1384 CALL build_quartet_data_screen(d, c, b, a, zeta_d, zeta_c, zeta_b, zeta_a, m_max, &
1385 potential_parameter, libint, r1, r2)
1386
1387 CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
1388 DO i = 1, mysize
1389 max_val = max(max_val, abs(p_work(i)))
1390 END DO
1391 END SELECT
1392 ELSE
1393 CALL build_quartet_data_screen(a, b, c, d, zeta_a, zeta_b, zeta_c, zeta_d, m_max, &
1394 potential_parameter, libint, r1, r2)
1395 max_val = abs(get_ssss_f_val(libint))
1396 END IF
1397
1398 END SUBROUTINE evaluate_eri_screen
1399
1400! **************************************************************************************************
1401!> \brief Evaluate electron repulsion integrals for a primitive quartet
1402!> \param libint ...
1403!> \param nproducts ...
1404!> \param pgf_product_list ...
1405!> \param n_a ...
1406!> \param n_b ...
1407!> \param n_c ...
1408!> \param n_d ...
1409!> \param ncoa ...
1410!> \param ncob ...
1411!> \param ncoc ...
1412!> \param ncod ...
1413!> \param nsgfa ...
1414!> \param nsgfb ...
1415!> \param nsgfc ...
1416!> \param nsgfd ...
1417!> \param primitives ...
1418!> \param max_contraction ...
1419!> \param tmp_max ...
1420!> \param eps_schwarz ...
1421!> \param neris ...
1422!> \param ZetaInv ...
1423!> \param EtaInv ...
1424!> \param s_offset_a ...
1425!> \param s_offset_b ...
1426!> \param s_offset_c ...
1427!> \param s_offset_d ...
1428!> \param nl_a ...
1429!> \param nl_b ...
1430!> \param nl_c ...
1431!> \param nl_d ...
1432!> \param nsoa ...
1433!> \param nsob ...
1434!> \param nsoc ...
1435!> \param nsod ...
1436!> \param sphi_a ...
1437!> \param sphi_b ...
1438!> \param sphi_c ...
1439!> \param sphi_d ...
1440!> \param work ...
1441!> \param work2 ...
1442!> \param buffer1 ...
1443!> \param buffer2 ...
1444!> \param primitives_tmp ...
1445!> \param p_work ...
1446!> \par History
1447!> 11.2006 created [Manuel Guidon]
1448!> 08.2007 refactured permutation part [Manuel Guidon]
1449!> \author Manuel Guidon
1450! **************************************************************************************************
1451 SUBROUTINE evaluate_eri(libint, nproducts, pgf_product_list, &
1452 n_a, n_b, n_c, n_d, &
1453 ncoa, ncob, ncoc, ncod, &
1454 nsgfa, nsgfb, nsgfc, nsgfd, &
1455 primitives, max_contraction, tmp_max, &
1456 eps_schwarz, neris, &
1457 ZetaInv, EtaInv, &
1458 s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
1459 nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, &
1460 sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, &
1461 primitives_tmp, p_work)
1462
1463 TYPE(cp_libint_t) :: libint
1464 INTEGER, INTENT(IN) :: nproducts
1465 TYPE(hfx_pgf_product_list), DIMENSION(nproducts) :: pgf_product_list
1466 INTEGER, INTENT(IN) :: n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, &
1467 ncod, nsgfa, nsgfb, nsgfc, nsgfd
1468 REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitives
1469 REAL(dp) :: max_contraction, tmp_max, eps_schwarz
1470 INTEGER(int_8) :: neris
1471 REAL(dp), INTENT(IN) :: zetainv, etainv
1472 INTEGER :: s_offset_a, s_offset_b, s_offset_c, &
1473 s_offset_d, nl_a, nl_b, nl_c, nl_d, &
1474 nsoa, nsob, nsoc, nsod
1475 REAL(dp), DIMENSION(ncoa, nsoa*nl_a), INTENT(IN) :: sphi_a
1476 REAL(dp), DIMENSION(ncob, nsob*nl_b), INTENT(IN) :: sphi_b
1477 REAL(dp), DIMENSION(ncoc, nsoc*nl_c), INTENT(IN) :: sphi_c
1478 REAL(dp), DIMENSION(ncod, nsod*nl_d), INTENT(IN) :: sphi_d
1479 REAL(dp) :: work(ncoa*ncob*ncoc*ncod), &
1480 work2(ncoa, ncob, ncoc, ncod)
1481 REAL(dp), DIMENSION(ncoa*ncob*ncoc*ncod) :: buffer1, buffer2
1482 REAL(dp), DIMENSION(nsoa*nl_a, nsob*nl_b, nsoc*& nl_c, nsod*nl_d) :: primitives_tmp
1483 REAL(dp), DIMENSION(:), POINTER :: p_work
1484
1485 INTEGER :: a_mysize(1), i, j, k, l, m_max, mysize, &
1486 p1, p2, p3, p4, perm_case
1487 REAL(dp) :: a(3), ab(3), b(3), c(3), cd(3), d(3), &
1488 p(3), q(3), rho, rhoinv, w(3), &
1489 zetapetainv
1490 REAL(kind=dp), DIMENSION(prim_data_f_size) :: f
1491
1492 m_max = n_a + n_b + n_c + n_d
1493 mysize = ncoa*ncob*ncoc*ncod
1494 a_mysize = mysize
1495
1496 work = 0.0_dp
1497 IF (m_max /= 0) THEN
1498 perm_case = 1
1499 IF (n_a < n_b) THEN
1500 perm_case = perm_case + 1
1501 END IF
1502 IF (n_c < n_d) THEN
1503 perm_case = perm_case + 2
1504 END IF
1505 IF (n_a + n_b > n_c + n_d) THEN
1506 perm_case = perm_case + 4
1507 END IF
1508 SELECT CASE (perm_case)
1509 CASE (1)
1510 DO i = 1, nproducts
1511 a = pgf_product_list(i)%ra
1512 b = pgf_product_list(i)%rb
1513 c = pgf_product_list(i)%rc
1514 d = pgf_product_list(i)%rd
1515 rho = pgf_product_list(i)%Rho
1516 rhoinv = pgf_product_list(i)%RhoInv
1517 p = pgf_product_list(i)%P
1518 q = pgf_product_list(i)%Q
1519 w = pgf_product_list(i)%W
1520 ab = pgf_product_list(i)%AB
1521 cd = pgf_product_list(i)%CD
1522 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1523 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1524
1525 CALL build_quartet_data(libint, a, b, c, d, zetainv, etainv, zetapetainv, rho, &
1526 p, q, w, m_max, f)
1527 CALL cp_libint_get_eris(n_d, n_c, n_b, n_a, libint, p_work, a_mysize)
1528 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1529 neris = neris + mysize
1530 END DO
1531 DO i = 1, mysize
1532 tmp_max = max(tmp_max, abs(work(i)))
1533 END DO
1534 tmp_max = tmp_max*max_contraction
1535 IF (tmp_max < eps_schwarz) THEN
1536 RETURN
1537 END IF
1538
1539 DO i = 1, ncoa
1540 p1 = (i - 1)*ncob
1541 DO j = 1, ncob
1542 p2 = (p1 + j - 1)*ncoc
1543 DO k = 1, ncoc
1544 p3 = (p2 + k - 1)*ncod
1545 DO l = 1, ncod
1546 p4 = p3 + l
1547 work2(i, j, k, l) = work(p4)
1548 END DO
1549 END DO
1550 END DO
1551 END DO
1552 CASE (2)
1553 DO i = 1, nproducts
1554 a = pgf_product_list(i)%ra
1555 b = pgf_product_list(i)%rb
1556 c = pgf_product_list(i)%rc
1557 d = pgf_product_list(i)%rd
1558 rho = pgf_product_list(i)%Rho
1559 rhoinv = pgf_product_list(i)%RhoInv
1560 p = pgf_product_list(i)%P
1561 q = pgf_product_list(i)%Q
1562 w = pgf_product_list(i)%W
1563 ab = pgf_product_list(i)%AB
1564 cd = pgf_product_list(i)%CD
1565 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1566 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1567
1568 CALL build_quartet_data(libint, b, a, c, d, &
1569 zetainv, etainv, zetapetainv, rho, &
1570 p, q, w, m_max, f)
1571
1572 CALL cp_libint_get_eris(n_d, n_c, n_a, n_b, libint, p_work, a_mysize)
1573 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1574 neris = neris + mysize
1575 END DO
1576 DO i = 1, mysize
1577 tmp_max = max(tmp_max, abs(work(i)))
1578 END DO
1579 tmp_max = tmp_max*max_contraction
1580 IF (tmp_max < eps_schwarz) THEN
1581 RETURN
1582 END IF
1583
1584 DO j = 1, ncob
1585 p1 = (j - 1)*ncoa
1586 DO i = 1, ncoa
1587 p2 = (p1 + i - 1)*ncoc
1588 DO k = 1, ncoc
1589 p3 = (p2 + k - 1)*ncod
1590 DO l = 1, ncod
1591 p4 = p3 + l
1592 work2(i, j, k, l) = work(p4)
1593 END DO
1594 END DO
1595 END DO
1596 END DO
1597 CASE (3)
1598 DO i = 1, nproducts
1599 a = pgf_product_list(i)%ra
1600 b = pgf_product_list(i)%rb
1601 c = pgf_product_list(i)%rc
1602 d = pgf_product_list(i)%rd
1603 rho = pgf_product_list(i)%Rho
1604 rhoinv = pgf_product_list(i)%RhoInv
1605 p = pgf_product_list(i)%P
1606 q = pgf_product_list(i)%Q
1607 w = pgf_product_list(i)%W
1608 ab = pgf_product_list(i)%AB
1609 cd = pgf_product_list(i)%CD
1610 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1611 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1612
1613 CALL build_quartet_data(libint, a, b, d, c, &
1614 zetainv, etainv, zetapetainv, rho, &
1615 p, q, w, m_max, f)
1616
1617 CALL cp_libint_get_eris(n_c, n_d, n_b, n_a, libint, p_work, a_mysize)
1618 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1619 neris = neris + mysize
1620 END DO
1621 DO i = 1, mysize
1622 tmp_max = max(tmp_max, abs(work(i)))
1623 END DO
1624 tmp_max = tmp_max*max_contraction
1625 IF (tmp_max < eps_schwarz) THEN
1626 RETURN
1627 END IF
1628
1629 DO i = 1, ncoa
1630 p1 = (i - 1)*ncob
1631 DO j = 1, ncob
1632 p2 = (p1 + j - 1)*ncod
1633 DO l = 1, ncod
1634 p3 = (p2 + l - 1)*ncoc
1635 DO k = 1, ncoc
1636 p4 = p3 + k
1637 work2(i, j, k, l) = work(p4)
1638 END DO
1639 END DO
1640 END DO
1641 END DO
1642 CASE (4)
1643 DO i = 1, nproducts
1644 a = pgf_product_list(i)%ra
1645 b = pgf_product_list(i)%rb
1646 c = pgf_product_list(i)%rc
1647 d = pgf_product_list(i)%rd
1648 rho = pgf_product_list(i)%Rho
1649 rhoinv = pgf_product_list(i)%RhoInv
1650 p = pgf_product_list(i)%P
1651 q = pgf_product_list(i)%Q
1652 w = pgf_product_list(i)%W
1653 ab = pgf_product_list(i)%AB
1654 cd = pgf_product_list(i)%CD
1655 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1656 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1657
1658 CALL build_quartet_data(libint, b, a, d, c, &
1659 zetainv, etainv, zetapetainv, rho, &
1660 p, q, w, m_max, f)
1661
1662 CALL cp_libint_get_eris(n_c, n_d, n_a, n_b, libint, p_work, a_mysize)
1663 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1664 neris = neris + mysize
1665 END DO
1666 DO i = 1, mysize
1667 tmp_max = max(tmp_max, abs(work(i)))
1668 END DO
1669 tmp_max = tmp_max*max_contraction
1670 IF (tmp_max < eps_schwarz) THEN
1671 RETURN
1672 END IF
1673
1674 DO j = 1, ncob
1675 p1 = (j - 1)*ncoa
1676 DO i = 1, ncoa
1677 p2 = (p1 + i - 1)*ncod
1678 DO l = 1, ncod
1679 p3 = (p2 + l - 1)*ncoc
1680 DO k = 1, ncoc
1681 p4 = p3 + k
1682 work2(i, j, k, l) = work(p4)
1683 END DO
1684 END DO
1685 END DO
1686 END DO
1687 CASE (5)
1688 DO i = 1, nproducts
1689 a = pgf_product_list(i)%ra
1690 b = pgf_product_list(i)%rb
1691 c = pgf_product_list(i)%rc
1692 d = pgf_product_list(i)%rd
1693 rho = pgf_product_list(i)%Rho
1694 rhoinv = pgf_product_list(i)%RhoInv
1695 p = pgf_product_list(i)%P
1696 q = pgf_product_list(i)%Q
1697 w = pgf_product_list(i)%W
1698 ab = pgf_product_list(i)%AB
1699 cd = pgf_product_list(i)%CD
1700 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1701 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1702
1703 CALL build_quartet_data(libint, c, d, a, b, &
1704 etainv, zetainv, zetapetainv, rho, &
1705 q, p, w, m_max, f)
1706
1707 CALL cp_libint_get_eris(n_b, n_a, n_d, n_c, libint, p_work, a_mysize)
1708 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1709 neris = neris + mysize
1710 END DO
1711 DO i = 1, mysize
1712 tmp_max = max(tmp_max, abs(work(i)))
1713 END DO
1714 tmp_max = tmp_max*max_contraction
1715 IF (tmp_max < eps_schwarz) THEN
1716 RETURN
1717 END IF
1718
1719 DO k = 1, ncoc
1720 p1 = (k - 1)*ncod
1721 DO l = 1, ncod
1722 p2 = (p1 + l - 1)*ncoa
1723 DO i = 1, ncoa
1724 p3 = (p2 + i - 1)*ncob
1725 DO j = 1, ncob
1726 p4 = p3 + j
1727 work2(i, j, k, l) = work(p4)
1728 END DO
1729 END DO
1730 END DO
1731 END DO
1732 CASE (6)
1733 DO i = 1, nproducts
1734 a = pgf_product_list(i)%ra
1735 b = pgf_product_list(i)%rb
1736 c = pgf_product_list(i)%rc
1737 d = pgf_product_list(i)%rd
1738 rho = pgf_product_list(i)%Rho
1739 rhoinv = pgf_product_list(i)%RhoInv
1740 p = pgf_product_list(i)%P
1741 q = pgf_product_list(i)%Q
1742 w = pgf_product_list(i)%W
1743 ab = pgf_product_list(i)%AB
1744 cd = pgf_product_list(i)%CD
1745 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1746 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1747
1748 CALL build_quartet_data(libint, c, d, b, a, &
1749 etainv, zetainv, zetapetainv, rho, &
1750 q, p, w, m_max, f)
1751
1752 CALL cp_libint_get_eris(n_a, n_b, n_d, n_c, libint, p_work, a_mysize)
1753 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1754 neris = neris + mysize
1755 END DO
1756 DO i = 1, mysize
1757 tmp_max = max(tmp_max, abs(work(i)))
1758 END DO
1759 tmp_max = tmp_max*max_contraction
1760 IF (tmp_max < eps_schwarz) THEN
1761 RETURN
1762 END IF
1763
1764 DO k = 1, ncoc
1765 p1 = (k - 1)*ncod
1766 DO l = 1, ncod
1767 p2 = (p1 + l - 1)*ncob
1768 DO j = 1, ncob
1769 p3 = (p2 + j - 1)*ncoa
1770 DO i = 1, ncoa
1771 p4 = p3 + i
1772 work2(i, j, k, l) = work(p4)
1773 END DO
1774 END DO
1775 END DO
1776 END DO
1777 CASE (7)
1778 DO i = 1, nproducts
1779 a = pgf_product_list(i)%ra
1780 b = pgf_product_list(i)%rb
1781 c = pgf_product_list(i)%rc
1782 d = pgf_product_list(i)%rd
1783 rho = pgf_product_list(i)%Rho
1784 rhoinv = pgf_product_list(i)%RhoInv
1785 p = pgf_product_list(i)%P
1786 q = pgf_product_list(i)%Q
1787 w = pgf_product_list(i)%W
1788 ab = pgf_product_list(i)%AB
1789 cd = pgf_product_list(i)%CD
1790 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1791 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1792
1793 CALL build_quartet_data(libint, d, c, a, b, &
1794 etainv, zetainv, zetapetainv, rho, &
1795 q, p, w, m_max, f)
1796
1797 CALL cp_libint_get_eris(n_b, n_a, n_c, n_d, libint, p_work, a_mysize)
1798 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1799 neris = neris + mysize
1800 END DO
1801 DO i = 1, mysize
1802 tmp_max = max(tmp_max, abs(work(i)))
1803 END DO
1804 tmp_max = tmp_max*max_contraction
1805 IF (tmp_max < eps_schwarz) THEN
1806 RETURN
1807 END IF
1808
1809 DO l = 1, ncod
1810 p1 = (l - 1)*ncoc
1811 DO k = 1, ncoc
1812 p2 = (p1 + k - 1)*ncoa
1813 DO i = 1, ncoa
1814 p3 = (p2 + i - 1)*ncob
1815 DO j = 1, ncob
1816 p4 = p3 + j
1817 work2(i, j, k, l) = work(p4)
1818 END DO
1819 END DO
1820 END DO
1821 END DO
1822 CASE (8)
1823 DO i = 1, nproducts
1824 a = pgf_product_list(i)%ra
1825 b = pgf_product_list(i)%rb
1826 c = pgf_product_list(i)%rc
1827 d = pgf_product_list(i)%rd
1828 rho = pgf_product_list(i)%Rho
1829 rhoinv = pgf_product_list(i)%RhoInv
1830 p = pgf_product_list(i)%P
1831 q = pgf_product_list(i)%Q
1832 w = pgf_product_list(i)%W
1833 ab = pgf_product_list(i)%AB
1834 cd = pgf_product_list(i)%CD
1835 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1836 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1837
1838 CALL build_quartet_data(libint, d, c, b, a, &
1839 etainv, zetainv, zetapetainv, rho, &
1840 q, p, w, m_max, f)
1841
1842 CALL cp_libint_get_eris(n_a, n_b, n_c, n_d, libint, p_work, a_mysize)
1843 work(1:mysize) = work(1:mysize) + p_work(1:mysize)
1844 neris = neris + mysize
1845 END DO
1846 DO i = 1, mysize
1847 tmp_max = max(tmp_max, abs(work(i)))
1848 END DO
1849 tmp_max = tmp_max*max_contraction
1850 IF (tmp_max < eps_schwarz) THEN
1851 RETURN
1852 END IF
1853
1854 DO l = 1, ncod
1855 p1 = (l - 1)*ncoc
1856 DO k = 1, ncoc
1857 p2 = (p1 + k - 1)*ncob
1858 DO j = 1, ncob
1859 p3 = (p2 + j - 1)*ncoa
1860 DO i = 1, ncoa
1861 p4 = p3 + i
1862 work2(i, j, k, l) = work(p4)
1863 END DO
1864 END DO
1865 END DO
1866 END DO
1867 END SELECT
1868 ELSE
1869 DO i = 1, nproducts
1870 a = pgf_product_list(i)%ra
1871 b = pgf_product_list(i)%rb
1872 c = pgf_product_list(i)%rc
1873 d = pgf_product_list(i)%rd
1874 rho = pgf_product_list(i)%Rho
1875 rhoinv = pgf_product_list(i)%RhoInv
1876 p = pgf_product_list(i)%P
1877 q = pgf_product_list(i)%Q
1878 w = pgf_product_list(i)%W
1879 zetapetainv = pgf_product_list(i)%ZetapEtaInv
1880 f(1:m_max + 1) = pgf_product_list(i)%Fm(1:m_max + 1)
1881
1882 CALL build_quartet_data(libint, a, b, c, d, & !todo: check
1883 zetainv, etainv, zetapetainv, rho, &
1884 p, q, w, m_max, f)
1885 work(1) = work(1) + f(1)
1886 neris = neris + mysize
1887 END DO
1888 work2(1, 1, 1, 1) = work(1)
1889 tmp_max = max_contraction*abs(work(1))
1890 IF (tmp_max < eps_schwarz) RETURN
1891 END IF
1892
1893 IF (tmp_max < eps_schwarz) RETURN
1894 primitives_tmp = 0.0_dp
1895
1896 CALL contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, &
1897 n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work2, &
1898 sphi_a, &
1899 sphi_b, &
1900 sphi_c, &
1901 sphi_d, &
1902 primitives_tmp, &
1903 buffer1, buffer2)
1904
1905 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
1906 s_offset_b + 1:s_offset_b + nsob*nl_b, &
1907 s_offset_c + 1:s_offset_c + nsoc*nl_c, &
1908 s_offset_d + 1:s_offset_d + nsod*nl_d) = &
1909 primitives(s_offset_a + 1:s_offset_a + nsoa*nl_a, &
1910 s_offset_b + 1:s_offset_b + nsob*nl_b, &
1911 s_offset_c + 1:s_offset_c + nsoc*nl_c, &
1912 s_offset_d + 1:s_offset_d + nsod*nl_d) + primitives_tmp(:, :, :, :)
1913
1914 END SUBROUTINE evaluate_eri
1915
1916! **************************************************************************************************
1917!> \brief Fill data structure used in libint
1918!> \param libint ...
1919!> \param A ...
1920!> \param B ...
1921!> \param C ...
1922!> \param D ...
1923!> \param ZetaInv ...
1924!> \param EtaInv ...
1925!> \param ZetapEtaInv ...
1926!> \param Rho ...
1927!> \param P ...
1928!> \param Q ...
1929!> \param W ...
1930!> \param m_max ...
1931!> \param F ...
1932!> \par History
1933!> 03.2007 created [Manuel Guidon]
1934!> \author Manuel Guidon
1935! **************************************************************************************************
1936 SUBROUTINE build_quartet_data(libint, A, B, C, D, &
1937 ZetaInv, EtaInv, ZetapEtaInv, Rho, &
1938 P, Q, W, m_max, F)
1939
1940 TYPE(cp_libint_t) :: libint
1941 REAL(kind=dp), INTENT(IN) :: a(3), b(3), c(3), d(3)
1942 REAL(dp), INTENT(IN) :: zetainv, etainv, zetapetainv, rho, p(3), &
1943 q(3), w(3)
1944 INTEGER, INTENT(IN) :: m_max
1945 REAL(kind=dp), DIMENSION(:) :: f
1946
1947 CALL cp_libint_set_params_eri(libint, a, b, c, d, zetainv, etainv, zetapetainv, rho, p, q, w, m_max, f)
1948 END SUBROUTINE build_quartet_data
1949
1950END MODULE hfx_libint_interface
1951
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:486
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
hfx_contraction_methods
Contains routines for contraction without dgemms. PLEASE DO NOT MODIFY. \notes Contains specific rout...
Definition hfx_contraction_methods.F:22
hfx_contraction_methods::contract
subroutine, public contract(ncoa, ncob, ncoc, ncod, nsoa, nsob, nsoc, nsod, n_a, n_b, n_c, n_d, nl_a, nl_b, nl_c, nl_d, work, sphi_a, sphi_b, sphi_c, sphi_d, primitives, buffer1, buffer2)
...
Definition hfx_contraction_methods.F:71
hfx_libint_interface
Interface to the Libint-Library.
Definition hfx_libint_interface.F:15
hfx_libint_interface::evaluate_eri
subroutine, public evaluate_eri(libint, nproducts, pgf_product_list, n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd, primitives, max_contraction, tmp_max, eps_schwarz, neris, zetainv, etainv, s_offset_a, s_offset_b, s_offset_c, s_offset_d, nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, primitives_tmp, p_work)
Evaluate electron repulsion integrals for a primitive quartet.
Definition hfx_libint_interface.F:1464
hfx_libint_interface::evaluate_eri_screen
subroutine, public evaluate_eri_screen(libint, a, b, c, d, zeta_a, zeta_b, zeta_c, zeta_d, n_a, n_b, n_c, n_d, max_val, potential_parameter, r1, r2, p_work)
Evaluates max-abs values of electron repulsion integrals for a primitive quartet.
Definition hfx_libint_interface.F:1300
hfx_libint_interface::evaluate_deriv_eri
subroutine, public evaluate_deriv_eri(deriv, nproducts, pgf_product_list, n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd, primitives, max_contraction, tmp_max_all, eps_schwarz, neris, zeta_a, zeta_b, zeta_c, zeta_d, zetainv, etainv, s_offset_a, s_offset_b, s_offset_c, s_offset_d, nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, sphi_a, sphi_b, sphi_c, sphi_d, work, work2, work_forces, buffer1, buffer2, primitives_tmp, use_virial, work_virial, work2_virial, primitives_tmp_virial, primitives_virial, cell, tmp_max_all_virial)
Evaluates derivatives of electron repulsion integrals for a primitive quartet.
Definition hfx_libint_interface.F:391
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
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_mix_cl
integer, parameter, public do_potential_mix_cl
Definition input_constants.F:801
input_constants::do_potential_gaussian
integer, parameter, public do_potential_gaussian
Definition input_constants.F:801
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:801
input_constants::do_potential_mix_lg
integer, parameter, public do_potential_mix_lg
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_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
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::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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_eris
subroutine, public cp_libint_get_eris(n_d, n_c, n_b, n_a, lib, p_work, a_mysize)
...
Definition libint_wrapper.F:844
libint_wrapper::get_ssss_f_val
real(kind=dp) function, public get_ssss_f_val(lib)
Definition libint_wrapper.F:1070
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::prim_data_f_size
integer, parameter, public prim_data_f_size
Definition libint_wrapper.F:63
libint_wrapper::cp_libint_set_params_eri_screen
subroutine, public cp_libint_set_params_eri_screen(libint, a, b, c, d, p, q, w, zetainv, etainv, zetapetainv, rho, m_max, f)
Definition libint_wrapper.F:81
libint_wrapper::cp_libint_get_derivs
subroutine, public cp_libint_get_derivs(n_d, n_c, n_b, n_a, lib, work_forces, a_mysize)
...
Definition libint_wrapper.F:1032
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
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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
hfx_types::hfx_pgf_product_list
Definition hfx_types.F:329
hfx_types::hfx_potential_type
Definition hfx_types.F:138
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:69