diff --git a/PyDuck.py b/PyDuck.py index c74dc22..a8e9033 100644 --- a/PyDuck.py +++ b/PyDuck.py @@ -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('-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('--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('-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.') @@ -32,6 +33,7 @@ frozen_core=args.frozen_core multiplicity=args.multiplicity xyz=args.xyz + '.xyz' cartesian=args.cartesian +print_2e=args.print_2e working_dir=args.working_dir #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 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']) f = open(working_dir+'/int/nBas.dat','w') -f.write(str(norb)) -f.write(' ') +f.write(" {} ".format(str(norb))) 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']) write_matrix_to_file(z,norb,working_dir+'/int/z.dat') -#Write two-electron integrals eri_ao = mol.intor('int2e') 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 l in range(j,size): 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('\n') f.close() -subprocess.call(['rm', working_dir + '/int/ERI.dat']) -write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat') +# Write two-electron integrals +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 subprocess.call(QuAcK_dir+'/bin/QuAcK') diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 6af981b..d1aa007 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -15,7 +15,7 @@ program QuAcK logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh - integer :: nNuc,nBas + integer :: nNuc,nBas_AOs integer :: nC(nspin) integer :: nO(nspin) integer :: nV(nspin) @@ -120,15 +120,15 @@ program QuAcK doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA) -!------------------------------------------------! -! Read input information ! -!------------------------------------------------! -! nC = number of core orbitals ! -! nO = number of occupied orbitals ! -! nV = number of virtual orbitals (see below) ! -! nR = number of Rydberg orbitals ! -! nBas = number of basis functions (see below) ! -!------------------------------------------------! +!---------------------------------------------------! +! Read input information ! +!---------------------------------------------------! +! nC = number of core orbitals ! +! nO = number of occupied orbitals ! +! nV = number of virtual orbitals (see below) ! +! nR = number of Rydberg orbitals ! +! nBas_AOs = number of basis functions in AOs ! +!---------------------------------------------------! call read_molecule(nNuc,nO,nC,nR) allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) @@ -141,7 +141,7 @@ program QuAcK ! 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 ! @@ -149,26 +149,31 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), & - ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart)) + allocate(S(nBas_AOs,nBas_AOs)) + 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 call wall_time(start_int) - call read_integrals(nBas,S,T,V,Hc,ERI_AO) - call read_dipole_integrals(nBas,dipole_int_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_AOs, dipole_int_AO) call wall_time(end_int) t_int = end_int - start_int 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(*,*) ! Compute orthogonalization matrix - call orthogonalization_matrix(nBas,S,X) + call orthogonalization_matrix(nBas_AOs, S, X) !---------------------! ! Choose QuAcK branch ! @@ -200,7 +205,7 @@ program QuAcK 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, & 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, & 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, & @@ -216,7 +221,7 @@ program QuAcK 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, & 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, & 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, & @@ -228,13 +233,13 @@ program QuAcK !--------------------------! if(doGQuAcK) & - call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & - 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, & - 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_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & + call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & + doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + 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_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, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) !-----------! @@ -256,7 +261,7 @@ program QuAcK call wall_time(end_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(*,*) end program diff --git a/src/utils/read_basis_pyscf.f90 b/src/utils/read_basis_pyscf.f90 index a677323..027fcac 100644 --- a/src/utils/read_basis_pyscf.f90 +++ b/src/utils/read_basis_pyscf.f90 @@ -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 @@ -14,23 +14,23 @@ subroutine read_basis_pyscf(nBas,nO,nV) ! Output variables integer,intent(out) :: nV(nspin) - integer,intent(out) :: nBas + integer,intent(out) :: nBas_AOs !------------------------------------------------------------------------ ! Primary basis set information !------------------------------------------------------------------------ open(unit=3,file='int/nBas.dat') - read(3,*) nBas + read(3, *) nBas_AOs close(unit=3) - write(*,'(A28)') '------------------' - write(*,'(A28,1X,I16)') 'Number of basis functions',nBas - write(*,'(A28)') '------------------' + write(*,'(A38)') '--------------------------------------' + write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs + write(*,'(A38)') '--------------------------------------' write(*,*) ! Number of virtual orbitals - nV(:) = nBas - nO(:) + nV(:) = nBas_AOs - nO(:) end subroutine diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index 9e215e4..ef8e1d9 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -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 @@ -7,7 +7,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs ! Local variables @@ -18,11 +18,11 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) ! Output variables - double precision,intent(out) :: S(nBas,nBas) - double precision,intent(out) :: T(nBas,nBas) - double precision,intent(out) :: V(nBas,nBas) - double precision,intent(out) :: Hc(nBas,nBas) - double precision,intent(out) :: G(nBas,nBas,nBas,nBas) + double precision,intent(out) :: S(nBas_AOs,nBas_AOs) + double precision,intent(out) :: T(nBas_AOs,nBas_AOs) + double precision,intent(out) :: V(nBas_AOs,nBas_AOs) + double precision,intent(out) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) ! 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=9 ,file='int/Kin.dat') open(unit=10,file='int/Nuc.dat') - open(unit=11,file='int/ERI.dat') open(unit=21,file='int/x.dat') open(unit=22,file='int/y.dat') @@ -75,31 +74,29 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) Hc(:,:) = T(:,:) + V(:,:) -! Read nuclear integrals +! Read 2e-integrals - G(:,:,:,:) = 0d0 - do - read(11,*,end=11) mu,nu,la,si,ERI +! ! formatted file +! open(unit=11, file='int/ERI.dat') +! 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 -! <12|34> - G(mu,nu,la,si) = ERI -! <32|14> - 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) + ! binary file + open(unit=11, file='int/ERI.bin', form='unformatted', access='stream') + read(11) G + close(11) ! Print results @@ -107,24 +104,24 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Overlap integrals' write(*,'(A28)') '----------------------' - call matout(nBas,nBas,S) + call matout(nBas_AOs,nBas_AOs,S) write(*,*) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Kinetic integrals' write(*,'(A28)') '----------------------' - call matout(nBas,nBas,T) + call matout(nBas_AOs,nBas_AOs,T) write(*,*) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Nuclear integrals' write(*,'(A28)') '----------------------' - call matout(nBas,nBas,V) + call matout(nBas_AOs,nBas_AOs,V) write(*,*) write(*,'(A28)') '----------------------' write(*,'(A28)') 'Electron repulsion integrals' write(*,'(A28)') '----------------------' - do la=1,nBas - do si=1,nBas - call matout(nBas,nBas,G(1,1,la,si)) + do la=1,nBas_AOs + do si=1,nBas_AOs + call matout(nBas_AOs, nBas_AOs, G(1,1,la,si)) end do end do write(*,*)