Dear wien2k users,
In the OMP version of lapw0 (with more than one thread in parallel) the
correlation energy and potential of the PW91 GGA-functional is wrong.
There are no problems with the standard functionals (LDA, PBE, PBESOL,
WC, mBJ) and even for PW91 the problem appears only for OMP_NUM_THREADS > 1.
The 2 attached subroutines fix the problem.
Regards
--
P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300 FAX: +43-1-58801-165982
Email: bl...@theochem.tuwien.ac.at WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------
SUBROUTINE VPW91(FU,GXU,GYU,GZU,GMAGU,G2U,GGXU,GGYU,GGZU, &
FD,GXD,GYD,GZD,GMAGD,G2D,GGXD,GGYD,GGZD,vxupw,vcupw, & !
vxdpw,vcdpw,exupw,ecupw,exdpw,ecdpw,GDGGT) !
! VXU,VXD, & !THORUH, 31.Jänner
! exu,exd,GDGGT) !
IMPLICIT REAL*8(A-H,O-Z)
!! save
! COMMON /GAS/ FK,SK,G,EC,ECRS,ECZET
DATA THRD,FTDS/0.3333333333333333D0,1.3333333333333333D0/
DATA TTDS/0.6666666666666667/
PI=4.D0*ATAN(1.D0)
PI2=PI*PI
! EXCHANGE - UP
D=2.D0*FU
GDGG=GXU*GGXU+GYU*GGYU+GZU*GGZU
TKF=2.D0*(3.D0*PI2*D)**THRD
S=GMAGU/(TKF*FU)
U=GDGG/(TKF*(FU*TKF)**2)
V=G2U/(FU*TKF**2)
CALL EXCH(D,S,U,V,EXUP,VXUP)
! EXCHANGE - DN
D=2.D0*FD
GDGG=GXD*GGXD+GYD*GGYD+GZD*GGZD
TKF=2.D0*(3.D0*PI2*D)**THRD
S=GMAGD/(TKF*FD)
U=GDGG/(TKF*(FD*TKF)**2)
V=G2D/(FD*TKF**2)
CALL EXCH(D,S,U,V,EXDN,VXDN)
! LOCAL CORRELATION
D=FU+FD
RS=(0.75D0/(PI*D))**THRD
ZET=(FU-FD)/D
CALL CORLSD(RS,ZET,ECLSD,VLSDU,VLSDD,ECRS,ECZET,ALFC)
EC=ECLSD
! GRADIENT CORRECTION TO CORRELATION
G=0.5D0*((1.D0+ZET)**TTDS+(1.D0-ZET)**TTDS)
GP=GXU*GXD+GYU*GYD+GZU*GZD
GDGZET=((1.D0-ZET)*GMAGU**2-(1.D0+ZET)*GMAGD**2-2.D0*ZET*GP)/D
FK=(3.D0*PI2*D)**THRD
SK=SQRT(4.D0*FK/PI)
TKSG=2.D0*G*SK
T=SQRT((GXU+GXD)**2+(GYU+GYD)**2+(GZU+GZD)**2)/(D*TKSG)
U=GDGGT/(TKSG*(TKSG*D)**2)
V=(G2U+G2D)/(D*TKSG**2)
W=GDGZET/(D*TKSG**2)
CALL CORGGA(RS,ZET,T,U,V,W,ECGGA,VGGAU,VGGAD,FK,SK,G,EC,ECRS,ECZET)
!
! VXU=2.D0*(VXUP+VLSDU+VGGAU) !
! VXD=2.D0*(VXDN+VLSDD+VGGAD) !
! exu=2.d0*(exup+eclsd+ecgga) !
! exd=2.d0*(exdn+eclsd+ecgga) !
exupw=2.d0*exup !
ecupw=2.d0*(eclsd+ecgga) !
exdpw=2.d0*exdn ! changed by THORUH, 31.Jänner
ecdpw=2.d0*(eclsd+ecgga) !
vxupw=2.d0*vxup !
vcupw=2.d0*(vlsdu+vggau) !
vxdpw=2.d0*vxdn !
vcdpw=2.d0*(vlsdd+vggad) !
! exu=2.d0*(exup+0.5d0*eclsd+0.5d0*ecgga)
! exd=2.d0*(exdn+0.5d0*eclsd+0.5d0*ecgga)
!
RETURN
END
SUBROUTINE CORGGA(RS,ZET,T,UU,VV,WW,H,DVCUP,DVCDN,fk,sk,g,ec,ecrs,eczet)
! GGA91 CORRELATION
! INPUT RS: SEITZ RADIUS
! INPUT ZET: RELATIVE SPIN POLARIZATION
! INPUT T: ABS(GRAD D)/(D*2.*KS*G)
! INPUT UU: (GRAD D)*GRAD(ABS(GRAD D))/(D**2 * (2*KS*G)**3)
! INPUT VV: (LAPLACIAN D)/(D * (2*KS*G)**2)
! INPUT WW: (GRAD D)*(GRAD ZET)/(D * (2*KS*G)**2
! OUTPUT H: NONLOCAL PART OF CORRELATION ENERGY PER ELECTRON
! OUTPUT DVCUP,DVCDN: NONLOCAL PARTS OF CORRELATION POTENTIALS
IMPLICIT REAL*8 (A-H,O-Z)
! COMMON/GAS/FK,SK,G,EC,ECRS,ECZET
DATA XNU,CC0,CX,ALF/15.75592D0,0.004235D0,-0.001667212D0,0.09D0/
DATA C1,C2,C3,C4/0.002568D0,0.023266D0,7.389D-6,8.723D0/
DATA C5,C6,A4/0.472D0,7.389D-2,100.D0/
DATA THRDM,THRD2/-0.333333333333D0,0.666666666667D0/
BET = XNU*CC0
DELT = 2.D0*ALF/BET
G3 = G**3
G4 = G3*G
PON = -DELT*EC/(G3*BET)
B = DELT/(DEXP(PON)-1.D0)
B2 = B*B
T2 = T*T
T4 = T2*T2
T6 = T4*T2
RS2 = RS*RS
RS3 = RS2*RS
Q4 = 1.D0+B*T2
Q5 = 1.D0+B*T2+B2*T4
Q6 = C1+C2*RS+C3*RS2
Q7 = 1.D0+C4*RS+C5*RS2+C6*RS3
CC = -CX + Q6/Q7
R0 = (SK/FK)**2
R1 = A4*R0*G4
COEFF = CC-CC0-3.D0*CX/7.D0
R2 = XNU*COEFF*G3
R3 = DEXP(-R1*T2)
H0 = G3*(BET/DELT)*DLOG(1.D0+DELT*Q4*T2/Q5)
H1 = R3*R2*T2
H = H0+H1
! LOCAL CORRELATION OPTION:
! H = 0.0D0
! ENERGY DONE. NOW THE POTENTIAL:
CCRS = (C2+2.*C3*RS)/Q7 - Q6*(C4+2.*C5*RS+3.*C6*RS2)/Q7**2
RSTHRD = RS/3.D0
R4 = RSTHRD*CCRS/COEFF
GZ = ((1.D0+ZET)**THRDM - (1.D0-ZET)**THRDM)/3.D0
FAC = DELT/B+1.D0
BG = -3.D0*B2*EC*FAC/(BET*G4)
BEC = B2*FAC/(BET*G3)
Q8 = Q5*Q5+DELT*Q4*Q5*T2
Q9 = 1.D0+2.D0*B*T2
Q82 = Q8*Q8
H0B = -BET*G3*B*T6*(2.D0+B*T2)/Q8
H0RS = -RSTHRD*H0B*BEC*ECRS
FACT0 = 2.D0*DELT-6.D0*B
FACT1 = Q5*Q9+Q4*Q9*Q9
H0BT = 2.D0*BET*G3*T4*(Q4*Q5*FACT0-DELT*FACT1)/Q82
H0RST = RSTHRD*T2*H0BT*BEC*ECRS
H0Z = 3.D0*GZ*H0/G + H0B*(BG*GZ+BEC*ECZET)
H0T = 2.*BET*G3*Q9/Q8
H0ZT = 3.D0*GZ*H0T/G+H0BT*(BG*GZ+BEC*ECZET)
FACT2 = Q4*Q5+B*T2*(Q4*Q9+Q5)
FACT3 = 2.D0*B*Q5*Q9+DELT*FACT2
H0TT = 4.D0*BET*G3*T*(2.D0*B/Q8-Q9*FACT3/Q82)
H1RS = R3*R2*T2*(-R4+R1*T2/3.D0)
FACT4 = 2.D0-R1*T2
H1RST = R3*R2*T2*(2.D0*R4*(1.D0-R1*T2)-THRD2*R1*T2*FACT4)
H1Z = GZ*R3*R2*T2*(3.D0-4.D0*R1*T2)/G
H1T = 2.D0*R3*R2*(1.D0-R1*T2)
H1ZT = 2.D0*GZ*R3*R2*(3.D0-11.D0*R1*T2+4.D0*R1*R1*T4)/G
H1TT = 4.D0*R3*R2*R1*T*(-2.D0+R1*T2)
HRS = H0RS+H1RS
HRST = H0RST+H1RST
HT = H0T+H1T
HTT = H0TT+H1TT
HZ = H0Z+H1Z
HZT = H0ZT+H1ZT
COMM = H+HRS+HRST+T2*HT/6.D0+7.D0*T2*T*HTT/6.D0
PREF = HZ-GZ*T2*HT/G
FACT5 = GZ*(2.D0*HT+T*HTT)/G
COMM = COMM-PREF*ZET-UU*HTT-VV*HT-WW*(HZT-FACT5)
DVCUP = COMM + PREF
DVCDN = COMM - PREF
! LOCAL CORRELATION OPTION:
! DVCUP = 0.0D0
! DVCDN = 0.0D0
RETURN
END
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
