4
1
mirror of https://github.com/pfloos/quack synced 2024-11-19 04:22:39 +01:00

level shift and proper guess

This commit is contained in:
Pierre-Francois Loos 2022-02-02 15:06:51 +01:00
parent a1e81ab7e5
commit 03f41ef20c
13 changed files with 180 additions and 177 deletions

View File

@ -19,11 +19,11 @@
# Number of states in ensemble (nEns)
2
# occupation numbers
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
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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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.1 0.0 0.0
0.0 0.0 0.0
# Ncentered ?
F
# Parameters for CC weight-dependent exchange functional

View File

@ -11,11 +11,11 @@
# RPA* RPAx* crRPA ppRPA
F F F F
# G0F2* evGF2* qsGF2* G0F3 evGF3
T F F F F
F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW
F F F F F
T F F F F
# G0T0 evGT qsGT
F F F
T F F
# MCMP2
F
# * unrestricted version available

View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability
1024 0.00001 T 5 1 1 T F
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability
256 0.0000001 T 5 2 1 T T F
# MP:
# CC: maxSCF thresh DIIS n_diis

View File

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

View File

@ -1,4 +1,5 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx)
subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx)
! Perform restricted Hartree-Fock calculation
@ -7,8 +8,11 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Input variables
integer,intent(in) :: maxSCF,max_diis,guess_type
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
logical,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nO
@ -46,7 +50,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: ON(:)
! Output variables
@ -71,30 +74,22 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), &
cp(nBas,nBas),Fp(nBas,nBas),ON(nBas), &
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),Fp(nBas,nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and eigenvalues
! Guess coefficients and density matrix
call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! ON(:) = 0d0
! do i=1,nO
! ON(i) = 1d0
! ON(i) = dble(2*i-1)
! end do
! call density_matrix(nBas,ON,c,P)
call mo_guess(nBas,guess_type,S,Hc,X,c)
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
Conv = 1d0
nSCF = 0
rcond = 0d0
Conv = 1d0
n_diis = 0
nSCF = 0
rcond = 0d0
!------------------------------------------------------------------------
! Main SCF loop
@ -123,17 +118,25 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Check convergence
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error))
Conv = maxval(abs(error))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
if(abs(rcond) > 1d-7) then
if(abs(rcond) > 1d-15) then
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
else
n_diis = 0
end if
! Level-shifting
if(level_shift .and. Conv > thresh) call level_shifting(nBas,nO,S,c,F)
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
@ -144,14 +147,12 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! call density_matrix(nBas,ON,c,P)
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.25d0*trace_matrix(nBas,matmul(P,K))
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.25d0*trace_matrix(nBas,matmul(P,K))
! Compute HOMO-LUMO gap

View File

@ -1,4 +1,5 @@
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx)
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx)
! Perform unrestricted Hartree-Fock calculation
@ -11,6 +12,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
logical,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
@ -79,22 +81,12 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
K(nBas,nBas,nspin),err(nBas,nBas,nspin),cp(nBas,nBas,nspin), &
err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin))
! Guess coefficients and eigenvalues
! Guess coefficients and demsity matrices
if(guess_type == 1) then
F(:,:,:) = 0d0
do ispin=1,nspin
F(:,:,ispin) = Hc(:,:)
end do
else if(guess_type == 2) then
do ispin=1,nspin
call random_number(F(:,:,ispin))
end do
end if
do ispin=1,nspin
call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do
! Initialization
@ -179,13 +171,28 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO
! 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
if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), &
F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin))
end do
else
n_diis = 0
end if
! Level-shifting
if(level_shift .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin))
end do
end if
!------------------------------------------------------------------------

33
src/HF/core_guess.f90 Normal file
View File

@ -0,0 +1,33 @@
subroutine core_guess(nBas,Hc,X,c)
! Core guess of the molecular orbitals for HF calculation
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
! Local variables
double precision,allocatable :: cp(:,:)
double precision,allocatable :: e(:)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
! Memory allocation
allocate(cp(nBas,nBas),e(nBas))
! Core guess
cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
end subroutine

View File

@ -1,47 +1,39 @@
subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
subroutine huckel_guess(nBas,S,Hc,X,c)
! Hickel guess of the molecular orbitals for HF calculation
! Hickel guess
implicit none
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Local variables
integer :: mu,nu
double precision :: a
double precision,allocatable :: F(:,:)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
! Memory allocation
allocate(F(nBas,nBas))
! Extended Huckel parameter
a = 1.75d0
Fp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
call Coulomb_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
! GWH approximation
do mu=1,nBas
F(mu,mu) = Hc(mu,mu)
do nu=mu+1,nBas
F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
@ -50,9 +42,6 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
enddo
enddo
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
call core_guess(nBas,F,X,c)
end subroutine

View File

@ -1,4 +1,4 @@
subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
! Guess of the molecular orbitals for HF calculation
@ -7,19 +7,10 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: guess_type
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: J(nBas,nBas)
double precision,intent(inout):: K(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(inout):: cp(nBas,nBas)
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: Fp(nBas,nBas)
double precision,intent(inout):: e(nBas)
double precision,intent(inout):: P(nBas,nBas)
! Local variables
@ -31,14 +22,11 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
if(guess_type == 1) then
Fp = matmul(transpose(X),matmul(Hc,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
c = matmul(X,cp)
call core_guess(nBas,Hc,X,c)
elseif(guess_type == 2) then
call huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
call huckel_guess(nBas,S,Hc,X,c)
elseif(guess_type == 3) then
@ -51,6 +39,4 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P)
endif
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
end subroutine

View File

@ -112,7 +112,7 @@ program QuAcK
integer :: maxSCF_HF,n_diis_HF
double precision :: thresh_HF
logical :: DIIS_HF,guess_type,ortho_type,mix
logical :: DIIS_HF,guess_type,ortho_type,mix,level_shift
integer :: maxSCF_CC,n_diis_CC
double precision :: thresh_CC
@ -180,15 +180,15 @@ program QuAcK
! Read options for methods
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Weird stuff
@ -292,7 +292,7 @@ program QuAcK
end if
call cpu_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, &
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
@ -312,7 +312,7 @@ program QuAcK
unrestricted = .true.
call cpu_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc, &
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,Vxc)
call cpu_time(end_HF)
@ -332,7 +332,7 @@ program QuAcK
unrestricted = .true.
call cpu_time(start_KS)
call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, &
call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, &
nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc)

View File

@ -1,12 +1,12 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, &
maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, &
maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, &
COHSEX,SOSEX,TDA_W,G0W,GW0, &
maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
BSE,dBSE,dTDA,evDyn, &
nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift)
! Read desired methods
@ -22,6 +22,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
integer,intent(out) :: guess_type
integer,intent(out) :: ortho_type
logical,intent(out) :: mix
logical,intent(out) :: level_shift
logical,intent(out) :: dostab
integer,intent(out) :: maxSCF_CC
@ -100,14 +101,16 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
guess_type = 1
ortho_type = 1
mix = .false.
level_shift = .false.
dostab = .false.
read(1,*)
read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3
read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3,answer4
if(answer1 == 'T') DIIS_HF = .true.
if(answer2 == 'T') mix = .true.
if(answer3 == 'T') dostab = .true.
if(answer1 == 'T') DIIS_HF = .true.
if(answer2 == 'T') mix = .true.
if(answer3 == 'T') level_shift = .true.
if(answer4 == 'T') dostab = .true.
if(.not.DIIS_HF) n_diis_HF = 1

View File

@ -1,5 +1,5 @@
subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,c,Pw,Vxc)
subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix,level_shift, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,c,Pw,Vxc)
! Perform unrestricted Kohn-Sham calculation for ensembles
@ -20,6 +20,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
logical,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
@ -45,12 +46,11 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
integer :: xc_rung
logical :: LDA_centered = .false.
logical :: do_level_shift = .false.
integer :: nSCF
integer :: nBasSq
integer :: n_diis
integer :: nO(nspin)
double precision :: conv
integer :: nO(nspin,nEns)
double precision :: Conv
double precision :: rcond(nspin)
double precision :: ET(nspin)
double precision :: EV(nspin)
@ -119,42 +119,23 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Guess coefficients and eigenvalues
nO(:) = 0
do ispin=1,nspin
nO(ispin) = int(sum(occnum(:,ispin,1)))
end do
do ispin=1,nspin
call mo_guess(nBas,nO(ispin),guess_type,S,Hc,ERI,J(:,:,ispin),Fx(:,:,ispin),X,cp(:,:,ispin),F(:,:,ispin), &
Fp(:,:,ispin),eKS(:,ispin),c(:,:,ispin),Pw(:,:,ispin))
nO(:,:) = 0
do iEns=1,nEns
do ispin=1,nspin
nO(ispin,iEns) = int(sum(occnum(:,ispin,iEns)))
end do
end do
if(mix) call mix_guess(nBas,nO,c)
do ispin=1,nspin
call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
end do
! if(guess_type == 1) then
! Mix guess for UHF solution in singlet states
! do ispin=1,nspin
! cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
! call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin))
! c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
! end do
! ! Mix guess to enforce symmetry breaking
! if(mix) call mix_guess(nBas,nO,c)
! else if(guess_type == 2) then
! do ispin=1,nspin
! call random_number(F(:,:,ispin))
! end do
! else
! print*,'Wrong guess option'
! stop
! end if
if(mix) then
write(*,*) '!! guess mixing disabled in UKS !!'
write(*,*)
end if
! Initialization
@ -184,7 +165,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
'|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|'
write(*,*)'------------------------------------------------------------------------------------------'
do while(conv > thresh .and. nSCF < maxSCF)
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
@ -273,18 +254,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin))
end do
if(nSCF > 1) conv = maxval(abs(err(:,:,:)))
if(nSCF > 1) Conv = maxval(abs(err(:,:,:)))
! Level-shifting
if(do_level_shift) then
do ispin=1,nspin
call level_shifting(nBas,nO(ispin),S,c,F(:,:,ispin))
end do
end if
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
@ -302,6 +273,16 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
end do
! Level-shifting
if(level_shift .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin))
end do
end if
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
@ -366,7 +347,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',conv,'|',sum(nEl(:)),'|'
'|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',Conv,'|',sum(nEl(:)),'|'
end do
write(*,*)'------------------------------------------------------------------------------------------'
@ -405,4 +386,4 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max
call individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, &
AO,dAO,T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew)
end subroutine eDFT_UKS
end subroutine UKS

View File

@ -1,4 +1,4 @@
subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, &
subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, &
nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS,Vxc)
@ -15,6 +15,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
logical,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nNuc
@ -165,8 +166,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
end do
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
@ -182,8 +184,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n
if(method == 'eDFT-UKS') then
call cpu_time(start_KS)
call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, &
nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight, &
maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc)
call cpu_time(end_KS)
t_KS = end_KS - start_KS