mirror of
https://github.com/pfloos/quack
synced 2024-11-03 20:53:53 +01:00
GHF search
This commit is contained in:
parent
67abf13740
commit
7847710059
@ -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)
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
210
src/HF/GHF_search.f90
Normal file
210
src/HF/GHF_search.f90
Normal file
@ -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
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
223
src/HF/UHF_search.f90
Normal file
223
src/HF/UHF_search.f90
Normal file
@ -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
|
||||
@ -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 !
|
||||
!-----------------------!
|
||||
|
||||
@ -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 !
|
||||
!-----------------------!
|
||||
|
||||
Loading…
Reference in New Issue
Block a user