From 03f41ef20cb229843e05c7ba46fc6722e0032a57 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 2 Feb 2022 15:06:51 +0100 Subject: [PATCH] level shift and proper guess --- input/dft | 10 ++-- input/methods | 6 +-- input/options | 4 +- mol/h2.xyz | 2 +- src/HF/RHF.f90 | 53 +++++++++---------- src/HF/UHF.f90 | 41 +++++++++------ src/HF/core_guess.f90 | 33 ++++++++++++ src/HF/huckel_guess.f90 | 37 +++++-------- src/HF/mo_guess.f90 | 20 ++----- src/QuAcK/QuAcK.f90 | 26 +++++----- src/QuAcK/read_options.f90 | 29 ++++++----- src/eDFT/{eDFT_UKS.f90 => UKS.f90} | 83 ++++++++++++------------------ src/eDFT/eDFT.f90 | 13 +++-- 13 files changed, 180 insertions(+), 177 deletions(-) create mode 100644 src/HF/core_guess.f90 rename src/eDFT/{eDFT_UKS.f90 => UKS.f90} (88%) diff --git a/input/dft b/input/dft index 24eda5d..b48d7c9 100644 --- a/input/dft +++ b/input/dft @@ -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 diff --git a/input/methods b/input/methods index 1ca495f..688d56a 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index a48a761..6d53fb5 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/mol/h2.xyz b/mol/h2.xyz index 3a1f2fe..a4e936a 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 0.741 +H 0. 0. 2.000000 diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 78e1dad..6cba349 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -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 diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index c6801c2..7917d00 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -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 !------------------------------------------------------------------------ diff --git a/src/HF/core_guess.f90 b/src/HF/core_guess.f90 new file mode 100644 index 0000000..7e45e8e --- /dev/null +++ b/src/HF/core_guess.f90 @@ -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 diff --git a/src/HF/huckel_guess.f90 b/src/HF/huckel_guess.f90 index 91da7cd..80efaa8 100644 --- a/src/HF/huckel_guess.f90 +++ b/src/HF/huckel_guess.f90 @@ -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 diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index d478ab6..387f371 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 67ff676..890b48f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 383aaee..6393922 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -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 diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/UKS.f90 similarity index 88% rename from src/eDFT/eDFT_UKS.f90 rename to src/eDFT/UKS.f90 index 8888243..3d55bf1 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/UKS.f90 @@ -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 diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 0ce947b..ffc0632 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -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