From ec797982deed1929e9b8246e9670b705ba7cef6e Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 10 Dec 2024 02:25:52 +0100 Subject: [PATCH] full Coulomb-Fock matrix with 8-fold symmetry --- src/AOtoMO/Hartree_matrix_AO_basis.f90 | 94 ++++++++++++-------------- src/HF/RHF_hpc.f90 | 4 +- 2 files changed, 47 insertions(+), 51 deletions(-) diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 11e4788..83a3e50 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -46,101 +46,97 @@ 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*8 :: nunununu, nunulala, nununula, nunulasi, lalanunu, lasinunu + 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 :: 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) - 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) + + 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 + 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 - 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 diff --git a/src/HF/RHF_hpc.f90 b/src/HF/RHF_hpc.f90 index ef6ea47..dc66299 100644 --- a/src/HF/RHF_hpc.f90 +++ b/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