""" This module defines class to compute
exchange-correlation potential and
exchange-correlation energy within LDA

class name:
    ExchangeCorrelation(type)  --     type can be 1,2,3,4

Be careful: All energies and potentials are in units of Hartree = 2*13.6058=2*Ry
If you work in Rydbergs, you have to multiply the potential by factor of 2.
"""

from math import *

class  ExchangeCorrelation(object):
    """******************************************************************************/
    Calculates Exchange&Correlation Energy and Potential                       */ 
    type=0 - due to U.von.Barth and L.Hedin, J.Phys.C5, 1629 (1972)            */
    type=1 - O.E.Gunnarsson and S.Lundqvist,  Phys.Rev.B                       */
    type=2 - V.L.Moruzzi, J.F.Janak, and A.R.Williams, Calculated              */
             Electronic Properties of Metals (New York, Pergamon Press, 1978)  */
    type=3 - S.H.Vosko, L.Wilk, and M.Nusair, Can.J.Phys.58, 1200 (1980)       */
    type=4 - Correlation of Perdew and Wang 1991                               */
    ******************************************************************************/
    """
    def __init__(self, type_=3):
        self.type = type_
        self.alphax = 0.610887057710857 #//(3/(2 Pi))^(2/3)
        self.Aw = 0.0311 
        self.Bw = -0.048 
        self.Cw = 0.002 
        self.D  = -0.0116 
        self.gamma  = -0.1423 
        self.beta1  =  1.0529 
        self.beta2  =  0.3334 
        self.Ap  =  0.0621814 
        self.xp0 = -0.10498 
        self.bp  =  3.72744 
        self.cp  =  12.9352 
        self.Qp  =  6.1519908 
        self.cp1 =  1.2117833 
        self.cp2 =  1.1435257 
        self.cp3 = -0.031167608 
        if self.type==0: self.C = 0.0504; self.A = 30
        if self.type==1: self.C = 0.0666; self.A = 11.4
        if self.type==2: self.C = 0.045;  self.A = 21

    def Vx(self, rs): # Vx
        return -self.alphax/rs
    
    def ExVx(self, rs): # Ex-Vx
        return 0.25*self.alphax/rs
    
    def Ex(self, rs):
        return -0.75*self.alphax/rs
    
    def Vc(self, rs): # Vc
        if (self.type<3):
            x = rs/self.A;
            return -0.5*self.C*log(1+1./x)
        elif(self.type<4): # type=3 WVN
            x=sqrt(rs)
            xpx= x*x + self.bp*x + self.cp
            atnp = atan(self.Qp/(2*x+self.bp))
            ecp = 0.5*self.Ap*(log(x*x/xpx)+self.cp1*atnp-self.cp3*(log((x-self.xp0)**2/xpx)+self.cp2*atnp))
            return ecp - self.Ap/6.*(self.cp*(x-self.xp0)-self.bp*x*self.xp0)/((x-self.xp0)*xpx)
        else:
            if rs>1:
                return self.gamma/(1+self.beta1*sqrt(rs)+self.beta2*rs)*(1+7/6.*self.beta1*sqrt(rs)+self.beta2*rs)/(1+self.beta1*sqrt(rs)+self.beta2*rs)
            else:
                return self.Aw*log(rs)+self.Bw-self.Aw/3.+2/3.*self.Cw*rs*log(rs)+(2*self.D-self.Cw)*rs/3.
            
    def EcVc(self, rs): # Ec-Vc
        if self.type<3 :
            x = rs/self.A
            epsilon = -0.5*self.C*((1+x*x*x)*log(1+1/x)+0.5*x-x*x-1/3.)
            return epsilon-Vc(rs)
        elif self.type<4: # type=3 WVN
            x = sqrt(rs)
            return self.Ap/6.*(self.cp*(x-self.xp0)-self.bp*x*self.xp0)/((x-self.xp0)*(x*x+self.bp*x+self.cp))
        else:
            if rs>1:
                return 2*self.gamma/(1+self.beta1*sqrt(rs)+self.beta2*rs)-Vc(rs)
            else:
                return self.Aw*log(rs)+self.Bw+self.Cw*rs*log(rs)+self.D*rs-Vc(rs)