From 7847710059ec3b6940ad4ff16aedbe7fc4878684 Mon Sep 17 00:00:00 2001 From: pfloos Date: Sat, 4 Nov 2023 15:17:39 +0100 Subject: [PATCH] GHF search --- input/methods | 2 +- input/options | 2 +- src/HF/GHF.f90 | 2 +- src/HF/GHF_search.f90 | 210 +++++++++++++++++++++++++++++++++++++++ src/HF/RHF_search.f90 | 18 +++- src/HF/UHF.f90 | 4 +- src/HF/UHF_search.f90 | 223 ++++++++++++++++++++++++++++++++++++++++++ src/QuAcK/GQuAcK.f90 | 13 +++ src/QuAcK/UQuAcK.f90 | 14 +++ 9 files changed, 479 insertions(+), 9 deletions(-) create mode 100644 src/HF/GHF_search.f90 create mode 100644 src/HF/UHF_search.f90 diff --git a/input/methods b/input/methods index 41e9098..faa85b7 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF GHF ROHF - T F F F + F T F F # MP2 MP3 F F # CCD pCCD DCD CCSD CCSD(T) diff --git a/input/options b/input/options index e30e10a..fa0c3ce 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS guess mix shift stab search - 1000 0.0000001 5 1 0.0 0.0 F T + 1000 0.00001 5 3 0.0 0.0 F T # MP: reg F # CC: maxSCF thresh DIIS diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index 19e0819..bef529e 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -66,7 +66,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision,intent(out) :: EHF double precision,intent(out) :: e(nBas2) - double precision,intent(out) :: C(nBas2,nBas2) + double precision,intent(inout):: C(nBas2,nBas2) double precision,intent(out) :: P(nBas2,nBas2) ! Hello world diff --git a/src/HF/GHF_search.f90 b/src/HF/GHF_search.f90 new file mode 100644 index 0000000..f342da1 --- /dev/null +++ b/src/HF/GHF_search.f90 @@ -0,0 +1,210 @@ +subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) + +! Search for GHF solutions + + implicit none + include 'parameters.h' + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + integer,intent(in) :: guess_type + double precision,intent(in) :: thresh + double precision,intent(inout):: mix + double precision,intent(in) :: level_shift + + integer,intent(in) :: nBas + integer,intent(in) :: nBas2 + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + 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,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + +! Local variables + + double precision :: start_HF ,end_HF ,t_HF + double precision :: start_stab ,end_stab ,t_stab + double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO + + logical :: unstab + integer :: guess + double precision,allocatable :: ERI_MO(:,:,:,:) + double precision,allocatable :: ERI_tmp(:,:,:,:) + double precision,allocatable :: Ca(:,:),Cb(:,:) + integer :: nS + + integer,parameter :: maxS = 20 + integer :: ia,i,a,mu + integer :: ispin + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: AB(:,:) + double precision,allocatable :: Om(:) + + integer :: eig + double precision :: kick,step + +! Output variables + + double precision,intent(out) :: EHF + double precision,intent(out) :: e(nBas) + double precision,intent(inout):: c(nBas,nBas) + double precision,intent(out) :: P(nBas,nBas) + +! Memory allocation + + write(*,*) + write(*,*) '****************************' + write(*,*) '* Search for GHF solutions *' + write(*,*) '****************************' + write(*,*) + +!-------------------! +! Memory allocation ! +!-------------------! + + nS = (nO - nC)*(nV - nR) + + allocate(ERI_MO(nBas2,nBas2,nBas2,nBas2),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS)) + +!------------------! +! Search algorithm ! +!------------------! + + unstab = .true. + guess = 0 + + do while(unstab) + +!---------------------! +! Hartree-Fock module ! +!---------------------! + + call wall_time(start_HF) + call GHF(maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nBas2,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) + call wall_time(end_HF) + + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GHF = ',t_HF,' seconds' + write(*,*) + +!----------------------------------! +! AO to MO integral transformation ! +!----------------------------------! + + call wall_time(start_AOtoMO) + write(*,*) + write(*,*) 'AO to MO transformation... Please be patient' + write(*,*) + + allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2)) + + Ca(:,:) = c(1:nBas,1:nBas2) + Cb(:,:) = c(nBas+1:nBas2,1:nBas2) + + ! 4-index transform + + call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Ca,Ca,Ca,ERI_AO,ERI_tmp) + ERI_MO(:,:,:,:) = ERI_tmp(:,:,:,:) + + call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Cb,Ca,Cb,ERI_AO,ERI_tmp) + ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) + + call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Ca,Cb,Ca,ERI_AO,ERI_tmp) + ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) + + call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp) + ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) + + deallocate(Ca,Cb,ERI_tmp) + + t_AOtoMO = end_AOtoMO - start_AOtoMO + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,*) + +!-------------------------------------------------------------! +! Stability analysis: Real GHF -> Real GHF +!-------------------------------------------------------------! + + ispin = 3 + + 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(:,:) = Aph(:,:) + Bph(:,:) + + call diagonalize_matrix(nS,AB,Om) + Om(:) = 0.5d0*Om(:) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real GHF -> Real GHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, GHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue:',Om(1),' au' + write(*,*) + write(*,'(1X,A40,1X,A10)') 'Which one would you like to follow?','[Exit:0]' + read(*,*) eig + + if(eig < 0 .or. eig > nS) then + write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...' + write(*,*) + stop + end if + + if(eig == 0) return + + step = 1d0 + + do mu=1,nBas + ia = 0 + do i=nC+1,nO + kick = 0d0 + do a=nO+1,nBas-nR + ia = ia + 1 + kick = kick + AB(ia,eig)*c(mu,a) + end do + c(mu,i) = c(mu,i) + step*kick + end do + end do + + else + + write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au' + + unstab = .false. + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +!---------------! +! End of Search ! +!---------------! + end do + +end subroutine diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index 968210c..938fb73 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -1,7 +1,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) -! Restricted branch of QuAcK +! Search for RHF solutions implicit none include 'parameters.h' @@ -48,7 +48,6 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,allocatable :: Bph(:,:) double precision,allocatable :: AB(:,:) double precision,allocatable :: Om(:) - double precision,allocatable :: R(:,:) integer :: eig double precision :: kick,step @@ -73,7 +72,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),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas,nBas)) + allocate(ERI_MO(nBas,nBas,nBas,nBas),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS)) !------------------! ! Search algorithm ! @@ -143,9 +142,17 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,'(1X,A40,1X)') 'Too bad, RHF solution is unstable!' write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue:',Om(1),' au' write(*,*) - write(*,'(1X,A40,1X)') 'Which one would you like to follow?' + write(*,'(1X,A40,1X,A10)') 'Which one would you like to follow?','[Exit:0]' read(*,*) eig + if(eig < 0 .or. eig > nS) then + write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...' + write(*,*) + stop + end if + + if(eig == 0) return + step = 1d0 do mu=1,nBas @@ -171,6 +178,9 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,*)'-------------------------------------------------------------' write(*,*) +!---------------! +! End of Search ! +!---------------! end do end subroutine diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 8435170..9899b9d 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -59,7 +59,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision,intent(out) :: EHF double precision,intent(out) :: e(nBas,nspin) - double precision,intent(out) :: c(nBas,nBas,nspin) + double precision,intent(inout):: c(nBas,nBas,nspin) double precision,intent(out) :: P(nBas,nBas,nspin) ! Hello world @@ -181,7 +181,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Mix guess for UHF solution in singlet states - if(nSCF == 1) call mix_guess(nBas,nO,mix,c) + if(nSCF == 1 .and. mix > 0d0) call mix_guess(nBas,nO,mix,c) ! Compute density matrix diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 new file mode 100644 index 0000000..57477ca --- /dev/null +++ b/src/HF/UHF_search.f90 @@ -0,0 +1,223 @@ +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,dipole_int_AO,X,EHF,e,c,P) + +! Search for UHF solutions + + implicit none + include 'parameters.h' + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + integer,intent(in) :: guess_type + double precision,intent(in) :: thresh + double precision,intent(inout):: mix + double precision,intent(in) :: level_shift + + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + 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,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + +! Local variables + + double precision :: start_HF ,end_HF ,t_HF + double precision :: start_stab ,end_stab ,t_stab + double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO + + logical :: unstab + integer :: guess + double precision,allocatable :: ERI_aaaa(:,:,:,:) + double precision,allocatable :: ERI_aabb(:,:,:,:) + double precision,allocatable :: ERI_bbbb(:,:,:,:) + + integer :: nS(nspin) + + integer,parameter :: maxS = 20 + integer :: ia,i,a,mu + integer :: ispin + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Om_sc(:) + double precision,allocatable :: A_sc(:,:) + double precision,allocatable :: B_sc(:,:) + double precision,allocatable :: AB_sc(:,:) + + integer :: eig + double precision :: kick,step + +! Output variables + + double precision,intent(out) :: EHF + double precision,intent(out) :: e(nBas,nspin) + double precision,intent(inout):: c(nBas,nBas,nspin) + double precision,intent(out) :: P(nBas,nBas,nspin) + +! Memory allocation + + write(*,*) + write(*,*) '****************************' + write(*,*) '* Search for UHF solutions *' + write(*,*) '****************************' + write(*,*) + +!-------------------! +! Memory allocation ! +!-------------------! + + nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:)) + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas),ERI_bbbb(nBas,nBas,nBas,nBas)) + allocate(Om_sc(nS_sc),A_sc(nS_sc,nS_sc),B_sc(nS_sc,nS_sc),AB_sc(nS_sc,nS_sc)) + +!------------------! +! Search algorithm ! +!------------------! + + unstab = .true. + guess = 0 + mix = 0d0 + + do while(unstab) + +!---------------------! +! Hartree-Fock module ! +!---------------------! + + call wall_time(start_HF) + call UHF(maxSCF,thresh,max_diis,guess,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,e,c,P) + call wall_time(end_HF) + + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds' + write(*,*) + +!----------------------------------! +! AO to MO integral transformation ! +!----------------------------------! + + call wall_time(start_AOtoMO) + write(*,*) + write(*,*) 'AO to MO transformation... Please be patient' + write(*,*) + + ! 4-index transform for (aa|aa) block + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_aaaa) + + ! 4-index transform for (aa|bb) block + call AOtoMO_integral_transform(1,1,2,2,nBas,c,ERI_AO,ERI_aabb) + + ! 4-index transform for (bb|bb) block + call AOtoMO_integral_transform(2,2,2,2,nBas,c,ERI_AO,ERI_bbbb) + + call wall_time(end_AOtoMO) + + t_AOtoMO = end_AOtoMO - start_AOtoMO + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,*) + +!-------------------------------------------------------------! +! Stability analysis: Real UHF -> Real UHF +!-------------------------------------------------------------! + + ispin = 1 + + call phULR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc) + call phULR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,ERI_aaaa,ERI_aabb,ERI_bbbb,B_sc) + + AB_sc(:,:) = A_sc(:,:) + B_sc(:,:) + + call diagonalize_matrix(nS_sc,AB_sc,Om_sc) + Om_sc(:) = 0.5d0*Om_sc(:) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real UHF -> Real UHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS_sc,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om_sc(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, UHF solution is unstable!' + write(*,'(1X,A40,1X,F10.6,A3)') 'Largest negative eigenvalue:',Om_sc(1),' au' + write(*,*) + write(*,'(1X,A40,1X,A10)') 'Which one would you like to follow?','[Exit:0]' + read(*,*) eig + + if(eig < 0 .or. eig > nS_sc) then + write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...' + write(*,*) + stop + end if + + if(eig == 0) return + + step = 1d0 + + ! Spin-up kick + + do mu=1,nBas + ia = 0 + do i=nC(1)+1,nO(1) + kick = 0d0 + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + kick = kick + AB_sc(ia,eig)*c(mu,a,1) + end do + c(mu,i,1) = c(mu,i,1) + step*kick + end do + end do + + ! Spin-down kick + + do mu=1,nBas + ia = nS_aa + do i=nC(2)+1,nO(2) + kick = 0d0 + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + kick = kick + AB_sc(ia,eig)*c(mu,a,2) + end do + c(mu,i,2) = c(mu,i,2) + step*kick + end do + end do + + else + + write(*,'(1X,A40,1X)') 'Well done, UHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sc(1),' au' + + unstab = .false. + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +!---------------! +! End of Search ! +!---------------! + end do + +end subroutine diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index b4bac6c..3eb5fe5 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -174,6 +174,19 @@ subroutine GQuAcK(doGHF,dostab,dosearch,doMP2,doMP3,dophRPA,dophRPAx,doppRPA, end if + if(dosearch) then + + call wall_time(start_stab) + call GHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nBas2,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + call wall_time(end_stab) + + t_stab = end_stab - start_stab + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds' + write(*,*) + + end if + !-----------------------! ! Moller-Plesset module ! !-----------------------! diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index 02e1b61..bdc3d7b 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -201,6 +201,20 @@ subroutine UQuAcK(doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,do end if + if(dosearch) then + + call wall_time(start_stab) + 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,dipole_int_AO,X,EHF,epsHF,cHF,PHF) + + call wall_time(end_stab) + + t_stab = end_stab - start_stab + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds' + write(*,*) + + end if + !-----------------------! ! Moller-Plesset module ! !-----------------------!