(git:1f285aa)
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 ! **************************************************************************************************
15 MODULE 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 
59 CONTAINS
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 
1950 END 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_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_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_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_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::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
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
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