From 3b958f56107ef9905c90f8f02692e04fcacafdb9 Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 3 Nov 2023 19:48:12 +0100 Subject: [PATCH] HF solution search --- input/options | 4 ++-- src/HF/RHF_search.f90 | 43 ++++++++++++++++++++++++++++++-------- src/HF/mo_guess.f90 | 3 ++- src/QuAcK/GQuAcK.f90 | 12 ++++++----- src/QuAcK/QuAcK.f90 | 33 +++++++++++++++-------------- src/QuAcK/RQuAcK.f90 | 5 +++-- src/QuAcK/UQuAcK.f90 | 3 ++- src/QuAcK/read_options.f90 | 23 +++++++++++--------- 8 files changed, 80 insertions(+), 46 deletions(-) diff --git a/input/options b/input/options index 7c56e8a..193f265 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index bc94d6a..bf7e5ec 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -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!' diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index 81c4cbc..6637651 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -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 diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 7ce71cc..b4bac6c 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 7401992..647a99a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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) !--------------! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index cc6bbc1..e42afbb 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -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, & diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index 303d5a1..02e1b61 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -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 diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 246514b..b9cdae1 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -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