Merge pull request #354 from AbdAmmar/dev-stable

Dev stable
QuantumPackage · Oct 18, 2024 · 340ed7a · 340ed7a
2 parents ed6e77d + 7dff773
commit 340ed7a
Show file tree

Hide file tree

Showing 12 changed files with 1,210 additions and 743 deletions.
diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg
@@ -73,7 +73,7 @@ doc: If true, use cgtos for AO integrals
 interface: ezfio
 default: False
 
-[ao_expo_im_cgtos]
+[ao_expo_im]
 type: double precision
 doc: imag part for Exponents for each primitive of each cGTOs |AO|
 size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
@@ -82,12 +82,12 @@ interface: ezfio, provider
 [ao_expo_pw]
 type: double precision
 doc: plane wave part for each primitive GTOs |AO|
-size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max)
+size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max)
 interface: ezfio, provider
 
 [ao_expo_phase]
 type: double precision
 doc: phase shift for each primitive GTOs |AO|
-size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max)
+size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max)
 interface: ezfio, provider
 
diff --git a/src/ao_basis/README.rst b/src/ao_basis/README.rst
@@ -40,6 +40,6 @@ Complex Gaussian-Type Orbitals (cGTOs) are also supported:
    \chi_i(\mathbf{r}) = x^a y^b z^c \sum_k c_{ki} \left( e^{-\alpha_{ki} \mathbf{r}^2 - \imath \mathbf{k}_{ki} \cdot \mathbf{r} - \imath \phi_{ki}} + \text{C.C.} \right)
 
 where:
-   - :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im_cgtos``),
+   - :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im``),
    - :math:`\mathbf{k} = (k_x, k_y, k_z) \in \mathbb{R}^3` (specified by ``ao_expo_pw``),
    - :math:`\phi = \phi_x + \phi_y + \phi_z \in \mathbb{R}` (specified by ``ao_expo_phase``).
diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f
@@ -75,10 +75,6 @@
       enddo
     endif
 
-    powA(1) = ao_power(i,1)
-    powA(2) = ao_power(i,2)
-    powA(3) = ao_power(i,3)
-
     ! Normalization of the contracted basis functions
     if (ao_normalized) then
       norm = 0.d0

diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f
@@ -4,6 +4,7 @@
 BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (ao_prim_num_max, ao_num)]
 
   implicit none
+
   integer :: i, j
 
   do j = 1, ao_num
@@ -17,25 +18,22 @@
 ! ---
 
  BEGIN_PROVIDER [complex*16, ao_expo_cgtos_ord_transp, (ao_prim_num_max, ao_num)]
-&BEGIN_PROVIDER [complex*16, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)]
-&BEGIN_PROVIDER [complex*16, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)]
+&BEGIN_PROVIDER [double precision, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)]
+&BEGIN_PROVIDER [double precision, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)]
 
   implicit none
+
   integer :: i, j, m
 
   do j = 1, ao_num
     do i = 1, ao_prim_num_max
+
       ao_expo_cgtos_ord_transp(i,j) = ao_expo_cgtos_ord(j,i)
-      do m = 1, 3
+
+      do m = 1, 4
         ao_expo_pw_ord_transp(m,i,j) = ao_expo_pw_ord(m,j,i)
         ao_expo_phase_ord_transp(m,i,j) = ao_expo_phase_ord(m,j,i)
       enddo
-      ao_expo_pw_ord_transp(4,i,j) = ao_expo_pw_ord_transp(1,i,j) * ao_expo_pw_ord_transp(1,i,j) &
-                                   + ao_expo_pw_ord_transp(2,i,j) * ao_expo_pw_ord_transp(2,i,j) &
-                                   + ao_expo_pw_ord_transp(3,i,j) * ao_expo_pw_ord_transp(3,i,j)
-      ao_expo_phase_ord_transp(4,i,j) = ao_expo_phase_ord_transp(1,j,i) &
-                                      + ao_expo_phase_ord_transp(2,j,i) &
-                                      + ao_expo_phase_ord_transp(3,j,i)
     enddo
   enddo
 
@@ -50,16 +48,12 @@
   integer          :: i, j, ii, m, powA(3), nz
   double precision :: norm
   double precision :: kA2, phiA
-  complex*16       :: expo, expo_inv, C_A(3)
+  complex*16       :: expo, expo_inv, C_Ae(3), C_Ap(3)
   complex*16       :: overlap_x, overlap_y, overlap_z
   complex*16       :: integ1, integ2, C1, C2
 
   nz = 100
 
-  C_A(1) = (0.d0, 0.d0)
-  C_A(2) = (0.d0, 0.d0)
-  C_A(3) = (0.d0, 0.d0)
-
   ao_coef_norm_cgtos = 0.d0
 
   do i = 1, ao_num
@@ -69,27 +63,34 @@
     powA(2) = ao_power(i,2)
     powA(3) = ao_power(i,3)
 
-    ! Normalization of the primitives
     if(primitives_normalized) then
 
+      ! Normalization of the primitives
       do j = 1, ao_prim_num(i)
 
-        expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im_cgtos(i,j)
+        expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im(i,j)
         expo_inv = (1.d0, 0.d0) / expo
         do m = 1, 3
-          C_A(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j)
+          C_Ap(m) = nucl_coord(ii,m)
+          C_Ae(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j)
         enddo
-        phiA = ao_expo_phase(4,i,j)
-        KA2 = ao_expo_pw(4,i,j)
+        phiA = ao_expo_phase(1,i,j) + ao_expo_phase(2,i,j) + ao_expo_phase(3,i,j)
+        KA2 = ao_expo_pw(1,i,j) * ao_expo_pw(1,i,j) &
+            + ao_expo_pw(2,i,j) * ao_expo_pw(2,i,j) &
+            + ao_expo_pw(3,i,j) * ao_expo_pw(3,i,j)
 
         C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2)
         C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2)
 
-        call overlap_cgaussian_xyz(C_A,        C_A,        expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz)
-        call overlap_cgaussian_xyz(conjg(C_A), C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz)
+        call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, &
+                                   C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz)
+
+        call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, &
+                                   conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz)
 
         norm = 2.d0 * real(C1 * integ1 + C2 * integ2)
 
+        !ao_coef_norm_cgtos(i,j) = 1.d0 / dsqrt(norm)
         ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
       enddo
 
@@ -99,7 +100,7 @@
         ao_coef_norm_cgtos(i,j) = ao_coef(i,j)
       enddo
 
-    endif
+    endif ! primitives_normalized
 
   enddo
 
@@ -126,12 +127,17 @@
       iorder(j) = j
       d(j,1) = ao_expo(i,j)
       d(j,2) = ao_coef_norm_cgtos(i,j)
-      d(j,3) = ao_expo_im_cgtos(i,j)
+      d(j,3) = ao_expo_im(i,j)
 
-      do m = 1, 4
+      do m = 1, 3
         d(j,3+m) = ao_expo_pw(m,i,j)
+      enddo
+      d(j,7) = d(j,4) * d(j,4) + d(j,5) * d(j,5) + d(j,6) * d(j,6)
+
+      do m = 1, 3
         d(j,7+m) = ao_expo_phase(m,i,j)
       enddo
+      d(j,11) = d(j,8) + d(j,9) + d(j,10)
     enddo
 
     call dsort(d(1,1), iorder, ao_prim_num(i))
@@ -165,8 +171,8 @@
   double precision :: c, overlap, overlap_x, overlap_y, overlap_z
   double precision :: KA2(3), phiA(3)
   double precision :: KB2(3), phiB(3)
-  complex*16       :: alpha, alpha_inv, A_center(3)
-  complex*16       :: beta, beta_inv, B_center(3)
+  complex*16       :: alpha, alpha_inv, Ae_center(3), Ap_center(3)
+  complex*16       :: beta, beta_inv, Be_center(3), Bp_center(3)
   complex*16       :: C1(1:4), C2(1:4)
   complex*16       :: overlap1, overlap_x1, overlap_y1, overlap_z1
   complex*16       :: overlap2, overlap_x2, overlap_y2, overlap_z2
@@ -178,18 +184,18 @@
 
   dim1 = 100
 
- !$OMP PARALLEL DO SCHEDULE(GUIDED)                                       &
- !$OMP DEFAULT(NONE)                                                      &
- !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2,                          &
- !$OMP         alpha, alpha_inv, A_center, power_A, KA2, phiA,            &
- !$OMP         beta, beta_inv, B_center, power_B, KB2, phiB,              &
- !$OMP         overlap_x , overlap_y , overlap_z , overlap,               &
- !$OMP         overlap_x1, overlap_y1, overlap_z1, overlap1,              &
- !$OMP         overlap_x2, overlap_y2, overlap_z2, overlap2)              &
- !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1,   &
- !$OMP        ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp,    &
- !$OMP        ao_expo_pw_ord_transp, ao_expo_phase_ord_transp,            &
- !$OMP        ao_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z, &
+ !$OMP PARALLEL DO SCHEDULE(GUIDED)                                        &
+ !$OMP DEFAULT(NONE)                                                       &
+ !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2,                           &
+ !$OMP         alpha, alpha_inv, Ae_center, Ap_center, power_A, KA2, phiA, &
+ !$OMP         beta, beta_inv, Be_center, Bp_center, power_B, KB2, phiB,   &
+ !$OMP         overlap_x , overlap_y , overlap_z , overlap,                &
+ !$OMP         overlap_x1, overlap_y1, overlap_z1, overlap1,               &
+ !$OMP         overlap_x2, overlap_y2, overlap_z2, overlap2)               &
+ !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1,    &
+ !$OMP        ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp,     &
+ !$OMP        ao_expo_pw_ord_transp, ao_expo_phase_ord_transp,             &
+ !$OMP        ao_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z,  &
  !$OMP        ao_overlap_cgtos)
 
   do j = 1, ao_num
@@ -212,7 +218,8 @@
         alpha_inv = (1.d0, 0.d0) / alpha
         do m = 1, 3
           phiA(m) = ao_expo_phase_ord_transp(m,n,j)
-          A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j)
+          Ap_center(m) = nucl_coord(jj,m)
+          Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j)
           KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j)
         enddo
 
@@ -222,7 +229,8 @@
           beta_inv = (1.d0, 0.d0) / beta
           do m = 1, 3
             phiB(m) = ao_expo_phase_ord_transp(m,l,i)
-            B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i)
+            Bp_center(m) = nucl_coord(ii,m)
+            Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i)
             KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i)
           enddo
 
@@ -238,11 +246,11 @@
           C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3)))
           C2(4) = C2(1) * C2(2) * C2(3)
 
-          call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
-                                     overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
+          call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
+                                     Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
 
-          call overlap_cgaussian_xyz(conjg(A_center), B_center, conjg(alpha), beta, power_A, power_B, &
-                                     overlap_x2, overlap_y2, overlap_z2, overlap2, dim1)
+          call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, &
+                                     conjg(Ap_center), Bp_center, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1)
 
           overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2)
           overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2)