diff --git a/input/dft b/input/dft index d6c9473..be7c0f3 100644 --- a/input/dft +++ b/input/dft @@ -1,12 +1,12 @@ # Restricted or unrestricted KS calculation - eDFT-UKS + GOK-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 1 UCC + 1 RCC # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index b8814b8..2533f99 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -1,4 +1,4 @@ -subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, & +subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c) ! Perform restricted Kohn-Sham calculation for ensembles @@ -14,6 +14,8 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: maxSCF @@ -234,7 +236,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri ! Compute exchange potential - call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,Pw(:,:),ERI(:,:,:,:), & + call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,Pw(:,:),ERI(:,:,:,:), & AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fx(:,:),FxHF(:,:)) ! Compute correlation potential @@ -292,7 +294,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri ! Exchange energy - call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas, & + call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex) ! Correlation energy @@ -340,7 +342,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & + call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), & nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), & J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:)) diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90 index 0ab3741..c0e9145 100644 --- a/src/eDFT/GOK_UKS.f90 +++ b/src/eDFT/GOK_UKS.f90 @@ -1,4 +1,4 @@ -subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh,max_diis,guess_type, & +subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew) ! Perform unrestricted Kohn-Sham calculation for ensembles @@ -14,6 +14,8 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) integer,intent(in) :: maxSCF,max_diis,guess_type double precision,intent(in) :: thresh integer,intent(in) :: nBas @@ -246,7 +248,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute exchange potential do ispin=1,nspin - call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,Pw(:,:,ispin),ERI(:,:,:,:), & + call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,nBas,Pw(:,:,ispin),ERI(:,:,:,:), & AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin),FxHF(:,:,ispin)) end do @@ -304,7 +306,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Exchange energy do ispin=1,nspin - call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & + call exchange_energy(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)) end do @@ -355,7 +357,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & + call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc, & Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), & Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:)) diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index 661a31f..a3ce34b 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -1,4 +1,4 @@ -subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & +subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,aCC_w1,aCC_w2,nGrid,weight, & maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) ! Perform restricted Kohn-Sham calculation for ensembles @@ -12,6 +12,8 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & character(len=12),intent(in) :: x_DFA,c_DFA logical,intent(in) :: LDA_centered integer,intent(in) :: nEns + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: maxSCF,max_diis,guess_type @@ -70,7 +72,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c) !------------------------------------------------------------------------ @@ -91,7 +93,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c) !------------------------------------------------------------------------ @@ -112,7 +114,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) -! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & +! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & ! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c) diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90 index 24b716d..01af545 100644 --- a/src/eDFT/MOM_RKS.f90 +++ b/src/eDFT/MOM_RKS.f90 @@ -1,5 +1,5 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & - maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) + aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) ! Perform restricted Kohn-Sham calculation for ensembles @@ -14,6 +14,8 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & integer,intent(in) :: nEns integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) integer,intent(in) :: maxSCF,max_diis,guess_type double precision,intent(in) :: thresh integer,intent(in) :: nBas @@ -70,7 +72,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c) !------------------------------------------------------------------------ @@ -91,7 +93,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) -! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & +! call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & ! max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c) !------------------------------------------------------------------------ @@ -112,7 +114,7 @@ subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c) !------------------------------------------------------------------------ diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 110bde0..759236c 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -164,7 +164,7 @@ program eDFT if(method == 'GOK-RKS') then call cpu_time(start_KS) - call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & + call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), & maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,c(:,:)) call cpu_time(end_KS) @@ -183,7 +183,7 @@ program eDFT call cpu_time(start_KS) call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & - maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & + aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:)) call cpu_time(end_KS) @@ -201,7 +201,7 @@ program eDFT call cpu_time(start_KS) call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & - maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & + aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:)) call cpu_time(end_KS) @@ -218,7 +218,7 @@ program eDFT if(method == 'GOK-UKS') then call cpu_time(start_KS) - call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & + call GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew) call cpu_time(end_KS) diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options.f90 index 34fa3a7..f133ba0 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options.f90 @@ -118,12 +118,6 @@ subroutine read_options(method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aC write(*,*)' parameters for w2-dependant exchange functional coefficient ' write(*,*)'----------------------------------------------------------' call matout(3,1,aCC_w2) - write(*,*) - - write(*,*)'----------------------------------------------------------' - write(*,*)' Ensemble weights ' - write(*,*)'----------------------------------------------------------' - call matout(nEns,1,wEns) write(*,*) ! Read KS options diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 29ce420..173376e 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas, & +subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om) ! Compute individual energies as well as excitation energies @@ -13,6 +13,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n logical,intent(in) :: LDA_centered integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: nBas @@ -89,7 +91,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n !------------------------------------------------------------------------ do iEns=1,nEns - call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), & + call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,ERI(:,:,:,:), & Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) end do @@ -112,7 +114,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n ! Compute derivative discontinuities !------------------------------------------------------------------------ - call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) + call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),EcDD(:))