LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-10-20 06:48:21 +02:00
Code Releases Activity

passing Fock matrix and others in U branch

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-02 22:35:46 +02:00
parent fda6a45276
commit a90e6203ec
5 changed files with 51 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
src/HF/RHF_search.f90
View File

@ -1,9 +1,6 @@
subroutine RHF_search(maxSCF,thresh,max_diis,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,e,c,P,F)
subroutine RHF_search(maxSCF, thresh, max_diis, 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, e, c, P)
! Search for RHF solutions ! Search for RHF solutions
@ -65,6 +62,7 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
double precision,intent(out) :: e(nOrb) double precision,intent(out) :: e(nOrb)
double precision,intent(inout):: c(nBas,nOrb) double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
! Memory allocation ! Memory allocation
@ -79,8 +77,8 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
!-------------------! !-------------------!
nS = (nO - nC)*(nV - nR) nS = (nO - nC)*(nV - nR)
allocate(Aph(nS,nS), Bph(nS,nS), AB(nS,nS), Om(nS)) allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS))
allocate(R(nOrb,nOrb), ExpR(nOrb,nOrb)) allocate(R(nOrb,nOrb),ExpR(nOrb,nOrb))
!------------------! !------------------!
! Search algorithm ! ! Search algorithm !
@ -96,8 +94,8 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
!---------------------! !---------------------!
call wall_time(start_HF) call wall_time(start_HF)
call RHF(.false., maxSCF, thresh, max_diis, guess, level_shift, nNuc, ZNuc, rNuc, ENuc, & call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc,&
nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, e, c, P) nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,e,c,P,F)
call wall_time(end_HF) call wall_time(end_HF)
t_HF = end_HF - start_HF t_HF = end_HF - start_HF
@ -112,10 +110,10 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
write(*,*) write(*,*)
write(*,*) 'AO to MO transformation... Please be patient' write(*,*) 'AO to MO transformation... Please be patient'
write(*,*) write(*,*)
do ixyz = 1, ncart do ixyz = 1,ncart
call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) call AOtoMO(nBas,nOrb,c,dipole_int_AO(1,1,ixyz),dipole_int_MO(1,1,ixyz))
end do end do
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) call AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
call wall_time(end_AOtoMO) call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO t_AOtoMO = end_AOtoMO - start_AOtoMO
@ -133,7 +131,7 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
AB(:,:) = Aph(:,:) + Bph(:,:) AB(:,:) = Aph(:,:) + Bph(:,:)
call diagonalize_matrix(nS, AB, Om) call diagonalize_matrix(nS,AB,Om)
Om(:) = 2d0*Om(:) Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------' write(*,*)'-------------------------------------------------------------'
@ -160,7 +158,7 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
if(eig < 0 .or. eig > nS) then if(eig < 0 .or. eig > nS) then
write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...' write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...'
write(*,*) write(*,*)
deallocate(Aph, Bph, AB, Om, R, ExpR) deallocate(Aph,Bph,AB,Om,R,ExpR)
stop stop
end if end if
@ -176,8 +174,8 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
end do end do
end do end do
call matrix_exponential(nOrb, R, ExpR) call matrix_exponential(nOrb,R,ExpR)
c = matmul(c, ExpR) c = matmul(c,ExpR)
else else
@ -196,6 +194,6 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z
!---------------! !---------------!
end do end do
deallocate(Aph, Bph, AB, Om, R, ExpR) deallocate(Aph,Bph,AB,Om,R,ExpR)
end subroutine end subroutine

8
src/HF/UHF.f90
View File

@ -1,5 +1,5 @@
subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,eHF,c,P) nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,eHF,c,P,F)
! Perform unrestricted Hartree-Fock calculation ! Perform unrestricted Hartree-Fock calculation
@ -47,7 +47,6 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,allocatable :: cp(:,:,:) double precision,allocatable :: cp(:,:,:)
double precision,allocatable :: J(:,:,:) double precision,allocatable :: J(:,:,:)
double precision,allocatable :: F(:,:,:)
double precision,allocatable :: Fp(:,:,:) double precision,allocatable :: Fp(:,:,:)
double precision,allocatable :: K(:,:,:) double precision,allocatable :: K(:,:,:)
double precision,allocatable :: err(:,:,:) double precision,allocatable :: err(:,:,:)
@ -63,6 +62,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,intent(out) :: eHF(nBas,nspin) double precision,intent(out) :: eHF(nBas,nspin)
double precision,intent(inout):: c(nBas,nBas,nspin) double precision,intent(inout):: c(nBas,nBas,nspin)
double precision,intent(out) :: P(nBas,nBas,nspin) double precision,intent(out) :: P(nBas,nBas,nspin)
double precision,intent(out) :: F(nBas,nBas,nspin)
! Hello world ! Hello world
@ -78,8 +78,8 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Memory allocation ! Memory allocation
allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & allocate(J(nBas,nBas,nspin),K(nBas,nBas,nspin),Fp(nBas,nBas,nspin), &
K(nBas,nBas,nspin),err(nBas,nBas,nspin),cp(nBas,nBas,nspin), & err(nBas,nBas,nspin),cp(nBas,nBas,nspin), &
err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
! Guess coefficients and demsity matrices ! Guess coefficients and demsity matrices

5
src/HF/UHF_search.f90
View File

@ -1,6 +1,6 @@
subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_AO,dipole_int_aa,dipole_int_bb,X,EUHF,e,c,P) dipole_int_AO,dipole_int_aa,dipole_int_bb,X,EUHF,e,c,P,F)
! Search for UHF solutions ! Search for UHF solutions
@ -68,6 +68,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,intent(out) :: e(nBas,nspin) double precision,intent(out) :: e(nBas,nspin)
double precision,intent(inout):: c(nBas,nBas,nspin) double precision,intent(inout):: c(nBas,nBas,nspin)
double precision,intent(out) :: P(nBas,nBas,nspin) double precision,intent(out) :: P(nBas,nBas,nspin)
double precision,intent(out) :: F(nBas,nBas,nspin)
! Memory allocation ! Memory allocation
@ -105,7 +106,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
call wall_time(start_HF) call wall_time(start_HF)
call UHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & call UHF(.false.,maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,e,c,P) nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,e,c,P,F)
call wall_time(end_HF) call wall_time(end_HF)
t_HF = end_HF - start_HF t_HF = end_HF - start_HF

44
src/QuAcK/RQuAcK.f90
View File

@ -184,7 +184,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(dostab) then if(dostab) then
call wall_time(start_stab) call wall_time(start_stab)
call RHF_stability(nOrb, nC, nO, nV, nR, nS, eHF, ERI_MO) call RHF_stability(nOrb,nC,nO,nV,nR,nS,eHF,ERI_MO)
call wall_time(end_stab) call wall_time(end_stab)
t_stab = end_stab - start_stab t_stab = end_stab - start_stab
@ -196,9 +196,9 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(dosearch) then if(dosearch) then
call wall_time(start_stab) call wall_time(start_stab)
call RHF_search(maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & 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, & nBas,nOrb,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X, &
dipole_int_MO, X, ERHF, eHF, cHF, PHF) ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_stab) call wall_time(end_stab)
t_stab = end_stab - start_stab t_stab = end_stab - start_stab
@ -216,7 +216,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doMP) then if(doMP) then
call wall_time(start_MP) call wall_time(start_MP)
call RMP(dotest, doMP2, doMP3, reg_MP, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, ERHF, eHF) call RMP(dotest,doMP2,doMP3,reg_MP,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_MP) call wall_time(end_MP)
t_MP = end_MP - start_MP t_MP = end_MP - start_MP
@ -255,8 +255,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doCI) then if(doCI) then
call wall_time(start_CI) call wall_time(start_CI)
call RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, & call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nOrb, &
nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, ERHF) nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,ERHF)
call wall_time(end_CI) call wall_time(end_CI)
t_CI = end_CI - start_CI t_CI = end_CI - start_CI
@ -274,8 +274,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRPA) then if(doRPA) then
call wall_time(start_RPA) call wall_time(start_RPA)
call RRPA(dotest, dophRPA, dophRPAx, docrRPA, doppRPA, TDA, doACFDT, exchange_kernel, singlet, triplet, & call RRPA(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) nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
@ -293,10 +293,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGF) then if(doGF) then
call wall_time(start_GF) call wall_time(start_GF)
call RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm_GF, maxSCF_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, & 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, & 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) S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
@ -314,10 +314,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGW) then if(doGW) then
call wall_time(start_GW) call wall_time(start_GW)
call RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF_GW, thresh_GW, max_diis_GW, & call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, &
doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, & 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, & 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) V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW
@ -335,11 +335,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGT) then if(doGT) then
call wall_time(start_GT) call wall_time(start_GT)
call RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
maxSCF_GT, thresh_GT, max_diis_GT, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, & 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, & 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, & 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) dipole_int_MO,PHF,cHF,eHF)
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT

12
src/QuAcK/UQuAcK.f90
View File

@ -90,7 +90,7 @@ 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,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:) 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(:,:,:,:)
@ -107,9 +107,9 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Memory allocation ! ! Memory allocation !
!-------------------! !-------------------!
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),PHF(nBas,nBas,nspin), & allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),PHF(nBas,nBas,nspin),FHF(nBas,nBas,nspin), &
dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart), & dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart), &
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))
!---------------------! !---------------------!
@ -120,7 +120,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
call wall_time(start_HF) call wall_time(start_HF)
call UHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & call UHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF) nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF) call wall_time(end_HF)
t_HF = end_HF - start_HF t_HF = end_HF - start_HF
@ -208,7 +208,7 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
call wall_time(start_stab) call wall_time(start_stab)
call UHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & call UHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO, & nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO, &
dipole_int_aa,dipole_int_bb,X,EUHF,eHF,cHF,PHF) dipole_int_aa,dipole_int_bb,X,EUHF,eHF,cHF,PHF,FHF)
call wall_time(end_stab) call wall_time(end_stab)