mirror of https://github.com/QMCPACK/qmcpack.git
Remove three-body bspline jastrow
This commit is contained in:
Jaron Krogel
parent
df44e6b2b4
commit
3f8983e8db
|
|
@ -1,618 +0,0 @@
|
|||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
// This file is distributed under the University of Illinois/NCSA Open Source License.
|
||||
// See LICENSE file in top directory for details.
|
||||
//
|
||||
// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
|
||||
//
|
||||
// File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_BSPLINE3D_FUNCTOR_H
|
||||
#define QMCPLUSPLUS_BSPLINE3D_FUNCTOR_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
#include "Numerics/HDFNumericAttrib.h"
|
||||
#include "Utilities/ProgressReportEngine.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
#include "Numerics/LinearFit.h"
|
||||
#include <cstdio>
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
struct BsplineFunctor3D: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
typedef real_type value_type;
|
||||
int NumParams_eI, NumParams_ee;
|
||||
int Dummy;
|
||||
const TinyVector<real_type,16> A, dA, d2A, d3A;
|
||||
//static const real_type A[16], dA[16], d2A[16];
|
||||
real_type Deltax, DeltaxInv;
|
||||
real_type DeltaR_eI, DeltaRInv_eI;
|
||||
real_type Y, dY, d2Y;
|
||||
Array<real_type,3> 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<real_type> Parameters;
|
||||
Array<real_type,3> ParamArray;
|
||||
std::vector<std::string> ParameterNames;
|
||||
std::string iSpecies, eSpecies1, eSpecies2;
|
||||
int ResetCount;
|
||||
|
||||
///constructor
|
||||
BsplineFunctor3D(real_type ecusp=0.0, real_type icusp=0.0) :
|
||||
NumParams_eI(0), NumParams_ee(0),
|
||||
A(-1.0/6.0, 3.0/6.0, -3.0/6.0, 1.0/6.0,
|
||||
3.0/6.0, -6.0/6.0, 0.0/6.0, 4.0/6.0,
|
||||
-3.0/6.0, 3.0/6.0, 3.0/6.0, 1.0/6.0,
|
||||
1.0/6.0, 0.0/6.0, 0.0/6.0, 0.0/6.0),
|
||||
dA(0.0, -0.5, 1.0, -0.5,
|
||||
0.0, 1.5, -2.0, 0.0,
|
||||
0.0, -1.5, 1.0, 0.5,
|
||||
0.0, 0.5, 0.0, 0.0),
|
||||
d2A(0.0, 0.0, -1.0, 1.0,
|
||||
0.0, 0.0, 3.0, -2.0,
|
||||
0.0, 0.0, -3.0, 1.0,
|
||||
0.0, 0.0, 1.0, 0.0),
|
||||
d3A(0.0, 0.0, 0.0, -1.0,
|
||||
0.0, 0.0, 0.0, 3.0,
|
||||
0.0, 0.0, 0.0, -3.0,
|
||||
0.0, 0.0, 0.0, 1.0),
|
||||
ResetCount(0)
|
||||
{
|
||||
cutoff_radius = 0.0;
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new BsplineFunctor3D(*this);
|
||||
}
|
||||
|
||||
void resize(int neI, int nee)
|
||||
{
|
||||
NumParams_eI = neI;
|
||||
NumParams_ee = nee;
|
||||
int numCoefs_eI = NumParams_eI + 4;
|
||||
int numCoefs_ee = NumParams_ee + 4;
|
||||
SplineCoefs.resize(numCoefs_ee, numCoefs_eI, numCoefs_eI);
|
||||
ParamArray.resize(NumParams_ee, NumParams_eI, NumParams_eI);
|
||||
int numParams = NumParams_eI*(NumParams_eI+1)/2 * NumParams_ee;
|
||||
Parameters.resize(numParams);
|
||||
int numKnots_eI = numCoefs_eI - 2;
|
||||
DeltaR_eI = 0.5*cutoff_radius / (real_type)(numKnots_eI - 1);
|
||||
DeltaRInv_eI = 1.0/DeltaR_eI;
|
||||
int numKnots_ee = numCoefs_ee - 2;
|
||||
Deltax = 1.0 / (real_type)(numKnots_ee - 1);
|
||||
DeltaxInv = 1.0/Deltax;
|
||||
}
|
||||
|
||||
inline int getNumParameters()
|
||||
{
|
||||
return Parameters.size();
|
||||
}
|
||||
|
||||
void reset()
|
||||
{
|
||||
int numCoefs_eI = NumParams_eI + 4;
|
||||
int numKnots_eI = numCoefs_eI - 2;
|
||||
int numCoefs_ee = NumParams_ee + 4;
|
||||
int numKnots_ee = numCoefs_ee - 2;
|
||||
DeltaR_eI = 0.5*cutoff_radius / (real_type)(numKnots_eI - 1);
|
||||
Deltax = 1.0 / (real_type)(numKnots_ee - 1);
|
||||
DeltaRInv_eI = 1.0/DeltaR_eI;
|
||||
DeltaxInv = 1.0/Deltax;
|
||||
// Zero out all coefficients
|
||||
for (int i=0; i<SplineCoefs.size(0); i++)
|
||||
for (int j=0; j<SplineCoefs.size(1); j++)
|
||||
for (int k=0; k<SplineCoefs.size(2); k++)
|
||||
SplineCoefs(i,j,k) = 0.0;
|
||||
// Set unconstrained coefficients
|
||||
for (int i=2; i<NumParams_ee; i++)
|
||||
for (int j=2; j<NumParams_eI; j++)
|
||||
for (int k=2; k<NumParams_eI; k++)
|
||||
SplineCoefs(i+1,j+1,k+1) = ParamArray(i,j,k);
|
||||
// j-k plane
|
||||
// Set e-I cusp parameters
|
||||
for (int j=2; j<NumParams_eI; j++)
|
||||
for (int k=2; k<NumParams_eI; k++)
|
||||
{
|
||||
SplineCoefs(1,j+1,k+1) = ParamArray(0,j,k);
|
||||
SplineCoefs(2,j+1,k+1) = ParamArray(1,j,k);
|
||||
SplineCoefs(0,k+1,j+1) = ParamArray(1,j,k);
|
||||
}
|
||||
// i-j plane
|
||||
// Set e-e cusp parameters
|
||||
for (int i=2; i<NumParams_ee; i++)
|
||||
for (int j=2; j<NumParams_eI; j++)
|
||||
{
|
||||
SplineCoefs(i+1,j+1,1) = ParamArray(i,j,0);
|
||||
SplineCoefs(i+1,j+1,2) = ParamArray(i,j,1);
|
||||
SplineCoefs(i+1,j+1,0) = ParamArray(i,j,1);
|
||||
}
|
||||
// i-k plane
|
||||
// Set e-e cusp parameters
|
||||
for (int i=2; i<NumParams_ee; i++)
|
||||
for (int k=2; k<NumParams_eI; k++)
|
||||
{
|
||||
SplineCoefs(i+1,1,k+1) = ParamArray(i,0,k);
|
||||
SplineCoefs(i+1,2,k+1) = ParamArray(i,1,k);
|
||||
SplineCoefs(i+1,0,k+1) = ParamArray(i,1,k);
|
||||
}
|
||||
// i edge
|
||||
for (int i=2; i<NumParams_ee; i++)
|
||||
{
|
||||
SplineCoefs(i+1,1,1) = ParamArray(i,0,0);
|
||||
SplineCoefs(i+1,2,1) = ParamArray(i,1,0);
|
||||
SplineCoefs(i+1,0,1) = ParamArray(i,1,0);
|
||||
SplineCoefs(i+1,1,2) = ParamArray(i,0,1);
|
||||
SplineCoefs(i+1,2,2) = ParamArray(i,1,1);
|
||||
SplineCoefs(i+1,0,2) = ParamArray(i,1,1);
|
||||
SplineCoefs(i+1,1,0) = ParamArray(i,0,1);
|
||||
SplineCoefs(i+1,2,0) = ParamArray(i,1,1);
|
||||
SplineCoefs(i+1,0,0) = ParamArray(i,1,1);
|
||||
}
|
||||
// j edge
|
||||
for (int j=2; j<NumParams_eI; j++)
|
||||
{
|
||||
SplineCoefs(1,j+1,1) = ParamArray(0,j,0);
|
||||
SplineCoefs(2,j+1,1) = ParamArray(1,j,0);
|
||||
SplineCoefs(0,j+1,1) = ParamArray(1,j,0);
|
||||
SplineCoefs(1,j+1,2) = ParamArray(0,j,1);
|
||||
SplineCoefs(2,j+1,2) = ParamArray(1,j,1);
|
||||
SplineCoefs(0,j+1,2) = ParamArray(1,j,1);
|
||||
SplineCoefs(1,j+1,0) = ParamArray(0,j,1);
|
||||
SplineCoefs(2,j+1,0) = ParamArray(1,j,1);
|
||||
SplineCoefs(0,j+1,0) = ParamArray(1,j,1);
|
||||
}
|
||||
// k edge
|
||||
for (int k=2; k<NumParams_eI; k++)
|
||||
{
|
||||
SplineCoefs(1,1,k+1) = ParamArray(0,0,k);
|
||||
SplineCoefs(2,1,k+1) = ParamArray(1,0,k);
|
||||
SplineCoefs(0,1,k+1) = ParamArray(1,0,k);
|
||||
SplineCoefs(1,2,k+1) = ParamArray(0,1,k);
|
||||
SplineCoefs(2,2,k+1) = ParamArray(1,1,k);
|
||||
SplineCoefs(0,2,k+1) = ParamArray(1,1,k);
|
||||
SplineCoefs(1,0,k+1) = ParamArray(0,1,k);
|
||||
SplineCoefs(2,0,k+1) = ParamArray(1,1,k);
|
||||
SplineCoefs(0,0,k+1) = ParamArray(1,1,k);
|
||||
}
|
||||
// Copy the 8 uniquely determined values
|
||||
SplineCoefs(1,1,1) = ParamArray(0,0,0);
|
||||
SplineCoefs(1,1,2) = ParamArray(0,0,1);
|
||||
SplineCoefs(1,2,1) = ParamArray(0,1,0);
|
||||
SplineCoefs(1,2,2) = ParamArray(0,1,1);
|
||||
SplineCoefs(2,1,1) = ParamArray(1,0,0);
|
||||
SplineCoefs(2,1,2) = ParamArray(1,0,1);
|
||||
SplineCoefs(2,2,1) = ParamArray(1,1,0);
|
||||
SplineCoefs(2,2,2) = ParamArray(1,1,1);
|
||||
// Now satisfy cusp constraints
|
||||
// ee
|
||||
SplineCoefs(1,1,0) = ParamArray(0,0,1);
|
||||
SplineCoefs(1,2,0) = ParamArray(0,1,1);
|
||||
SplineCoefs(2,1,0) = ParamArray(1,0,1);
|
||||
SplineCoefs(2,2,0) = ParamArray(1,1,1);
|
||||
SplineCoefs(1,0,1) = ParamArray(0,1,0);
|
||||
SplineCoefs(1,0,2) = ParamArray(0,1,1);
|
||||
SplineCoefs(2,0,1) = ParamArray(1,1,0);
|
||||
SplineCoefs(2,0,2) = ParamArray(1,1,1);
|
||||
// eI
|
||||
SplineCoefs(0,1,1) = ParamArray(1,0,0);
|
||||
SplineCoefs(0,1,2) = ParamArray(1,0,1);
|
||||
SplineCoefs(0,2,1) = ParamArray(1,1,0);
|
||||
SplineCoefs(0,2,2) = ParamArray(1,1,1);
|
||||
// More than one cusp constraint
|
||||
SplineCoefs(0,0,1) = ParamArray(1,1,0);
|
||||
SplineCoefs(0,1,0) = ParamArray(1,0,1);
|
||||
SplineCoefs(1,0,0) = ParamArray(0,1,1);
|
||||
SplineCoefs(0,0,2) = ParamArray(1,1,1);
|
||||
SplineCoefs(0,2,0) = ParamArray(1,1,1);
|
||||
SplineCoefs(2,0,0) = ParamArray(1,1,1);
|
||||
SplineCoefs(0,0,0) = ParamArray(1,1,1);
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r_12,
|
||||
real_type r_1I,
|
||||
real_type r_2I) const
|
||||
{
|
||||
if (r_12 >= cutoff_radius || r_1I >= 0.5*cutoff_radius ||
|
||||
r_2I >= 0.5*cutoff_radius)
|
||||
return 0.0;
|
||||
real_type x = r_12 / (r_1I + r_2I);
|
||||
x *= DeltaxInv;
|
||||
r_1I *= DeltaRInv_eI;
|
||||
r_2I *= DeltaRInv_eI;
|
||||
real_type ipart, t, u, v;
|
||||
int i, j, k;
|
||||
t = std::modf(x, &ipart);
|
||||
i = (int) ipart;
|
||||
u = std::modf(r_1I, &ipart);
|
||||
j = (int) ipart;
|
||||
v = std::modf(r_2I, &ipart);
|
||||
k = (int) ipart;
|
||||
real_type tp[4], up[4], vp[4], a[4], b[4], c[4];
|
||||
tp[0] = t*t*t;
|
||||
tp[1] = t*t;
|
||||
tp[2] = t;
|
||||
tp[3] = 1.0;
|
||||
up[0] = u*u*u;
|
||||
up[1] = u*u;
|
||||
up[2] = u;
|
||||
up[3] = 1.0;
|
||||
vp[0] = v*v*v;
|
||||
vp[1] = v*v;
|
||||
vp[2] = v;
|
||||
vp[3] = 1.0;
|
||||
int index=0;
|
||||
for (int m=0; m<4; m++)
|
||||
{
|
||||
a[m] = b[m] = c[m] = 0.0;
|
||||
for (int n=0; n<4; n++)
|
||||
{
|
||||
a[m] += A[index] * tp[n];
|
||||
b[m] += A[index] * up[n];
|
||||
c[m] += A[index] * vp[n];
|
||||
index++;
|
||||
}
|
||||
}
|
||||
real_type val = 0.0;
|
||||
for (int ia=0; ia<4; ia++)
|
||||
for (int ib=0; ib<4; ib++)
|
||||
for (int ic=0; ic<4; ic++)
|
||||
val += SplineCoefs(i+ia, j+ib, k+ic)*a[ia]*b[ib]*c[ic];
|
||||
return val;
|
||||
}
|
||||
|
||||
inline real_type evaluateV(int Nptcl,
|
||||
const real_type* restrict r_12_array,
|
||||
const real_type* restrict r_1I_array,
|
||||
const real_type* restrict r_2I_array) const
|
||||
{
|
||||
real_type val_tot(0);
|
||||
for(int ptcl=0; ptcl<Nptcl; ptcl++)
|
||||
val_tot+=evaluate(r_12_array[ptcl],r_1I_array[ptcl],r_2I_array[ptcl]);
|
||||
return val_tot;
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r_12, real_type r_1I, real_type r_2I,
|
||||
TinyVector<real_type,3> &grad,
|
||||
Tensor<real_type,3> &hess) const
|
||||
{
|
||||
if (r_12 >= cutoff_radius || r_1I >= 0.5*cutoff_radius ||
|
||||
r_2I >= 0.5*cutoff_radius)
|
||||
{
|
||||
grad = 0.0;
|
||||
hess = 0.0;
|
||||
return 0.0;
|
||||
}
|
||||
// double eps = 1.0e-6;
|
||||
// grad[0] = (evaluate (r_12+eps, r_1I, r_2I) -evaluate (r_12-eps, r_1I, r_2I))/(2.0*eps);
|
||||
// grad[1] = (evaluate (r_12, r_1I+eps, r_2I) -evaluate (r_12, r_1I-eps, r_2I))/(2.0*eps);
|
||||
// grad[2] = (evaluate (r_12, r_1I, r_2I+eps) -evaluate (r_12, r_1I, r_2I-eps))/(2.0*eps);
|
||||
real_type qInv = 1.0/(r_1I + r_2I);
|
||||
real_type x = r_12 * qInv;
|
||||
x *= DeltaxInv;
|
||||
r_1I *= DeltaRInv_eI;
|
||||
r_2I *= DeltaRInv_eI;
|
||||
real_type ipart, t, u, v;
|
||||
int i, j, k;
|
||||
t = std::modf(x, &ipart);
|
||||
i = (int) ipart;
|
||||
u = std::modf(r_1I, &ipart);
|
||||
j = (int) ipart;
|
||||
v = std::modf(r_2I, &ipart);
|
||||
k = (int) ipart;
|
||||
real_type tp[4], up[4], vp[4], a[4], b[4], c[4],
|
||||
da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4];
|
||||
tp[0] = t*t*t;
|
||||
tp[1] = t*t;
|
||||
tp[2] = t;
|
||||
tp[3] = 1.0;
|
||||
up[0] = u*u*u;
|
||||
up[1] = u*u;
|
||||
up[2] = u;
|
||||
up[3] = 1.0;
|
||||
vp[0] = v*v*v;
|
||||
vp[1] = v*v;
|
||||
vp[2] = v;
|
||||
vp[3] = 1.0;
|
||||
int index=0;
|
||||
for (int m=0; m<4; m++)
|
||||
{
|
||||
a[m]=b[m]=c[m]=da[m]=db[m]=dc[m]=d2a[m]=d2b[m]=d2c[m]=0.0;
|
||||
for (int n=0; n<4; n++)
|
||||
{
|
||||
a[m]+=A[index]*tp[n];
|
||||
da[m]+=dA[index]*tp[n];
|
||||
d2a[m]+=d2A[index]*tp[n];
|
||||
b[m]+=A[index]*up[n];
|
||||
db[m]+=dA[index]*up[n];
|
||||
d2b[m]+=d2A[index]*up[n];
|
||||
c[m]+=A[index]*vp[n];
|
||||
dc[m]+=dA[index]*vp[n];
|
||||
d2c[m]+=d2A[index]*vp[n];
|
||||
index++;
|
||||
}
|
||||
}
|
||||
real_type val = 0.0;
|
||||
TinyVector<real_type,3> gradF;
|
||||
Tensor<real_type,3> hessF;
|
||||
gradF = 0.0;
|
||||
hessF = 0.0;
|
||||
for (int ia=0; ia<4; ia++)
|
||||
for (int ib=0; ib<4; ib++)
|
||||
for (int ic=0; ic<4; ic++)
|
||||
{
|
||||
real_type coef = SplineCoefs(i+ia, j+ib, k+ic);
|
||||
val += coef * a[ia] * b[ib] * c[ic];
|
||||
gradF[0] += coef * da[ia] * b[ib] * c[ic];
|
||||
gradF[1] += coef * a[ia] * db[ib] * c[ic];
|
||||
gradF[2] += coef * a[ia] * b[ib] * dc[ic];
|
||||
hessF(0,0) += coef *d2a[ia] * b[ib] * c[ic];
|
||||
hessF(0,1) += coef * da[ia] * db[ib] * c[ic];
|
||||
hessF(0,2) += coef * da[ia] * b[ib] * dc[ic];
|
||||
hessF(1,1) += coef * a[ia] *d2b[ib] * c[ic];
|
||||
hessF(1,2) += coef * a[ia] * db[ib] * dc[ic];
|
||||
hessF(2,2) += coef * a[ia] * b[ib] *d2c[ic];
|
||||
}
|
||||
gradF[0] *= DeltaxInv;
|
||||
gradF[1] *= DeltaRInv_eI;
|
||||
gradF[2] *= DeltaRInv_eI;
|
||||
grad[0] = qInv*gradF[0];
|
||||
grad[1] = gradF[1] - r_12*qInv*qInv*gradF[0];
|
||||
grad[2] = gradF[2] - r_12*qInv*qInv*gradF[0];
|
||||
hessF(0,0) *= DeltaxInv * DeltaxInv;
|
||||
hessF(0,1) *= DeltaxInv * DeltaRInv_eI;
|
||||
hessF(0,2) *= DeltaxInv * DeltaRInv_eI;
|
||||
hessF(1,1) *= DeltaRInv_eI * DeltaRInv_eI;
|
||||
hessF(1,2) *= DeltaRInv_eI * DeltaRInv_eI;
|
||||
hessF(2,2) *= DeltaRInv_eI * DeltaRInv_eI;
|
||||
hess(0,0) = qInv*qInv*hessF(0,0);
|
||||
hess(0,1) = qInv*hessF(0,1) - gradF[0]*qInv*qInv - hessF(0,0)*r_12*qInv*qInv*qInv;
|
||||
hess(0,2) = qInv*hessF(0,2) - gradF[0]*qInv*qInv - hessF(0,0)*r_12*qInv*qInv*qInv;
|
||||
hess(1,1) = hessF(1,1) + 2.0*gradF[0]*r_12 * qInv*qInv*qInv - 2.0*hessF(0,1) *r_12*qInv*qInv + hessF(0,0)*r_12*r_12*qInv*qInv*qInv*qInv;
|
||||
hess(1,2) = hessF(1,2) + 2.0*gradF[0]*r_12 * qInv*qInv*qInv - (hessF(0,2)+hessF(0,1))*r_12*qInv*qInv + hessF(0,0)*r_12*r_12*qInv*qInv*qInv*qInv;
|
||||
hess(2,2) = hessF(2,2) + 2.0*gradF[0]*r_12 * qInv*qInv*qInv - 2.0*hessF(0,2) *r_12*qInv*qInv + hessF(0,0)*r_12*r_12*qInv*qInv*qInv*qInv;
|
||||
hess(1,0) = hess(0,1);
|
||||
hess(2,0) = hess(0,2);
|
||||
hess(2,1) = hess(1,2);
|
||||
return val;
|
||||
}
|
||||
|
||||
// assume r_1I < L && r_2I < L, compression and screening is handled outside
|
||||
inline void evaluateVGL(int Nptcl, const real_type* restrict r_12_array,
|
||||
const real_type* restrict r_1I_array,
|
||||
const real_type* restrict r_2I_array,
|
||||
real_type* restrict val_array,
|
||||
real_type* restrict grad0_array,
|
||||
real_type* restrict grad1_array,
|
||||
real_type* restrict grad2_array,
|
||||
real_type* restrict hess00_array,
|
||||
real_type* restrict hess11_array,
|
||||
real_type* restrict hess22_array,
|
||||
real_type* restrict hess01_array,
|
||||
real_type* restrict hess02_array) const
|
||||
{
|
||||
APP_ABORT("BsplineFunctor3D::evaluateVGL not implemented yet!");
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r_12, real_type r_1I, real_type r_2I,
|
||||
TinyVector<real_type,3> &grad,
|
||||
Tensor<real_type,3> &hess,
|
||||
TinyVector<Tensor<real_type,3>,3> &d3)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
|
||||
inline real_type evaluate(real_type r, real_type rinv)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
|
||||
inline bool
|
||||
evaluateDerivatives(real_type r, std::vector<TinyVector<real_type,3> >& derivs)
|
||||
{
|
||||
//what is this?
|
||||
return false;
|
||||
}
|
||||
|
||||
inline bool
|
||||
evaluateDerivatives (real_type r_12, real_type r_1I, real_type r_2I,
|
||||
std::vector<real_type> &d_vals,
|
||||
std::vector<TinyVector<real_type,3> >& d_grads,
|
||||
std::vector<Tensor<real_type,3> > &d_hess)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
inline real_type f(real_type r)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
inline real_type df(real_type r)
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
ReportEngine PRE("BsplineFunctor3D","put(xmlNodePtr)");
|
||||
//CuspValue = -1.0e10;
|
||||
NumParams_eI = NumParams_ee = 0;
|
||||
cutoff_radius = 0.0;
|
||||
OhmmsAttributeSet rAttrib;
|
||||
rAttrib.add(NumParams_ee, "esize");
|
||||
rAttrib.add(NumParams_eI, "isize");
|
||||
rAttrib.add(cutoff_radius, "rcut");
|
||||
rAttrib.put(cur);
|
||||
if (NumParams_eI == 0)
|
||||
PRE.error("You must specify a positive number for \"isize\"",true);
|
||||
if (NumParams_ee == 0)
|
||||
PRE.error("You must specify a positive number for \"esize\"",true);
|
||||
app_log() << " esize = " << NumParams_ee << " parameters " << std::endl;
|
||||
app_log() << " isize = " << NumParams_eI << " parameters " << std::endl;
|
||||
app_log() << " rcut = " << cutoff_radius << std::endl;
|
||||
resize(NumParams_eI, NumParams_ee);
|
||||
// Now read coefficents
|
||||
xmlNodePtr xmlCoefs = cur->xmlChildrenNode;
|
||||
while (xmlCoefs != NULL)
|
||||
{
|
||||
std::string cname((const char*)xmlCoefs->name);
|
||||
if (cname == "coefficients")
|
||||
{
|
||||
std::string type("0"), id("0");
|
||||
OhmmsAttributeSet cAttrib;
|
||||
cAttrib.add(id, "id");
|
||||
cAttrib.add(type, "type");
|
||||
cAttrib.put(xmlCoefs);
|
||||
if (type != "Array")
|
||||
{
|
||||
PRE.error("Unknown correlation type " + type +
|
||||
" in BsplineFunctor3D." + "Resetting to \"Array\"");
|
||||
xmlNewProp(xmlCoefs, (const xmlChar*) "type", (const xmlChar*) "Array");
|
||||
}
|
||||
std::vector<real_type> params;
|
||||
putContent(params, xmlCoefs);
|
||||
if (params.size() == Parameters.size())
|
||||
Parameters = params;
|
||||
else
|
||||
if (params.size() == 0)
|
||||
{
|
||||
app_log()<<" Initializing all parameters to zero"<< std::endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
app_error() << "Expected " << Parameters.size() << " parameters,"
|
||||
<< " but found only " << params.size()
|
||||
<< " in BsplineFunctor3D.\n";
|
||||
abort();
|
||||
}
|
||||
// Setup parameter names
|
||||
int index=0;
|
||||
for (int i=0; i< NumParams_ee; i++)
|
||||
for (int j=0; j < NumParams_eI; j++)
|
||||
for (int k=0; k<=j; k++)
|
||||
{
|
||||
std::stringstream sstr;
|
||||
sstr << id << "_" << i << "_" << j << "_" << k;
|
||||
myVars.insert(sstr.str(),Parameters[index],true,optimize::LOGLINEAR_P);
|
||||
ParamArray(i,j,k) = ParamArray(i,k,j) = Parameters[index];
|
||||
index++;
|
||||
}
|
||||
app_log() << "Parameter Name Value\n";
|
||||
myVars.print(app_log());
|
||||
}
|
||||
xmlCoefs = xmlCoefs->next;
|
||||
}
|
||||
reset();
|
||||
print();
|
||||
return true;
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int iparam = 0;
|
||||
for (int i=0; i<NumParams_ee; i++)
|
||||
for (int j=0; j<NumParams_eI; j++)
|
||||
for (int k=0; k<=j; k++)
|
||||
{
|
||||
int loc = myVars.where(iparam);
|
||||
if (loc >=0)
|
||||
Parameters[iparam] = myVars[iparam] = active[loc];
|
||||
ParamArray(i,j,k) = Parameters[iparam];
|
||||
ParamArray(i,k,j) = Parameters[iparam];
|
||||
iparam++;
|
||||
}
|
||||
reset();
|
||||
// for(int i=0; i<Parameters.size(); ++i) {
|
||||
// int loc=myVars.where(i);
|
||||
// if(loc>=0) Parameters[i]=myVars[i]=active[loc];
|
||||
// }
|
||||
if (ResetCount++ == 100)
|
||||
{
|
||||
ResetCount = 0;
|
||||
//print();
|
||||
}
|
||||
reset();
|
||||
}
|
||||
|
||||
|
||||
void print()
|
||||
{
|
||||
const int N = 100;
|
||||
std::string fname = iSpecies + ".J3.h5";
|
||||
hid_t hid = H5Fcreate(fname.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
|
||||
Array<real_type,3> val(N,N,N);
|
||||
for (int i=0; i<N; i++)
|
||||
{
|
||||
double r_12 = (real_type)i/(real_type)(N-1);
|
||||
for (int j=0; j<N; j++)
|
||||
{
|
||||
double r_1I = (real_type)j/(real_type)(N-1) * 0.5*cutoff_radius;
|
||||
for (int k=0; k<N; k++)
|
||||
{
|
||||
double r_2I = (real_type)k/(real_type)(N-1) * 0.5*cutoff_radius;
|
||||
val(i,j,k) = evaluate(r_12*(r_1I+r_2I), r_1I, r_2I);
|
||||
}
|
||||
}
|
||||
}
|
||||
Array<double,3> SplineCoefs_DP(NumParams_ee + 4, NumParams_eI + 4, NumParams_eI + 4);
|
||||
Array<double,3> ParamArray_DP(NumParams_ee, NumParams_eI, NumParams_eI);
|
||||
Array<double,3> val_DP(N,N,N);
|
||||
HDFAttribIO<Array<double,3> > coefs_attrib(SplineCoefs_DP);
|
||||
HDFAttribIO<Array<double,3> > param_attrib(ParamArray_DP);
|
||||
HDFAttribIO<Array<double,3> > val_attrib(val_DP);
|
||||
SplineCoefs_DP = SplineCoefs;
|
||||
ParamArray_DP = ParamArray;
|
||||
val_DP = val;
|
||||
val_attrib.write(hid, "val");
|
||||
coefs_attrib.write(hid, "coefs");
|
||||
param_attrib.write(hid, "params");
|
||||
H5Fclose(hid);
|
||||
// std::string fname = (elementType != "") ? elementType : pairType;
|
||||
// fname = fname + ".dat";
|
||||
// //cerr << "Writing " << fname << " file.\n";
|
||||
// FILE *fout = fopen (fname.c_str(), "w");
|
||||
// for (double r=0.0; r<cutoff_radius; r+=0.001)
|
||||
// fprintf (fout, "%8.3f %16.10f\n", r, evaluate(r));
|
||||
// fclose(fout);
|
||||
}
|
||||
|
||||
|
||||
void print(std::ostream& os)
|
||||
{
|
||||
/* no longer correct. Ye Luo
|
||||
int n=100;
|
||||
real_type d=cutoff_radius/100.,r=0;
|
||||
real_type u,du,d2du;
|
||||
for (int i=0; i<n; ++i)
|
||||
{
|
||||
u=evaluate(r,du,d2du);
|
||||
os << std::setw(22) << r << std::setw(22) << u << std::setw(22) << du
|
||||
<< std::setw(22) << d2du << std::endl;
|
||||
r+=d;
|
||||
}
|
||||
*/
|
||||
}
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -17,14 +17,12 @@
|
|||
#include "Particle/DistanceTableData.h"
|
||||
#include "Particle/DistanceTable.h"
|
||||
#include "QMCWaveFunctions/Jastrow/eeI_JastrowBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/BsplineFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/eeI_JastrowOrbital.h"
|
||||
#include "QMCWaveFunctions/Jastrow/JeeIOrbitalSoA.h"
|
||||
#include "QMCWaveFunctions/Jastrow/DiffOneBodyJastrowOrbital.h"
|
||||
#include "QMCWaveFunctions/Jastrow/TwoBodyJastrowOrbital.h"
|
||||
#include "QMCWaveFunctions/Jastrow/DiffTwoBodyJastrowOrbital.h"
|
||||
#include "Utilities/ProgressReportEngine.h"
|
||||
#include "QMCWaveFunctions/Jastrow/BsplineFunctor3D.h"
|
||||
#include "QMCWaveFunctions/Jastrow/PolynomialFunctor3D.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
|
|
@ -109,29 +107,17 @@ bool eeI_JastrowBuilder::put(xmlNodePtr cur)
|
|||
ReportEngine PRE(ClassName,"put(xmlNodePtr)");
|
||||
bool PrintTables=true;
|
||||
xmlNodePtr kids = cur->xmlChildrenNode;
|
||||
typedef BsplineFunctor3D FuncType;
|
||||
// Create a three-body Jastrow
|
||||
if (sourcePtcl)
|
||||
{
|
||||
// std::cerr << "sourcePtcl = " << sourcePtcl << std::endl;
|
||||
std::string ftype("Bspline");
|
||||
std::string ftype("polynomial");
|
||||
OhmmsAttributeSet tAttrib;
|
||||
tAttrib.add(ftype,"function");
|
||||
tAttrib.put (cur);
|
||||
SpeciesSet &iSet = sourcePtcl->getSpeciesSet();
|
||||
SpeciesSet &eSet = targetPtcl.getSpeciesSet();
|
||||
int numiSpecies = iSet.getTotalNum();
|
||||
if (ftype == "Bspline")
|
||||
{
|
||||
#ifdef ENABLE_SOA
|
||||
typedef JeeIOrbitalSoA<BsplineFunctor3D> J3Type;
|
||||
#else
|
||||
typedef eeI_JastrowOrbital<BsplineFunctor3D> J3Type;
|
||||
#endif
|
||||
J3Type &J3 = *(new J3Type(*sourcePtcl, targetPtcl, true));
|
||||
putkids (kids, J3);
|
||||
}
|
||||
else if (ftype == "polynomial")
|
||||
if (ftype == "polynomial")
|
||||
{
|
||||
#ifdef ENABLE_SOA
|
||||
typedef JeeIOrbitalSoA<PolynomialFunctor3D> J3Type;
|
||||
|
|
@ -147,179 +133,9 @@ bool eeI_JastrowBuilder::put(xmlNodePtr cur)
|
|||
<< " eeI_JastrowBuilder. Aborting.\n";
|
||||
abort();
|
||||
}
|
||||
// // Find the number of the source species
|
||||
// bool success=false;
|
||||
// while (kids != NULL) {
|
||||
// std::string kidsname = (char*)kids->name;
|
||||
// if (kidsname == "correlation") {
|
||||
// RealType ee_cusp=0.0;
|
||||
// RealType eI_cusp=0.0;
|
||||
// std::string iSpecies, eSpecies1("u"), eSpecies2("u");
|
||||
// OhmmsAttributeSet rAttrib;
|
||||
// rAttrib.add(iSpecies,"ispecies");
|
||||
// rAttrib.add(eSpecies1,"especies1");
|
||||
// rAttrib.add(eSpecies2,"especies2");
|
||||
// rAttrib.add(ee_cusp,"ecusp");
|
||||
// rAttrib.add(eI_cusp,"icusp");
|
||||
// rAttrib.put(kids);
|
||||
// BsplineFunctor3D *functor = new BsplineFunctor3D(ee_cusp, eI_cusp);
|
||||
// functor->iSpecies = iSpecies;
|
||||
// functor->eSpecies1 = eSpecies1;
|
||||
// functor->eSpecies2 = eSpecies2;
|
||||
// int iNum = iSet.findSpecies (iSpecies);
|
||||
// int eNum1 = eSet.findSpecies (eSpecies1);
|
||||
// int eNum2 = eSet.findSpecies (eSpecies2);
|
||||
// functor->put (kids);
|
||||
// std::strstream aname;
|
||||
// aname << iSpecies << "_" << eSpecies1 << "_" << eSpecies2;
|
||||
// J3->addFunc(aname.str(), iNum, eNum1, eNum2, functor);
|
||||
// }
|
||||
// kids = kids->next;
|
||||
// }
|
||||
// targetPsi.addOrbital(J3,"eeI_bspline");
|
||||
// J3->setOptimizable(true);
|
||||
// }
|
||||
}
|
||||
else
|
||||
app_error() << "You must specify the \"source\" particleset for a three-body Jastrow.\n";
|
||||
return true;
|
||||
// // Find the number of the source species
|
||||
// SpeciesSet &sSet = sourcePtcl->getSpeciesSet();
|
||||
// int numSpecies = sSet.getTotalNum();
|
||||
// bool success=false;
|
||||
// while (kids != NULL)
|
||||
// {
|
||||
// std::string kidsname = (char*)kids->name;
|
||||
// if (kidsname == "correlation")
|
||||
// {
|
||||
// RealType cusp=0.0;
|
||||
// std::string elementType;
|
||||
// OhmmsAttributeSet rAttrib;
|
||||
// rAttrib.add(elementType,"elementType");
|
||||
// rAttrib.add(cusp,"cusp");
|
||||
// rAttrib.put(kids);
|
||||
// BsplineFunctor<double> *functor = new BsplineFunctor<double>(cusp);
|
||||
// functor->elementType = elementType;
|
||||
// int ig = sSet.findSpecies (elementType);
|
||||
// if (ig < numSpecies)
|
||||
// {//ignore
|
||||
// functor->put (kids);
|
||||
// if (functor->cutoff_radius < 1.0e-6) {
|
||||
// app_log() << " BsplineFunction rcut is currently zero.\n"
|
||||
// << " Setting to Wigner-Seitz radius = "
|
||||
// << sourcePtcl->Lattice.WignerSeitzRadius << std::endl;
|
||||
// functor->cutoff_radius = sourcePtcl->Lattice.WignerSeitzRadius;
|
||||
// functor->reset();
|
||||
// }
|
||||
// J1->addFunc (ig,functor);
|
||||
// success = true;
|
||||
// dJ1->addFunc(ig,functor);
|
||||
// if(ReportLevel)
|
||||
// {
|
||||
// std::string fname="J1."+elementType+".dat";
|
||||
// std::ofstream fout(fname.c_str());
|
||||
// fout.setf(std::ios::scientific, std::ios::floatfield);
|
||||
// fout << "# One-body Jastrow generated by BsplineJastrowBuilder" << std::endl;
|
||||
// functor->print(fout);
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// kids = kids->next;
|
||||
// }
|
||||
// if(success)
|
||||
// {
|
||||
// //assign derivatives to J1
|
||||
// //dJ1->initialize();
|
||||
// //J1->setDiffOrbital(dJ1);
|
||||
// J1->dPsi=dJ1;
|
||||
// targetPsi.addOrbital(J1,"J1_bspline");
|
||||
// J1->setOptimizable(true);
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// PRE.warning("eeI_JastrowBuilder failed to add an One-Body Jastrow.");
|
||||
// delete J1;
|
||||
// delete dJ1;
|
||||
// }
|
||||
// }
|
||||
// // Create a two-body Jastrow
|
||||
// else
|
||||
// {
|
||||
// typedef TwoBodyJastrowOrbital<BsplineFunctor<RealType> > J2Type;
|
||||
// typedef DiffTwoBodyJastrowOrbital<BsplineFunctor<RealType> > dJ2Type;
|
||||
// J2Type *J2 = new J2Type(targetPtcl,targetPsi.is_manager());
|
||||
// dJ2Type *dJ2 = new dJ2Type(targetPtcl);
|
||||
// SpeciesSet& species(targetPtcl.getSpeciesSet());
|
||||
// RealType q=species(0,species.addAttribute("charge"));
|
||||
// //std::map<std::string,RadFuncType*> functorMap;
|
||||
// while (kids != NULL)
|
||||
// {
|
||||
// std::string kidsname((const char*)kids->name);
|
||||
// if (kidsname == "correlation")
|
||||
// {
|
||||
// OhmmsAttributeSet rAttrib;
|
||||
// RealType cusp=-1e10;
|
||||
// std::string pairType("0");
|
||||
// std::string spA(species.speciesName[0]);
|
||||
// std::string spB(species.speciesName[0]);
|
||||
// rAttrib.add(spA,"speciesA");
|
||||
// rAttrib.add(spB,"speciesB");
|
||||
// rAttrib.add(pairType,"pairType");
|
||||
// rAttrib.add(cusp,"cusp");
|
||||
// rAttrib.put(kids);
|
||||
// if(pairType[0]=='0')
|
||||
// {
|
||||
// pairType=spA+spB;
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// PRE.warning("pairType is deprecated. Use speciesA/speciesB");
|
||||
// //overwrite the species
|
||||
// spA=pairType[0]; spB=pairType[1];
|
||||
// }
|
||||
// int ia = species.findSpecies(spA);
|
||||
// int ib = species.findSpecies(spB);
|
||||
// if(ia==species.size() || ib == species.size())
|
||||
// {
|
||||
// PRE.error("Failed. Species are incorrect.",true);
|
||||
// }
|
||||
// if(cusp<-1e6)
|
||||
// {
|
||||
// if(ia==ib)
|
||||
// cusp=-0.25*q*q;
|
||||
// else
|
||||
// cusp=-0.5*q*q;
|
||||
// }
|
||||
// app_log() << " eeI_JastrowBuilder adds a functor with cusp = " << cusp << std::endl;
|
||||
// RadFuncType *functor = new RadFuncType(cusp);
|
||||
// functor->put (kids);
|
||||
// functor->elementType=pairType;
|
||||
// if (functor->cutoff_radius < 1.0e-6) {
|
||||
// app_log() << " BsplineFunction rcut is currently zero.\n"
|
||||
// << " Setting to Wigner-Seitz radius = "
|
||||
// << targetPtcl.Lattice.WignerSeitzRadius << std::endl;
|
||||
// functor->cutoff_radius = targetPtcl.Lattice.WignerSeitzRadius;
|
||||
// functor->reset();
|
||||
// }
|
||||
// J2->addFunc(pairType,ia,ib,functor);
|
||||
// dJ2->addFunc(pairType,ia,ib,functor);
|
||||
// if(ReportLevel)
|
||||
// {
|
||||
// std::string fname="J2."+pairType+".dat";
|
||||
// std::ofstream fout(fname.c_str());
|
||||
// fout.setf(std::ios::scientific, std::ios::floatfield);
|
||||
// fout << "# Two-body Jastrow generated by eeI_JastrowBuilder" << std::endl;
|
||||
// functor->print(fout);
|
||||
// }
|
||||
// }
|
||||
// kids = kids->next;
|
||||
// }
|
||||
// //dJ2->initialize();
|
||||
// //J2->setDiffOrbital(dJ2);
|
||||
// J2->dPsi=dJ2;
|
||||
// targetPsi.addOrbital(J2,"J2_bspline");
|
||||
// J2->setOptimizable(true);
|
||||
// }
|
||||
// return true;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue