diff --git a/input/basis b/input/basis index a1f7a8d..dc1936c 100644 --- a/input/basis +++ b/input/basis @@ -1,27 +1,43 @@ -1 5 +1 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +2 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 diff --git a/input/dft b/input/dft index 01b4f72..883e724 100644 --- a/input/dft +++ b/input/dft @@ -6,7 +6,7 @@ # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 1 RGIC + 666 HF # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 @@ -19,6 +19,6 @@ # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.5 + 1.0 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.00001 T 5 1 1 diff --git a/input/weight b/input/weight index a1f7a8d..dc1936c 100644 --- a/input/weight +++ b/input/weight @@ -1,27 +1,43 @@ -1 5 +1 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 -2 5 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 +2 9 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 33.8700000 0.0060680 + 2 5.0950000 0.0453080 + 3 1.1590000 0.2028220 S 1 - 1 0.1220000 1.0000000 + 1 0.3258000 1.0000000 S 1 - 1 0.0297400 1.0000000 + 1 0.1027000 1.0000000 +S 1 + 1 0.0252600 1.0000000 P 1 - 1 0.7270000 1.0000000 + 1 1.4070000 1.0000000 P 1 - 1 0.1410000 1.0000000 + 1 0.3880000 1.0000000 +P 1 + 1 0.1020000 1.0000000 +D 1 + 1 1.0570000 1.0000000 +D 1 + 1 0.2470000 1.0000000 diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index 5bd67f9..7df8b16 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -1,5 +1,5 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,EwGIC,c) + 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 @@ -80,7 +80,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri ! Output variables double precision,intent(out) :: Ew - double precision,intent(out) :: EwGIC ! Hello world @@ -343,6 +342,6 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGri call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), & - J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:)) + J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:)) end subroutine GOK_RKS diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index 7708b1c..7a9f80a 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -35,10 +35,10 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight ! Local variables integer :: iEns - double precision :: EwZW,EwGICZW - double precision :: EwEW,EwGICEW + double precision :: EwZW + double precision :: EwEW double precision :: wLIM(nEns) - double precision :: Om(nEns),OmGIC(nEns) + double precision :: Om(nEns) ! Hello world @@ -67,7 +67,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight write(*,*) call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,EwGICZW,c) + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,c) !------------------------------------------------------------------------ ! Equiensemble calculation @@ -84,20 +84,15 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight write(*,*) call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW,c) + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,c) !------------------------------------------------------------------------ ! LIM excitation energies !------------------------------------------------------------------------ - Om(1) = 0d0 - OmGIC(1) = 0d0 - - Om(2) = 0d0 - OmGIC(2) = 0d0 + Om(:) = 0d0 if(wEns(2) > 10d-3) then - Om(2) = (EwEW - EwZW)/wEns(2) - OmGIC(2) = (EwGICEW - EwGICZW)/wEns(2) + Om(2) = (EwEW - EwZW)/wEns(2) end if write(*,'(A60)') '-------------------------------------------------' @@ -105,20 +100,12 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight write(*,'(A60)') '-------------------------------------------------' write(*,'(A44,F16.10,A3)') ' Zero-weight ensemble energy',EwZW, ' au' - write(*,'(A44,F16.10,A3)') ' Zero-weight GIC ensemble energy',EwGICZW,' au' - write(*,*) write(*,'(A44,F16.10,A3)') ' Equi-weight ensemble energy',EwEW, ' au' - write(*,'(A44,F16.10,A3)') ' Equi-weight GIC ensemble energy',EwGICEW, ' au' write(*,'(A60)') '-------------------------------------------------' do iEns=2,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' end do - write(*,*) - do iEns=2,nEns - write(*,'(A40,I2,A2,F16.10,A3)') ' GIC Excitation energy 1 ->',iEns,': ',OmGIC(iEns), ' au' - write(*,'(A40,I2,A2,F16.10,A3)') ' GIC Excitation energy 1 ->',iEns,': ',OmGIC(iEns)*HaToeV,' eV' - end do write(*,'(A60)') '-------------------------------------------------' diff --git a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 index bc96f43..5f8e758 100644 --- a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 @@ -43,15 +43,15 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho ! Parameters for stretch H2 - a = + 0.01922622507087411d0 - b = - 0.01799647558018601d0 - c = - 0.022945430666782573d0 +! a = + 0.01922622507087411d0 +! b = - 0.01799647558018601d0 +! c = - 0.022945430666782573d0 ! Parameters for He -! a = 1.9125735895875828d0 -! b = 2.715266992840757d0 -! c = 2.1634223380633086d0 + a = 1.9125735895875828d0 + b = 2.715266992840757d0 + c = 2.1634223380633086d0 w = wEns(2) dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) diff --git a/src/eDFT/RGIC_lda_exchange_energy.f90 b/src/eDFT/RGIC_lda_exchange_energy.f90 index 4c16b47..fb60b5d 100644 --- a/src/eDFT/RGIC_lda_exchange_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_energy.f90 @@ -35,15 +35,15 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Parameters for stretch H2 - a = + 0.01922622507087411d0 - b = - 0.01799647558018601d0 - c = - 0.022945430666782573d0 +! a = + 0.01922622507087411d0 +! b = - 0.01799647558018601d0 +! c = - 0.022945430666782573d0 ! Parameters for He -! a = 1.9125735895875828d0 -! b = 2.715266992840757d0 -! c = 2.1634223380633086d0 + a = 1.9125735895875828d0 + b = 2.715266992840757d0 + c = 2.1634223380633086d0 w = wEns(2) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) diff --git a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 index b1730bf..e2b7b42 100644 --- a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 @@ -38,15 +38,15 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E ! Parameters for stretch H2 - a = + 0.01922622507087411d0 - b = - 0.01799647558018601d0 - c = - 0.022945430666782573d0 +! a = + 0.01922622507087411d0 +! b = - 0.01799647558018601d0 +! c = - 0.022945430666782573d0 ! Parameters for He -! a = 1.9125735895875828d0 -! b = 2.715266992840757d0 -! c = 2.1634223380633086d0 + a = 1.9125735895875828d0 + b = 2.715266992840757d0 + c = 2.1634223380633086d0 w = wEns(2) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) diff --git a/src/eDFT/RGIC_lda_exchange_potential.f90 b/src/eDFT/RGIC_lda_exchange_potential.f90 index 468ef27..ad6e66f 100644 --- a/src/eDFT/RGIC_lda_exchange_potential.f90 +++ b/src/eDFT/RGIC_lda_exchange_potential.f90 @@ -38,15 +38,15 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Parameters for stretch H2 - a = + 0.01922622507087411d0 - b = - 0.01799647558018601d0 - c = - 0.022945430666782573d0 +! a = + 0.01922622507087411d0 +! b = - 0.01799647558018601d0 +! c = - 0.022945430666782573d0 ! Parameters for He -! a = 1.9125735895875828d0 -! b = 2.715266992840757d0 -! c = 2.1634223380633086d0 + a = 1.9125735895875828d0 + b = 2.715266992840757d0 + c = 2.1634223380633086d0 w = wEns(2) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 993e8c8..4a6aec6 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,Ew,EwGIC + double precision :: ENuc,Ew double precision,allocatable :: ZNuc(:),rNuc(:,:) @@ -164,7 +164,7 @@ program eDFT call cpu_time(start_KS) call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & - S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,EwGIC,c(:,:)) + S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,Ew,c(:,:)) 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 376a94d..295fde8 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,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & +subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD) ! Print individual energies for eDFT calculation @@ -11,7 +11,6 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec, integer,intent(in) :: nEns double precision,intent(in) :: ENuc double precision,intent(in) :: Ew - double precision,intent(in) :: EwGIC double precision,intent(in) :: ET(nEns) double precision,intent(in) :: EV(nEns) double precision,intent(in) :: EJ(nEns) @@ -36,7 +35,6 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET,EV,EJ,Ex,Ec, write(*,'(A60)') ' ENSEMBLE ENERGIES' write(*,'(A60)') '-------------------------------------------------' write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au' - write(*,'(A44,F16.10,A3)') ' GIC Ensemble energy: ',EwGIC + ENuc,' au' write(*,'(A60)') '-------------------------------------------------' write(*,*) diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index e501ce3..48c6bfd 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,LDA_centered,nEns,wEns,nGrid,weight,nBas, & - nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,EwGIC,E,Om) + 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 @@ -55,7 +55,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n ! Output variables - double precision,intent(out) :: EwGIC double precision,intent(out) :: E(nEns) double precision,intent(out) :: Om(nEns) @@ -129,15 +128,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n + Ex(iEns) + Ec(iEns) + ExcDD(iEns) end do -!------------------------------------------------------------------------ -! Total energy with ghost-interaction correction -!------------------------------------------------------------------------ - - EwGIC = 0d0 - do iEns=1,nEns - EwGIC = EwGIC + wEns(iEns)*E(iEns) - end do - !------------------------------------------------------------------------ ! Excitation energies !------------------------------------------------------------------------ @@ -162,7 +152,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n ! Dump results !------------------------------------------------------------------------ - call print_restricted_individual_energy(nEns,ENuc,Ew,EwGIC,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & + call print_restricted_individual_energy(nEns,ENuc,Ew,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), & Eaux(:),ExDD(:),EcDD(:),ExcDD(:),E(:), & Om(:),Omx(:),Omc(:),Omxc(:),Omaux,OmxDD(:),OmcDD(:),OmxcDD(:))