1) Minor refactoring of libint 2- and 3-centers integrals 2) Fixed

4-center identity potential if HFX
cp2k · Nov 26, 2019 · c38db1d · c38db1d
1 parent 54b179d
commit c38db1d
Show file tree

Hide file tree

Showing 5 changed files with 148 additions and 102 deletions.
diff --git a/src/hfx_libint_interface.F b/src/hfx_libint_interface.F
@@ -7,6 +7,7 @@
 !> \brief Interface to the Libint-Library
 !> \par History
 !>      11.2006 created [Manuel Guidon]
+!>      11.2019 Fixed potential_id initial values (A. Bussy)
 !> \author Manuel Guidon
 ! **************************************************************************************************
 MODULE hfx_libint_interface
@@ -86,9 +87,9 @@ SUBROUTINE build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D,
 
       INTEGER                                            :: i
       LOGICAL                                            :: use_gamma
-      REAL(KIND=dp) :: AB(3), AB2, CD(3), CD2, Eta, EtaInv, factor, omega2, omega_corr, &
-         omega_corr2, P(3), PQ(3), PQ2, Q(3), R, R1, R2, Rho, RhoInv, S1234, T, tmp, W(3), Zeta, &
-         ZetaInv, ZetapEtaInv
+      REAL(KIND=dp) :: AB(3), AB2, CD(3), CD2, den, Eta, EtaInv, factor, G(3), num, omega2, &
+         omega_corr, omega_corr2, P(3), PQ(3), PQ2, Q(3), R, R1, R2, Rho, RhoInv, S1234, ssss, T, &
+         tmp, W(3), Zeta, ZetaInv, ZetapEtaInv
       REAL(KIND=dp), DIMENSION(prim_data_f_size)         :: F, Fm
 
       Zeta = Zeta_A + Zeta_B
@@ -236,9 +237,20 @@ SUBROUTINE build_quartet_data_screen(A, B, C, D, Zeta_A, Zeta_B, Zeta_C, Zeta_D,
          END DO
          factor = 1.0_dp
       CASE (do_potential_id)
-         F(1) = (Pi*RhoInv)**(1.5_dp)
-         F(2:m_max + 1) = 0.0_dp
+         ssss = -Zeta_A*Zeta_B*ZetaInv*AB2
+
+         num = (Zeta_A+Zeta_B)*Zeta_C
+         den = Zeta_A+Zeta_B+Zeta_C
+         ssss = ssss-num/den*SUM((P-C)**2)
+
+         G(:) = (Zeta_A*A(:)+Zeta_B*B(:)+Zeta_C*C(:))/den
+         num = den*Zeta_D
+         den = den+Zeta_D
+         ssss = ssss-num/den*SUM((G-D)**2)
+
+         F(:) = EXP(ssss)
          factor = 1.0_dp
+         IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
       END SELECT
 
       tmp = (Pi*ZetapEtaInv)**3

diff --git a/src/hfx_pair_list_methods.F b/src/hfx_pair_list_methods.F
@@ -7,6 +7,7 @@
 !> \brief Rountines for optimizing load balance between processes in HFX calculations
 !> \par History
 !>      04.2008 created [Manuel Guidon]
+!>      11.2019 fixed initial value for potential_id (A. Bussy)
 !> \author Manuel Guidon
 ! **************************************************************************************************
 MODULE hfx_pair_list_methods
@@ -77,10 +78,10 @@ SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
 
       INTEGER                                            :: i, j, k, l, nimages1, nimages2, tmp_i4
       LOGICAL                                            :: use_gamma
-      REAL(dp) :: C11(3), Eta, EtaInv, factor, Fm(prim_data_f_size), omega2, omega_corr, &
-         omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), Q(3), R, R1, R2, ra(3), rb(3), rc(3), &
-         rd(3), Rho, RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), T, temp(3), temp_CC(3), &
-         temp_DD(3), tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv
+      REAL(dp) :: C11(3), den, Eta, EtaInv, factor, Fm(prim_data_f_size), G(3), num, omega2, &
+         omega_corr, omega_corr2, P(3), pgf_max_1, pgf_max_2, PQ(3), Q(3), R, R1, R2, ra(3), &
+         rb(3), rc(3), rd(3), Rho, RhoInv, rpq2, S1234, S1234a, S1234b, shift(3), ssss, T, &
+         temp(3), temp_CC(3), temp_DD(3), tmp, tmp_D(3), W(3), Zeta1, Zeta_C, Zeta_D, ZetapEtaInv
       TYPE(hfx_pgf_product_list), ALLOCATABLE, &
          DIMENSION(:)                                    :: tmp_product_list
 
@@ -270,11 +271,22 @@ SUBROUTINE build_pgf_product_list(list1, list2, product_list, nproducts, &
                   END DO
                   factor = 1.0_dp
                CASE (do_potential_id)
-                  product_list(nproducts)%Fm(1) = (Pi*RhoInv)**(1.5_dp)
-                  DO l = 2, m_max + 1
-                     product_list(nproducts)%Fm(l) = 0.0_dp
-                  END DO
+                  num = list1%zeta*list1%zetb
+                  den = list1%zeta+list1%zetb
+                  ssss = -num/den*SUM((ra-rb)**2)
+
+                  num = den*Zeta_C
+                  den = den+Zeta_C
+                  ssss = ssss-num/den*SUM((P-rc)**2)
+
+                  G(:) = (list1%zeta*ra(:)+list1%zetb*rb(:)+Zeta_C*rc(:))/den
+                  num = den*Zeta_D
+                  den = den+Zeta_D
+                  ssss = ssss-num/den*SUM((G-rd)**2)
+
+                  product_list(nproducts)%Fm(:) = EXP(ssss)
                   factor = 1.0_dp
+                  IF (S1234 > EPSILON(0.0_dp)) factor = 1.0_dp/S1234
                END SELECT
 
                tmp = (Pi*ZetapEtaInv)**3