From d28c0339ff9190164631d580825d65ac2bb164e7 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 12 Dec 2024 18:13:47 +0100 Subject: [PATCH] // H calc --- src/AOtoMO/Hartree_matrix_AO_basis.f90 | 31 ++++++- src/AOtoMO/exchange_matrix_AO_basis.f90 | 103 ++++++++++++++---------- src/HF/RHF_hpc.f90 | 28 +++---- 3 files changed, 101 insertions(+), 61 deletions(-) diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 4d12a5b..010e5bb 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -53,7 +53,32 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H) integer*8 :: numunula, numulasi, lasinumu, nununumu integer*8 :: nunununu0, numunumu0 +! integer*8 :: munusila +! integer*8, external :: Yoshimine_4ind +! +! do nu = 1, nBas +! do mu = 1, nu +! H(mu,nu) = 0.d0 +! do la = 1, nBas +! do si = 1, nBas +! munusila = Yoshimine_4ind(int(mu, kind=8), & +! int(nu, kind=8), & +! int(si, kind=8), & +! int(la, kind=8)) +! H(mu,nu) = H(mu,nu) + P(si,la) * ERI_chem(munusila) +! enddo +! enddo +! enddo +! enddo + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (nu, la, si, mu, & + !$OMP nunu0, nunu, nula, lala0, lala, lasi, numu, & + !$OMP nunununu0, nunununu, nununula, numulala, numunula, & + !$OMP nunulala, lalanunu, lalanumu, nunulasi, lasinunu, & + !$OMP numunumu0, nununumu, numulasi, lasinumu) & + !$OMP SHARED (nBas, H, P, ERI_chem) + !$OMP DO do nu = 1, nBas nunu0 = shiftr(nu * (nu - 1), 1) @@ -63,7 +88,7 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H) nunununu = nunununu0 + nunu H(nu,nu) = P(nu,nu) * ERI_chem(nunununu) - do la = 1, nu-1 + do la = 1, nu - 1 lala0 = shiftr(la * (la - 1), 1) @@ -150,7 +175,8 @@ subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H) enddo ! mu enddo ! nu - + !$OMP END DO + !$OMP END PARALLEL do nu = 1, nBas do mu = nu+1, nBas @@ -163,3 +189,4 @@ end subroutine ! --- + diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 4781436..581e8db 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/src/AOtoMO/exchange_matrix_AO_basis.f90 @@ -47,22 +47,46 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) integer :: mu, nu, la, si integer :: nunu, nula, lanu, lasi, nusi, sinu integer :: numu, mumu, mula, lamu, musi, simu + integer :: nunu0, lala0, mumu0 integer*8 :: nunununu, nulanula, lanulanu, nulanusi integer*8 :: munulasi, lanunusi, lanusinu, numumumu integer*8 :: nulamula, nulalamu, lanulamu, nulamusi integer*8 :: nulasimu, lanumusi, lanusimu, simunula - integer*8 :: simulanu + integer*8 :: simulanu, nulanula0, lanulanu0 + +! integer*8 :: munusila +! integer*8, external :: Yoshimine_4ind +! +! do nu = 1, nBas +! do mu = 1, nu +! K(mu,nu) = 0.d0 +! do la = 1, nBas +! do si = 1, nBas +! munusila = Yoshimine_4ind(int(mu, kind=8), & +! int(si, kind=8), & +! int(la, kind=8), & +! int(nu, kind=8)) +! K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(munusila) +! enddo +! enddo +! enddo +! enddo +! !$OMP PARALLEL SHARED (NONE) +! !$OMP PRIVATE () +! !$OMP SHARED () +! !$OMP DO do nu = 1, nBas - nunu = shiftr(nu * (nu - 1), 1) + nu + nunu0 = shiftr(nu * (nu - 1), 1) + nunu = nunu0 + nu nunununu = shiftr(nunu * (nunu - 1), 1) + nunu K(nu,nu) = -P(nu,nu) * ERI_chem(nunununu) do la = 1, nu - 1 - nula = shiftr(nu * (nu - 1), 1) + la + nula = nunu0 + la nulanula = shiftr(nula * (nula - 1), 1) + nula K(nu,nu) = K(nu,nu) - P(la,la) * ERI_chem(nulanula) enddo @@ -74,24 +98,23 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) enddo do la = 1, nu - nula = shiftr(nu * (nu - 1), 1) + la + nula = nunu0 + la + nulanula0 = shiftr(nula * (nula - 1), 1) do si = 1, la - 1 - nusi = shiftr(nu * (nu - 1), 1) + si - nulanusi = shiftr(nula * (nula - 1), 1) + nusi + nulanusi = nulanula0 + nunu0 + si K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(nulanusi) enddo enddo do la = nu + 1, nBas lanu = shiftr(la * (la - 1), 1) + nu + lanulanu0 = shiftr(lanu * (lanu - 1), 1) do si = 1, nu - nusi = shiftr(nu * (nu - 1), 1) + si - lanunusi = shiftr(lanu * (lanu - 1), 1) + nusi + lanunusi = lanulanu0 + nunu0 + si K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(lanunusi) enddo do si = nu + 1, la - 1 - sinu = shiftr(si * (si - 1), 1) + nu - lanusinu = shiftr(lanu * (lanu - 1), 1) + sinu + lanusinu = lanulanu0 + shiftr(si * (si - 1), 1) + nu K(nu,nu) = K(nu,nu) - 2.d0 * P(si,la) * ERI_chem(lanusinu) enddo enddo @@ -99,75 +122,71 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) do mu = 1, nu - 1 - numu = shiftr(nu * (nu - 1), 1) + mu - mumu = shiftr(mu * (mu - 1), 1) + mu + numu = nunu0 + mu + mumu0 = shiftr(mu * (mu - 1), 1) + mumu = mumu0 + mu numumumu = shiftr(numu * (numu - 1), 1) + mumu K(mu,nu) = - P(mu,mu) * ERI_chem(numumumu) do la = 1, mu - 1 - mula = shiftr(mu * (mu - 1), 1) + la - nula = shiftr(nu * (nu - 1), 1) + la - nulamula = shiftr(nula * (nula - 1), 1) + mula + nula = nunu0 + la + nulamula = shiftr(nula * (nula - 1), 1) + mumu0 + la K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(nulamula) enddo do la = mu + 1, nu - lamu = shiftr(la * (la - 1), 1) + mu - nula = shiftr(nu * (nu - 1), 1) + la - nulalamu = shiftr(nula * (nula - 1), 1) + lamu + nula = nunu0 + la + nulalamu = shiftr(nula * (nula - 1), 1) + shiftr(la * (la - 1), 1) + mu K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(nulalamu) enddo do la = nu + 1, nBas - lamu = shiftr(la * (la - 1), 1) + mu - lanu = shiftr(la * (la - 1), 1) + nu - lanulamu = shiftr(lanu * (lanu - 1), 1) + lamu + lala0 = shiftr(la * (la - 1), 1) + lanu = lala0 + nu + lanulamu = shiftr(lanu * (lanu - 1), 1) + lala0 + mu K(mu,nu) = K(mu,nu) - P(la,la) * ERI_chem(lanulamu) enddo do la = 1, mu - nula = shiftr(nu * (nu - 1), 1) + la + nula = nunu0 + la + nulanula0 = shiftr(nula * (nula - 1), 1) do si = 1, la - 1 - musi = shiftr(mu * (mu - 1), 1) + si - nulamusi = shiftr(nula * (nula - 1), 1) + musi + nulamusi = nulanula0 + mumu0 + si K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi) enddo enddo do la = mu + 1, nu - nula = shiftr(nu * (nu - 1), 1) + la + nula = nunu0 + la + nulanula0 = shiftr(nula * (nula - 1), 1) do si = 1, mu - musi = shiftr(mu * (mu - 1), 1) + si - nulamusi = shiftr(nula * (nula - 1), 1) + musi + nulamusi = nulanula0 + mumu0 + si K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulamusi) enddo do si = mu + 1, la - 1 - simu = shiftr(si * (si - 1), 1) + mu - nulasimu = shiftr(nula * (nula - 1), 1) + simu + nulasimu = nulanula0 + shiftr(si * (si - 1), 1) + mu K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) enddo enddo do la = nu + 1, nBas lanu = shiftr(la * (la - 1), 1) + nu + lanulanu0 = shiftr(lanu * (lanu - 1), 1) do si = 1, mu - musi = shiftr(mu * (mu - 1), 1) + si - lanumusi = shiftr(lanu * (lanu - 1), 1) + musi + lanumusi = lanulanu0 + mumu0 + si K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(lanumusi) enddo do si = mu + 1, la - 1 - simu = shiftr(si * (si - 1), 1) + mu - lanusimu = shiftr(lanu * (lanu - 1), 1) + simu + lanusimu = lanulanu0 + shiftr(si * (si - 1), 1) + mu K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(lanusimu) enddo enddo do la = 1, mu - nula = shiftr(nu * (nu - 1), 1) + la + nula = nunu0 + la + nulanula0 = shiftr(nula * (nula - 1), 1) do si = la + 1, mu - musi = shiftr(mu * (mu - 1) , 1) + si - nulamusi = shiftr(nula * (nula - 1), 1) + musi + nulamusi = nulanula0 + mumu0 + si 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 + nulasimu = nulanula0 + shiftr(si * (si - 1), 1) + mu K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) enddo do si = nu, nBas @@ -177,10 +196,10 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) enddo enddo do la = mu + 1, nu - nula = shiftr(nu * (nu - 1), 1) + la + nula = nunu0 + la + nulanula0 = shiftr(nula * (nula - 1), 1) do si = la + 1, nu - simu = shiftr(si * (si - 1), 1) + mu - nulasimu = shiftr(nula * (nula - 1), 1) + simu + nulasimu = nulanula0 + shiftr(si * (si - 1), 1) + mu K(mu,nu) = K(mu,nu) - P(si,la) * ERI_chem(nulasimu) enddo do si = nu + 1, nBas @@ -200,6 +219,8 @@ subroutine exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) enddo ! mu enddo ! nu +! !$OMP END DO +! !$OMP END PARALLEL do nu = 1, nBas diff --git a/src/HF/RHF_hpc.f90 b/src/HF/RHF_hpc.f90 index 44095d3..740bc94 100644 --- a/src/HF/RHF_hpc.f90 +++ b/src/HF/RHF_hpc.f90 @@ -107,18 +107,14 @@ 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 wall_time(t1) - do ii = 1, 5 - call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J) - enddo + call Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, J) call wall_time(t2) - print*, " J built in (sec):", (t2-t1) / 5.d0 + print*, " J built in (sec):", (t2-t1) call wall_time(t1) - do ii = 1, 5 - call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) - enddo + call exchange_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, K) call wall_time(t2) - print*, " K built in (sec):", (t2-t1) / 5.d0 + print*, " K built in (sec):", (t2-t1) allocate(ERI_phys(nBas,nBas,nBas,nBas)) @@ -128,25 +124,21 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh call read_2e_integrals(working_dir, nBas, ERI_phys) call wall_time(t1) - do ii = 1, 5 - call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb) - enddo + call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb) call wall_time(t2) - print*, " J_deb built in (sec):", (t2-t1) / 5.d0 + print*, " J_deb built in (sec):", (t2-t1) call wall_time(t1) - do ii = 1, 5 - call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb) - enddo + call exchange_matrix_AO_basis(nBas, P, ERI_phys, K_deb) call wall_time(t2) - print*, " K_deb built in (sec):", (t2-t1) / 5.d0 + print*, " K_deb built in (sec):", (t2-t1) 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 + if(diff_loc .gt. 1d-10) then print*, 'error in J on: ', jj, ii print*, J(jj,ii), J_deb(jj,ii) stop @@ -161,7 +153,7 @@ subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_sh 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 + if(diff_loc .gt. 1d-10) then print*, 'error in K on: ', jj, ii print*, K(jj,ii), K_deb(jj,ii) stop