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

ERI dat -> binary

Browse Source
This commit is contained in:
Abdallah Ammar 2024-08-28 15:07:20 +02:00
parent 1de213dc89
commit 567e77dd4a
4 changed files with 90 additions and 77 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
PyDuck.py
View File

@ -18,6 +18,7 @@ parser.add_argument('-b', '--basis', type=str, required=True, help='Name of the
parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.') parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.')
parser.add_argument('-c', '--charge', type=int, default=0, help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0') parser.add_argument('-c', '--charge', type=int, default=0, help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0')
parser.add_argument('--cartesian', default=False, action='store_true', help='Add this option if you want to use cartesian basis functions.') parser.add_argument('--cartesian', default=False, action='store_true', help='Add this option if you want to use cartesian basis functions.')
parser.add_argument('--print_2e', default=False, action='store_true', help='Add this option if you want to print 2e-integrals.')
parser.add_argument('-fc', '--frozen_core', type=bool, default=False, help='Freeze core MOs. Default is false') parser.add_argument('-fc', '--frozen_core', type=bool, default=False, help='Freeze core MOs. Default is false')
parser.add_argument('-m', '--multiplicity', type=int, default=1, help='Spin multiplicity. Default is 1 therefore singlet') parser.add_argument('-m', '--multiplicity', type=int, default=1, help='Spin multiplicity. Default is 1 therefore singlet')
parser.add_argument('--working_dir', type=str, default=QuAcK_dir, help='Set a working directory to run the calculation.') parser.add_argument('--working_dir', type=str, default=QuAcK_dir, help='Set a working directory to run the calculation.')
@ -32,6 +33,7 @@ frozen_core=args.frozen_core
multiplicity=args.multiplicity multiplicity=args.multiplicity
xyz=args.xyz + '.xyz' xyz=args.xyz + '.xyz'
cartesian=args.cartesian cartesian=args.cartesian
print_2e=args.print_2e
working_dir=args.working_dir working_dir=args.working_dir
#Read molecule #Read molecule
@ -90,11 +92,10 @@ t1e = mol.intor('int1e_kin') #Kinetic energy matrix elements
dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators
x,y,z = dipole[0],dipole[1],dipole[2] x,y,z = dipole[0],dipole[1],dipole[2]
norb = len(ovlp) norb = len(ovlp) # nBAS_AOs
subprocess.call(['rm', working_dir + '/int/nBas.dat']) subprocess.call(['rm', working_dir + '/int/nBas.dat'])
f = open(working_dir+'/int/nBas.dat','w') f = open(working_dir+'/int/nBas.dat','w')
f.write(str(norb)) f.write(" {} ".format(str(norb)))
f.write(' ')
f.close() f.close()
@ -122,7 +123,6 @@ write_matrix_to_file(y,norb,working_dir+'/int/y.dat')
subprocess.call(['rm', working_dir + '/int/z.dat']) subprocess.call(['rm', working_dir + '/int/z.dat'])
write_matrix_to_file(z,norb,working_dir+'/int/z.dat') write_matrix_to_file(z,norb,working_dir+'/int/z.dat')
#Write two-electron integrals
eri_ao = mol.intor('int2e') eri_ao = mol.intor('int2e')
def write_tensor_to_file(tensor,size,file,cutoff=1e-15): def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
@ -132,12 +132,23 @@ def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
for k in range(i,size): for k in range(i,size):
for l in range(j,size): for l in range(j,size):
if abs(tensor[i][k][j][l]) > cutoff: if abs(tensor[i][k][j][l]) > cutoff:
#f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l])) f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
f.write('\n') f.write('\n')
f.close() f.close()
subprocess.call(['rm', working_dir + '/int/ERI.dat']) # Write two-electron integrals
write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat') if print_2e:
# (formatted)
subprocess.call(['rm', working_dir + '/int/ERI.dat'])
write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat')
else:
# (binary)
subprocess.call(['rm', working_dir + '/int/ERI.bin'])
# chem -> phys notation
eri_ao = eri_ao.transpose(0, 2, 1, 3)
eri_ao.tofile('int/ERI.bin')
#Execute the QuAcK fortran program #Execute the QuAcK fortran program
subprocess.call(QuAcK_dir+'/bin/QuAcK') subprocess.call(QuAcK_dir+'/bin/QuAcK')

59
src/QuAcK/QuAcK.f90
View File

@ -15,7 +15,7 @@ program QuAcK
logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh
integer :: nNuc,nBas integer :: nNuc,nBas_AOs
integer :: nC(nspin) integer :: nC(nspin)
integer :: nO(nspin) integer :: nO(nspin)
integer :: nV(nspin) integer :: nV(nspin)
@ -120,15 +120,15 @@ program QuAcK
doACFDT,exchange_kernel,doXBS, & doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA) dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
!------------------------------------------------! !---------------------------------------------------!
! Read input information ! ! Read input information !
!------------------------------------------------! !---------------------------------------------------!
! nC = number of core orbitals ! ! nC = number of core orbitals !
! nO = number of occupied orbitals ! ! nO = number of occupied orbitals !
! nV = number of virtual orbitals (see below) ! ! nV = number of virtual orbitals (see below) !
! nR = number of Rydberg orbitals ! ! nR = number of Rydberg orbitals !
! nBas = number of basis functions (see below) ! ! nBas_AOs = number of basis functions in AOs !
!------------------------------------------------! !---------------------------------------------------!
call read_molecule(nNuc,nO,nC,nR) call read_molecule(nNuc,nO,nC,nR)
allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
@ -141,7 +141,7 @@ program QuAcK
! Read basis set information from PySCF ! ! Read basis set information from PySCF !
!---------------------------------------! !---------------------------------------!
call read_basis_pyscf(nBas,nO,nV) call read_basis_pyscf(nBas_AOs, nO, nV)
!--------------------------------------! !--------------------------------------!
! Read one- and two-electron integrals ! ! Read one- and two-electron integrals !
@ -149,26 +149,31 @@ program QuAcK
! Memory allocation for one- and two-electron integrals ! Memory allocation for one- and two-electron integrals
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), & allocate(S(nBas_AOs,nBas_AOs))
ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart)) allocate(T(nBas_AOs,nBas_AOs))
allocate(V(nBas_AOs,nBas_AOs))
allocate(Hc(nBas_AOs,nBas_AOs))
allocate(X(nBas_AOs,nBas_AOs))
allocate(ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs))
allocate(dipole_int_AO(nBas_AOs,nBas_AOs,ncart))
! Read integrals ! Read integrals
call wall_time(start_int) call wall_time(start_int)
call read_integrals(nBas,S,T,V,Hc,ERI_AO) call read_integrals(nBas_AOs, S(1,1), T(1,1), V(1,1), Hc(1,1), ERI_AO(1,1,1,1))
call read_dipole_integrals(nBas,dipole_int_AO) call read_dipole_integrals(nBas_AOs, dipole_int_AO)
call wall_time(end_int) call wall_time(end_int)
t_int = end_int - start_int t_int = end_int - start_int
write(*,*) write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds'
write(*,*) write(*,*)
! Compute orthogonalization matrix ! Compute orthogonalization matrix
call orthogonalization_matrix(nBas,S,X) call orthogonalization_matrix(nBas_AOs, S, X)
!---------------------! !---------------------!
! Choose QuAcK branch ! ! Choose QuAcK branch !
@ -200,7 +205,7 @@ program QuAcK
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & nNuc,nBas_AOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
@ -216,7 +221,7 @@ program QuAcK
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & nNuc,nBas_AOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
@ -228,13 +233,13 @@ program QuAcK
!--------------------------! !--------------------------!
if(doGQuAcK) & if(doGQuAcK) &
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & nNuc,nBas_AOs,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
!-----------! !-----------!
@ -256,7 +261,7 @@ program QuAcK
call wall_time(end_QuAcK) call wall_time(end_QuAcK)
t_QuAcK = end_QuAcK - start_QuAcK t_QuAcK = end_QuAcK - start_QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for QuAcK = ',t_QuAcK,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds'
write(*,*) write(*,*)
end program end program

14
src/utils/read_basis_pyscf.f90
View File

@ -1,4 +1,4 @@
subroutine read_basis_pyscf(nBas,nO,nV) subroutine read_basis_pyscf(nBas_AOs, nO, nV)
! Read basis set information from PySCF ! Read basis set information from PySCF
@ -14,23 +14,23 @@ subroutine read_basis_pyscf(nBas,nO,nV)
! Output variables ! Output variables
integer,intent(out) :: nV(nspin) integer,intent(out) :: nV(nspin)
integer,intent(out) :: nBas integer,intent(out) :: nBas_AOs
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Primary basis set information ! Primary basis set information
!------------------------------------------------------------------------ !------------------------------------------------------------------------
open(unit=3,file='int/nBas.dat') open(unit=3,file='int/nBas.dat')
read(3,*) nBas read(3, *) nBas_AOs
close(unit=3) close(unit=3)
write(*,'(A28)') '------------------' write(*,'(A38)') '--------------------------------------'
write(*,'(A28,1X,I16)') 'Number of basis functions',nBas write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs
write(*,'(A28)') '------------------' write(*,'(A38)') '--------------------------------------'
write(*,*) write(*,*)
! Number of virtual orbitals ! Number of virtual orbitals
nV(:) = nBas - nO(:) nV(:) = nBas_AOs - nO(:)
end subroutine end subroutine

71
src/utils/read_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine read_integrals(nBas,S,T,V,Hc,G) subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
! Read one- and two-electron integrals from files ! Read one- and two-electron integrals from files
@ -7,7 +7,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
! Input variables ! Input variables
integer,intent(in) :: nBas integer,intent(in) :: nBas_AOs
! Local variables ! Local variables
@ -18,11 +18,11 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
! Output variables ! Output variables
double precision,intent(out) :: S(nBas,nBas) double precision,intent(out) :: S(nBas_AOs,nBas_AOs)
double precision,intent(out) :: T(nBas,nBas) double precision,intent(out) :: T(nBas_AOs,nBas_AOs)
double precision,intent(out) :: V(nBas,nBas) double precision,intent(out) :: V(nBas_AOs,nBas_AOs)
double precision,intent(out) :: Hc(nBas,nBas) double precision,intent(out) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(out) :: G(nBas,nBas,nBas,nBas) double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
! Open file with integrals ! Open file with integrals
@ -35,7 +35,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
open(unit=8 ,file='int/Ov.dat') open(unit=8 ,file='int/Ov.dat')
open(unit=9 ,file='int/Kin.dat') open(unit=9 ,file='int/Kin.dat')
open(unit=10,file='int/Nuc.dat') open(unit=10,file='int/Nuc.dat')
open(unit=11,file='int/ERI.dat')
open(unit=21,file='int/x.dat') open(unit=21,file='int/x.dat')
open(unit=22,file='int/y.dat') open(unit=22,file='int/y.dat')
@ -75,31 +74,29 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
Hc(:,:) = T(:,:) + V(:,:) Hc(:,:) = T(:,:) + V(:,:)
! Read nuclear integrals ! Read 2e-integrals
G(:,:,:,:) = 0d0 ! ! formatted file
do ! open(unit=11, file='int/ERI.dat')
read(11,*,end=11) mu,nu,la,si,ERI ! G(:,:,:,:) = 0d0
! do
! read(11,*,end=11) mu, nu, la, si, ERI
! ERI = lambda*ERI
! G(mu,nu,la,si) = ERI ! <12|34>
! G(la,nu,mu,si) = ERI ! <32|14>
! G(mu,si,la,nu) = ERI ! <14|32>
! G(la,si,mu,nu) = ERI ! <34|12>
! G(si,mu,nu,la) = ERI ! <41|23>
! G(nu,la,si,mu) = ERI ! <23|41>
! G(nu,mu,si,la) = ERI ! <21|43>
! G(si,la,nu,mu) = ERI ! <43|21>
! end do
! 11 close(unit=11)
ERI = lambda*ERI ! binary file
! <12|34> open(unit=11, file='int/ERI.bin', form='unformatted', access='stream')
G(mu,nu,la,si) = ERI read(11) G
! <32|14> close(11)
G(la,nu,mu,si) = ERI
! <14|32>
G(mu,si,la,nu) = ERI
! <34|12>
G(la,si,mu,nu) = ERI
! <41|23>
G(si,mu,nu,la) = ERI
! <23|41>
G(nu,la,si,mu) = ERI
! <21|43>
G(nu,mu,si,la) = ERI
! <43|21>
G(si,la,nu,mu) = ERI
end do
11 close(unit=11)
! Print results ! Print results
@ -107,24 +104,24 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Overlap integrals' write(*,'(A28)') 'Overlap integrals'
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
call matout(nBas,nBas,S) call matout(nBas_AOs,nBas_AOs,S)
write(*,*) write(*,*)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Kinetic integrals' write(*,'(A28)') 'Kinetic integrals'
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
call matout(nBas,nBas,T) call matout(nBas_AOs,nBas_AOs,T)
write(*,*) write(*,*)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Nuclear integrals' write(*,'(A28)') 'Nuclear integrals'
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
call matout(nBas,nBas,V) call matout(nBas_AOs,nBas_AOs,V)
write(*,*) write(*,*)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals' write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
do la=1,nBas do la=1,nBas_AOs
do si=1,nBas do si=1,nBas_AOs
call matout(nBas,nBas,G(1,1,la,si)) call matout(nBas_AOs, nBas_AOs, G(1,1,la,si))
end do end do
end do end do
write(*,*) write(*,*)