diff --git a/input/dft b/input/dft index 4743c0f..19bf050 100644 --- a/input/dft +++ b/input/dft @@ -6,7 +6,7 @@ # GGA = 2: B88,G96,PBE # Hybrid = 4: B3LYP,PBE0 # Hartree-Fock = 666 - 2 B88 + 1 S51 # correlation rung: # Hartree = 0: H # LDA = 1: VWN5,eVWN5 @@ -19,8 +19,8 @@ # Number of states in ensemble (nEns) 3 # occupation numbers - 1 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 + 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 + 1 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 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/src/eDFT/ULYP_gga_correlation_energy.f90 b/src/eDFT/ULYP_gga_correlation_energy.f90 index 974d860..e0788b8 100644 --- a/src/eDFT/ULYP_gga_correlation_energy.f90 +++ b/src/eDFT/ULYP_gga_correlation_energy.f90 @@ -47,10 +47,10 @@ subroutine ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) if(r > threshold) then - ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**2 - gb = drho(1,iG,2)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 + ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1) + gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2) gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2) - g = ga + gab + gb + g = ga + 2d0*gab + gb omega = exp(-c*r**(-1d0/3d0))/(1d0 + d*r**(-1d0/3d0))*r**(-11d0/3d0) delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0)) diff --git a/src/eDFT/ULYP_gga_correlation_potential.f90 b/src/eDFT/ULYP_gga_correlation_potential.f90 index 4eac161..702aba3 100644 --- a/src/eDFT/ULYP_gga_correlation_potential.f90 +++ b/src/eDFT/ULYP_gga_correlation_potential.f90 @@ -55,27 +55,27 @@ subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) if(r > threshold) then - ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**2 - gb = drho(1,iG,2)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 + ga = drho(1,iG,1)*drho(1,iG,1) + drho(2,iG,1)*drho(2,iG,1) + drho(3,iG,1)*drho(3,iG,1) + gb = drho(1,iG,2)*drho(1,iG,2) + drho(2,iG,2)*drho(2,iG,2) + drho(3,iG,2)*drho(3,iG,2) gab = drho(1,iG,1)*drho(1,iG,2) + drho(2,iG,1)*drho(2,iG,2) + drho(3,iG,1)*drho(3,iG,2) - g = ga + gab + gb + g = ga + 2d0*gab + gb omega = exp(-c*r**(-1d0/3d0))/(1d0 + d*r**(-1d0/3d0))*r**(-11d0/3d0) delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0)) vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) - dodra = (d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0)) + c/(3d0*r**(4d0/3d0)) - 11d0/(3d0*r)))*omega + dodra = (d/(3d0*r**(4d0/3d0)*(1d0 + d*r**(-1d0/3d0))) + c/(3d0*r**(4d0/3d0)) - 11d0/(3d0*r))*omega dodrb = dodra dddra = - c/3d0*r**(-4d0/3d0) & + d**2/(3d0*(1d0 + d*r**(-1d0/3d0))**2)*r**(-5d0/3d0) & - - d/(3d0*(1d0 + d*r**(-1d0/3d0))*r**(-4d0/3d0)) + - d/(3d0*(1d0 + d*r**(-1d0/3d0)))*r**(-4d0/3d0) dddrb = dddra dfdra = - 4d0*a/(1d0 + d*r**(-1d0/3d0))*rb/r & - 4d0/3d0*a*d/(1d0 + d*r**(-1d0/3d0))**2*ra*rb/r**(7d0/3d0) & - + 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r & + + 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r**2 & - a*b*omega*rb*( & + 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) & + (47d0/18d0 - 7d0*delta/18d0)*g & @@ -102,14 +102,14 @@ subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) dfdrb = - 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra/r & - 4d0/3d0*a*d/(1d0 + d*r**(-1d0/3d0))**2*ra*rb/r**(7d0/3d0) & - + 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r & + + 4d0*a/(1d0 + d*r**(-1d0/3d0))*ra*rb/r**2 & - a*b*omega*ra*( & + 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) & + (47d0/18d0 - 7d0*delta/18d0)*g & - (5d0/2d0 - delta/18d0)*(ga + gb) & - (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) & - 4d0/3d0*r/ra*g & - + (4d0/3d0*r/ra - 2d0*ra/rb)*ga & + + (4d0/3d0*r/ra - 2d0*rb/ra)*ga & + 4d0/3d0*r/ra*gb ) & - a*b*omega*ra*rb*( & + 8d0/3d0*2d0**(11d0/3d0)*Cf*rb**(5d0/3d0) & @@ -128,7 +128,6 @@ subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + (2d0*r**2/3d0 - rb**2)*ga ) Fc(mu,nu,1) = Fc(mu,nu,1) + vAO*dfdra - Fc(mu,nu,2) = Fc(mu,nu,2) + vAO*dfdrb gaAO = drho(1,iG,1)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & @@ -142,7 +141,7 @@ subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) gbAO = weight(iG)*gbAO dfdga = -a*b*omega*(-rb**2 + ra*rb*(1d0/9d0 - (delta-11d0)/9d0*ra/r - delta/3d0)) - dfdgab = -a*b*omega*(-2d0/3d0*r**2 + ra*rb*(47d0/18d0 - 7d0*delta/18d0)) + dfdgab = -a*b*omega*(-4d0/3d0*r**2 + 2d0*ra*rb*(47d0/18d0 - 7d0*delta/18d0)) dfdgb = -a*b*omega*(-ra**2 + ra*rb*(1d0/9d0 - (delta-11d0)/9d0*rb/r - delta/3d0)) Fc(mu,nu,1) = Fc(mu,nu,1) + 2d0*gaAO*dfdga + gbAO*dfdgab