LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

Working on s8 rep of 2e-integrals

Browse Source
This commit is contained in:
Abdallah Ammar 2024-12-07 02:20:05 +01:00
parent ce1882cd4e
commit b7af468f11
14 changed files with 1019 additions and 172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
PyDuck.py
View File

@ -24,8 +24,10 @@ 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('--print_2e', default=True, action='store_true', help='If True, print 2e-integrals to disk.')
parser.add_argument('--formatted_2e', default=False, action='store_true', help='Add this option if you want to print formatted 2e-integrals.')
parser.add_argument('--mmap_2e', default=False, action='store_true', help='If True, avoid using DRAM when generating 2e-integrals.') parser.add_argument('--mmap_2e', default=False, action='store_true', help='If True, avoid using DRAM when generating 2e-integrals.')
parser.add_argument('--aosym_2e', default=False, action='store_true', help='If True, use 8-fold symmetry 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.')
@ -41,7 +43,9 @@ multiplicity=args.multiplicity
xyz=args.xyz + '.xyz' xyz=args.xyz + '.xyz'
cartesian=args.cartesian cartesian=args.cartesian
print_2e=args.print_2e print_2e=args.print_2e
formatted_2e=args.formatted_2e
mmap_2e=args.mmap_2e mmap_2e=args.mmap_2e
aosym_2e=args.aosym_2e
working_dir=args.working_dir working_dir=args.working_dir
#Read molecule #Read molecule
@ -63,6 +67,7 @@ mol = gto.M(
basis = input_basis, basis = input_basis,
charge = charge, charge = charge,
spin = multiplicity - 1 spin = multiplicity - 1
# symmetry = True # Enable symmetry
) )
#Fix the unit for the lengths #Fix the unit for the lengths
@ -144,16 +149,25 @@ def write_tensor_to_file(tensor,size,file_name,cutoff=1e-15):
f.write('\n') f.write('\n')
f.close() f.close()
# Write two-electron integrals to HD
ti_2e = time.time()
if print_2e: if print_2e:
# (formatted) # Write two-electron integrals to HD
ti_2e = time.time()
if formatted_2e:
output_file_path = working_dir + '/int/ERI.dat' output_file_path = working_dir + '/int/ERI.dat'
subprocess.call(['rm', '-f', output_file_path]) subprocess.call(['rm', '-f', output_file_path])
eri_ao = mol.intor('int2e') eri_ao = mol.intor('int2e')
write_tensor_to_file(eri_ao, norb, output_file_path) write_tensor_to_file(eri_ao, norb, output_file_path)
else:
# (binary) if aosym_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()
else:
output_file_path = working_dir + '/int/ERI.bin' output_file_path = working_dir + '/int/ERI.bin'
subprocess.call(['rm', '-f', output_file_path]) subprocess.call(['rm', '-f', output_file_path])
if(mmap_2e): if(mmap_2e):
@ -170,8 +184,11 @@ else:
f = open(output_file_path, 'w') f = open(output_file_path, 'w')
eri_ao.tofile(output_file_path) eri_ao.tofile(output_file_path)
f.close() f.close()
te_2e = time.time()
print("Wall time for writing 2e-integrals (physicist notation) to disk: {:.3f} seconds".format(te_2e - ti_2e)) te_2e = time.time()
print("Wall time for writing 2e-integrals to disk: {:.3f} seconds".format(te_2e - ti_2e))
sys.stdout.flush()

2
input/hardware
View File

@ -1,2 +0,0 @@
# if True (T), use GPU
F

4
input/hpc_flags Normal file
View File

@ -0,0 +1,4 @@
# if True (T), switch to HPC mode
F
# if True (T), use GPU
F

103
src/AOtoMO/Hartree_matrix_AO_basis.f90
View File

@ -32,3 +32,106 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H)
end do end do
end subroutine end subroutine
! ---
subroutine Hartree_matrix_AO_basis_hpc(nBas, ERI_size, P, ERI_chem, H)
implicit none
integer, intent(in) :: nBas
integer*8, intent(in) :: ERI_size
double precision, intent(in) :: P(nBas,nBas)
double precision, intent(in) :: ERI_chem(ERI_size)
double precision, intent(out) :: H(nBas,nBas)
integer :: mu, nu, la, si
integer*8 :: munu0, munu
integer*8 :: sila0, sila
integer*8 :: munulasi0, munulasi
integer*8, external :: Yoshimine_ind
do nu = 1, nBas
do mu = 1, nBas
H(mu,nu) = 0.d0
do si = 1, nBas
do la = 1, nBas
munulasi = Yoshimine_ind(mu, nu, la, si)
H(mu,nu) = H(mu,nu) + P(la,si) * ERI_chem(munulasi)
enddo
enddo
enddo
enddo
! do nu = 1, nBas
! munu0 = (nu * (nu + 1)) / 2
!
! do mu = 1, nu
! munu = munu0 + mu
! munulasi0 = (munu * (munu + 1)) / 2
!
! H(mu,nu) = 0.d0
!
! do si = 1, nu
! sila0 = (si * (si + 1)) / 2
!
! do la = 1, si
! sila = sila0 + la
!
! if(nu == si .and. mu < la) cycle
!
! munulasi = munulasi0 + sila
!
! H(mu,nu) = H(mu,nu) + 4.d0 * P(la,si) * ERI_chem(munulasi)
! enddo
! enddo
! enddo
! enddo
!
!
! do nu = 1, nBas
! do mu = nu+1, nBas
! H(mu,nu) = H(nu,mu)
! enddo
! enddo
return
end subroutine
! ---
integer*8 function Yoshimine_ind(a, b, c, d)
implicit none
integer, intent(in) :: a, b, c, d
integer*8 :: ab, cd, abcd
if(a > b) then
ab = (a * (a - 1)) / 2 + b
else
ab = (b * (b - 1)) / 2 + a
endif
if(c > d) then
cd = (c * (c - 1)) / 2 + d
else
cd = (d * (d - 1)) / 2 + c
endif
if(ab > cd) then
abcd = (ab * (ab - 1)) / 2 + cd
else
abcd = (cd * (cd - 1)) / 2 + ab
endif
Yoshimine_ind = abcd
return
end
! ---

247
src/HF/RHF_hpc.f90 Normal file
View File

@ -0,0 +1,247 @@
subroutine RHF_hpc(working_dir,dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,dipole_int,X,ERHF,eHF,c,P,F)
! Perform restricted Hartree-Fock calculation
implicit none
include 'parameters.h'
! Input variables
character(len=256),intent(in) :: working_dir
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: ii, jj
integer :: nSCF
integer :: nBas_Sq
integer :: n_diis
integer*8 :: ERI_size
double precision :: diff, diff_loc
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: dipole(ncart)
double precision :: Conv
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: err(:,:)
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ERI_chem(:)
double precision,allocatable :: ERI_phys(:,:,:,:), J_deb(:,:)
! Output variables
double precision,intent(out) :: ERHF
double precision,intent(out) :: eHF(nOrb)
double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'****************************************'
write(*,*)'* Restricted HF Calculation (HPC mode) *'
write(*,*)'****************************************'
write(*,*)
! Useful quantities
nBas_Sq = nBas*nBas
! Memory allocation
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(err(nBas,nBas))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Guess coefficients and density matrix
call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
c(1,1), nBas, c(1,1), nBas, &
0.d0, P(1,1), nBas)
ERI_size = (nBas * (nBas + 1)) / 2
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)
allocate(J_deb(nBas,nBas))
allocate(ERI_phys(nBas,nBas,nBas,nBas))
call read_2e_integrals(working_dir, nBas, ERI_phys)
call Hartree_matrix_AO_basis(nBas, P, ERI_phys, J_deb)
print*, 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-13) then
print*, 'error on: ', jj, ii
print*, J(jj,ii), J_deb(jj,ii)
stop
endif
diff = diff + diff_loc
enddo
enddo
print*, 'total diff = ', diff
stop
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'-----------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(RHF)','|','EJ(RHF)','|','EK(RHF)','|','Conv','|'
write(*,*)'-----------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Build Fock matrix
call Hartree_matrix_AO_basis(nBas,P,ERI_phys,J)
call exchange_matrix_AO_basis(nBas,P,ERI_phys,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
! Check convergence
err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
if(nSCF > 1) Conv = maxval(abs(err))
! Kinetic energy
ET = trace_matrix(nBas, matmul(P, T))
! Potential energy
EV = trace_matrix(nBas, matmul(P, V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas, matmul(P, K))
! Total energy
ERHF = ET + EV + EJ + EK
! DIIS extrapolation
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
endif
! Level shift
if(level_shift > 0d0 .and. Conv > thresh) then
call level_shifting(level_shift,nBas,nOrb,nO,S,c,F)
endif
! Diagonalize Fock matrix
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X,cp)
! Density matrix
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
c(1,1), nBas, c(1,1), nBas, &
0.d0, P(1,1), nBas)
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',ERHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
end do
write(*,*)'-----------------------------------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
stop
end if
! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_RHF(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Testing zone
if(dotest) then
call dump_test_value('R','RHF energy',ERHF)
call dump_test_value('R','RHF HOMO energy',eHF(nO))
call dump_test_value('R','RHF LUMO energy',eHF(nO+1))
call dump_test_value('R','RHF dipole moment',norm2(dipole))
end if
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
end subroutine

20
src/QuAcK/GQuAcK.f90
View File

@ -1,7 +1,7 @@
subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & subroutine GQuAcK(working_dir,dotest,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,doG0T0pp,doevGTpp,doqsGTpp, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_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, & maxSCF_CC,max_diis_CC,thresh_CC, &
TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
@ -12,6 +12,8 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
implicit none implicit none
include 'parameters.h' include 'parameters.h'
character(len=256),intent(in) :: working_dir
logical,intent(in) :: dotest logical,intent(in) :: dotest
logical,intent(in) :: doGHF logical,intent(in) :: doGHF
@ -41,7 +43,6 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
integer,intent(in) :: maxSCF_HF,max_diis_HF integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix double precision,intent(in) :: thresh_HF,level_shift,mix
@ -86,9 +87,11 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
double precision :: start_GW ,end_GW ,t_GW double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT double precision :: start_GT ,end_GT ,t_GT
double precision :: start_int, end_int, t_int
double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:),FHF(:,:) double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:),FHF(:,:)
double precision :: EGHF double precision :: EGHF
double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:)
double precision,allocatable :: ERI_tmp(:,:,:,:) double precision,allocatable :: ERI_tmp(:,:,:,:)
double precision,allocatable :: Ca(:,:),Cb(:,:) double precision,allocatable :: Ca(:,:),Cb(:,:)
@ -112,6 +115,17 @@ subroutine GQuAcK(dotest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
allocate(cHF(nBas2,nBas2),eHF(nBas2),PHF(nBas2,nBas2),FHF(nBas2,nBas2), & allocate(cHF(nBas2,nBas2),eHF(nBas2),PHF(nBas2,nBas2),FHF(nBas2,nBas2), &
dipole_int_MO(nBas2,nBas2,ncart),ERI_MO(nBas2,nBas2,nBas2,nBas2)) dipole_int_MO(nBas2,nBas2,ncart),ERI_MO(nBas2,nBas2,nBas2,nBas2))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds'
write(*,*)
!---------------------! !---------------------!
! Hartree-Fock module ! ! Hartree-Fock module !
!---------------------! !---------------------!

44
src/QuAcK/QuAcK.f90
View File

@ -32,7 +32,6 @@ program QuAcK
double precision,allocatable :: Hc(:,:) double precision,allocatable :: Hc(:,:)
double precision,allocatable :: X(:,:),X_tmp(:,:) double precision,allocatable :: X(:,:),X_tmp(:,:)
double precision,allocatable :: dipole_int_AO(:,:,:) double precision,allocatable :: dipole_int_AO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: Uvec(:,:), Uval(:) double precision,allocatable :: Uvec(:,:), Uval(:)
double precision :: start_QuAcK,end_QuAcK,t_QuAcK double precision :: start_QuAcK,end_QuAcK,t_QuAcK
@ -44,6 +43,7 @@ program QuAcK
logical :: reg_MP logical :: reg_MP
logical :: switch_hpc
logical :: use_gpu logical :: use_gpu
integer :: maxSCF_CC,max_diis_CC integer :: maxSCF_CC,max_diis_CC
@ -140,7 +140,7 @@ program QuAcK
! Hardware ! ! Hardware !
!------------------! !------------------!
call read_hardware(working_dir,use_gpu) call read_hpc_flags(working_dir,switch_hpc,use_gpu)
!------------------------------------! !------------------------------------!
! Read input information ! ! Read input information !
@ -176,20 +176,19 @@ program QuAcK
allocate(T(nBas,nBas)) allocate(T(nBas,nBas))
allocate(V(nBas,nBas)) allocate(V(nBas,nBas))
allocate(Hc(nBas,nBas)) allocate(Hc(nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
allocate(dipole_int_AO(nBas,nBas,ncart)) allocate(dipole_int_AO(nBas,nBas,ncart))
! Read integrals ! Read integrals
call wall_time(start_int) call wall_time(start_int)
call read_integrals(working_dir,nBas,S,T,V,Hc,ERI_AO) call read_1e_integrals(working_dir,nBas,S,T,V,Hc)
call read_dipole_integrals(working_dir,nBas,dipole_int_AO) call read_dipole_integrals(working_dir,nBas,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 wall time for reading integrals = ',t_int,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 1e-integrals = ',t_int,' seconds'
write(*,*) write(*,*)
! Compute orthogonalization matrix ! Compute orthogonalization matrix
@ -225,29 +224,44 @@ program QuAcK
! Restricted QuAcK branch ! ! Restricted QuAcK branch !
!-------------------------! !-------------------------!
if(doRQuAcK) & if(doRQuAcK) then
call RQuAcK(use_gpu,doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
if(switch_hpc) then
call RQuAcK_hpc(working_dir,use_gpu,doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
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, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & nNuc,nBas,nOrb,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,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, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
else
call RQuAcK(working_dir,use_gpu,doRtest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
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, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_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, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
endif
endif
!---------------------------! !---------------------------!
! Unrestricted QuAcK branch ! ! Unrestricted QuAcK branch !
!---------------------------! !---------------------------!
if(doUQuAcK) & if(doUQuAcK) &
call UQuAcK(doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & call UQuAcK(working_dir,doUtest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
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, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
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,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,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, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
@ -257,10 +271,10 @@ program QuAcK
! Generalized QuAcK branch ! ! Generalized QuAcK branch !
!--------------------------! !--------------------------!
if(doGQuAcK) & if(doGQuAcK) &
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & call GQuAcK(working_dir,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,doG0T0pp,doevGTpp,doqsGTpp, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,doG0T0pp,doevGTpp,doqsGTpp, &
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,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_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, &
@ -289,4 +303,10 @@ program QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds' write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds'
write(*,*) write(*,*)
deallocate(S)
deallocate(T)
deallocate(V)
deallocate(Hc)
deallocate(dipole_int_AO)
end program end program

27
src/QuAcK/RQuAcK.f90
View File

@ -1,9 +1,9 @@
subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
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, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & nNuc,nBas,nOrb,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,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,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, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
@ -14,6 +14,8 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,
implicit none implicit none
include 'parameters.h' include 'parameters.h'
character(len=256),intent(in) :: working_dir
logical,intent(in) :: use_gpu logical,intent(in) :: use_gpu
logical,intent(in) :: dotest logical,intent(in) :: dotest
@ -46,7 +48,6 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb) double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
integer,intent(in) :: maxSCF_HF,max_diis_HF integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix double precision,intent(in) :: thresh_HF,level_shift,mix
@ -94,12 +95,14 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,
double precision :: start_GW ,end_GW ,t_GW double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT double precision :: start_GT ,end_GT ,t_GT
double precision :: start_int, end_int, t_int
double precision,allocatable :: eHF(:) double precision,allocatable :: eHF(:)
double precision,allocatable :: cHF(:,:) double precision,allocatable :: cHF(:,:)
double precision,allocatable :: PHF(:,:) double precision,allocatable :: PHF(:,:)
double precision,allocatable :: FHF(:,:) double precision,allocatable :: FHF(:,:)
double precision :: ERHF double precision :: ERHF
double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: ixyz integer :: ixyz
integer :: nS integer :: nS
@ -121,6 +124,15 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,
allocate(dipole_int_MO(nOrb,nOrb,ncart)) allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds'
write(*,*)
!---------------------! !---------------------!
! Hartree-Fock module ! ! Hartree-Fock module !
!---------------------! !---------------------!
@ -350,4 +362,13 @@ subroutine RQuAcK(use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,
end if end if
deallocate(eHF)
deallocate(cHF)
deallocate(PHF)
deallocate(FHF)
deallocate(dipole_int_MO)
deallocate(ERI_MO)
deallocate(ERI_AO)
end subroutine end subroutine

353
src/QuAcK/RQuAcK_hpc.f90 Normal file
View File

@ -0,0 +1,353 @@
subroutine RQuAcK_hpc(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
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, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, &
maxSCF_GF,max_diis_GF,renorm_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_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
! Restricted branch of QuAcK
implicit none
include 'parameters.h'
character(len=256),intent(in) :: working_dir
logical,intent(in) :: use_gpu
logical,intent(in) :: dotest
logical,intent(in) :: doRHF,doROHF
logical,intent(in) :: dostab
logical,intent(in) :: dosearch
logical,intent(in) :: doMP2,doMP3
logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD
logical,intent(in) :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA
logical,intent(in) :: doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3
logical,intent(in) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
integer,intent(in) :: guess_type
logical,intent(in) :: reg_MP
integer,intent(in) :: maxSCF_CC,max_diis_CC
double precision,intent(in) :: thresh_CC
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: TDA
integer,intent(in) :: maxSCF_GF,max_diis_GF,renorm_GF
double precision,intent(in) :: thresh_GF
logical,intent(in) :: lin_GF,reg_GF
double precision,intent(in) :: eta_GF
integer,intent(in) :: maxSCF_GW,max_diis_GW
double precision,intent(in) :: thresh_GW
logical,intent(in) :: TDA_W,lin_GW,reg_GW
double precision,intent(in) :: eta_GW
integer,intent(in) :: maxSCF_GT,max_diis_GT
double precision,intent(in) :: thresh_GT
logical,intent(in) :: TDA_T,lin_GT,reg_GT
double precision,intent(in) :: eta_GT
logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA
logical,intent(in) :: doACFDT,exchange_kernel,doXBS
! Local variables
logical :: doMP,doCC,doCI,doRPA,doGF,doGW,doGT
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_stab ,end_stab ,t_stab
double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO
double precision :: start_MP ,end_MP ,t_MP
double precision :: start_CC ,end_CC ,t_CC
double precision :: start_CI ,end_CI ,t_CI
double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_GF ,end_GF ,t_GF
double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT
double precision :: start_int, end_int, t_int
double precision,allocatable :: eHF(:)
double precision,allocatable :: cHF(:,:)
double precision,allocatable :: PHF(:,:)
double precision,allocatable :: FHF(:,:)
double precision :: ERHF
integer :: ixyz
integer :: nS
write(*,*)
write(*,*) '*****************************************'
write(*,*) '* Restricted Branch of QuAcK (HPC mode) *'
write(*,*) '*****************************************'
write(*,*)
!-------------------!
! Memory allocation !
!-------------------!
allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(PHF(nBas,nBas))
allocate(FHF(nBas,nBas))
!---------------------!
! Hartree-Fock module !
!---------------------!
if(doRHF) then
call wall_time(start_HF)
call RHF_hpc(working_dir,dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*)
end if
! if(doROHF) then
!
! call wall_time(start_HF)
! call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
! call wall_time(end_HF)
!
! t_HF = end_HF - start_HF
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ROHF = ',t_HF,' seconds'
! write(*,*)
!
! end if
!
!!----------------------------------!
!! AO to MO integral transformation !
!!----------------------------------!
!
! call wall_time(start_AOtoMO)
!
! write(*,*)
! write(*,*) 'AO to MO transformation... Please be patient'
! write(*,*)
!
! ! Read and transform dipole-related integrals
!
! do ixyz=1,ncart
! call AOtoMO(nBas,nOrb,cHF,dipole_int_AO(1,1,ixyz),dipole_int_MO(1,1,ixyz))
! end do
!
! ! 4-index transform
!
! call AOtoMO_ERI_RHF(nBas,nOrb,cHF,ERI_AO,ERI_MO)
!
! call wall_time(end_AOtoMO)
!
! t_AOtoMO = end_AOtoMO - start_AOtoMO
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds'
! write(*,*)
!
!!-----------------------------------!
!! Stability analysis of HF solution !
!!-----------------------------------!
!
! nS = (nO - nC)*(nV - nR)
!
! if(dostab) then
!
! call wall_time(start_stab)
! call RHF_stability(nOrb,nC,nO,nV,nR,nS,eHF,ERI_MO)
! call wall_time(end_stab)
!
! t_stab = end_stab - start_stab
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds'
! write(*,*)
!
! end if
!
! if(dosearch) then
!
! call wall_time(start_stab)
! call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nOrb,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X, &
! ERHF,eHF,cHF,PHF,FHF)
! call wall_time(end_stab)
!
! t_stab = end_stab - start_stab
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds'
! write(*,*)
!
! end if
!
!!-----------------------!
!! Moller-Plesset module !
!!-----------------------!
!
! doMP = doMP2 .or. doMP3
!
! if(doMP) then
!
! call wall_time(start_MP)
! call RMP(dotest,doMP2,doMP3,reg_MP,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
! call wall_time(end_MP)
!
! t_MP = end_MP - start_MP
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP = ',t_MP,' seconds'
! write(*,*)
!
! end if
!
!!------------------------!
!! Coupled-cluster module !
!!------------------------!
!
! doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT .or. &
! dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD
!
! if(doCC) then
!
! call wall_time(start_CC)
! call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
! maxSCF_CC,thresh_CC,max_diis_CC,nBas,nOrb,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO, &
! ENuc,ERHF,eHF,cHF,PHF,FHF)
! call wall_time(end_CC)
!
! t_CC = end_CC - start_CC
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CC = ',t_CC,' seconds'
! write(*,*)
!
! end if
!
!!----------------------------------!
!! Configuration interaction module !
!!----------------------------------!
!
! doCI = doCIS .or. doCID .or. doCISD .or. doFCI
!
! if(doCI) then
!
! call wall_time(start_CI)
! call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nOrb, &
! nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,ERHF)
! call wall_time(end_CI)
!
! t_CI = end_CI - start_CI
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CI = ',t_CI,' seconds'
! write(*,*)
!
! end if
!
!!-----------------------------------!
!! Random-phase approximation module !
!!-----------------------------------!
!
! doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA
!
! if(doRPA) then
!
! call wall_time(start_RPA)
! call RRPA(use_gpu,dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
! nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
! call wall_time(end_RPA)
!
! t_RPA = end_RPA - start_RPA
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds'
! write(*,*)
!
! end if
!
!!-------------------------!
!! Green's function module !
!!-------------------------!
!
! doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3
!
! if(doGF) then
!
! call wall_time(start_GF)
! call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF, &
! thresh_GF,max_diis_GF,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF, &
! eta_GF,reg_GF,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, &
! S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! call wall_time(end_GF)
!
! t_GF = end_GF - start_GF
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
! write(*,*)
!
! end if
!
!!-----------!
!! GW module !
!!-----------!
!
! doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW
!
! if(doGW) then
!
! call wall_time(start_GW)
! call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF_GW,thresh_GW,max_diis_GW, &
! doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
! lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, &
! V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! call wall_time(end_GW)
!
! t_GW = end_GW - start_GW
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GW = ',t_GW,' seconds'
! write(*,*)
!
! end if
!
!!-----------------!
!! T-matrix module !
!!-----------------!
!
! doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp .or. doufG0T0pp .or. doG0T0eh .or. doevGTeh .or. doqsGTeh
!
! if(doGT) then
!
! call wall_time(start_GT)
! call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
! maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, &
! TDA_T,TDA,dBSE,dTDA,singlet,triplet,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO, &
! dipole_int_MO,PHF,cHF,eHF)
! call wall_time(end_GT)
!
! t_GT = end_GT - start_GT
! write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GT = ',t_GT,' seconds'
! write(*,*)
!
! end if
return
end subroutine

18
src/QuAcK/UQuAcK.f90
View File

@ -1,9 +1,9 @@
subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & subroutine UQuAcK(working_dir,dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
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, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW, &
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,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,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, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, & TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
@ -12,6 +12,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
implicit none implicit none
include 'parameters.h' include 'parameters.h'
character(len=256),intent(in) :: working_dir
logical,intent(in) :: dotest logical,intent(in) :: dotest
logical,intent(in) :: doUHF logical,intent(in) :: doUHF
@ -42,7 +44,6 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
integer,intent(in) :: maxSCF_HF,max_diis_HF integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix double precision,intent(in) :: thresh_HF,level_shift,mix
@ -90,10 +91,12 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
double precision :: start_GW ,end_GW ,t_GW double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT double precision :: start_GT ,end_GT ,t_GT
double precision :: start_int, end_int, t_int
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:),FHF(:,:,:) double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:),FHF(:,:,:)
double precision :: EUHF double precision :: EUHF
double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:) double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:)
double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:) double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
integer :: ixyz integer :: ixyz
integer :: nS(nspin) integer :: nS(nspin)
@ -112,6 +115,15 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas), & ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas), &
ERI_bbbb(nBas,nBas,nBas,nBas)) ERI_bbbb(nBas,nBas,nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds'
write(*,*)
!---------------------! !---------------------!
! Hartree-Fock module ! ! Hartree-Fock module !
!---------------------! !---------------------!

45
src/QuAcK/read_hardware.f90
View File

@ -1,45 +0,0 @@
subroutine read_hardware(working_dir,use_gpu)
! Read desired methods
implicit none
! Input variables
character(len=256),intent(in) :: working_dir
! Output variables
logical,intent(out) :: use_gpu
! Local variables
character(len=1) :: ans
integer :: ios
character(len=256) :: file_path
! Open file with method specification
file_path = trim(working_dir) // '/input/hardware'
open(unit=1, file=file_path, status='old', action='read', iostat=ios)
if(ios /= 0) then
use_gpu = .False.
else
read(1,*)
read(1,*) ans
if(ans == 'T') then
use_gpu = .true.
else
use_gpu = .False.
endif
endif
! Close file with options
close(unit=1)
end subroutine

39
src/QuAcK/read_hpc_flags.f90 Normal file
View File

@ -0,0 +1,39 @@
subroutine read_hpc_flags(working_dir, switch_hpc, use_gpu)
implicit none
character(len=256), intent(in) :: working_dir
logical, intent(out) :: switch_hpc
logical, intent(out) :: use_gpu
character(len=1) :: ans
integer :: ios
character(len=256) :: file_path
file_path = trim(working_dir) // '/input/hpc_flags'
open(unit=1, file=file_path, status='old', action='read', iostat=ios)
if(ios /= 0) then
switch_hpc = .False.
use_gpu = .False.
else
switch_hpc = .False.
use_gpu = .False.
read(1,*)
read(1,*) ans
if(ans == 'T') switch_hpc = .true.
read(1,*)
read(1,*) ans
if(ans == 'T') use_gpu = .true.
endif
close(unit=1)
end subroutine

68
src/utils/read_integrals.f90 → src/utils/read_1e_integrals.f90
View File

@ -1,6 +1,6 @@
subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G) subroutine read_1e_integrals(working_dir,nBas_AOs,S,T,V,Hc)
! Read one- and two-electron integrals from files ! Read one-electron integrals from files
implicit none implicit none
include 'parameters.h' include 'parameters.h'
@ -13,9 +13,8 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
! Local variables ! Local variables
logical :: debug logical :: debug
integer :: mu,nu,la,si integer :: mu,nu
double precision :: Ov,Kin,Nuc,ERI double precision :: Ov,Kin,Nuc
double precision :: lambda
! Output variables ! Output variables
@ -23,26 +22,21 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
double precision,intent(out) :: T(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) :: V(nBas_AOs,nBas_AOs)
double precision,intent(out) :: Hc(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)
integer :: status, ios integer :: ios
character(len=256) :: file_path character(len=256) :: file_path
! Open file with integrals ! Open file with integrals
debug = .false. debug = .false.
lambda = 1d0
print*, 'Scaling integrals by ',lambda
! --- ! ---
! Read overlap integrals ! Read overlap integrals
file_path = trim(working_dir) // '/int/Ov.dat' file_path = trim(working_dir) // '/int/Ov.dat'
open(unit=8, file=file_path, status='old', action='read', iostat=status) open(unit=8, file=file_path, status='old', action='read', iostat=ios)
if(status /= 0) then if(ios /= 0) then
print *, "Error opening file: ", file_path print *, "Error opening file: ", file_path
stop stop
else else
@ -60,8 +54,8 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
! Read kinetic integrals ! Read kinetic integrals
file_path = trim(working_dir) // '/int/Kin.dat' file_path = trim(working_dir) // '/int/Kin.dat'
open(unit=9, file=file_path, status='old', action='read', iostat=status) open(unit=9, file=file_path, status='old', action='read', iostat=ios)
if(status /= 0) then if(ios /= 0) then
print *, "Error opening file: ", file_path print *, "Error opening file: ", file_path
stop stop
else else
@ -79,8 +73,8 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
! Read nuclear integrals ! Read nuclear integrals
file_path = trim(working_dir) // '/int/Nuc.dat' file_path = trim(working_dir) // '/int/Nuc.dat'
open(unit=10, file=file_path, status='old', action='read', iostat=status) open(unit=10, file=file_path, status='old', action='read', iostat=ios)
if(status /= 0) then if(ios /= 0) then
print *, "Error opening file: ", file_path print *, "Error opening file: ", file_path
stop stop
else else
@ -99,37 +93,6 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
! Define core Hamiltonian ! Define core Hamiltonian
Hc(:,:) = T(:,:) + V(:,:) Hc(:,:) = T(:,:) + V(:,:)
! Read 2e-integrals
! ! 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)
! binary file
file_path = trim(working_dir) // '/int/ERI.bin'
open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=status)
if(status /= 0) then
print *, "Error opening file: ", file_path
stop
else
read(11) G
endif
close(unit=11)
! Print results ! Print results
if(debug) then if(debug) then
@ -148,15 +111,6 @@ subroutine read_integrals(working_dir,nBas_AOs,S,T,V,Hc,G)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
call matout(nBas_AOs,nBas_AOs,V) call matout(nBas_AOs,nBas_AOs,V)
write(*,*) write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
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(*,*)
end if end if
end subroutine end subroutine

110
src/utils/read_2e_integrals.f90 Normal file
View File

@ -0,0 +1,110 @@
subroutine read_2e_integrals(working_dir,nBas_AOs,G)
! Read two-electron integrals from files
implicit none
! Input variables
integer,intent(in) :: nBas_AOs
character(len=256),intent(in) :: working_dir
! Local variables
logical :: debug
integer :: mu,nu,la,si
double precision :: ERI
double precision :: lambda
! Output variables
double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
integer :: ios
character(len=256) :: file_path
! Open file with integrals
debug = .false.
lambda = 1d0
print*, 'Scaling integrals by ',lambda
! Read 2e-integrals
! ! 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)
! binary file
file_path = trim(working_dir) // '/int/ERI.bin'
open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=ios)
if(ios /= 0) then
print *, "Error opening file: ", file_path
stop
else
read(11) G
endif
close(unit=11)
G = G * lambda
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
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(*,*)
end if
end subroutine
! ---
subroutine read_2e_integrals_hpc(working_dir, ERI_size, ERI_chem)
implicit none
character(len=256), intent(in) :: working_dir
integer*8, intent(in) :: ERI_size
double precision, intent(out) :: ERI_chem(ERI_size)
integer :: ios
character(len=256) :: file_path
file_path = trim(working_dir) // '/int/ERI_chem.bin'
open(unit=11, file=file_path, status='old', action='read', form='unformatted', access='stream', iostat=ios)
if(ios /= 0) then
print *, "Error opening file: ", file_path
stop
else
read(11) ERI_chem
endif
close(unit=11)
return
end subroutine
! ---