From 0d60666ab2c0295e35d37faab0efd98cada7bb2a Mon Sep 17 00:00:00 2001
From: Ken Esler <kpesler@gmail.com>
Date: Tue, 25 Mar 2008 21:08:02 +0000
Subject: [PATCH] Hopefully fixed parameter derivative evaluation for
 BsplineFunctor.

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@2545 e5b18d87-469d-4833-9cc0-8cdfa06e9491
---
 src/QMCWaveFunctions/Jastrow/BsplineFunctor.h | 50 +++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/src/QMCWaveFunctions/Jastrow/BsplineFunctor.h b/src/QMCWaveFunctions/Jastrow/BsplineFunctor.h
index 1f78d5472..f8a978909 100644
--- a/src/QMCWaveFunctions/Jastrow/BsplineFunctor.h
+++ b/src/QMCWaveFunctions/Jastrow/BsplineFunctor.h
@@ -30,6 +30,9 @@ namespace qmcplusplus {
     
     std::string elementType, pairType;
     std::vector<double> SplineCoefs;
+    // Stores the derivatives w.r.t. SplineCoefs 
+    // of the u, du/dr, and d2u/dr2
+    std::vector<TinyVector<real_type,3> > SplineDerivs;
     std::vector<double> Parameters;
     std::vector<std::string> ParameterNames;
     double Rcut, DeltaR, DeltaRInv;
@@ -52,6 +55,7 @@ namespace qmcplusplus {
 
       Parameters.resize(n);
       SplineCoefs.resize(numCoefs);
+      SplineDerivs.resize(numCoefs);
     }
     
     void reset() 
@@ -151,6 +155,52 @@ namespace qmcplusplus {
 
     }
 
+    inline bool
+    evaluate (real_type r, vector<TinyVector<real_type,3> >& derivs)
+    {
+      if (r >= Rcut) 
+	return false;
+
+      r *= DeltaRInv;
+      double ipart, t;
+      t = modf (r, &ipart);
+      int i = (int) ipart;
+      
+      double tp[4];
+      tp[0] = t*t*t;  tp[1] = t*t;  tp[2] = t;  tp[3] = 1.0;
+
+      // d/dp_i u(r)
+      SplineDerivs[i+0][0] = A[ 0]*tp[0] + A[ 1]*tp[1] + A[ 2]*tp[2] + A[ 3]*tp[3];
+      SplineDerivs[i+1][0] = A[ 4]*tp[0] + A[ 5]*tp[1] + A[ 6]*tp[2] + A[ 7]*tp[3];
+      SplineDerivs[i+2][0] = A[ 8]*tp[0] + A[ 9]*tp[1] + A[10]*tp[2] + A[11]*tp[3];
+      SplineDerivs[i+3][0] = A[12]*tp[0] + A[13]*tp[1] + A[14]*tp[2] + A[15]*tp[3];
+
+      // d/dp_i du/dr
+      SplineDerivs[i+0][1] = DeltaRInv * (dA[ 1]*tp[1] + dA[ 2]*tp[2] + dA[ 3]*tp[3]);
+      SplineDerivs[i+1][1] = DeltaRInv * (dA[ 5]*tp[1] + dA[ 6]*tp[2] + dA[ 7]*tp[3]);
+      SplineDerivs[i+2][1] = DeltaRInv * (dA[ 9]*tp[1] + dA[10]*tp[2] + dA[11]*tp[3]);
+      SplineDerivs[i+3][1] = DeltaRInv * (dA[13]*tp[1] + dA[14]*tp[2] + dA[15]*tp[3]);
+
+      // d/dp_i d2u/dr2
+      SplineDerivs[i+0][2] = DeltaRInv * DeltaRInv * (d2A[ 2]*tp[2] + d2A[ 3]*tp[3]);
+      SplineDerivs[i+1][2] = DeltaRInv * DeltaRInv * (d2A[ 6]*tp[2] + d2A[ 7]*tp[3])
+      SplineDerivs[i+2][2] = DeltaRInv * DeltaRInv * (d2A[10]*tp[2] + d2A[11]*tp[3]);
+      SplineDerivs[i+3][2] = DeltaRInv * DeltaRInv * (d2A[14]*tp[2] + d2A[15]*tp[3]);
+      
+      if (i == 0) 
+	for (int j=0; j<3; j++) {
+	  derivs[0][j] = SplineDerivs[1][j];
+	  derivs[1][j] = SplineDerivs[0][j] + SplineDerivs[2][j];
+	  derivs[2][j] = SplineDerivs[3][j];
+	  derivs[3][j] = SplineDerivs[4][j];
+	}
+      else 
+	for (int n=1; n<i+4; n++) 
+	  for (int j=0; j<3; j++) 
+	    derivs[n-1][j] = SplineDerivs[n][j];
+      return true;
+    }
+
     inline real_type f(real_type r) {
       if(r>Rcut) return 0.0;
       return evaluate (r);