diff --git a/input/dft b/input/dft index 2b261ca..f6badeb 100644 --- a/input/dft +++ b/input/dft @@ -1,7 +1,7 @@ # Restricted or unrestricted KS calculation GOK-RKS # exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666 - 1 RS51 + 666 HF # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666 0 HF # quadrature grid SG-n diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index 3d3f331..2238f8c 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) +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) ! Perform restricted Kohn-Sham calculation for ensembles @@ -206,15 +206,15 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres rhow(:) = rhow(:) + wEns(iEns)*rho(:,iEns) end do - if(xc_rung > 1 .and. xc_rung /= 666) then + if(xc_rung > 1) then -! Ground state density +! Compute gradient of the one-electron density do iEns=1,nEns call gradient_density(nGrid,nBas,P(:,:,iEns),AO(:,:),dAO(:,:,:),drho(:,:,iEns)) end do -! Weight-dependent one-electron density +! Weight-dependent one-electron density gradient drhow(:,:) = 0d0 do iEns=1,nEns @@ -234,11 +234,12 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute correlation potential - call correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:)) + call restricted_correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), & + nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:)) ! Build Fock operator - F(:,:) = Hc(:,:) + J(:,:) + J(:,:) + Fx(:,:) + Fc(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*Fx(:,:) + Fc(:,:) ! Check convergence @@ -275,10 +276,11 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex) + Ex = 0.5d0*Ex ! Correlation energy - call correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec) + call restricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec) ! Total energy @@ -321,9 +323,9 @@ subroutine GOK_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call 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(:,:,:),E(:),Om(:)) + 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(:,:,:),E(:),Om(:)) end subroutine GOK_RKS diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 136d660..3312ce8 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -138,7 +138,7 @@ program eDFT 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(:),nV(:),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS) + nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1),S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,EKS) call cpu_time(end_KS) t_KS = end_KS - start_KS diff --git a/src/eDFT/gradient_density.f90 b/src/eDFT/gradient_density.f90 index 7d9e0cd..0e8c8f2 100644 --- a/src/eDFT/gradient_density.f90 +++ b/src/eDFT/gradient_density.f90 @@ -13,7 +13,7 @@ subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho) integer,intent(in) :: nBas double precision,intent(in) :: P(nBas,nBas) double precision,intent(in) :: AO(nBas,nGrid) - double precision,intent(in) :: dAO(3,nBas,nGrid) + double precision,intent(in) :: dAO(ncart,nBas,nGrid) ! Local variables @@ -22,13 +22,13 @@ subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho) ! Output variables - double precision,intent(out) :: drho(3,nGrid) + double precision,intent(out) :: drho(ncart,nGrid) drho(:,:) = 0d0 do iG=1,nGrid do mu=1,nBas do nu=1,nBas - do ixyz=1,3 + do ixyz=1,ncart drho(ixyz,iG) = drho(ixyz,iG) & + P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG)) enddo @@ -36,10 +36,10 @@ subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho) enddo enddo - do iG=1,nGrid - do ixyz=1,3 - if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh - enddo - enddo +! do iG=1,nGrid +! do ixyz=1,ncart +! if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh +! enddo +! enddo end subroutine gradient_density diff --git a/src/eDFT/print_RKS.f90 b/src/eDFT/print_RKS.f90 index 1cc19c1..346d8b6 100644 --- a/src/eDFT/print_RKS.f90 +++ b/src/eDFT/print_RKS.f90 @@ -43,7 +43,6 @@ subroutine print_RKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) write(*,'(A40,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex + Ec,' au' write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',EJ,' au' write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au' - write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au' write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',Ec,' au' write(*,'(A60)') '-------------------------------------------------' write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',Ew,' au' diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90 index 0f1ea16..a5add09 100644 --- a/src/eDFT/restricted_density_matrix.f90 +++ b/src/eDFT/restricted_density_matrix.f90 @@ -29,7 +29,12 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P) ! Doubly-excited state density matrix iEns = 2 - P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) & - + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1))) + if(nO > 1) then + + P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) + + end if + + P(:,:,iEns) = P(:,:,iEns) + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1))) end subroutine restricted_density_matrix