``excVWN.m''

% VWN parameterization of the exchange correlation energy
function out=excVWN(n)
  %# Constants
  X1 = 0.75*(3.0/(2.0*pi))^(2.0/3.0);
  A  =  0.0310907;
  x0 = -0.10498;
  b  = 3.72744;
  c  = 12.9352;
  Q  = sqrt(4*c-b*b);
  X0 = x0*x0+b*x0+c;

  rs=(4*pi/3*n).^(-1/3); %# Added internal conversion to rs

  x=sqrt(rs); X=x.*x+b*x+c;

  out=-X1./rs ...
    + A*( ...
        +log(x.*x./X)+2*b/Q*atan(Q./(2*x+b)) ...
        -(b*x0)/X0*( ...
                    log((x-x0).*(x-x0)./X)+2*(2*x0+b)/Q*atan(Q./(2*x+b)) ...
                    ) ...
        );
endfunction