full Coulomb-Fock matrix with 8-fold symmetry

src/AOtoMO/Hartree_matrix_AO_basis.f90

@@ -46,102 +46,98 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
   double precision, intent(out) :: H(nBas,nBas)
 
   integer                       :: mu, nu, la, si
-  integer                       :: nunu, lala, nula, lasi
+  integer                       :: nunu, lala, nula, lasi, numu
-  integer*8                     :: nunununu, nunulala, nununula, nunulasi, lalanunu, lasinunu
+  integer*8                     :: nunununu, nunulala, nununula, nunulasi
+  integer*8                     :: lalanunu, lasinunu, numulala, lalanumu
+  integer*8                     :: numunula, numulasi, lasinumu
   integer*8                     :: munu0, munu
   integer*8                     :: sila0, sila
   integer*8                     :: munulasi0, munulasi
 
   integer*8, external           :: Yoshimine_4ind
 
-!  do nu = 1, nBas
-!    do mu = 1, nBas
-!      H(mu,nu) = 0.d0
-!      do si = 1, nBas
-!        do la = 1, nBas
-!          munulasi = Yoshimine_4ind(mu, nu, la, si)
-!          H(mu,nu) = H(mu,nu) + P(la,si) * ERI_chem(munulasi)
-!        enddo
-!      enddo
-!    enddo
-!  enddo
-
-!  do nu = 1, nBas
-!    do mu = 1, nu
-!      H(mu,nu) = 0.d0
-!      do si = 1, nBas
-!        munulasi = Yoshimine_4ind(mu, nu, si, si)
-!        H(mu,nu) = H(mu,nu) + P(si,si) * ERI_chem(munulasi)
-!        do la = 1, si-1
-!          munulasi = Yoshimine_4ind(mu, nu, la, si)
-!          H(mu,nu) = H(mu,nu) + 2.d0 * P(la,si) * ERI_chem(munulasi)
-!        enddo
-!      enddo
-!    enddo
-!  enddo
-
-
 
   do nu = 1, nBas
 
     nunu = (nu * (nu - 1)) / 2 + nu
     nunununu = (nunu * (nunu - 1)) / 2 + nunu
-    !nunununu = Yoshimine_4ind(nu, nu, nu, nu)
     H(nu,nu) = P(nu,nu) * ERI_chem(nunununu)
 
     do la = 1, nu-1
 
-      ! la < nu
       lala = (la * (la - 1)) / 2 + la
       nunulala = (nunu * (nunu - 1)) / 2 + lala
-      !nunulala = Yoshimine_4ind(nu, nu, la, la)
       H(nu,nu) = H(nu,nu) + P(la,la) * ERI_chem(nunulala)
 
-      ! la < nu
       nula = (nu * (nu - 1)) / 2 + la
       nununula = (nunu * (nunu - 1)) / 2 + nula
-      !nununula = Yoshimine_4ind(nu, nu, la, nu)
       H(nu,nu) = H(nu,nu) + 2.d0 * P(la,nu) * ERI_chem(nununula)
 
       do si = 1, la - 1
-        ! lasi < nunu
         lasi = (la * (la - 1)) / 2 + si
         nunulasi = (nunu * (nunu - 1)) / 2 + lasi
-        !nunulasi = Yoshimine_4ind(nu, nu, si, la)
         H(nu,nu) = H(nu,nu) + 2.d0 * P(si,la) * ERI_chem(nunulasi)
       enddo
     enddo
 
     do la = nu+1, nBas
 
-      ! nu < la
       lala = (la * (la - 1)) / 2 + la
       lalanunu = (lala * (lala - 1)) / 2 + nunu
-      !lalanunu = Yoshimine_4ind(nu, nu, la, la)
       H(nu,nu) = H(nu,nu) + P(la,la) * ERI_chem(lalanunu)
 
       do si = 1, la - 1
-        ! nunu < lasi
         lasi = (la * (la - 1)) / 2 + si
         lasinunu = (lasi * (lasi - 1)) / 2 + nunu
-        !lasinunu = Yoshimine_4ind(nu, nu, si, la)
         H(nu,nu) = H(nu,nu) + 2.d0 * P(si,la) * ERI_chem(lasinunu)
       enddo
     enddo
 
-    do mu = 1, nu-1
+    do mu = 1, nu - 1
+
+      numu = (nu * (nu - 1)) / 2 + mu
+
       H(mu,nu) = 0.d0
-      do si = 1, nBas
+
-        munulasi = Yoshimine_4ind(mu, nu, si, si)
+      do la = 1, nu - 1
-        H(mu,nu) = H(mu,nu) + P(si,si) * ERI_chem(munulasi)
+        lala = (la * (la - 1)) / 2 + la
-        do la = 1, si-1
+        numulala = (numu * (numu - 1)) / 2 + lala
-          munulasi = Yoshimine_4ind(mu, nu, la, si)
+        H(mu,nu) = H(mu,nu) + p(la,la) * ERI_chem(numulala)
-          H(mu,nu) = H(mu,nu) + 2.d0 * P(la,si) * ERI_chem(munulasi)
+      enddo
+      do la = nu, 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)
+      enddo
+
+      do la = 1, mu
+        nula = (nu * (nu - 1)) / 2 + la
+        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
+          numulasi = (numu * (numu - 1)) / 2 + lasi
+          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
+          lasinumu = (lasi * (lasi - 1)) / 2 + numu
+          H(mu,nu) = H(mu,nu) + 2.d0 * P(si,la) * ERI_chem(lasinumu)
         enddo
       enddo
 
+    enddo ! mu
+  enddo ! nu
+
 
   do nu = 1, nBas
     do mu = nu+1, nBas

src/HF/RHF_hpc.f90

@@ -111,12 +111,12 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh
   call read_2e_integrals(working_dir, nBas, ERI_phys)
   call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
 
-  print*, 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-13) then
+      if(diff_loc .gt. 1d-12) then
         print*, 'error on: ', jj, ii
         print*, J(jj,ii), J_deb(jj,ii)
         stop