diff --git a/input/basis b/input/basis index 0e9b57e..d92a3df 100644 --- a/input/basis +++ b/input/basis @@ -1,13 +1,26 @@ -1 3 +1 8 +S 6 + 1 1264.5857000 0.0019448 + 2 189.9368100 0.0148351 + 3 43.1590890 0.0720906 + 4 12.0986630 0.2371542 + 5 3.8063232 0.4691987 + 6 1.2728903 0.3565202 S 3 - 1 30.1678710 0.15432897 - 2 5.4951153 0.53532814 - 3 1.4871927 0.44463454 -S 3 - 1 1.3148331 -0.0999672 - 2 0.3055389 0.3995128 - 3 0.0993707 0.7001155 + 1 3.1964631 -0.1126487 + 2 0.7478133 -0.2295064 + 3 0.2199663 1.1869167 P 3 - 1 1.3148331 0.1559163 - 2 0.3055389 0.6076837 - 3 0.0993707 0.3919574 + 1 3.1964631 0.0559802 + 2 0.7478133 0.2615506 + 3 0.2199663 0.7939723 +S 1 + 1 0.0823099 1.0000000 +P 1 + 1 0.0823099 1.0000000 +S 1 + 1 0.0207000 1.0000000 +P 1 + 1 0.0207000 1.0000000 +D 1 + 1 0.4000000 1.0000000 diff --git a/input/dft b/input/dft index 53337d0..ef6a189 100644 --- a/input/dft +++ b/input/dft @@ -1,19 +1,19 @@ # Restricted or unrestricted KS calculation - GOK-RKS + LIM-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,S51,RMFL20 # GGA = 2: G96,B88 # Hybrid = 4 # Hartree-Fock = 666 - 1 RMFL20 + 1 RS51 # correlation rung: # Hartree = 0 # LDA = 1: W38,VWN5,C16,RMFL20 # GGA = 2: LYP # Hybrid = 4: B3LYP # Hartree-Fock = 666 - 1 RMFL20 + 1 RVWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) diff --git a/input/weight b/input/weight index 0e9b57e..d92a3df 100644 --- a/input/weight +++ b/input/weight @@ -1,13 +1,26 @@ -1 3 +1 8 +S 6 + 1 1264.5857000 0.0019448 + 2 189.9368100 0.0148351 + 3 43.1590890 0.0720906 + 4 12.0986630 0.2371542 + 5 3.8063232 0.4691987 + 6 1.2728903 0.3565202 S 3 - 1 30.1678710 0.15432897 - 2 5.4951153 0.53532814 - 3 1.4871927 0.44463454 -S 3 - 1 1.3148331 -0.0999672 - 2 0.3055389 0.3995128 - 3 0.0993707 0.7001155 + 1 3.1964631 -0.1126487 + 2 0.7478133 -0.2295064 + 3 0.2199663 1.1869167 P 3 - 1 1.3148331 0.1559163 - 2 0.3055389 0.6076837 - 3 0.0993707 0.3919574 + 1 3.1964631 0.0559802 + 2 0.7478133 0.2615506 + 3 0.2199663 0.7939723 +S 1 + 1 0.0823099 1.0000000 +P 1 + 1 0.0823099 1.0000000 +S 1 + 1 0.0207000 1.0000000 +P 1 + 1 0.0207000 1.0000000 +D 1 + 1 0.4000000 1.0000000 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index 00299cc..076ee58 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -1,5 +1,5 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew) + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,EwGIC) ! Perform restricted Kohn-Sham calculation for ensembles @@ -41,7 +41,6 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres double precision :: EJ double precision :: Ex double precision :: Ec - double precision :: Ew double precision,allocatable :: eps(:) double precision,allocatable :: c(:,:) @@ -72,6 +71,11 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres integer :: iEns +! Output variables + + double precision,intent(out) :: Ew + double precision,intent(out) :: EwGIC + ! Hello world write(*,*) @@ -326,6 +330,6 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & AO(:,:),dAO(:,:,:),nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc, & Pw(:,:),rhow(:),drhow(:,:),J(:,:),Fx(:,:),FxHF(:,:), & - Fc(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:)) + Fc(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:)) end subroutine GOK_RKS diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 3312ce8..3414e1f 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -7,7 +7,7 @@ program eDFT integer :: nNuc,nBas integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin) - double precision :: ENuc,EKS + double precision :: ENuc,Ew,EwGIC double precision,allocatable :: ZNuc(:),rNuc(:,:) @@ -131,31 +131,48 @@ program eDFT call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO) !------------------------------------------------------------------------ -! Compute RKS energy +! Compute GOK-RKS energy !------------------------------------------------------------------------ if(method == 'GOK-RKS') then call cpu_time(start_KS) call GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & - nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS) + nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,EwGIC) call cpu_time(end_KS) t_KS = end_KS - start_KS - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RKS = ',t_KS,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GOC-RKS = ',t_KS,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Compute RKS energy +!------------------------------------------------------------------------ + + if(method == 'LIM-RKS') then + + call cpu_time(start_KS) + call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & + nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc) + call cpu_time(end_KS) + + t_KS = end_KS - start_KS + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for LIM-RKS = ',t_KS,' seconds' write(*,*) end if !------------------------------------------------------------------------ -! Compute UKS energy +! Compute GOK-UKS energy (BROKEN) !------------------------------------------------------------------------ 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, & - nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS) + nBas,AO(:,:),dAO(:,:,:),nO(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew) call cpu_time(end_KS) t_KS = end_KS - start_KS diff --git a/src/eDFT/print_restricted_individual_energy.f90 b/src/eDFT/print_restricted_individual_energy.f90 index 8fedd75..c560b18 100644 --- a/src/eDFT/print_restricted_individual_energy.f90 +++ b/src/eDFT/print_restricted_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine print_restricted_individual_energy(nEns,Ew,EwGOC,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & +subroutine print_restricted_individual_energy(nEns,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,ExDD,EcDD,ExcDD,E, & Om,Omx,Omc,Omxc,OmxDD,OmcDD,OmxcDD) ! Print individual energies for eDFT calculation @@ -10,7 +10,7 @@ subroutine print_restricted_individual_energy(nEns,Ew,EwGOC,ET,EV,EJ,Ex,Ec,Exc,E integer,intent(in) :: nEns double precision,intent(in) :: Ew - double precision,intent(in) :: EwGOC + double precision,intent(in) :: EwGIC double precision,intent(in) :: ET(nEns) double precision,intent(in) :: EV(nEns) double precision,intent(in) :: EJ(nEns) @@ -33,7 +33,7 @@ subroutine print_restricted_individual_energy(nEns,Ew,EwGOC,ET,EV,EJ,Ex,Ec,Exc,E write(*,'(A60)') ' ENSEMBLE ENERGIES' write(*,'(A60)') '-------------------------------------------------' write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew, ' au' - write(*,'(A44,F16.10,A3)') ' GOC Ensemble energy: ',EwGOC,' au' + write(*,'(A44,F16.10,A3)') ' GIC Ensemble energy: ',EwGIC,' au' write(*,'(A60)') '-------------------------------------------------' write(*,*) diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index eaae936..69d26cb 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -1,5 +1,5 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO, & - nO,nV,T,V,ERI,ENuc,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,E,Om) + nO,nV,T,V,ERI,ENuc,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,EwGIC,E,Om) ! Compute individual energies as well as excitation energies @@ -41,7 +41,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Local variables - double precision :: EwGOC double precision :: ET(nEns) double precision :: EV(nEns) double precision :: EJ(nEns) @@ -56,6 +55,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Output variables + double precision,intent(out) :: EwGIC double precision,intent(out) :: E(nEns) double precision,intent(out) :: Om(nEns) @@ -130,9 +130,9 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Total energy with ghost-interaction correction !------------------------------------------------------------------------ - EwGOC = 0d0 + EwGIC = 0d0 do iEns=1,nEns - EwGOC = EwGOC + wEns(iEns)*E(iEns) + EwGIC = EwGIC + wEns(iEns)*E(iEns) end do !------------------------------------------------------------------------ @@ -157,7 +157,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri ! Dump results !------------------------------------------------------------------------ - call print_restricted_individual_energy(nEns,Ew,EwGOC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & + call print_restricted_individual_energy(nEns,Ew,EwGIC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & ExDD(:),EcDD(:),ExcDD(:),E(:), & Om(:),Omx(:),Omc(:),Omxc(:),OmxDD(:),OmcDD(:),OmxcDD(:))