PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 04:13:52 +01:00
Code Releases Activity

added Hartree and EX terms with s8-sym

Browse Source
This commit is contained in:
Abdallah Ammar 2024-12-10 19:14:33 +01:00
parent 671c6f4f12
commit 37420d1207
4 changed files with 115 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
PyDuck.py
View File

@ -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()

165
src/AOtoMO/exchange_matrix_AO_basis.f90
View File

@ -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

35
src/HF/RHF_hpc.f90
View File

@ -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

6
src/utils/utils.f90
View File

@ -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