diff --git a/input/dft b/input/dft index b48d7c9..aaee5fd 100644 --- a/input/dft +++ b/input/dft @@ -6,24 +6,24 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4 HF,B3LYP,PBE -1 S51 +1 S51 # correlation rung: # Hartree = 0: H # LDA = 1: PW92,VWN3,VWN5,eVWN5 # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE -1 VWN5 +1 VWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) -2 +4 # occupation numbers -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 +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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 @@ -36,8 +36,8 @@ F # Parameters for CC weight-dependent exchange functional 4 -0.642674 -0.07818 -0.0280307 0.00144198 -0.254939 -0.0893405 0.00765581 0. +-0.718713,-0.133321,0.226288,-0.250718 +-0.525899,0.687216,-0.13866,-0.0226579 0.0 0.0 0.0 0.0 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2 1 diff --git a/mol/h2.xyz b/mol/h2.xyz index 3a1f2fe..9e53116 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 - + H 0. 0. 0. H 0. 0. 0.741 diff --git a/src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 index 96e8f7c..e7998d2 100644 --- a/src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 @@ -1,4 +1,6 @@ -subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,ExDD) +subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,& + doNcentered,kappa,ExDD) + ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative @@ -16,6 +18,7 @@ subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weig double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) ! Local variables @@ -153,7 +156,13 @@ subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weig do iEns=1,nEns do jEns=2,nEns - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) + if(doNcentered) then + + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - kappa(iEns)*wEns(jEns))*dExdw(jEns) + else + + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) + end if end do end do diff --git a/src/eDFT/UKS.f90 b/src/eDFT/UKS.f90 index 152b35f..afa4914 100644 --- a/src/eDFT/UKS.f90 +++ b/src/eDFT/UKS.f90 @@ -352,6 +352,9 @@ subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,t end do write(*,*)'------------------------------------------------------------------------------------------' +! print*,'Ensemble energy:',Ew + ENuc,'au' + + !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ diff --git a/src/eDFT/exchange_derivative_discontinuity.f90 b/src/eDFT/exchange_derivative_discontinuity.f90 index 68d400a..bc8485e 100644 --- a/src/eDFT/exchange_derivative_discontinuity.f90 +++ b/src/eDFT/exchange_derivative_discontinuity.f90 @@ -1,5 +1,5 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,drhow,& - Cx_choice,doNcentered,ExDD) + Cx_choice,doNcentered,kappa,ExDD) ! Compute the exchange part of the derivative discontinuity @@ -20,11 +20,12 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,we double precision,intent(in) :: drhow(ncart,nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) -! Local variables +!Local variables -! Output variables +!Output variables double precision,intent(out) :: ExDD(nEns) @@ -41,7 +42,7 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,we case(1) call lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& - rhow(:),Cx_choice,doNcentered,ExDD(:)) + rhow(:),Cx_choice,doNcentered,kappa,ExDD(:)) ! GGA functionals case(2) diff --git a/src/eDFT/exchange_energy.f90 b/src/eDFT/exchange_energy.f90 index 558a4d9..a787910 100644 --- a/src/eDFT/exchange_energy.f90 +++ b/src/eDFT/exchange_energy.f90 @@ -49,7 +49,7 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight, case(2) - call gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + call gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex) ! MGGA functionals diff --git a/src/eDFT/exchange_potential.f90 b/src/eDFT/exchange_potential.f90 index c37ace9..3322ab0 100644 --- a/src/eDFT/exchange_potential.f90 +++ b/src/eDFT/exchange_potential.f90 @@ -59,7 +59,8 @@ subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weig case(2) - call gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,& + Cx_choice,Fx) ! MGGA functionals diff --git a/src/eDFT/gga_exchange_energy.f90 b/src/eDFT/gga_exchange_energy.f90 index 3f0d99b..fba9cd0 100644 --- a/src/eDFT/gga_exchange_energy.f90 +++ b/src/eDFT/gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) +subroutine gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex) ! Select GGA exchange functional for energy calculation @@ -11,11 +11,15 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) integer,intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nCC + double precision,intent(in) :: aCC(nCC,nEns-1) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) + integer,intent(in) :: Cx_choice double precision,intent(in) :: drho(ncart,nGrid) + ! Output variables double precision :: Ex @@ -34,6 +38,11 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) call PBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + case (4) + + call CC_B88_gga_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,& + Cx_choice,Ex) + case default call print_warning('!!! GGA exchange energy not available !!!') diff --git a/src/eDFT/gga_exchange_potential.f90 b/src/eDFT/gga_exchange_potential.f90 index fb2b60c..81f05a3 100644 --- a/src/eDFT/gga_exchange_potential.f90 +++ b/src/eDFT/gga_exchange_potential.f90 @@ -1,4 +1,5 @@ -subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,& + rho,drho,Cx_choice,Fx) ! Select GGA exchange functional for potential calculation @@ -10,6 +11,8 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh integer,intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nCC + double precision,intent(in) :: aCC(nCC,nEns-1) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: nBas @@ -17,6 +20,7 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh double precision,intent(in) :: dAO(3,nBas,nGrid) double precision,intent(in) :: rho(nGrid) double precision,intent(in) :: drho(3,nGrid) + integer,intent(in) :: Cx_choice ! Output variables @@ -38,6 +42,11 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh call PBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + case (4) + + call CC_B88_gga_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,& + Cx_choice,Fx) + case default call print_warning('!!! GGA exchange potential not available !!!') diff --git a/src/eDFT/individual_energy.f90 b/src/eDFT/individual_energy.f90 index 47baed8..b546bab 100644 --- a/src/eDFT/individual_energy.f90 +++ b/src/eDFT/individual_energy.f90 @@ -150,21 +150,22 @@ subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nC do ispin=1,nspin call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nCC,aCC,nGrid,weight, & - rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,ExDD(ispin,:)) + rhow(:,ispin),drhow(:,:,ispin),Cx_choice,& + doNcentered,kappa,ExDD(ispin,:)) end do call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD) ! Scaling derivative discontinuity for N-centered ensembles - if(doNcentered) then +! if(doNcentered) then - do iEns=1,nEns - ExDD(:,iEns) = (1d0 - kappa(iEns))*ExDD(:,iEns) - EcDD(:,iEns) = (1d0 - kappa(iEns))*EcDD(:,iEns) - end do +! do iEns=1,nEns +! ExDD(:,iEns) = (1d0 - kappa(iEns))*ExDD(:,iEns) +! EcDD(:,iEns) = (1d0 - kappa(iEns))*EcDD(:,iEns) +! end do - end if +! end if !------------------------------------------------------------------------ ! Total energy @@ -187,6 +188,8 @@ subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nC end do end if + +print*,'LZshift=',sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) ! do iEns=1,nEns ! E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) & diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index 2a32f53..5065d69 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -1,5 +1,5 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& - Cx_choice,doNcentered,ExDD) + Cx_choice,doNcentered,kappa,ExDD) ! Compute the exchange LDA part of the derivative discontinuity @@ -19,6 +19,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,wei double precision,intent(in) :: rhow(nGrid) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered + double precision,intent(in) :: kappa(nEns) ! Local variables @@ -38,7 +39,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,wei case (2) call CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,& - Cx_choice,doNcentered,ExDD) + Cx_choice,doNcentered,kappa,ExDD) case default diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90 index d718477..ea972a5 100644 --- a/src/eDFT/read_options_dft.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -117,7 +117,11 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, case ('PBE') x_DFA = 3 - + + case ('CC-B88') + + x_DFA = 4 + case default call print_warning('!!! GGA exchange functional not available !!!')