diff --git a/input/dft b/input/dft index 8059451..b528057 100644 --- a/input/dft +++ b/input/dft @@ -1,5 +1,5 @@ # Restricted or unrestricted KS calculation - eDFT-UKS + UKS # exchange rung: # Hartree = 0: H # LDA = 1: S51,CC-S51 @@ -28,7 +28,7 @@ 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 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 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.5 0.0 + 0.25 0.0 # N-centered? F # Parameters for CC weight-dependent exchange functional diff --git a/input/methods b/input/methods index 24288c1..768c20d 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F T F F + F F T F # MP2* MP3 MP2-F12 F F F # CCD DCD CCSD CCSD(T) @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* - F F T + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/mol/h2.xyz b/mol/h2.xyz index 85810e4..f8e2dab 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0.0 0.0 0.0 -H 0.0 0.0 2.0 +H 0.0 0.0 0.71 diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 387e44a..f4a5b40 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -226,8 +226,8 @@ 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,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type,mix, & - nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) + call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,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 cpu_time(end_KS) t_KS = end_KS - start_KS @@ -243,8 +243,8 @@ 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,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type,mix, & - nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) + call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,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 cpu_time(end_KS) t_KS = end_KS - start_KS diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index d0bf4e0..03a9212 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -1,5 +1,5 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & - nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw,Vxc) + nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eps,c,Pw,Vxc) ! Perform unrestricted Kohn-Sham calculation for ensembles @@ -23,13 +23,18 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig double precision,intent(in) :: AO(nBas,nGrid) double precision,intent(in) :: dAO(ncart,nBas,nGrid) + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ENuc + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: occnum(nBas,nspin,nEns) integer,intent(in) :: Cx_choice logical,intent(in) :: doNcentered @@ -48,6 +53,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig double precision :: EJ(nsp) double precision :: Ex(nspin) double precision :: Ec(nsp) + double precision :: dipole(ncart) double precision,allocatable :: cp(:,:,:) double precision,allocatable :: J(:,:,:) @@ -132,6 +138,11 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Guess coefficients and eigenvalues + nO(:) = 0 + do ispin=1,nspin + nO(ispin) = int(sum(occnum(:,ispin,1))) + end do + if(guess_type == 1) then do ispin=1,nspin @@ -142,11 +153,6 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Mix guess to enforce symmetry breaking - nO(:) = 0 - do ispin=1,nspin - nO(ispin) = int(sum(occnum(:,ispin,1))) - end do - if(mix) call mix_guess(nBas,nO,c) else if(guess_type == 2) then @@ -385,7 +391,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Compute final KS energy - call print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) + call dipole_moment(nBas,Pw(:,:,1)+Pw(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_UKS(nBas,nEns,nO,S,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew,dipole) ! Compute Vxc for post-HF calculations diff --git a/src/eDFT/print_UKS.f90 b/src/eDFT/print_UKS.f90 index f063849..be30958 100644 --- a/src/eDFT/print_UKS.f90 +++ b/src/eDFT/print_UKS.f90 @@ -1,4 +1,4 @@ -subroutine print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) +subroutine print_UKS(nBas,nEns,nO,Ov,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew,dipole) ! Print one- and two-electron energies and other stuff for KS calculation @@ -9,7 +9,8 @@ subroutine print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) integer,intent(in) :: nBas integer,intent(in) :: nEns - double precision,intent(in) :: occnum(nBas,nspin,nEns) + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: Ov(nBas,nBas) double precision,intent(in) :: wEns(nEns) double precision,intent(in) :: eps(nBas,nspin) double precision,intent(in) :: c(nBas,nBas,nspin) @@ -20,40 +21,39 @@ subroutine print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: Ec(nsp) double precision,intent(in) :: Ew + double precision,intent(in) :: dipole(ncart) ! Local variables + integer :: ixyz integer :: ispin integer :: iEns integer :: iBas integer :: HOMO(nspin) integer :: LUMO(nspin) double precision :: Gap(nspin) - double precision :: nEl(nspin) - -! Number of electrons in the ensemble - - nEl(:) = 0d0 - do ispin=1,nspin - do iEns=1,nEns - do iBas=1,nBas - nEl(ispin) = nEl(ispin) + wEns(iEns)*occnum(iBas,ispin,iEns) - end do - end do - end do + double precision :: S_exact,S2_exact + double precision :: S,S2 ! HOMO and LUMO do ispin=1,nspin - HOMO(ispin) = ceiling(nEl(ispin)) + HOMO(ispin) = nO(ispin) LUMO(ispin) = HOMO(ispin) + 1 Gap(ispin) = eps(LUMO(ispin),ispin) - eps(HOMO(ispin),ispin) end do -! Dump results +! Spin comtamination + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2) + + S_exact = 0.5d0*dble(nO(1) - nO(2)) + S = -0.5d0 + 0.5d0*sqrt(1d0 + 4d0*S2) + +! Dump results write(*,*) write(*,'(A60)') '-------------------------------------------------' @@ -97,6 +97,15 @@ subroutine print_UKS(nBas,nEns,occnum,wEns,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) write(*,'(A40,F13.6,A3)') ' KS LUMO b energy:',eps(LUMO(2),2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' KS HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.6)') ' S (exact) :',2d0*S_exact + 1d0 + write(*,'(A40,1X,F16.6)') ' S :',2d0*S + 1d0 + write(*,'(A40,1X,F16.6)') ' (exact) :',S2_exact + write(*,'(A40,1X,F16.6)') ' :',S2 + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A45)') ' Dipole moment (Debye) ' + write(*,'(19X,4A10)') 'X','Y','Z','Tot.' + write(*,'(19X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A60)') '-------------------------------------------------' write(*,*) ! Print results