PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-31 08:35:51 +01:00
Code Releases Activity

fix bug in GHF search

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-07 16:11:01 +01:00
parent 66f6890082
commit b794290c96
9 changed files with 54 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF GHF ROHF
T F F F
F F T F
# MP2 MP3
F F
# CCD pCCD DCD CCSD CCSD(T)

2
input/options
View File

@ -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.0000001 5 1 0.0 0.0 T T
# MP: reg
F
# CC: maxSCF thresh DIIS

49
src/HF/GHF_search.f90
View File

@ -52,6 +52,9 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision,allocatable :: Bph(:,:)
double precision,allocatable :: AB(:,:)
double precision,allocatable :: Om(:)
double precision,allocatable :: R(:,:)
double precision,allocatable :: ExpR(:,:)
integer :: eig
double precision :: kick,step
@ -59,9 +62,9 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! 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)
double precision,intent(out) :: e(nBas2)
double precision,intent(inout):: c(nBas2,nBas2)
double precision,intent(out) :: P(nBas2,nBas2)
! Memory allocation
@ -77,7 +80,8 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
nS = (nO - nC)*(nV - nR)
allocate(ERI_MO(nBas2,nBas2,nBas2,nBas2),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS))
allocate(ERI_MO(nBas2,nBas2,nBas2,nBas2),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS), &
R(nBas2,nBas2),ExpR(nBas2,nBas2))
!------------------!
! Search algorithm !
@ -145,11 +149,11 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
call phLR_A(ispin,.false.,nBas2,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phLR_B(ispin,.false.,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
AB(:,:) = Aph(:,:) - Bph(:,:)
AB(:,:) = Aph(:,:) + Bph(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real GHF -> Real GHF |'
@ -182,18 +186,31 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
step = 1d0
do mu=1,nBas2
ia = 0
do i=nC+1,nO
kick = 0d0
do a=nO+1,nBas2-nR
ia = ia + 1
kick = kick + AB(ia,eig)*c(mu,a)
end do
c(mu,i) = c(mu,i) + step*kick
! do mu=1,nBas2
! ia = 0
! do i=nC+1,nO
! kick = 0d0
! do a=nO+1,nBas2-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
R(:,:) = 0d0
ia = 0
do i=nC+1,nO
do a=nO+1,nBas2-nR
ia = ia + 1
R(a,i) = +AB(ia,eig)
R(i,a) = -AB(ia,eig)
end do
end do
call matrix_exponential(nBas2,R,ExpR)
c = matmul(c,ExpR)
else
write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!'

4
src/HF/GHF_stability.f90
View File

@ -44,7 +44,7 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
AB(:,:) = A(:,:) + B(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real GHF -> Real GHF |'
@ -79,7 +79,7 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
AB(:,:) = A(:,:) - B(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real GHF -> Complex GHF |'

4
src/HF/RHF_search.f90
View File

@ -123,10 +123,10 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
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(:,:)
AB(:,:) = Aph(:,:) - Bph(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real RHF -> Real RHF |'

6
src/HF/RHF_stability.f90
View File

@ -44,7 +44,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
AB(:,:) = A(:,:) + B(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real RHF -> Real RHF |'
@ -81,7 +81,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
AB(:,:) = A(:,:) - B(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real RHF -> Complex RHF |'
@ -121,7 +121,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
AB(:,:) = A(:,:) + B(:,:)
call diagonalize_matrix(nS,AB,Om)
Om(:) = 0.5d0*Om(:)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real RHF -> Real UHF |'

2
src/HF/UHF_search.f90
View File

@ -144,7 +144,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
AB_sc(:,:) = A_sc(:,:) + B_sc(:,:)
call diagonalize_matrix(nS_sc,AB_sc,Om_sc)
Om_sc(:) = 0.5d0*Om_sc(:)
Om_sc(:) = 2d0*Om_sc(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real UHF -> Real UHF |'

6
src/HF/UHF_stability.f90
View File

@ -57,7 +57,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
AB_sc(:,:) = A_sc(:,:) + B_sc(:,:)
call diagonalize_matrix(nS_sc,AB_sc,Om_sc)
Om_sc(:) = 0.5d0*Om_sc(:)
Om_sc(:) = 2d0*Om_sc(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real UHF -> Real UHF |'
@ -92,7 +92,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
AB_sc(:,:) = A_sc(:,:) - B_sc(:,:)
call diagonalize_matrix(nS_sc,AB_sc,Om_sc)
Om_sc(:) = 0.5d0*Om_sc(:)
Om_sc(:) = 2d0*Om_sc(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real UHF -> Complex UHF |'
@ -142,7 +142,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
AB_sf(:,:) = A_sf(:,:) + B_sf(:,:)
call diagonalize_matrix(nS_sf,AB_sf,Om_sf)
Om_sf(:) = 0.5d0*Om_sf(:)
Om_sf(:) = 2d0*Om_sf(:)
write(*,*)'-------------------------------------------------------------'
write(*,*)'| Stability analysis: Real UHF -> Real GHF |'

9
src/utils/utils.f90
View File

@ -72,6 +72,7 @@ subroutine matrix_exponential(N,A,ExpA)
implicit none
integer,intent(in) :: N
integer :: i
double precision,intent(in) :: A(N,N)
double precision,allocatable :: W(:,:)
double precision,allocatable :: tau(:)
@ -90,7 +91,13 @@ subroutine matrix_exponential(N,A,ExpA)
W(:,:) = - matmul(A,A)
call diagonalize_matrix(N,W,tau)
tau(:) = sqrt(tau(:))
call matout(N,1,tau)
! do i=1,N
! tau(i) = max(abs(tau(i)),1d-14)
! end do
tau(:) = sqrt(abs(tau(:)))
call matout(N,1,tau)
! Construct cos part