Changeset b70503 for pcp/src/excor.c

Timestamp:

Apr 21, 2008, 2:19:24 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

f7754e

Parents:

6c8ee7

git-author:

Frederik Heber <heber@…> (04/18/08 15:29:28)

git-committer:

Frederik Heber <heber@…> (04/21/08 14:19:24)

Message:

sinh_inverse(), CalcSE_GC() and DensityGradient() added for Gradient Correction XC

File:

• pcp/src/excor.c (modified) (1 diff)

pcp/src/excor.c

@@ -787 +787 @@
   E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
   E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*(SumEx); // - SumEx_GC);
+}
+
+/** Returns \f$y = \sinh^{-1}(x)\f$.
+ * Recall that \f$ x = sinh (y) =  \frac{\exp(y)-\exp(-y)}{2}\f$. Then, solve the in \f$\exp(y)\f$ quadratic
+ * equation: \f$\exp(2y) - 1 - 2x\exp(y) = 0\f$ by pq-formula.
+ * \param arg argument \f$x\f$
+ * \return \f$\sinh^{-1}(x) = \ln(x+\sqrt{x^2+1})\f$
+ */
+inline static double sinh_inverse(double arg)
+{
+  return log(arg+sqrt(arg*arg+1));
+}
+
+/** Gradient correction according to Becke '92.
+ * \param EC ExCor structure
+ * \param p \f$\rho\f$ local density
+ * \param Dp \f$|\nabla \rho|\f$ local magnitude of density gradient
+ * \return \f$b \rho^{4/3} \frac{x^2}{(1+6bx \sinh^-1 x}\f$ with \f$x = \frac{|\nabla \rho|}{\rho^{4/3}}\f$
+ */
+double CalcSE_GC(struct ExCor *EC, double p, double Dp)
+{
+  double res = 0.;
+  double b = 0.0042;  // in atomic units, empirical constant fitted to noble gas atoms
+  double x = Dp/pow(p, 4./3.);
+  res = b * pow(p, 4./3.) * x/(1+6.*b*x*sinh_inverse(x));
+  
+  return res;
+}
+
+/** Evaluates magnitude of gradient by simple 3-star.
+ * \param *density density array
+ * \param i current node
+ * \param *Lev LatticeLevel structure for number of nodes per axis
+ * \param *Lat Lattice structure for axis lengths
+ * \return \f$|\nabla\rho|\f$
+ */
+double DensityGradient(fftw_real *density, int i, struct LatticeLevel *Lev, struct Lattice *Lat)
+{
+  double res=0., right_diff, left_diff;
+  int neighbour[NDIM], nodes[NDIM];
+  nodes[0] = Lev->Plan0.plan->local_nx;
+  nodes[1] = Lev->Plan0.plan->N[1];
+  nodes[2] = Lev->Plan0.plan->N[2];  
+  neighbour[0] = Lev->Plan0.plan->N[1] * Lev->Plan0.plan->N[2];
+  neighbour[1] = Lev->Plan0.plan->N[2];
+  neighbour[2] = 1.;
+  int k, l, nr, i_check = i;
+  double h;
+
+  //iS = n[2] + N[2]*(n[1] + N[1]*n0);  // howto access the array ...
+
+  for (k=NDIM-1;k>=0;k--) { // for each axis
+    h = 0.;
+    for (l=0;l<NDIM;l++)
+      h += Lat->RealBasis[k*NDIM+l]/Lev->Plan0.plan->N[k];  // finite distance
+    h = sqrt(h);
+    // check which limit exists: right, left, both?
+    right_diff = 0.;
+    left_diff = 0.;
+    nr = 0; // neighbour counter
+    if (i_check % nodes[k] != nodes[k]-1) {// right neighbour?
+      right_diff = density[i] - density[i+neighbour[k]];
+      nr++;
+    }
+    if (i_check % nodes[k] != 0) { // left neighbour?
+      left_diff = density[i] - density[i-neighbour[k]];
+      nr++;
+    }
+    res += (left_diff + right_diff)/(h*nr) * (left_diff + right_diff)/(h*nr);
+    i_check /= nodes[k];  // remove axis from i_check
+  }
+  
+  return sqrt(res);
 }