From 995bcfd2362eff3ecff11470d8450d53dbbfd678 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 10 Dec 2024 11:29:31 +0100 Subject: [PATCH] working on X-Fock matrix (with s8 sym) --- src/AOtoMO/Hartree_matrix_AO_basis.f90 | 16 +-- src/AOtoMO/exchange_matrix_AO_basis.f90 | 149 ++++++++++++++++++++++++ src/HF/RHF_hpc.f90 | 30 ++++- src/utils/utils.f90 | 6 +- 4 files changed, 186 insertions(+), 15 deletions(-) diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 83a3e50..0d0810d 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/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 - - H(mu,nu) = 0.d0 + nununumu = (nunu * (nunu - 1)) / 2 + numu + 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 diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 1a29038..d1c96b4 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/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 + +! --- + diff --git a/src/HF/RHF_hpc.f90 b/src/HF/RHF_hpc.f90 index dc66299..d9bf77d 100644 --- a/src/HF/RHF_hpc.f90 +++ b/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 = ", maxval(dabs(J - J_deb)) + print*, "max error on J = ", 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 diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index 79d6d82..ad528c9 100644 --- a/src/utils/utils.f90 +++ b/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), & - Yoshimine_2ind(c, d)) + ab = Yoshimine_2ind(a, b) + cd = Yoshimine_2ind(c, d) + Yoshimine_4ind = Yoshimine_2ind(ab, cd) return end