PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

LYP done

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-02-13 22:31:46 +01:00
parent e505d6d1a2
commit 337e801ecd
3 changed files with 15 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/dft
View File

@ -6,7 +6,7 @@
# GGA = 2: B88,G96,PBE # GGA = 2: B88,G96,PBE
# Hybrid = 4: B3LYP,PBE0 # Hybrid = 4: B3LYP,PBE0
# Hartree-Fock = 666 # Hartree-Fock = 666
2 B88 1 S51
# correlation rung: # correlation rung:
# Hartree = 0: H # Hartree = 0: H
# LDA = 1: VWN5,eVWN5 # LDA = 1: VWN5,eVWN5
@ -19,8 +19,8 @@
# Number of states in ensemble (nEns) # Number of states in ensemble (nEns)
3 3
# occupation numbers # occupation numbers
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 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
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 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

6
src/eDFT/ULYP_gga_correlation_energy.f90
View File

@ -47,10 +47,10 @@ subroutine ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec)
if(r > threshold) then if(r > threshold) then
ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**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)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 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) 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) 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)) delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0))

19
src/eDFT/ULYP_gga_correlation_potential.f90
View File

@ -55,27 +55,27 @@ subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
if(r > threshold) then if(r > threshold) then
ga = drho(1,iG,1)**2 + drho(2,iG,1)**2 + drho(3,iG,1)**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)**2 + drho(2,iG,2)**2 + drho(3,iG,2)**2 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) 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) 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)) delta = c*r**(-1d0/3d0) + d*r**(-1d0/3d0)/(1d0 + d*r**(-1d0/3d0))
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) 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 dodrb = dodra
dddra = - c/3d0*r**(-4d0/3d0) & dddra = - c/3d0*r**(-4d0/3d0) &
+ d**2/(3d0*(1d0 + d*r**(-1d0/3d0))**2)*r**(-5d0/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 dddrb = dddra
dfdra = - 4d0*a/(1d0 + d*r**(-1d0/3d0))*rb/r & 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/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*( & - a*b*omega*rb*( &
+ 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) & + 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g & + (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 & 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/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*( & - a*b*omega*ra*( &
+ 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) & + 2d0**(11d0/3d0)*Cf*(ra**(8d0/3d0) + rb**(8d0/3d0)) &
+ (47d0/18d0 - 7d0*delta/18d0)*g & + (47d0/18d0 - 7d0*delta/18d0)*g &
- (5d0/2d0 - delta/18d0)*(ga + gb) & - (5d0/2d0 - delta/18d0)*(ga + gb) &
- (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) & - (delta - 11d0)/9d0*(ra/r*ga + rb/r*gb) &
- 4d0/3d0*r/ra*g & - 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 ) & + 4d0/3d0*r/ra*gb ) &
- a*b*omega*ra*rb*( & - a*b*omega*ra*rb*( &
+ 8d0/3d0*2d0**(11d0/3d0)*Cf*rb**(5d0/3d0) & + 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 ) + (2d0*r**2/3d0 - rb**2)*ga )
Fc(mu,nu,1) = Fc(mu,nu,1) + vAO*dfdra Fc(mu,nu,1) = Fc(mu,nu,1) + vAO*dfdra
Fc(mu,nu,2) = Fc(mu,nu,2) + vAO*dfdrb 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)) & 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 gbAO = weight(iG)*gbAO
dfdga = -a*b*omega*(-rb**2 + ra*rb*(1d0/9d0 - (delta-11d0)/9d0*ra/r - delta/3d0)) 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)) 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 Fc(mu,nu,1) = Fc(mu,nu,1) + 2d0*gaAO*dfdga + gbAO*dfdgab