4
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00

fix bug in eDFT rcond not initialized

This commit is contained in:
Pierre-Francois Loos 2022-01-24 10:29:27 +01:00
parent 2fa1010fd9
commit 0f7cde9766
11 changed files with 25 additions and 24 deletions

View File

@ -19,11 +19,11 @@
# Number of states in ensemble (nEns)
2
# occupation numbers
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
@ -31,7 +31,7 @@
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.6 0.0 0.0
0.5 0.0 0.0
# Ncentered ?
F
# Parameters for CC weight-dependent exchange functional

View File

@ -1,5 +1,5 @@
# RHF UHF KS MOM
F T F F
F F T F
# MP2* MP3 MP2-F12
F F F
# CCD pCCD DCD CCSD CCSD(T)

View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis
64 0.00001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip
F T F T T
F T T T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
256 0.00001 T 5 T 0.0 3 F
# GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg

View File

@ -1,4 +1,4 @@
2
H 0. 0. 0.
H 0. 0. 1.500000
H 0. 0. 2.000000

View File

@ -86,7 +86,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
write(*,*) 'nH = ',nH
write(*,*)
maxH = min(nH,21)
maxH = min(nH,41)
! Memory allocation

View File

@ -96,13 +96,13 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO
! First-order correction
if(ispin == 1) then
ZDyn(ia) = dot_product(X,matmul(ZAt+ZAs,X))
OmDyn(ia) = dot_product(X,matmul(dTAt+dTAs,X)) - dot_product(X,matmul(TAt+TAs,X))
ZDyn(ia) = - dot_product(X,matmul(ZAt+ZAs,X))
OmDyn(ia) = - dot_product(X,matmul(dTAt+dTAs,X)) + dot_product(X,matmul(TAt+TAs,X))
end if
if(ispin == 2) then
ZDyn(ia) = dot_product(X,matmul(ZAt-ZAs,X))
OmDyn(ia) = dot_product(X,matmul(dTAt-dTAs,X)) - dot_product(X,matmul(TAt-TAs,X))
ZDyn(ia) = - dot_product(X,matmul(ZAt-ZAs,X))
OmDyn(ia) = - dot_product(X,matmul(dTAt-dTAs,X)) + dot_product(X,matmul(TAt-TAs,X))
end if
ZDyn(ia) = 1d0/(1d0 - ZDyn(ia))

View File

@ -67,7 +67,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
end do
TA(ia,jb) = TA(ia,jb) + lambda*chi
TA(ia,jb) = TA(ia,jb) - lambda*chi
chi = 0d0
@ -81,7 +81,7 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2
end do
ZA(ia,jb) = ZA(ia,jb) - lambda*chi
ZA(ia,jb) = ZA(ia,jb) + lambda*chi
end do
end do

View File

@ -46,13 +46,13 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
chi = 0d0
do cd=1,nVV
eps = - Omega1(cd)
eps = + Omega1(cd)
! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2)
enddo
do kl=1,nOO
eps = + Omega2(kl)
eps = - Omega2(kl)
! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2)
enddo

View File

@ -46,13 +46,13 @@ subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O
chi = 0d0
do cd=1,nVV
eps = - Omega1(cd)
eps = + Omega1(cd)
! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2
chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2)
enddo
do kl=1,nOO
eps = + Omega2(kl)
eps = - Omega2(kl)
! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2
chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2)
enddo

View File

@ -47,7 +47,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
! print*,'TA'
! call matout(nS,nS,A_BSE)
A(:,:) = A(:,:) + A_BSE(:,:)
A(:,:) = A(:,:) - A_BSE(:,:)
! Tamm-Dancoff approximation
@ -68,7 +68,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda
! print*,'TB'
! call matout(nS,nS,B_BSE)
B(:,:) = B(:,:) + B_BSE(:,:)
B(:,:) = B(:,:) - B_BSE(:,:)
! Build A + B and A - B matrices

View File

@ -163,6 +163,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
n_diis = 0
F_diis(:,:,:) = 0d0
err_diis(:,:,:) = 0d0
rcond(:) = 1d0
!------------------------------------------------------------------------
! Main SCF loop
@ -268,7 +269,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
if(minval(rcond(:)) > 1d-7) then
if(minval(rcond(:)) > 1d-15) then
do ispin=1,nspin
call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, &
err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin))