working on X-Fock matrix (with s8 sym)

src/AOtoMO/Hartree_matrix_AO_basis.f90

@@ -49,13 +49,11 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
   integer                       :: nunu, lala, nula, lasi, numu
   integer*8                     :: nunununu, nunulala, nununula, nunulasi
   integer*8                     :: lalanunu, lasinunu, numulala, lalanumu
-  integer*8                     :: numunula, numulasi, lasinumu
+  integer*8                     :: numunula, numulasi, lasinumu, nununumu
   integer*8                     :: munu0, munu
   integer*8                     :: sila0, sila
   integer*8                     :: munulasi0, munulasi
 
-  integer*8, external           :: Yoshimine_4ind
-
 
   do nu = 1, nBas
 
@@ -80,7 +78,7 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
       enddo
     enddo
 
-    do la = nu+1, nBas
+    do la = nu + 1, nBas
 
       lala = (la * (la - 1)) / 2 + la
       lalanunu = (lala * (lala - 1)) / 2 + nunu
@@ -96,15 +94,16 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
     do mu = 1, nu - 1
 
       numu = (nu * (nu - 1)) / 2 + mu
-
+      nununumu = (nunu * (nunu - 1)) / 2 + numu
-      H(mu,nu) = 0.d0
+      H(mu,nu) = p(nu,nu) * ERI_chem(nununumu)
 
       do la = 1, nu - 1
         lala = (la * (la - 1)) / 2 + la
         numulala = (numu * (numu - 1)) / 2 + lala
         H(mu,nu) = H(mu,nu) + p(la,la) * ERI_chem(numulala)
       enddo
-      do la = nu, nBas
+
+      do la = nu + 1, nBas
         lala = (la * (la - 1)) / 2 + la
         lalanumu = (lala * (lala - 1)) / 2 + numu
         H(mu,nu) = H(mu,nu) + p(la,la) * ERI_chem(lalanumu)
@@ -115,11 +114,13 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
         numunula = (numu * (numu - 1)) / 2 + nula
         H(mu,nu) = H(mu,nu) + 2.d0 * P(la,nu) * ERI_chem(numunula)
       enddo
+
       do la = mu + 1, nu - 1
         nula = (nu * (nu - 1)) / 2 + la
         numunula = (nula * (nula - 1)) / 2 + numu
         H(mu,nu) = H(mu,nu) + 2.d0 * P(la,nu) * ERI_chem(numunula)
       enddo
+
       do la = 2, nu - 1
         do si = 1, la - 1
           lasi = (la * (la - 1)) / 2 + si
@@ -127,6 +128,7 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
           H(mu,nu) = H(mu,nu) + 2.d0 * P(si,la) * ERI_chem(numulasi)
         enddo
       enddo
+
       do la = nu + 1, nBas
         do si = 1, la - 1
           lasi = (la * (la - 1)) / 2 + si

src/AOtoMO/exchange_matrix_AO_basis.f90

@@ -31,3 +31,152 @@ subroutine exchange_matrix_AO_basis(nBas,P,ERI,K)
   end do
 
 end subroutine 
+
+! ---
+
+subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
+
+  implicit none
+
+  integer,          intent(in)  :: nBas
+  integer*8,        intent(in)  :: ERI_size
+  double precision, intent(in)  :: P(nBas,nBas)
+  double precision, intent(in)  :: ERI_chem(ERI_size)
+  double precision, intent(out) :: K(nBas,nBas)
+
+  integer                       :: mu, nu, la, si
+  integer                       :: nunu, lala, nula, lasi, numu
+  integer*8                     :: nunununu, nunulala, nununula, nunulasi
+  integer*8                     :: lalanunu, lasinunu, numulala, lalanumu
+  integer*8                     :: numunula, numulasi, lasinumu, nununumu
+  integer*8                     :: munu0, munu
+  integer*8                     :: sila0, sila
+  integer*8                     :: munulasi0, munulasi
+
+  integer*8, external           :: Yoshimine_4ind
+
+
+  do nu = 1, nBas
+
+    munulasi = Yoshimine_4ind(nu, nu, nu, nu)
+    K(nu,nu) = -P(nu,nu) * ERI_chem(munulasi)
+
+    do la = 1, nu - 1
+      munulasi = Yoshimine_4ind(nu, la, nu, la)
+      K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(munulasi)
+    enddo
+
+    do la = nu + 1, nBas
+      munulasi = Yoshimine_4ind(nu, la, nu, la)
+      K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(munulasi)
+    enddo
+
+    do la = 1, nu
+      do si = 1, la - 1
+        munulasi = Yoshimine_4ind(nu, la, nu, si)
+        K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(munulasi)
+      enddo
+    enddo
+
+    do la = nu + 1, nBas
+      do si = 1, nu
+        munulasi = Yoshimine_4ind(nu, la, nu, si)
+        K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(munulasi)
+      enddo
+    enddo
+
+    do la = nu + 1, nBas
+      do si = nu + 1, la - 1
+        munulasi = Yoshimine_4ind(nu, la, nu, si)
+        K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(munulasi)
+      enddo
+    enddo
+
+
+    do mu = 1, nu - 1
+
+      K(mu,nu) = 0.d0
+      do la = 1, mu
+        munulasi = Yoshimine_4ind(nu, la, mu, la)
+        K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(munulasi)
+      enddo
+      do la = mu + 1, nu
+        munulasi = Yoshimine_4ind(nu, la, mu, la)
+        K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(munulasi)
+      enddo
+      do la = nu + 1, nBas
+        munulasi = Yoshimine_4ind(nu, la, mu, la)
+        K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(munulasi)
+      enddo
+
+      do la = 1, mu
+        do si = 1, la - 1
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+      enddo
+      do la = mu+1, nu
+        do si = 1, mu
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+        do si = mu + 1, la - 1
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+      enddo
+      do la = nu + 1, nBas
+        do si = 1, mu
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+        do si = mu + 1, la - 1
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+      enddo
+
+      do la = 1, mu
+        do si = la + 1, mu
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+        do si = mu + 1, nBas
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+      enddo
+      do la = mu + 1, nu
+        do si = la + 1, nBas
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+      enddo
+      do la = nu + 1, nBas
+        do si = la + 1, nBas
+          munulasi = Yoshimine_4ind(nu, la, si, mu)
+          K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+        enddo
+      enddo
+
+      !do la = 1, nBas
+      !  do si = la + 1, nBas
+      !    munulasi = Yoshimine_4ind(nu, la, mu, si)
+      !    K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munulasi)
+      !  enddo
+      !enddo
+    enddo ! mu
+  enddo ! nu
+
+
+  do nu = 1, nBas
+    do mu = nu+1, nBas
+      K(mu,nu) = K(nu,mu)
+    enddo
+  enddo
+
+  return
+end subroutine 
+
+! ---
+

src/HF/RHF_hpc.f90

@@ -56,7 +56,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
   double precision,allocatable  :: cp(:,:)
   double precision,allocatable  :: Fp(:,:)
   double precision,allocatable  :: ERI_chem(:)
-  double precision,allocatable  :: ERI_phys(:,:,:,:), J_deb(:,:)
+  double precision,allocatable  :: ERI_phys(:,:,:,:), J_deb(:,:), K_deb(:,:)
 
 
 ! Output variables
@@ -106,25 +106,43 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
   call read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem)
   call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J)
   
-  allocate(J_deb(nBas,nBas))
   allocate(ERI_phys(nBas,nBas,nBas,nBas))
+  allocate(J_deb(nBas,nBas))
+  allocate(K_deb(nBas,nBas))
+
   call read_2e_integrals(working_dir, nBas, ERI_phys)
   call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
-
+  print*, "max error on J = ", maxval(dabs(J - J_deb))
-  print*, "max error = ", maxval(dabs(J - J_deb))
   diff = 0.d0
   do ii = 1, nBas
     do jj = 1, nBas
       diff_loc = dabs(J(jj,ii) - J_deb(jj,ii))
       if(diff_loc .gt. 1d-12) then
-        print*, 'error on: ', jj, ii
+        print*, 'error in J on: ', jj, ii
         print*, J(jj,ii), J_deb(jj,ii)
         stop
       endif
       diff = diff + diff_loc
     enddo
   enddo
-  print*, 'total diff = ', diff
+  print*, 'total diff on J = ', diff
+
+  call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K)
+  call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb)
+  print*, "max error on K = ", maxval(dabs(K - K_deb))
+  diff = 0.d0
+  do ii = 1, nBas
+    do jj = 1, nBas
+      diff_loc = dabs(K(jj,ii) - K_deb(jj,ii))
+      if(diff_loc .gt. 1d-12) then
+        print*, 'error in K on: ', jj, ii
+        print*, K(jj,ii), K_deb(jj,ii)
+        stop
+      endif
+      diff = diff + diff_loc
+    enddo
+  enddo
+  print*, 'total diff on K = ', diff
 
   stop
 

src/utils/utils.f90

@@ -914,9 +914,11 @@ integer*8 function Yoshimine_4ind(a, b, c, d)
   implicit none
   integer, intent(in) :: a, b, c, d
   integer*8, external :: Yoshimine_2ind
+  integer*8           :: ab, cd
 
-  Yoshimine_4ind = Yoshimine_2ind(Yoshimine_2ind(a, b), &
+  ab = Yoshimine_2ind(a, b)
-                                  Yoshimine_2ind(c, d))
+  cd = Yoshimine_2ind(c, d)
+  Yoshimine_4ind = Yoshimine_2ind(ab, cd)
 
   return
 end