From 17ce993aa2269fe4626a181fd26e7a36bae7d596 Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 3 Nov 2023 15:56:18 +0100 Subject: [PATCH] RHF search --- src/HF/RHF.f90 | 18 ++--- src/HF/RHF_search.f90 | 166 +++++++++++++++++++++++++++++++++++++++ src/HF/RHF_stability.f90 | 2 + src/HF/mo_guess.f90 | 9 ++- src/QuAcK/GQuAcK.f90 | 2 +- src/QuAcK/QuAcK.f90 | 2 +- src/QuAcK/RQuAcK.f90 | 15 +++- src/QuAcK/UQuAcK.f90 | 2 +- 8 files changed, 202 insertions(+), 14 deletions(-) create mode 100644 src/HF/RHF_search.f90 diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 4fae804..499a144 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,5 +1,5 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,eps,c,P) + nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P) ! Perform restricted Hartree-Fock calculation @@ -55,8 +55,8 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc ! Output variables double precision,intent(out) :: EHF - double precision,intent(out) :: eps(nBas) - double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: e(nBas) + double precision,intent(inout):: c(nBas,nBas) double precision,intent(out) :: P(nBas,nBas) ! Hello world @@ -136,7 +136,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eps) + call diagonalize_matrix(nBas,cp,e) c = matmul(X,cp) ! Density matrix @@ -153,20 +153,20 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc if(nBas > nO) then - Gap = eps(nO+1) - eps(nO) + Gap = e(nO+1) - e(nO) else Gap = 0d0 - endif + end if ! Dump results write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & '|',nSCF,'|',EHF+ENuc,'|',Conv,'|',Gap,'|' - enddo + end do write(*,*)'----------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop @@ -184,7 +184,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc stop - endif + end if ! Compute HF energy @@ -197,6 +197,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc ! Compute dipole moments call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_RHF(nBas,nO,eps,C,ENuc,ET,EV,EJ,EK,EHF,dipole) + call print_RHF(nBas,nO,e,C,ENuc,ET,EV,EJ,EK,EHF,dipole) end subroutine diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 new file mode 100644 index 0000000..bc94d6a --- /dev/null +++ b/src/HF/RHF_search.f90 @@ -0,0 +1,166 @@ +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 + + 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(in) :: level_shift + + integer,intent(in) :: nBas + 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(:,:,:,:) + integer :: nS + + integer,parameter :: maxS = 20 + integer :: ia + integer :: ispin + + double precision,allocatable :: A(:,:) + double precision,allocatable :: B(:,:) + double precision,allocatable :: AB(:,:) + double precision,allocatable :: Om(:) + + integer :: eig + double precision :: eigval + double precision,allocatable :: eigvec(:) + +! 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 RHF solutions *' + write(*,*) '****************************' + write(*,*) + +!-------------------! +! Memory allocation ! +!-------------------! + + 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)) + +!------------------! +! Search algorithm ! +!------------------! + + unstab = .true. + guess = 0 + + do while(unstab) + +!---------------------! +! Hartree-Fock module ! +!---------------------! + + call wall_time(start_HF) + call RHF(maxSCF,thresh,max_diis,guess,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 RHF = ',t_HF,' seconds' + write(*,*) + +!----------------------------------! +! AO to MO integral transformation ! +!----------------------------------! + + call wall_time(start_AOtoMO) + write(*,*) + write(*,*) 'AO to MO transformation... Please be patient' + write(*,*) + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) + 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 RHF -> Real RHF +!-------------------------------------------------------------! + + 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) + + AB(:,:) = A(:,:) + B(:,:) + + call diagonalize_matrix(nS,AB,Om) + Om(:) = 0.5d0*Om(:) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real RHF -> Real RHF |' + 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, 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?' + read(*,*) eig + + eigval = Om(eig) + eigvec = AB(:,eig) + + else + + write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au' + + unstab = .false. + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + + end do + +end subroutine diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 index 72a2099..a96ca7c 100644 --- a/src/HF/RHF_stability.f90 +++ b/src/HF/RHF_stability.f90 @@ -72,6 +72,8 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) write(*,*)'-------------------------------------------------------------' write(*,*) + + !-------------------------------------------------------------! ! Stability analysis: Real RHF -> Complex RHF !-------------------------------------------------------------! diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index 91e6e88..81c4cbc 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -16,18 +16,25 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c) double precision,intent(out) :: c(nBas,nBas) - if(guess_type == 1) then + if(guess_type == 0) then + write(*,*) 'Reading HF coefficients...' + + elseif(guess_type == 1) then + + write(*,*) 'Core guess...' call core_guess(nBas,Hc,X,c) elseif(guess_type == 2) then + write(*,*) 'Huckel guess...' call huckel_guess(nBas,S,Hc,X,c) elseif(guess_type == 3) then call random_number(c) + write(*,*) 'Random guess...' c(:,:) = 2d0*c(:,:) - 1d0 else diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 844f593..7ce71cc 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -35,7 +35,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,doMP3,dophRPA,dophRPAx,doppRPA,doG0W0,doevG integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix - logical,intent(in) :: guess_type + integer,intent(in) :: guess_type logical,intent(in) :: reg_MP diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index acc97e9..7401992 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -37,7 +37,7 @@ program QuAcK integer :: maxSCF_HF,max_diis_HF double precision :: thresh_HF,level_shift,mix - logical :: guess_type + integer :: guess_type logical :: reg_MP diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index a50dab6..cc6bbc1 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -43,7 +43,7 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix - logical,intent(in) :: guess_type + integer,intent(in) :: guess_type logical,intent(in) :: reg_MP @@ -182,6 +182,19 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC end if + if(dostab) 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,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 5eea1e8..303d5a1 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -41,7 +41,7 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix - logical,intent(in) :: guess_type + integer,intent(in) :: guess_type logical,intent(in) :: reg_MP