10
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:42 +01:00

fix bug for zero electrons

This commit is contained in:
Pierre-Francois Loos 2020-10-06 15:24:24 +02:00
parent e203b3c711
commit 9212d69aa2
4 changed files with 13 additions and 9 deletions

View File

@ -7,13 +7,13 @@
# drCCD rCCD lCCD pCCD
F F F F
# CIS* CIS(D) CID CISD
F F F F
T F F F
# RPA* RPAx* ppRPA
F F F
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0* evGW* qsGW
T F F
F F F
# G0T0 evGT qsGT
F F F
# MCMP2

View File

@ -80,7 +80,8 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Guess coefficients and eigenvalues
if(guess_type == 1) then
F(:,:,:) = 0d0
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:)
end do
@ -154,8 +155,9 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
do ispin=1,nspin
call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin))
end do
! Build Fock operator
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin)
end do
@ -166,7 +168,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin))
end do
conv = maxval(abs(err(:,:,:)))
if(nSCF > 1) conv = maxval(abs(err(:,:,:)))
! DIIS extrapolation

View File

@ -49,7 +49,7 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
! Compute <S**2>
call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2)
call unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
! Print details about spin-conserved excitations
@ -60,7 +60,6 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:))
Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:))
print*,'-------------------------------------------------------------'
write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' <S**2> = ',S2(ia)

View File

@ -1,4 +1,4 @@
subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,XpY,XmY,S2)
subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S,Omega,XpY,XmY,S2)
! Compute <S**2> for linear response excited states
@ -21,6 +21,7 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S
integer,intent(in) :: maxS
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Omega(nSt)
double precision,intent(in) :: XpY(nSt,nSt)
double precision,intent(in) :: XmY(nSt,nSt)
@ -152,7 +153,9 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S
Yat(:,:) = transpose(Ya(:,:))
Ybt(:,:) = transpose(Yb(:,:))
S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0 &
S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0
S2(m) = S2(m) &
+ trace_matrix(nV(1),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
- trace_matrix(nO(2),matmul(Xb,matmul(VO,matmul(VOt,Xbt)))) &