File: XC_PW92C.c

package info (click to toggle)
openmx 3.7.6-1
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166` ``````/********************************************************************** XC_PW92C.c: XC_PW92C.c is a subroutine to calculate the correlation energy density and potential developed by Perdew and Wang (Ref: J.P.Perdew and Y.Wang, PRB, 45, 13244 (1992)) for given up (dens[0]) and down (dens[1]) densities. This routine was written by T.Ozaki, based on the original fortran code provided by the SIESTA group through their website. Thanks to them. Log of PW92C.c: 22/Nov/2001 Released by T.Ozaki ***********************************************************************/ #include #include #include #include "openmx_common.h" #define den_min 1.0e-14 #define den_min_half 0.5*1.0e-14 void XC_PW92C(double dens[2], double Ec[1], double Vc[2]) { int i; double dtot,rs,srs,zeta,coe; double dum,dum1,dum2,b,c,dbdrs,dcdrs; double tmp0,tmp1,tmp2,tmp12,tmp14,tmp22,tmp24; double fpp0,f,dfdz; double G[3],dGdrs[3]; double dEcdd[2]; double dEcdrs,dEcdz; double drsdd,dzdd[2]; /**************************************************** parameters from Table I of Perdew and Wang, PRB, 45, 13244 (92) ****************************************************/ double p[3] = {1.0000000, 1.0000000, 1.0000000}; double A[3] = {0.0310910, 0.0155450, 0.0168870}; double alpha1[3] = {0.2137000, 0.2054800, 0.1112500}; double beta1[3] = {7.5957000, 14.1189000, 10.3570000}; double beta2[3] = {3.5876000, 6.1977000, 3.6231000}; double beta3[3] = {1.6382000, 3.3662000, 0.8802600}; double beta4[3] = {0.4929400, 0.6251700, 0.4967100}; /**************************************************** zeta and rs ****************************************************/ coe = 0.6203504908994; /* pow(3.0/4.0/PI,1.0/3.0); */ dtot = dens[0] + dens[1]; if (dtot<=den_min){ rs = 6203.504908994; /* coe*pow(1.0e-12,-1.0/3.0); */ dens[0] = den_min_half; dens[1] = den_min_half; dtot = den_min; } else rs = coe*pow(dtot,-0.33333333333333333333333); tmp0 = 1.0/dtot; zeta = tmp0*(dens[0] - dens[1]); if (0.99