From 37420d12078abdf512530d330a1d70b094c6b15f Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 10 Dec 2024 19:14:33 +0100 Subject: [PATCH] added Hartree and EX terms with s8-sym --- PyDuck.py | 1 - src/AOtoMO/exchange_matrix_AO_basis.f90 | 165 ++++++++++++------------ src/HF/RHF_hpc.f90 | 35 ++++- src/utils/utils.f90 | 6 +- 4 files changed, 115 insertions(+), 92 deletions(-) diff --git a/PyDuck.py b/PyDuck.py index 45113de..4779566 100644 --- a/PyDuck.py +++ b/PyDuck.py @@ -163,7 +163,6 @@ if print_2e: output_file_path = working_dir + '/int/ERI_chem.bin' subprocess.call(['rm', '-f', output_file_path]) eri_ao = mol.intor('int2e', aosym='s8') - print(eri_ao.shape) f = open(output_file_path, 'w') eri_ao.tofile(output_file_path) f.close() diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 47ffbe6..4781436 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/src/AOtoMO/exchange_matrix_AO_basis.f90 @@ -50,49 +50,48 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) integer*8 :: nunununu, nulanula, lanulanu, nulanusi integer*8 :: munulasi, lanunusi, lanusinu, numumumu integer*8 :: nulamula, nulalamu, lanulamu, nulamusi - integer*8 :: nulasimu, lanumusi, lanusimu + integer*8 :: nulasimu, lanumusi, lanusimu, simunula + integer*8 :: simulanu - integer*8, external :: Yoshimine_4ind - do nu = 1, nBas - nunu = (nu * (nu - 1)) / 2 + nu - nunununu = (nunu * (nunu - 1)) / 2 + nunu + nunu = shiftr(nu * (nu - 1), 1) + nu + nunununu = shiftr(nunu * (nunu - 1), 1) + nunu K(nu,nu) = -P(nu,nu) * ERI_chem(nunununu) do la = 1, nu - 1 - nula = (nu * (nu - 1)) / 2 + la - nulanula = (nula * (nula - 1)) / 2 + nula + nula = shiftr(nu * (nu - 1), 1) + la + nulanula = shiftr(nula * (nula - 1), 1) + nula K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(nulanula) enddo do la = nu + 1, nBas - lanu = (la * (la - 1)) / 2 + nu - lanulanu = (lanu * (lanu - 1)) / 2 + lanu + lanu = shiftr(la * (la - 1), 1) + nu + lanulanu = shiftr(lanu * (lanu - 1), 1) + lanu K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(lanulanu) enddo do la = 1, nu - nula = (nu * (nu - 1)) / 2 + la + nula = shiftr(nu * (nu - 1), 1) + la do si = 1, la - 1 - nusi = (nu * (nu - 1)) / 2 + si - nulanusi = (nula * (nula - 1)) / 2 + nusi + nusi = shiftr(nu * (nu - 1), 1) + si + nulanusi = shiftr(nula * (nula - 1), 1) + nusi K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(nulanusi) enddo enddo do la = nu + 1, nBas - lanu = (la * (la - 1)) / 2 + nu + lanu = shiftr(la * (la - 1), 1) + nu do si = 1, nu - nusi = (nu * (nu - 1)) / 2 + si - lanunusi = (lanu * (lanu - 1)) / 2 + nusi + nusi = shiftr(nu * (nu - 1), 1) + si + lanunusi = shiftr(lanu * (lanu - 1), 1) + nusi K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(lanunusi) enddo do si = nu + 1, la - 1 - sinu = (si * (si - 1)) / 2 + nu - lanusinu = (lanu * (lanu - 1)) / 2 + sinu + sinu = shiftr(si * (si - 1), 1) + nu + lanusinu = shiftr(lanu * (lanu - 1), 1) + sinu K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(lanusinu) enddo enddo @@ -100,108 +99,104 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) do mu = 1, nu - 1 - numu = (nu * (nu - 1)) / 2 + mu - mumu = (mu * (mu - 1)) / 2 + mu - numumumu = (numu * (numu - 1)) / 2 + mumu + numu = shiftr(nu * (nu - 1), 1) + mu + mumu = shiftr(mu * (mu - 1), 1) + mu + numumumu = shiftr(numu * (numu - 1), 1) + mumu K(mu,nu) = - P(mu,mu) * ERI_chem(numumumu) do la = 1, mu - 1 - mula = (mu * (mu - 1)) / 2 + la - nula = (nu * (nu - 1)) / 2 + la - nulamula = (nula * (nula - 1)) / 2 + mula + mula = shiftr(mu * (mu - 1), 1) + la + nula = shiftr(nu * (nu - 1), 1) + la + nulamula = shiftr(nula * (nula - 1), 1) + mula K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(nulamula) enddo do la = mu + 1, nu - lamu = (la * (la - 1)) / 2 + mu - nula = (nu * (nu - 1)) / 2 + la - nulalamu = (nula * (nula - 1)) / 2 + lamu + lamu = shiftr(la * (la - 1), 1) + mu + nula = shiftr(nu * (nu - 1), 1) + la + nulalamu = shiftr(nula * (nula - 1), 1) + lamu K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(nulalamu) enddo do la = nu + 1, nBas - lamu = (la * (la - 1)) / 2 + mu - lanu = (la * (la - 1)) / 2 + nu - lanulamu = (lanu * (lanu - 1)) / 2 + lamu + lamu = shiftr(la * (la - 1), 1) + mu + lanu = shiftr(la * (la - 1), 1) + nu + lanulamu = shiftr(lanu * (lanu - 1), 1) + lamu K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(lanulamu) enddo do la = 1, mu - nula = (nu * (nu - 1)) / 2 + la + nula = shiftr(nu * (nu - 1), 1) + la do si = 1, la - 1 - musi = (mu * (mu - 1)) / 2 + si - nulamusi = (nula * (nula - 1)) / 2 + musi + musi = shiftr(mu * (mu - 1), 1) + si + nulamusi = shiftr(nula * (nula - 1), 1) + musi K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi) enddo enddo do la = mu + 1, nu - nula = (nu * (nu - 1)) / 2 + la + nula = shiftr(nu * (nu - 1), 1) + la do si = 1, mu - musi = (mu * (mu - 1)) / 2 + si - nulamusi = (nula * (nula - 1)) / 2 + musi + musi = shiftr(mu * (mu - 1), 1) + si + nulamusi = shiftr(nula * (nula - 1), 1) + musi K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi) enddo do si = mu + 1, la - 1 - simu = (si * (si - 1)) / 2 + mu - nulasimu = (nula * (nula - 1)) / 2 + simu + simu = shiftr(si * (si - 1), 1) + mu + nulasimu = shiftr(nula * (nula - 1), 1) + simu K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) enddo enddo do la = nu + 1, nBas - lanu = (la * (la - 1)) / 2 + nu + lanu = shiftr(la * (la - 1), 1) + nu do si = 1, mu - musi = (mu * (mu - 1)) / 2 + si - lanumusi = (lanu * (lanu - 1)) / 2 + musi + musi = shiftr(mu * (mu - 1), 1) + si + lanumusi = shiftr(lanu * (lanu - 1), 1) + musi K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(lanumusi) enddo do si = mu + 1, la - 1 - simu = (si * (si - 1)) / 2 + mu - lanusimu = (lanu * (lanu - 1)) / 2 + simu + simu = shiftr(si * (si - 1), 1) + mu + lanusimu = shiftr(lanu * (lanu - 1), 1) + simu K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(lanusimu) enddo enddo -!TODO -! do la = 1, mu -! nula = (nu * (nu - 1)) / 2 + la -! do si = la + 1, mu -! musi = (mu * (mu - 1)) / 2 + si -! nulamusi = (nula * (nula - 1)) / 2 + musi -! !nulamusi = Yoshimine_4ind(nu, la, si, mu) -! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi) -! enddo -! do si = mu + 1, nBas -! simu = (si * (si - 1)) / 2 + mu -! nulasimu = (nula * (nula - 1)) / 2 + simu -! !nulasimu = Yoshimine_4ind(nu, la, si, mu) -! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) -! enddo -! enddo -! do la = mu + 1, nu -! nula = (nu * (nu - 1)) / 2 + la -! do si = la + 1, nu -! simu = (si * (si - 1)) / 2 + mu -! nulasimu = (nula * (nula - 1)) / 2 + simu -! !nulasimu = Yoshimine_4ind(nu, la, si, mu) -! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) -! enddo -! do si = nu + 1, nBas -! simu = (si * (si - 1)) / 2 + mu -! 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 -! lanu = (la * (la - 1)) / 2 + nu -! do si = la + 1, mu -! simu = (si * (si - 1)) / 2 + 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 -! musi = (mu * (mu - 1)) / 2 + si -! 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 + nula = shiftr(nu * (nu - 1), 1) + la + do si = la + 1, mu + musi = shiftr(mu * (mu - 1) , 1) + si + nulamusi = shiftr(nula * (nula - 1), 1) + musi + K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi) + enddo + do si = mu + 1, nu - 1 + simu = shiftr(si * (si - 1), 1) + mu + nulasimu = shiftr(nula * (nula - 1), 1) + simu + K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) + enddo + do si = nu, nBas + simu = shiftr(si * (si - 1), 1) + mu + simunula = shiftr(simu * (simu - 1), 1) + nula + K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(simunula) + enddo + enddo + do la = mu + 1, nu + nula = shiftr(nu * (nu - 1), 1) + la + do si = la + 1, nu + simu = shiftr(si * (si - 1), 1) + mu + nulasimu = shiftr(nula * (nula - 1), 1) + simu + K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) + enddo + do si = nu + 1, nBas + simu = shiftr(si * (si - 1), 1) + mu + simunula = shiftr(simu * (simu - 1), 1) + nula + K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(simunula) + enddo + enddo + do la = nu + 1, nBas + lanu = shiftr(la * (la - 1), 1) + nu + do si = la + 1, nBas + simu = shiftr(si * (si - 1), 1) + mu + simulanu = shiftr(simu * (simu - 1), 1) + lanu + K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(simulanu) + enddo + enddo enddo ! mu enddo ! nu diff --git a/src/HF/RHF_hpc.f90 b/src/HF/RHF_hpc.f90 index d9bf77d..44095d3 100644 --- a/src/HF/RHF_hpc.f90 +++ b/src/HF/RHF_hpc.f90 @@ -39,6 +39,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh integer :: nBas_Sq integer :: n_diis integer*8 :: ERI_size + double precision :: t1, t2 double precision :: diff, diff_loc double precision :: ET double precision :: EV @@ -104,14 +105,42 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh ERI_size = (ERI_size * (ERI_size + 1)) / 2 allocate(ERI_chem(ERI_size)) call read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem) - call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J) + + call wall_time(t1) + do ii = 1, 5 + call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J) + enddo + call wall_time(t2) + print*, " J built in (sec):", (t2-t1) / 5.d0 + + call wall_time(t1) + do ii = 1, 5 + call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) + enddo + call wall_time(t2) + print*, " K built in (sec):", (t2-t1) / 5.d0 + 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) + + call wall_time(t1) + do ii = 1, 5 + call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb) + enddo + call wall_time(t2) + print*, " J_deb built in (sec):", (t2-t1) / 5.d0 + + call wall_time(t1) + do ii = 1, 5 + call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb) + enddo + call wall_time(t2) + print*, " K_deb built in (sec):", (t2-t1) / 5.d0 + print*, "max error on J = ", maxval(dabs(J - J_deb)) diff = 0.d0 do ii = 1, nBas @@ -127,8 +156,6 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh enddo 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 diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index c6eaac0..cc7998a 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -929,9 +929,11 @@ integer*8 function Yoshimine_2ind(a, b) integer*8, intent(in) :: a, b if(a > b) then - Yoshimine_2ind = (a * (a - 1)) / 2 + b + !Yoshimine_2ind = (a * (a - 1)) / 2 + b + Yoshimine_2ind = shiftr(a * (a - 1), 1) + b else - Yoshimine_2ind = (b * (b - 1)) / 2 + a + !Yoshimine_2ind = (b * (b - 1)) / 2 + a + Yoshimine_2ind = shiftr(b * (b - 1), 1) + a endif return