4
1
mirror of https://github.com/pfloos/quack synced 2024-11-19 04:22:39 +01:00

HF solution search

This commit is contained in:
Pierre-Francois Loos 2023-11-03 19:48:12 +01:00
parent 17ce993aa2
commit 3b958f5610
8 changed files with 80 additions and 46 deletions

View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS guess mix_guess level_shift stability
1000 0.00001 5 1 0.0 0.0 T
# HF: maxSCF thresh DIIS guess mix_guess level_shift stab search
1000 0.00001 5 1 0.0 0.0 T T
# MP: reg
F
# CC: maxSCF thresh DIIS

View File

@ -41,13 +41,14 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
integer :: nS
integer,parameter :: maxS = 20
integer :: ia
integer :: ia,i,a
integer :: ispin
double precision,allocatable :: A(:,:)
double precision,allocatable :: B(:,:)
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: AB(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: R(:,:)
integer :: eig
double precision :: eigval
@ -74,7 +75,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
nS = (nO - nC)*(nV - nR)
allocate(ERI_MO(nBas,nBas,nBas,nBas),eigvec(nS))
allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS))
allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas,nBas))
!------------------!
! Search algorithm !
@ -119,10 +120,10 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
ispin = 1
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO,A)
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,B)
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
AB(:,:) = A(:,:) + B(:,:)
AB(:,:) = Aph(:,:) + Bph(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
@ -148,8 +149,32 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
read(*,*) eig
eigval = Om(eig)
eigvec = AB(:,eig)
eigvec(:) = AB(:,eig)
print*,eigval
call matout(nS,1,eigvec)
R(:,:) = 0d0
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
R(i,a) = +eigvec(ia)
R(a,i) = -eigvec(ia)
end do
end do
print*,'rotation matrix'
call matout(nBas,nBas,R)
print*,'old coefficients'
call matout(nBas,nBas,c)
c = c - 0.1d0*matmul(c,R)
print*,'new coefficients'
call matout(nBas,nBas,c)
else
write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!'

View File

@ -14,11 +14,12 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(inout) :: c(nBas,nBas)
if(guess_type == 0) then
write(*,*) 'Reading HF coefficients...'
c(:,:) = c(:,:)
elseif(guess_type == 1) then

View File

@ -1,8 +1,9 @@
subroutine GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
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,guess_type,mix,reg_MP, &
spin_conserved,spin_flip,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, &
subroutine GQuAcK(doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
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,guess_type,mix,reg_MP, &
spin_conserved,spin_flip,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)
implicit none
@ -10,6 +11,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevG
logical,intent(in) :: doGHF
logical,intent(in) :: dostab
logical,intent(in) :: dosearch
logical,intent(in) :: doMP2
logical,intent(in) :: doMP3
logical,intent(in) :: dophRPA,dophRPAx,doppRPA

View File

@ -5,7 +5,7 @@ program QuAcK
logical :: doRQuAcK,doUQuAcK,doGQuAcK
logical :: doRHF,doUHF,doGHF,doROHF
logical :: dostab
logical :: dostab,dosearch
logical :: doMP2,doMP3
logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical :: dodrCCD,dorCCD,docrCCD,dolCCD
@ -107,14 +107,14 @@ program QuAcK
! Read options for methods !
!--------------------------!
call read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab, &
reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, &
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
call read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, &
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
!------------------------------------------------!
@ -185,7 +185,7 @@ program QuAcK
!-------------------------!
if(doRQuAcK) &
call RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
call RQuAcK(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,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
@ -200,7 +200,7 @@ program QuAcK
!---------------------------!
if(doUQuAcK) &
call UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
call UQuAcK(doUHF,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,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
@ -215,11 +215,12 @@ program QuAcK
!--------------------------!
if(doGQuAcK) &
call GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,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, &
spin_conserved,spin_flip,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(doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,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, &
spin_conserved,spin_flip,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)
!--------------!

View File

@ -1,4 +1,4 @@
subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
subroutine RQuAcK(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,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
@ -15,6 +15,7 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC
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
@ -182,7 +183,7 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC
end if
if(dostab) then
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, &

View File

@ -1,4 +1,4 @@
subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
subroutine UQuAcK(doUHF,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,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
@ -13,6 +13,7 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,
logical,intent(in) :: doUHF
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

View File

@ -1,11 +1,11 @@
subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab, &
reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, &
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab,dosearch, &
reg_MP, &
maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, &
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
! Read desired methods
@ -21,6 +21,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
double precision,intent(out) :: mix
double precision,intent(out) :: level_shift
logical,intent(out) :: dostab
logical,intent(out) :: dosearch
logical,intent(out) :: reg_MP
@ -83,11 +84,13 @@ subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shi
mix = 0d0
level_shift = 0d0
dostab = .false.
dosearch = .false.
read(1,*)
read(1,*) maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,ans1
read(1,*) maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,ans1,ans2
if(ans1 == 'T') dostab = .true.
if(ans1 == 'T') dostab = .true.
if(ans2 == 'T') dosearch = .true.
! Read MPn options