RHF search

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)  :: e(nBas)
-  double precision,intent(out)  :: c(nBas,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 

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

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  
 !-------------------------------------------------------------!

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

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
 

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
 

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 !
 !-----------------------!

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