mirror of https://github.com/QMCPACK/qmcpack.git
Remove old/unused geminal and model Jastrows
This commit is contained in:
Jaron Krogel
parent
8605709864
commit
df44e6b2b4
|
|
@ -1,54 +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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_BLOCKMATRIXFUNCTIONS_H
|
||||
#define QMCPLUSPLUS_BLOCKMATRIXFUNCTIONS_H
|
||||
namespace qmcplusplus
|
||||
{
|
||||
template<typename T, unsigned D>
|
||||
struct GEMV
|
||||
{
|
||||
static inline
|
||||
void apply(const T* restrict a, const T* restrict x, T* restrict y, int n, int m)
|
||||
{
|
||||
for(int i=0; i<n; i++)
|
||||
{
|
||||
T tmp = 0;
|
||||
const T* restrict xcopy=x;
|
||||
for(int j=0; j<m; j++ )
|
||||
{
|
||||
tmp += (*a++)*(*xcopy++);
|
||||
}
|
||||
(*y++) = tmp;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
/** specialization for D=3
|
||||
*/
|
||||
template<typename T>
|
||||
struct GEMV<T,3>
|
||||
{
|
||||
static inline
|
||||
void apply(const T* restrict a, const T* restrict x, T* restrict y, int n, int m)
|
||||
{
|
||||
y[0]=a[0]*x[0]+a[1]*x[1]+a[2]*x[2];
|
||||
y[1]=a[3]*x[0]+a[4]*x[1]+a[5]*x[2];
|
||||
y[2]=a[6]*x[0]+a[7]*x[1]+a[8]*x[2];
|
||||
}
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -54,7 +54,6 @@ SET(JASTROW_SRCS
|
|||
Jastrow/kSpaceJastrowBuilder.cpp
|
||||
Jastrow/RPAJastrow.cpp
|
||||
Jastrow/singleRPAJastrowBuilder.cpp
|
||||
Jastrow/JAABuilder.cpp
|
||||
Jastrow/JABBuilder.cpp
|
||||
IonOrbital.cpp
|
||||
IonOrbitalBuilder.cpp
|
||||
|
|
@ -81,17 +80,12 @@ IF(OHMMS_DIM MATCHES 3)
|
|||
|
||||
SET(JASTROW_SRCS ${JASTROW_SRCS}
|
||||
Jastrow/eeI_JastrowBuilder.cpp
|
||||
Jastrow/ThreeBodyGeminal.cpp
|
||||
Jastrow/ThreeBodyBlockSparse.cpp
|
||||
Jastrow/JastrowBasisBuilder.cpp
|
||||
Jastrow/CBSOBuilder.cpp
|
||||
)
|
||||
|
||||
SET(FERMION_SRCS ${FERMION_SRCS}
|
||||
MolecularOrbitals/STOBuilder.cpp
|
||||
MolecularOrbitals/GTOBuilder.cpp
|
||||
MolecularOrbitals/NGOBuilder.cpp
|
||||
MolecularOrbitals/BsplineAOBuilder.cpp
|
||||
)
|
||||
|
||||
IF(NOT QMC_COMPLEX AND ENABLE_SOA)
|
||||
|
|
|
|||
|
|
@ -1,274 +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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_BESSELZERO_FUNCTOR_H
|
||||
#define QMCPLUSPLUS_BESSELZERO_FUNCTOR_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
#include "Utilities/ProgressReportEngine.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
#include <cstdio>
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
template<class T>
|
||||
struct BesselZero: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
typedef real_type value_type;
|
||||
int NumParams;
|
||||
std::vector<real_type> Parameters;
|
||||
real_type R_B,R_Binv;
|
||||
real_type C_0,C_0inv;
|
||||
std::string pr;
|
||||
|
||||
///constructor
|
||||
BesselZero():R_B(1.0)
|
||||
{
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new BesselZero(*this);
|
||||
}
|
||||
|
||||
void resize(int n)
|
||||
{
|
||||
NumParams = n;
|
||||
Parameters.resize(n,0);
|
||||
}
|
||||
|
||||
void reset()
|
||||
{
|
||||
if (pr=="yes")
|
||||
print();
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
real_type u=0;
|
||||
real_type A(C_0);
|
||||
real_type Ainv(C_0inv);
|
||||
real_type AX(A*r);
|
||||
real_type Xinv(1.0/r);
|
||||
real_type AXinv(Ainv*Xinv);
|
||||
for(int i=0; i<NumParams; i++)
|
||||
{
|
||||
real_type j(i+1);
|
||||
real_type siniAX = std::sin(j*AX);
|
||||
real_type iinv(1.0/j);
|
||||
u += Parameters[i]*AXinv*iinv*siniAX;
|
||||
}
|
||||
return u;
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r, real_type rinv)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
inline void evaluateAll(real_type r, real_type rinv)
|
||||
{
|
||||
real_type du(0),d2u(0);
|
||||
evaluate(r,du,d2u);
|
||||
}
|
||||
|
||||
inline real_type
|
||||
evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
real_type u(0);
|
||||
sumB(r,u,dudr,d2udr2);
|
||||
return u;
|
||||
}
|
||||
|
||||
|
||||
inline real_type
|
||||
evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type &d3udr3)
|
||||
{
|
||||
real_type u(0);
|
||||
sumB(r,u,dudr,d2udr2);
|
||||
return u;
|
||||
}
|
||||
|
||||
inline void sumB(real_type X, real_type& u, real_type& du, real_type& d2u)
|
||||
{
|
||||
u=0;
|
||||
du=0;
|
||||
d2u=0;
|
||||
real_type A(C_0);
|
||||
real_type Ainv(C_0inv);
|
||||
real_type AX(A*X);
|
||||
real_type AX2(AX*AX);
|
||||
real_type Xinv(1.0/X);
|
||||
real_type AXinv(Ainv*Xinv);
|
||||
real_type AX2inv(Ainv*Xinv*Xinv);
|
||||
real_type AX3inv(Ainv*Xinv*Xinv*Xinv);
|
||||
for(int i=0; i<NumParams; i++)
|
||||
{
|
||||
real_type j(i+1.0);
|
||||
real_type cosjAX = std::cos(j*AX);
|
||||
real_type sinjAX = std::sin(j*AX);
|
||||
real_type jinv(1.0/j);
|
||||
u += Parameters[i]*AXinv*jinv*sinjAX;
|
||||
du += Parameters[i]*(j*AX*cosjAX - sinjAX)*AX2inv*jinv;
|
||||
d2u += Parameters[i]*(-2.0*j*AX*cosjAX + (2.0-j*j*AX2)*sinjAX)*AX3inv*jinv;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
inline bool
|
||||
evaluateDerivatives (real_type r, std::vector<TinyVector<real_type,3> >& derivs)
|
||||
{
|
||||
real_type A(C_0);
|
||||
real_type Ainv(C_0inv);
|
||||
real_type X(r);
|
||||
real_type Xinv(1.0/r);
|
||||
real_type AX(A*X);
|
||||
real_type AX2(AX*AX);
|
||||
real_type AXinv(Ainv*Xinv);
|
||||
real_type AX2inv(Ainv*Xinv*Xinv);
|
||||
real_type AX3inv(Ainv*Xinv*Xinv*Xinv);
|
||||
for(int i=0; i<NumParams; i++)
|
||||
{
|
||||
real_type j(i+1.0);
|
||||
real_type cosjAX = std::cos(j*AX);
|
||||
real_type sinjAX = std::sin(j*AX);
|
||||
real_type jinv(1.0/j);
|
||||
derivs[i][0] += AXinv*jinv*sinjAX;
|
||||
derivs[i][1] += (j*AX*cosjAX - sinjAX)*AX2inv*jinv;
|
||||
derivs[i][2] += (-2.0*j*AX*cosjAX + (2.0-j*j*AX2)*sinjAX)*AX3inv*jinv;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
inline real_type f(real_type r)
|
||||
{
|
||||
real_type du, d2u;
|
||||
return evaluate (r, du, d2u);
|
||||
}
|
||||
|
||||
inline real_type df(real_type r)
|
||||
{
|
||||
real_type du, d2u;
|
||||
evaluate (r, du, d2u);
|
||||
return du;
|
||||
}
|
||||
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
ReportEngine PRE("BesselZero","put(xmlNodePtr)");
|
||||
//CuspValue = -1.0e10;
|
||||
NumParams = 0;
|
||||
pr = "no";
|
||||
OhmmsAttributeSet rAttrib;
|
||||
rAttrib.add(NumParams, "size");
|
||||
rAttrib.add(R_B, "RB");
|
||||
rAttrib.add(pr, "print");
|
||||
rAttrib.put(cur);
|
||||
app_log()<<" R_B is set to: "<<R_B<< std::endl;
|
||||
R_Binv = 1.0/R_B;
|
||||
C_0 = 3.1415926535897932384626433832795028841968*R_Binv;
|
||||
C_0inv = 1.0/C_0;
|
||||
if (NumParams == 0)
|
||||
{
|
||||
PRE.error("You must specify a positive number of parameters for the BesselZero jastrow function.",true);
|
||||
}
|
||||
app_log() << " size = " << NumParams << " parameters " << std::endl;
|
||||
resize (NumParams);
|
||||
// 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 BesselZero." + "Resetting to \"Array\"");
|
||||
xmlNewProp (xmlCoefs, (const xmlChar*) "type", (const xmlChar*) "Array");
|
||||
}
|
||||
std::vector<real_type> params;
|
||||
putContent(params, xmlCoefs);
|
||||
if (params.size() == NumParams)
|
||||
Parameters = params;
|
||||
// Setup parameter names
|
||||
for (int i=0; i< NumParams; i++)
|
||||
{
|
||||
std::stringstream sstr;
|
||||
sstr << id << "_" << i;
|
||||
myVars.insert(sstr.str(),Parameters[i],true);
|
||||
}
|
||||
app_log() << "Parameter Name Value\n";
|
||||
myVars.print(app_log());
|
||||
}
|
||||
xmlCoefs = xmlCoefs->next;
|
||||
}
|
||||
reset();
|
||||
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)
|
||||
{
|
||||
for(int i=0; i<Parameters.size(); ++i)
|
||||
{
|
||||
int loc=myVars.where(i);
|
||||
if(loc>=0)
|
||||
Parameters[i]=myVars[i]=active[loc];
|
||||
}
|
||||
reset();
|
||||
}
|
||||
|
||||
|
||||
void print()
|
||||
{
|
||||
// std::string fname = (elementType != "") ? elementType : pairType;
|
||||
std::string fname = "BesselZero.dat";
|
||||
// //cerr << "Writing " << fname << " file.\n";
|
||||
FILE *fout = fopen (fname.c_str(), "w");
|
||||
for (double r=0.0; r<R_B*4; r+=0.01)
|
||||
fprintf (fout, "%8.3f %16.10f\n", r, evaluate(r));
|
||||
fclose(fout);
|
||||
}
|
||||
|
||||
|
||||
void print(std::ostream& os)
|
||||
{
|
||||
// int n=100;
|
||||
// T d=cutoff_radius/100.,r=0;
|
||||
// T 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
|
||||
|
|
@ -1,148 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#include "Numerics/LibxmlNumericIO.h"
|
||||
#include "Numerics/GaussianBasisSet.h"
|
||||
#include "Numerics/SlaterBasisSet.h"
|
||||
#include "QMCWaveFunctions/Jastrow/CBSOBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/WMFunctor.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
CBSOBuilder::CBSOBuilder(xmlNodePtr cur):
|
||||
Normalized(true),m_rcut(5.0),NumGridPoints(51)
|
||||
{
|
||||
if(cur != NULL)
|
||||
{
|
||||
putCommon(cur);
|
||||
}
|
||||
}
|
||||
|
||||
bool
|
||||
CBSOBuilder::putCommon(xmlNodePtr cur)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
void
|
||||
CBSOBuilder::setOrbitalSet(CenteredOrbitalType* oset, const std::string& acenter)
|
||||
{
|
||||
m_orbitals = oset;
|
||||
m_species = acenter;
|
||||
}
|
||||
|
||||
bool
|
||||
CBSOBuilder::addGrid(xmlNodePtr cur)
|
||||
{
|
||||
if(!m_orbitals)
|
||||
{
|
||||
ERRORMSG("m_orbitals, SphericalOrbitals<ROT,GT>*, is not initialized")
|
||||
return false;
|
||||
}
|
||||
RealType delta=-1.0;
|
||||
OhmmsAttributeSet attrib;
|
||||
attrib.add(m_rcut,"rf");
|
||||
attrib.add(NumGridPoints,"npts");
|
||||
attrib.add(delta,"step");
|
||||
attrib.put(cur);
|
||||
if(delta>0)
|
||||
NumGridPoints = static_cast<int>(m_rcut/delta)+1;
|
||||
app_log() << " CBSOBuilder::addGrid Rcut = " << m_rcut << " NumGridPoints = " << NumGridPoints << std::endl;
|
||||
return true;
|
||||
}
|
||||
|
||||
/** Add a new Slater Type Orbital with quantum numbers \f$(n,l,m,s)\f$
|
||||
* \param cur the current xmlNode to be processed
|
||||
* \param nlms a vector containing the quantum numbers \f$(n,l,m,s)\f$
|
||||
* \return true is succeeds
|
||||
*
|
||||
This function puts the STO on a logarithmic grid and calculates the boundary
|
||||
conditions for the 1D Cubic Spline. The derivates at the endpoint
|
||||
are assumed to be all zero. Note: for the radial orbital we use
|
||||
\f[ f(r) = \frac{R(r)}{r^l}, \f] where \f$ R(r) \f$ is the usual
|
||||
radial orbital and \f$ l \f$ is the angular momentum.
|
||||
*/
|
||||
bool
|
||||
CBSOBuilder::addRadialOrbital(xmlNodePtr cur, const QuantumNumberType& nlms)
|
||||
{
|
||||
if(!m_orbitals)
|
||||
{
|
||||
ERRORMSG("m_orbitals, SphericalOrbitals<ROT,GT>*, is not initialized")
|
||||
return false;
|
||||
}
|
||||
std::string radtype("Gaussian");
|
||||
OhmmsAttributeSet attrib;
|
||||
attrib.add(radtype,"type");
|
||||
attrib.put(cur);
|
||||
int lastRnl = m_orbitals->Rnl.size();
|
||||
m_nlms = nlms;
|
||||
int L= m_nlms[1];
|
||||
typedef RadialOrbitalType::FNIN InputFuncType;
|
||||
InputFuncType* infunc=0;
|
||||
//create one
|
||||
if(radtype == "Gaussian" || radtype == "GTO")
|
||||
{
|
||||
GaussianCombo<RealType> *gset = new GaussianCombo<RealType>(L,Normalized);
|
||||
gset->putBasisGroup(cur);
|
||||
app_log() << " " << L << " basisGroup contains " << gset->size() << " radial functors." << std::endl;
|
||||
infunc=gset;
|
||||
}
|
||||
else
|
||||
if(radtype == "Slater")
|
||||
{
|
||||
SlaterCombo<RealType> *sto_set=new SlaterCombo<RealType>(L,Normalized);
|
||||
sto_set->putBasisGroup(cur);
|
||||
infunc=sto_set;
|
||||
}
|
||||
else
|
||||
if(radtype == "WM")
|
||||
{
|
||||
xmlNodePtr tcur=cur->children;
|
||||
int nr=0;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if(cname == "radfunc")
|
||||
{
|
||||
OhmmsAttributeSet rAttrib;
|
||||
RealType rcut=m_rcut;//use the global cutoff but can overwrite
|
||||
rAttrib.add(rcut,"rcut");
|
||||
rAttrib.put(tcur);
|
||||
infunc = new WMFunctor<RealType>(1.0,rcut);
|
||||
infunc->put(tcur);
|
||||
nr++;
|
||||
}
|
||||
tcur=tcur->next;
|
||||
}
|
||||
}
|
||||
if(infunc)
|
||||
{
|
||||
app_log()
|
||||
<< " CBSOBuilder::addRadialOrbital input " << radtype
|
||||
<< " output = CubicBspline " << std::endl;
|
||||
RadialOrbitalType *radorb=new RadialOrbitalType;
|
||||
radorb->initialize(infunc,m_rcut,NumGridPoints);
|
||||
m_orbitals->Rnl.push_back(radorb);
|
||||
m_orbitals->RnlID.push_back(m_nlms);
|
||||
//if(lastRnl && m_orbitals->Rnl.size()> lastRnl) {
|
||||
// //LOGMSG("\tSetting GridManager of " << lastRnl << " radial orbital to false")
|
||||
// m_orbitals->Rnl[lastRnl]->setGridManager(false);
|
||||
//}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -1,80 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_SPHERICALBASISBUILDER_H
|
||||
#define QMCPLUSPLUS_SPHERICALBASISBUILDER_H
|
||||
|
||||
#include "Configuration.h"
|
||||
#include "QMCWaveFunctions/Jastrow/SplineFunctors.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/SphericalBasisSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** Builder of cubic-bspline radial orbitals associated with a center
|
||||
*
|
||||
* CBSO stands for CubicBSplineOrbitals
|
||||
* For a center,
|
||||
* - only one linear grid is used
|
||||
* - any number of radial orbitals of CubicBsplineSingle
|
||||
* Will test the group function later using CubicBsplineGroup<T,GRIDTYPE>
|
||||
*/
|
||||
class CBSOBuilder: public QMCTraits
|
||||
{
|
||||
//typedef CubicBsplineGroup<RealType,LINEAR_1DGRID> RadialOrbitalType;
|
||||
public:
|
||||
///spline engine
|
||||
typedef CubicBspline<RealType,LINEAR_1DGRID,FIRSTDERIV_CONSTRAINTS> SplineEngineType;
|
||||
///numerical functor
|
||||
typedef CubicSplineSingle<RealType,SplineEngineType> RadialOrbitalType;
|
||||
///spherical basis set using splines
|
||||
typedef SphericalBasisSet<RadialOrbitalType> CenteredOrbitalType;
|
||||
|
||||
///true, if the RadialOrbitalType is normalized
|
||||
bool Normalized;
|
||||
///number of grid points [0,Npts]
|
||||
int NumGridPoints;
|
||||
///maximum cutoff
|
||||
RealType m_rcut;
|
||||
///the quantum number of this node
|
||||
QuantumNumberType m_nlms;
|
||||
///the radial orbitals
|
||||
CenteredOrbitalType* m_orbitals;
|
||||
///the species
|
||||
std::string m_species;
|
||||
|
||||
///constructor
|
||||
CBSOBuilder(xmlNodePtr cur=NULL);
|
||||
|
||||
///assign a CenteredOrbitalType to work on
|
||||
void setOrbitalSet(CenteredOrbitalType* oset, const std::string& acenter);
|
||||
|
||||
///add a grid
|
||||
bool addGrid(xmlNodePtr cur);
|
||||
|
||||
/** add a radial functor
|
||||
* @param cur xml element
|
||||
* @param nlms quantum number
|
||||
*/
|
||||
bool addRadialOrbital(xmlNodePtr cur, const QuantumNumberType& nlms);
|
||||
|
||||
/** put common element
|
||||
* @param cur xml element
|
||||
*/
|
||||
bool putCommon(xmlNodePtr cur);
|
||||
|
||||
private:
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,480 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_GAUSSIAN_H
|
||||
#define QMCPLUSPLUS_GAUSSIAN_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
template<class T>
|
||||
struct GaussianFunctor: public OptimizableFunctorBase
|
||||
{
|
||||
/// f(r) = A exp(-(r)^2/C^2)
|
||||
///for smooth truncation f(r)->f(r)+f(2rc-r)-2f(rc)
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
|
||||
///coefficients
|
||||
real_type A;
|
||||
real_type C;
|
||||
real_type RC;
|
||||
real_type c0,c1,c2,c3,c4;
|
||||
|
||||
std::string ID_A;
|
||||
std::string ID_C;
|
||||
std::string ID_RC;
|
||||
bool optimizable;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
GaussianFunctor(real_type a=1.0, real_type c=1.0,real_type rc=1.0):ID_A("G_A"), ID_C("G_C"),ID_RC("G_RC"), optimizable(true)
|
||||
{
|
||||
A=a;
|
||||
C=c;
|
||||
RC=rc;
|
||||
reset();
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new GaussianFunctor<T>(*this);
|
||||
}
|
||||
|
||||
inline void reset()
|
||||
{
|
||||
c0 = -1.0/(C*C);
|
||||
c1 = 2*RC;
|
||||
c2 = 2*A*std::exp(-(RC*RC)/(C*C));
|
||||
c3 = -2.0*A/(C*C);
|
||||
c4 = 4.0*A/(C*C*C*C);
|
||||
}
|
||||
|
||||
/** reset the internal variables.
|
||||
*
|
||||
* USE_resetParameters
|
||||
*/
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
if (!optimizable)
|
||||
return;
|
||||
int ia=myVars.where(0);
|
||||
if (ia>-1)
|
||||
A=myVars[0]=active[ia];
|
||||
int ic=myVars.where(1);
|
||||
if (ic>-1)
|
||||
C=myVars[1]=active[ic];
|
||||
// int id=myVars.where(2);
|
||||
// if (id>-1) RC=myVars[2]=active[id];
|
||||
reset();
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
if (optimizable)
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
if (optimizable)
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
real_type expart = std::exp(c0*r*r);
|
||||
real_type r1=c1-r;
|
||||
real_type expartRC = std::exp(c0*r1*r1);
|
||||
return A*(expart+expartRC)+c2;
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@param d3udr3 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr,
|
||||
real_type& d2udr2, real_type& d3udr3)
|
||||
{
|
||||
std::cerr << "Third derivative not implemented for GaussianFunctor.\n";
|
||||
if (r<RC)
|
||||
{
|
||||
real_type expart = std::exp(c0*r*r);
|
||||
real_type r1=c1-r;
|
||||
real_type expartRC = std::exp(c0*r1*r1);
|
||||
dudr = c3*(r*expart - r1*expartRC);
|
||||
d2udr2 = (c4*r*r+c3)*expart + (c4*r1*r1+c3)*expartRC;
|
||||
return A*(expart+expartRC)+c2;
|
||||
}
|
||||
else
|
||||
{
|
||||
dudr = 0.0;
|
||||
d2udr2 = 0.0;
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
real_type expart = std::exp(c0*r*r);
|
||||
real_type r1=c1-r;
|
||||
real_type expartRC = std::exp(c0*r1*r1);
|
||||
dudr = c3*(r*expart - r1*expartRC);
|
||||
d2udr2 = (c4*r*r+c3)*expart + (c4*r1*r1+c3)*expartRC;
|
||||
return A*(expart+expartRC)+c2;
|
||||
}
|
||||
else
|
||||
{
|
||||
dudr = 0.0;
|
||||
d2udr2 = 0.0;
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
real_type r1=(c1-r);
|
||||
return c3*(r*std::exp(c0*r*r) - r1*std::exp(c0*r1*r1));
|
||||
}
|
||||
else
|
||||
return 0.0;
|
||||
}
|
||||
inline bool
|
||||
evaluateDerivatives(real_type r, std::vector<TinyVector<real_type,3> >& derivs)
|
||||
{
|
||||
if ((r >= RC)||(!optimizable))
|
||||
return false;
|
||||
else
|
||||
{
|
||||
real_type rcmr =(c1-r);
|
||||
real_type rcmr2 =rcmr*rcmr;
|
||||
real_type r2 =r*r;
|
||||
real_type expart = std::exp(c0*r2);
|
||||
real_type expartRC = std::exp(c0*rcmr2);
|
||||
real_type Am1 = -2.0/(C*C);
|
||||
real_type cm2 = 1.0/(C*C);
|
||||
real_type cm3 = 1.0/(C*C*C);
|
||||
real_type cm4 = cm2*cm2;
|
||||
real_type cm5 = cm2*cm3;
|
||||
//Derivs of function
|
||||
//
|
||||
// df_dA
|
||||
derivs[0][0]= (expart+expartRC)+c2*Am1 ;
|
||||
// df_dC
|
||||
derivs[1][0]= 2*A*( expart*r2 - 2.0*expartRC*rcmr2 + expartRC*rcmr2)*cm3 ;
|
||||
// df_dRC
|
||||
// derivs[2][0]= (2.0*c2*RC - 2.0*A*expartRC*rcmr)*cm2 ;
|
||||
//grad derivs
|
||||
// dr_df_dA
|
||||
derivs[0][1]= Am1*(r*expart - rcmr*expartRC) ;
|
||||
// dr_df_dC
|
||||
derivs[1][1]= 4.0*A*(expart*r*cm3*(1.0-r2*cm2)-expartRC*rcmr*cm3*(1.0-rcmr2*cm2)) ;
|
||||
// dr_df_dRC
|
||||
// derivs[2][1]= 4.0*A*expartRC*cm2*(1.0-2.0*rcmr2*cm2) ;
|
||||
//lap derivs
|
||||
// dr2_df_dA
|
||||
derivs[0][2]= Am1*((Am1*r2+1)*expart + (Am1*rcmr2 + 1)*expartRC);
|
||||
// dr2_df_dC
|
||||
derivs[1][2]= 4.0*A*cm3*(expart*(1.0-5.0*r2*cm2+2.0*r2*r2*cm4)+expartRC*(1.0-5.0*rcmr2*cm2+2.0*rcmr2*rcmr2*cm4)) ;
|
||||
// dr2_df_dRC
|
||||
// derivs[2][2]= 24*A*expartRC*rcmr*cm4 - 16.0*A*rcmr*rcmr2*expartRC*cm3*cm3;
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
RC = cutoff_radius;
|
||||
std::string opt("yes");
|
||||
OhmmsAttributeSet Tattrib;
|
||||
Tattrib.add(opt,"optimize");
|
||||
Tattrib.put(cur);
|
||||
optimizable = (opt=="yes");
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while (tcur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if (cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
// std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if (aname == "A")
|
||||
{
|
||||
ID_A = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if (aname == "C")
|
||||
{
|
||||
ID_C = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(C,tcur);
|
||||
}
|
||||
else
|
||||
if (aname == "RC")
|
||||
{
|
||||
ID_RC = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(RC,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
if (optimizable)
|
||||
{
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_C,C);
|
||||
// myVars.insert(ID_RC,RC);
|
||||
}
|
||||
reset();
|
||||
return true;
|
||||
}
|
||||
};
|
||||
|
||||
template<class T>
|
||||
struct TruncatedShiftedGaussianFunctor: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
///coefficients
|
||||
real_type A, B, C, RC;
|
||||
real_type c0,c1,c2,c3,c4,c5,c6;
|
||||
///id of A
|
||||
std::string ID_A;
|
||||
///id of B
|
||||
std::string ID_B;
|
||||
///id of rc
|
||||
std::string ID_RC;
|
||||
std::string ID_C;
|
||||
///constructor
|
||||
TruncatedShiftedGaussianFunctor(real_type a=1.0 , real_type b=1.0 , real_type c=1.0):ID_A("SG_A"), ID_B("SG_B"),ID_C("SG_C")
|
||||
{
|
||||
A=a;
|
||||
B=b;
|
||||
C=c;
|
||||
reset();
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new TruncatedShiftedGaussianFunctor(*this);
|
||||
}
|
||||
|
||||
/**
|
||||
*@brief reset the internal variables.
|
||||
*/
|
||||
void reset()
|
||||
{
|
||||
c0 = -1.0/(C*C);
|
||||
c1 = std::exp(c0*B*B);
|
||||
c2 = A*c1/(C*C);
|
||||
c3 = -B*B*c2;
|
||||
c4 = 2.0*A*c0;
|
||||
c5 = 2.0*c2;
|
||||
c6 = 4.0*A*c0*c0;
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@return \f$ u(r) = a/(1+br^2) \f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
r-=B;
|
||||
real_type r2=r*r;
|
||||
return A*(std::exp(c0*r2)-c1)+c2*r2+c3;
|
||||
}
|
||||
else
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@param d3udr3 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr,
|
||||
real_type& d2udr2, real_type& d3udr3)
|
||||
{
|
||||
std::cerr << "Third derivative not implemented for GaussianFunctor.\n";
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr return value \f$ du/dr = -2abr/(1+br^2)^2 \f$
|
||||
@param d2udr2 return value \f$ d^2u/dr^2 =
|
||||
-2ab(1-3br^2)/(1+br^2)^3 \f$
|
||||
@return \f$ u(r) = a/(1+br^2) \f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
r-=B;
|
||||
real_type r2=r*r;
|
||||
real_type expart = std::exp(c0*r2);
|
||||
dudr = c4*r*expart+c5*r;
|
||||
d2udr2 = (c4+c6*r2)*expart+c5;
|
||||
return A*(expart-c1)+c2*r2+c3;
|
||||
}
|
||||
else
|
||||
{
|
||||
dudr=0.0;
|
||||
d2udr2 = 0.0;
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
real_type df(real_type r)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
r-=B;
|
||||
return c4*r*std::exp(c0*r*r)+c5*r;
|
||||
}
|
||||
else
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/** implements virtual function
|
||||
* @param cur xml node
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
RC = cutoff_radius;
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while (tcur != NULL)
|
||||
{
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if (cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if (aname == "A")
|
||||
{
|
||||
ID_A = idname;
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if (aname == "B")
|
||||
{
|
||||
ID_B = idname;
|
||||
putContent(B,tcur);
|
||||
}
|
||||
else
|
||||
if (aname == "C")
|
||||
{
|
||||
ID_C = idname;
|
||||
putContent(C,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
///Sim cell too small
|
||||
if (B<0.5*cutoff_radius)
|
||||
RC=2.0*B;
|
||||
else
|
||||
return false;
|
||||
reset();
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_B,B);
|
||||
myVars.insert(ID_C,C);
|
||||
LOGMSG(" TruncatedShiftedGaussianFunctor Parameters ")
|
||||
LOGMSG(" A (" << ID_A << ") = " << A << " B (" << ID_B << ") = " << B<<" C (" << ID_C << ") = " << C << " R_c (" << ID_RC << ") = " << RC)
|
||||
return true;
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
/** reset the internal variables.
|
||||
*
|
||||
* USE_resetParameters
|
||||
*/
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ia=myVars.where(0);
|
||||
if (ia>-1)
|
||||
A=active[ia];
|
||||
int ib=myVars.where(1);
|
||||
if (ib>-1)
|
||||
B=active[ib];
|
||||
int ic=myVars.where(2);
|
||||
if (ic>-1)
|
||||
C=active[ic];
|
||||
if (B>0.5*cutoff_radius)
|
||||
{
|
||||
B=0.5*cutoff_radius;
|
||||
RC=cutoff_radius;
|
||||
}
|
||||
else
|
||||
RC=2.0*B;
|
||||
reset();
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,190 +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: D.C. Yang, University of Illinois at Urbana-Champaign
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#include "Particle/DistanceTableData.h"
|
||||
#include "Particle/DistanceTable.h"
|
||||
#include "QMCWaveFunctions/Jastrow/JAABuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ModPadeFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/PadeFunctors.h"
|
||||
#include "QMCWaveFunctions/Jastrow/McMillanJ2Functor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/McMillanJ2GFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/GaussianFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/TwoBodyJastrowOrbital.h"
|
||||
#include "QMCWaveFunctions/Jastrow/DiffTwoBodyJastrowOrbital.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
JAABuilder::JAABuilder(ParticleSet& p, TrialWaveFunction& psi):OrbitalBuilderBase(p,psi),
|
||||
IgnoreSpin(true)
|
||||
{
|
||||
}
|
||||
/** Create a two-body Jatrow function with a template
|
||||
*@param cur the current xmlNode
|
||||
*@param dummy null pointer used to identify FN
|
||||
*
|
||||
*The template class JeeType is a functor which handles the
|
||||
*evaluation of the function value, gradients and laplacians using
|
||||
*distance tables. This is a specialized builder function for
|
||||
*spin-dependent Jastrow function,e.g., for electrons, two functions
|
||||
*are created for uu(dd) and ud(du).
|
||||
*/
|
||||
template <class FN> TwoBodyJastrowOrbital<FN>* JAABuilder::createJAA(xmlNodePtr cur, const std::string& jname)
|
||||
{
|
||||
std::string corr_tag("correlation");
|
||||
int ng = targetPtcl.groups();
|
||||
int ia=0, ib=0, iab=0;
|
||||
xmlNodePtr gridPtr=NULL;
|
||||
cur = cur->children;
|
||||
const SpeciesSet& species(targetPtcl.getSpeciesSet());
|
||||
typedef TwoBodyJastrowOrbital<FN> JeeType;
|
||||
int taskid=(targetPsi.is_manager())?targetPsi.getGroupID():-1;
|
||||
JeeType *J2 = new JeeType(targetPtcl,taskid);
|
||||
typedef DiffTwoBodyJastrowOrbital<FN> dJ2Type;
|
||||
dJ2Type *dJ2 = new dJ2Type(targetPtcl);
|
||||
RealType rc=targetPtcl.Lattice.WignerSeitzRadius;
|
||||
int pairs=0;
|
||||
while (cur != NULL)
|
||||
{
|
||||
std::string cname((const char*)(cur->name));
|
||||
if (cname == corr_tag)
|
||||
{
|
||||
std::string spA("u");
|
||||
std::string spB("u");
|
||||
OhmmsAttributeSet rAttrib;
|
||||
rAttrib.add(spA, "speciesA");
|
||||
rAttrib.add(spA, "species1");
|
||||
rAttrib.add(spB, "speciesB");
|
||||
rAttrib.add(spB, "species2");
|
||||
rAttrib.put(cur);
|
||||
if (spA==targetPsi.getName()) //could have used the particle name
|
||||
{
|
||||
spA=species.speciesName[0];
|
||||
spB=species.speciesName[0];
|
||||
}
|
||||
int ia = species.findSpecies(spA);
|
||||
int ib = species.findSpecies(spB);
|
||||
if (ia==species.size() || ia == species.size())
|
||||
{
|
||||
APP_ABORT("JAABuilder::createJAA is trying to use invalid species");
|
||||
}
|
||||
std::string pairID=spA+spB;
|
||||
FN *j= new FN;
|
||||
j->cutoff_radius=rc;
|
||||
j->put(cur);
|
||||
J2->addFunc( ia,ib,j);
|
||||
dJ2->addFunc( ia,ib,j);
|
||||
++pairs;
|
||||
}
|
||||
cur = cur->next;
|
||||
} // while cur
|
||||
if (pairs)
|
||||
{
|
||||
J2->dPsi=dJ2;
|
||||
std::string j2name="J2_"+jname;
|
||||
targetPsi.addOrbital(J2,j2name);
|
||||
return J2;
|
||||
}
|
||||
else
|
||||
{
|
||||
//clean up and delete the twobody orbitals
|
||||
APP_ABORT("JAABuilder::put Failed to create Two-Body with "+jname);
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
bool JAABuilder::put(xmlNodePtr cur)
|
||||
{
|
||||
std::string spinOpt("no");
|
||||
std::string typeOpt("Two-Body");
|
||||
std::string jastfunction("pade");
|
||||
OhmmsAttributeSet aAttrib;
|
||||
aAttrib.add(spinOpt,"spin");
|
||||
aAttrib.add(typeOpt,"type");
|
||||
aAttrib.add(jastfunction,"function");
|
||||
aAttrib.put(cur);
|
||||
IgnoreSpin=(spinOpt=="no");
|
||||
OrbitalBase* newJ = 0;
|
||||
//if(jastfunction == "pade") {
|
||||
// app_log() << " Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
// PadeJastrow<RealType> *dummy = 0;
|
||||
// success = createJAA(cur,dummy);
|
||||
//} else
|
||||
if (jastfunction == "short")
|
||||
{
|
||||
app_log() << " Modified Jastrow function Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
//ModPadeFunctor<RealType> *dummy = 0;
|
||||
newJ = createJAA<ModPadeFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "modmcmillan")
|
||||
{
|
||||
app_log() << " Modified McMillan Jastrow function Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
newJ = createJAA<ModMcMillanJ2Functor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "combomcmillan")
|
||||
{
|
||||
app_log() << " Combo McMillan Jastrow function Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
newJ = createJAA<comboMcMillanJ2Functor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "mcmillan")
|
||||
{
|
||||
app_log() << " McMillan (LONG RANGE!) Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
SpeciesSet& species(targetPtcl.getSpeciesSet());
|
||||
TwoBodyJastrowOrbital<McMillanJ2Functor<RealType> >* a = createJAA<McMillanJ2Functor<RealType> >(cur,jastfunction);
|
||||
species(species.addAttribute("J2_A"),species.addSpecies(species.speciesName[targetPtcl.GroupID[0]])) = (a->F[a->F.size()-1])->A;
|
||||
species(species.addAttribute("J2_B"),species.addSpecies(species.speciesName[targetPtcl.GroupID[0]])) = (a->F[a->F.size()-1])->B;
|
||||
newJ = a;
|
||||
}
|
||||
else if (jastfunction == "mcmillanj2g")
|
||||
{
|
||||
app_log() << " McMillan Two-Body Jastrow Function (Gaussian for r < 2.5) = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
SpeciesSet& species(targetPtcl.getSpeciesSet());
|
||||
TwoBodyJastrowOrbital<McMillanJ2GFunctor<RealType> >* a = createJAA<McMillanJ2GFunctor<RealType> >(cur,jastfunction);
|
||||
species(species.addAttribute("J2_A"),species.addSpecies(species.speciesName[targetPtcl.GroupID[0]])) = (a->F[a->F.size()-1])->A;
|
||||
species(species.addAttribute("J2_B"),species.addSpecies(species.speciesName[targetPtcl.GroupID[0]])) = (a->F[a->F.size()-1])->B;
|
||||
newJ = a;
|
||||
}
|
||||
else if (jastfunction == "gaussian")
|
||||
{
|
||||
app_log() << " Gaussian function Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
newJ = createJAA<GaussianFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "shiftedgaussian")
|
||||
{
|
||||
app_log() << " Gaussian function Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
IgnoreSpin=true;
|
||||
newJ = createJAA<TruncatedShiftedGaussianFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "padetwo2ndorderfunctor")
|
||||
{
|
||||
app_log() << " PadeTwo2ndOrderFunctor Jastrow function Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
//IgnoreSpin=true;
|
||||
newJ = createJAA<PadeTwo2ndOrderFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
//} else if(jastfunction == "rpa") {
|
||||
// app_log() << " Two-Body Jastrow Function = " << jastfunction << std::endl;
|
||||
// RPAJastrow<RealType> *dummy = 0;
|
||||
// success = createJAA(cur,dummy);
|
||||
//}
|
||||
return (newJ != 0);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,46 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// D.C. Yang, 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_ORIGINAL_JASTROW_AA_BUILDER_H
|
||||
#define QMCPLUSPLUS_ORIGINAL_JASTROW_AA_BUILDER_H
|
||||
#include "QMCWaveFunctions/OrbitalBuilderBase.h"
|
||||
#include "QMCWaveFunctions/Jastrow/TwoBodyJastrowOrbital.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
//forward declaration
|
||||
class ParticleSet;
|
||||
|
||||
/** Generic two-body Jastrow builder
|
||||
*
|
||||
* Replacement of JastrowBuilder::createTwoBodySpin and JastrowBuilder::createTwoBodyNoSpin
|
||||
*/
|
||||
struct JAABuilder: public OrbitalBuilderBase
|
||||
{
|
||||
|
||||
JAABuilder(ParticleSet& p, TrialWaveFunction& psi);
|
||||
|
||||
bool put(xmlNodePtr cur);
|
||||
|
||||
template <class FN> TwoBodyJastrowOrbital<FN>* createJAA(xmlNodePtr cur, const std::string& jname);
|
||||
|
||||
bool IgnoreSpin;
|
||||
|
||||
DistanceTableData* d_table;
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
|
@ -19,11 +19,7 @@
|
|||
#include "QMCWaveFunctions/Jastrow/PadeFunctors.h"
|
||||
#if QMC_BUILD_LEVEL>2
|
||||
#include "QMCWaveFunctions/Jastrow/BsplineJastrowBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/GaussianFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ModPadeFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/BesselZeroFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/BsplineFunctor.h"
|
||||
#include "QMCWaveFunctions/Jastrow/OpenGaussianSlaterFunctor.h"
|
||||
#endif
|
||||
#include "QMCWaveFunctions/Jastrow/OneBodyJastrowOrbital.h"
|
||||
#include "QMCWaveFunctions/Jastrow/DiffOneBodyJastrowOrbital.h"
|
||||
|
|
@ -121,48 +117,20 @@ bool JABBuilder::put(xmlNodePtr cur)
|
|||
{
|
||||
success = createJAB<PadeFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else
|
||||
if (jastfunction == "pade2")
|
||||
{
|
||||
success = createJAB<Pade2ndOrderFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "pade2")
|
||||
{
|
||||
success = createJAB<Pade2ndOrderFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
#if QMC_BUILD_LEVEL>2
|
||||
else
|
||||
if (jastfunction == "short")
|
||||
{
|
||||
success = createJAB<ModPadeFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else
|
||||
if (jastfunction == "gaussian")
|
||||
{
|
||||
success = createJAB<GaussianFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else
|
||||
if (jastfunction == "opengaussianslater")
|
||||
{
|
||||
app_log()<<jastfunction<< std::endl;
|
||||
success = createJAB<OpenGaussianSlaterFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else
|
||||
if (jastfunction == "besselzero")
|
||||
{
|
||||
success = createJAB<BesselZero<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else
|
||||
if (jastfunction == "shiftedgaussian")
|
||||
{
|
||||
success = createJAB<TruncatedShiftedGaussianFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else
|
||||
if (jastfunction == "bspline")
|
||||
{
|
||||
success = createJAB<BsplineFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
else if (jastfunction == "bspline")
|
||||
{
|
||||
success = createJAB<BsplineFunctor<RealType> >(cur,jastfunction);
|
||||
}
|
||||
#endif
|
||||
else
|
||||
{
|
||||
app_error() << "Unknown one body function: " << jastfunction << ".\n";
|
||||
}
|
||||
else
|
||||
{
|
||||
app_error() << "Unknown one body function: " << jastfunction << ".\n";
|
||||
}
|
||||
return success;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,128 +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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#include "QMCWaveFunctions/Jastrow/JastrowBasisBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/AtomicBasisBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/GTOBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/BsplineAOBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/CBSOBuilder.h"
|
||||
#include "QMCWaveFunctions/SparseLocalizedBasisSet.h"
|
||||
#include "Message/Communicate.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** constructor
|
||||
* \param els reference to the electrons
|
||||
* \param ions reference to the ions
|
||||
*/
|
||||
JastrowBasisBuilder::JastrowBasisBuilder(ParticleSet& els, ParticleSet& ions,
|
||||
const std::string& functype, bool usespline):
|
||||
targetPtcl(els), sourcePtcl(ions), UseSpline(usespline),FuncType(functype), myBasisSet(0)
|
||||
{
|
||||
}
|
||||
|
||||
JastrowBasisBuilder::~JastrowBasisBuilder()
|
||||
{
|
||||
}
|
||||
|
||||
template<typename RFBUILDER>
|
||||
void JastrowBasisBuilder::createLocalizedBasisSet(xmlNodePtr cur)
|
||||
{
|
||||
typedef typename RFBUILDER::CenteredOrbitalType COT;
|
||||
typedef SparseLocalizedBasisSet<COT> ThisBasisSetType;
|
||||
ThisBasisSetType* curBasis= new ThisBasisSetType(sourcePtcl,targetPtcl);
|
||||
//create the size vector with zeros
|
||||
SizeOfBasisPerCenter.resize(sourcePtcl.getSpeciesSet().getTotalNum(),0);
|
||||
//create the basis set
|
||||
cur = cur->xmlChildrenNode;
|
||||
while(cur!=NULL)
|
||||
{
|
||||
std::string cname((const char*)(cur->name));
|
||||
if(cname == "atomicBasisSet")
|
||||
{
|
||||
const xmlChar* eptr=xmlGetProp(cur,(const xmlChar*)"elementType");
|
||||
if(eptr == NULL)
|
||||
{
|
||||
app_error() << " Missing elementType attribute of atomicBasisSet.\n Abort at MOBasisBuilder::put " << std::endl;
|
||||
OHMMS::Controller->abort();
|
||||
}
|
||||
std::string elementType((const char*)eptr);
|
||||
std::map<std::string,BasisSetBuilder*>::iterator it = aoBuilders.find(elementType);
|
||||
if(it == aoBuilders.end())
|
||||
{
|
||||
AtomicBasisBuilder<RFBUILDER>* any = new AtomicBasisBuilder<RFBUILDER>(elementType);
|
||||
any->put(cur);
|
||||
COT* aoBasis= any->createAOSet(cur);
|
||||
if(aoBasis)
|
||||
//add the new atomic basis to the basis set
|
||||
{
|
||||
int activeCenter =sourcePtcl.getSpeciesSet().findSpecies(elementType);
|
||||
curBasis->add(activeCenter, aoBasis);
|
||||
aoBuilders[elementType]=any;
|
||||
printRadialFunctors(elementType,aoBasis);
|
||||
SizeOfBasisPerCenter[activeCenter]=aoBasis->getBasisSetSize();
|
||||
}
|
||||
}
|
||||
}
|
||||
cur = cur->next;
|
||||
}
|
||||
//resize the basis set
|
||||
curBasis->setBasisSetSize(-1);
|
||||
//clean up
|
||||
if(myBasisSet) delete myBasisSet;
|
||||
myBasisSet=curBasis;
|
||||
}
|
||||
|
||||
bool JastrowBasisBuilder::put(xmlNodePtr cur)
|
||||
{
|
||||
if(myBasisSet)
|
||||
return true;
|
||||
if(UseSpline)
|
||||
{
|
||||
createLocalizedBasisSet<CBSOBuilder>(cur);
|
||||
}
|
||||
else
|
||||
{
|
||||
if(FuncType == "Bspline")
|
||||
{
|
||||
createLocalizedBasisSet<BsplineAOBuilder>(cur);
|
||||
}
|
||||
if(FuncType == "gto" || FuncType == "GTO")
|
||||
createLocalizedBasisSet<GTOBuilder>(cur);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
template<typename COT>
|
||||
void JastrowBasisBuilder::printRadialFunctors(const std::string& elementType, COT* aoBasis)
|
||||
{
|
||||
#if !defined(HAVE_MPI)
|
||||
std::string fname(elementType);
|
||||
fname.append(".j3.dat");
|
||||
std::ofstream fout(fname.c_str());
|
||||
int nr=aoBasis->Rnl.size();
|
||||
fout << "# number of radial functors = " << nr << std::endl;
|
||||
double r=0.0;
|
||||
while(r<20)
|
||||
{
|
||||
fout << r ;
|
||||
for(int i=0; i<nr; i++)
|
||||
fout << " " << aoBasis->Rnl[i]->evaluate(r,1.0/r);
|
||||
fout << std::endl;
|
||||
r += 0.013;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
|
@ -1,81 +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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
|
||||
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_JASTROWBASISBUILDER_H
|
||||
#define QMCPLUSPLUS_JASTROWBASISBUILDER_H
|
||||
|
||||
#include "QMCWaveFunctions/BasisSetBase.h"
|
||||
//#include "QMCWaveFunctions/Jastrow/CBSOBuilder.h"
|
||||
//#include "QMCWaveFunctions/LocalizedBasisSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
|
||||
/** derived class from BasisSetBuilder
|
||||
*
|
||||
* Create a basis set of molecular orbital types as defined in MolecularOrbitalBasis
|
||||
* with radial wave functions on the radial grids.
|
||||
*/
|
||||
class JastrowBasisBuilder: public BasisSetBuilder
|
||||
{
|
||||
public:
|
||||
///keep the current basis set
|
||||
BasisSetBase<RealType>* myBasisSet;
|
||||
/** constructor
|
||||
* \param els reference to the electrons
|
||||
* \param ions reference to the ions
|
||||
*/
|
||||
JastrowBasisBuilder(ParticleSet& els, ParticleSet& ions, const std::string& functype, bool usespline=false);
|
||||
///destructor
|
||||
~JastrowBasisBuilder();
|
||||
|
||||
bool put(xmlNodePtr cur);
|
||||
|
||||
/** do not create anything. */
|
||||
SPOSetBase* createSPOSetFromXML(xmlNodePtr cur)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
///size of blocks
|
||||
std::vector<int> SizeOfBasisPerCenter;
|
||||
private:
|
||||
///target ParticleSet
|
||||
ParticleSet& targetPtcl;
|
||||
///source ParticleSet
|
||||
ParticleSet& sourcePtcl;
|
||||
///boolean to choose numerical or analytic form
|
||||
bool UseSpline;
|
||||
///function name
|
||||
std::string FuncType;
|
||||
///save AtomiBasisBuilder<RFB>*
|
||||
std::map<std::string,BasisSetBuilder*> aoBuilders;
|
||||
/** create a localized basis set
|
||||
*
|
||||
* The template parameter RFBUILDER is a builder class for radial
|
||||
* functors.
|
||||
*/
|
||||
template<typename RFBUILDER>
|
||||
void createLocalizedBasisSet(xmlNodePtr cur);
|
||||
|
||||
/** print the basis set
|
||||
* @param elementType element name
|
||||
* @param aoBasis atomic orbtial basis associated with elementType
|
||||
*/
|
||||
template<typename COT>
|
||||
void printRadialFunctors(const std::string& elementType, COT* aoBasis);
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -14,7 +14,6 @@
|
|||
|
||||
|
||||
#include "QMCWaveFunctions/Jastrow/JastrowBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/JastrowBasisBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/PadeJastrowBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/RPAJastrow.h"
|
||||
#include "QMCWaveFunctions/Jastrow/BsplineJastrowBuilder.h"
|
||||
|
|
@ -22,11 +21,8 @@
|
|||
#include "QMCWaveFunctions/Jastrow/singleRPAJastrowBuilder.h"
|
||||
#if OHMMS_DIM ==3
|
||||
#include "QMCWaveFunctions/Jastrow/eeI_JastrowBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ThreeBodyGeminal.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ThreeBodyBlockSparse.h"
|
||||
#endif
|
||||
#include "QMCWaveFunctions/Jastrow/JABBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/JAABuilder.h"
|
||||
#include "Utilities/ProgressReportEngine.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
|
|
@ -72,8 +68,6 @@ bool JastrowBuilder::put(xmlNodePtr cur)
|
|||
return addOneBody(cur);
|
||||
if(typeOpt.find("Two") < typeOpt.size())
|
||||
return addTwoBody(cur);
|
||||
if(typeOpt.find("Three") < typeOpt.size())
|
||||
return addThreeBody(cur);
|
||||
if(typeOpt.find("eeI") < typeOpt.size())
|
||||
return add_eeI(cur);
|
||||
if(typeOpt.find("kSpace") < typeOpt.size())
|
||||
|
|
@ -173,7 +167,6 @@ bool JastrowBuilder::addTwoBody(xmlNodePtr cur)
|
|||
bool ignoreSpin = (spinOpt == "no");
|
||||
//convert to lowercase
|
||||
tolower(funcOpt);
|
||||
// OrbitalConstraintsBase* control=0;
|
||||
if(funcOpt == "pade")
|
||||
{
|
||||
if (targetPtcl.Lattice.SuperCellEnum != SUPERCELL_OPEN)
|
||||
|
|
@ -184,176 +177,31 @@ bool JastrowBuilder::addTwoBody(xmlNodePtr cur)
|
|||
PadeJastrowBuilder pbuilder(targetPtcl,targetPsi,ptclPool);
|
||||
return pbuilder.put(cur);
|
||||
}
|
||||
else
|
||||
if((funcOpt == "yukawa") || (funcOpt == "rpa"))
|
||||
else if((funcOpt == "yukawa") || (funcOpt == "rpa"))
|
||||
{
|
||||
if(targetPtcl.Lattice.SuperCellEnum == SUPERCELL_OPEN)
|
||||
{
|
||||
if(targetPtcl.Lattice.SuperCellEnum == SUPERCELL_OPEN)
|
||||
{
|
||||
PRE.warning("RPA is requested for an open system. Please choose other functors.");
|
||||
return false;
|
||||
}
|
||||
// control = new RPAPBCConstraints(targetPtcl,targetPsi,ignoreSpin);
|
||||
RPAJastrow* rpajastrow = new RPAJastrow(targetPtcl,targetPsi.is_manager());
|
||||
rpajastrow->put(cur);
|
||||
targetPsi.addOrbital(rpajastrow,nameOpt);
|
||||
return true;
|
||||
PRE.warning("RPA is requested for an open system. Please choose other functors.");
|
||||
return false;
|
||||
}
|
||||
else
|
||||
if (funcOpt == "bspline" )
|
||||
{
|
||||
BsplineJastrowBuilder bbuilder(targetPtcl,targetPsi);
|
||||
return bbuilder.put(cur);
|
||||
}
|
||||
#if QMC_BUILD_LEVEL>2
|
||||
else //try other special analytic functions
|
||||
{
|
||||
app_log() << "\n Using JAABuilder for two-body jastrow with analytic functions" << std::endl;
|
||||
JAABuilder jb(targetPtcl,targetPsi);
|
||||
return jb.put(cur);
|
||||
}
|
||||
#endif
|
||||
// control->setReportLevel(ReportLevel);
|
||||
// success=control->put(cur);
|
||||
// if(!success) {
|
||||
// PRE.error("Failed to build "+nameOpt+" Jastrow function");
|
||||
// delete control;
|
||||
// return false;
|
||||
// }
|
||||
//
|
||||
// OrbitalBase* j2=control->createTwoBody();
|
||||
// if(j2== 0)
|
||||
// {
|
||||
// APP_ABORT("JastrowBuild::addTwoBody");
|
||||
// delete control;
|
||||
// return false;
|
||||
// }
|
||||
//
|
||||
// enum {MULTIPLE=0, LONGRANGE, ONEBODY, TWOBODY, THREEBODY, FOURBODY};
|
||||
//
|
||||
// if(control->JComponent[MULTIPLE])
|
||||
// {
|
||||
// //create a combo orbital
|
||||
// ComboOrbital* jcombo=new ComboOrbital(control);
|
||||
// jcombo->Psi.push_back(j2);
|
||||
// if(control->JComponent[ONEBODY] && sourceOpt != targetPtcl.getName())
|
||||
// {
|
||||
// app_log() << ("Adding one-body Jastrow function dependent upon two-body " + funcOpt + " Jastrow function\n");
|
||||
// std::map<std::string,ParticleSet*>::iterator pa_it(ptclPool.find(sourceOpt));
|
||||
// if(pa_it != ptclPool.end())
|
||||
// {
|
||||
// OrbitalBase* j1=control->createOneBody(*((*pa_it).second));
|
||||
// if(j1) jcombo->Psi.push_back(j1);
|
||||
// }
|
||||
// }
|
||||
// if(control->JComponent[LONGRANGE])
|
||||
// {
|
||||
// app_log() << ("Adding long-range component of " + funcOpt + " Jastrow function\n");
|
||||
// control->addExtra2ComboOrbital(jcombo);
|
||||
// }
|
||||
// targetPsi.addOrbital(jcombo,"J2_LR");
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// std::string j2name="J2_"+funcOpt;
|
||||
// targetPsi.addOrbital(j2,j2name);
|
||||
// }
|
||||
// Children.push_back(control);
|
||||
RPAJastrow* rpajastrow = new RPAJastrow(targetPtcl,targetPsi.is_manager());
|
||||
rpajastrow->put(cur);
|
||||
targetPsi.addOrbital(rpajastrow,nameOpt);
|
||||
return true;
|
||||
}
|
||||
else if (funcOpt == "bspline" )
|
||||
{
|
||||
BsplineJastrowBuilder bbuilder(targetPtcl,targetPsi);
|
||||
return bbuilder.put(cur);
|
||||
}
|
||||
else
|
||||
{
|
||||
app_error() << "Unknown two body function: " << funcOpt << ".\n";
|
||||
}
|
||||
|
||||
return success;
|
||||
}
|
||||
|
||||
bool JastrowBuilder::addThreeBody(xmlNodePtr cur)
|
||||
{
|
||||
#if OHMMS_DIM==3
|
||||
if(sourceOpt == targetPtcl.getName())
|
||||
{
|
||||
app_warning() << " Three-Body Jastrow Function needs a source different from " << targetPtcl.getName() << std::endl;
|
||||
APP_ABORT(" JastrowBuilder::addThreeBody");
|
||||
return false;
|
||||
}
|
||||
PtclPoolType::iterator pit(ptclPool.find(sourceOpt));
|
||||
if(pit == ptclPool.end())
|
||||
{
|
||||
app_error() << " JastrowBuilder::addThreeBody requires a center. " << sourceOpt << " is invalid " << std::endl;
|
||||
APP_ABORT(" JastrowBuilder::addThreeBody");
|
||||
return false;
|
||||
}
|
||||
app_log() << " JastrowBuilder::addThreeBody source="<< sourceOpt << std::endl;
|
||||
xmlNodePtr basisPtr=NULL;
|
||||
xmlNodePtr coeffPtr=NULL;
|
||||
std::string diagOnly("no");//default:use diagonal blocks only
|
||||
std::string sameblocks("yes");
|
||||
OhmmsAttributeSet tAttrib;
|
||||
tAttrib.add(diagOnly,"diagonal");
|
||||
tAttrib.add(sameblocks,"sameBlocksForGroup");
|
||||
tAttrib.put(cur);
|
||||
cur = cur->xmlChildrenNode;
|
||||
while(cur != NULL)
|
||||
{
|
||||
std::string cname((const char*)(cur->name));
|
||||
if(cname == basisset_tag)
|
||||
{
|
||||
basisPtr=cur;
|
||||
}
|
||||
else
|
||||
if(cname.find("oeff")<cname.size())
|
||||
{
|
||||
coeffPtr=cur;
|
||||
OhmmsAttributeSet oAttrib;
|
||||
oAttrib.add(diagOnly,"diagonal");
|
||||
oAttrib.put(cur);
|
||||
}
|
||||
cur=cur->next;
|
||||
}
|
||||
if(basisPtr == NULL)
|
||||
{
|
||||
app_error() << " JastrowBuilder::addThreeBody exit. Missing <basisset/>"<< std::endl;
|
||||
return false;
|
||||
}
|
||||
ParticleSet* sourcePtcl=(*pit).second;
|
||||
JastrowBasisBuilder* basisBuilder =
|
||||
new JastrowBasisBuilder(targetPtcl,*sourcePtcl,funcOpt,transformOpt == "yes");
|
||||
basisBuilder->put(basisPtr);
|
||||
if(diagOnly == "yes")
|
||||
{
|
||||
app_log() << "\n creating Three-Body Jastrow function using only diagnoal blocks." << std::endl;
|
||||
ThreeBodyBlockSparse* J3 = new ThreeBodyBlockSparse(*sourcePtcl, targetPtcl);
|
||||
J3->setBasisSet(basisBuilder->myBasisSet);
|
||||
J3->put(coeffPtr);
|
||||
J3->SameBlocksForGroup=(sameblocks == "yes");
|
||||
J3->setBlocks(basisBuilder->SizeOfBasisPerCenter);
|
||||
targetPsi.addOrbital(J3,"J3_block");
|
||||
}
|
||||
else
|
||||
{
|
||||
app_log() << "\n creating Three-Body Jastrow function using a complete Geminal matrix." << std::endl;
|
||||
ThreeBodyGeminal* J3 = new ThreeBodyGeminal(*sourcePtcl, targetPtcl);
|
||||
J3->setBasisSet(basisBuilder->myBasisSet);
|
||||
J3->put(coeffPtr);
|
||||
targetPsi.addOrbital(J3,"J3_full");
|
||||
}
|
||||
delete basisBuilder;
|
||||
#else
|
||||
app_error() << " Three-body Jastrow function is not supported for DIM != 3." << std::endl;
|
||||
//#error " Three-body Jastrow is disabled for QMC_DIM != 3\n "
|
||||
#endif
|
||||
//if(jbuilder)
|
||||
//{
|
||||
// jbuilder->put(cur);
|
||||
// Children.push_back(jbuilder);
|
||||
// return true;
|
||||
//}
|
||||
// } else if (jasttype == "Three-Body-Pade") {
|
||||
// app_log() << "\n creating Three-Body-Pade Jastrow function " << std::endl;
|
||||
// std::string source_name("i");
|
||||
// const xmlChar* iptr = xmlGetProp(cur, (const xmlChar *)"source");
|
||||
// //if(iptr != NULL) source_name=(const char*)iptr;
|
||||
// PtclPoolType::iterator pit(ptclPool.find(source_name));
|
||||
// if(pit != ptclPool.end()) {
|
||||
// jbuilder = new ThreeBodyPadeBuilder(*targetPtcl,*targetPsi,*((*pit).second));
|
||||
// }
|
||||
// }
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -56,8 +56,6 @@ private:
|
|||
bool addOneBody(xmlNodePtr cur);
|
||||
///add two-body term
|
||||
bool addTwoBody(xmlNodePtr cur);
|
||||
///add three-body term
|
||||
bool addThreeBody(xmlNodePtr cur);
|
||||
/// add electron-electron ion term
|
||||
bool add_eeI (xmlNodePtr cur);
|
||||
///add k-Space term
|
||||
|
|
|
|||
|
|
@ -1,752 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_MCMILLANJ2_H
|
||||
#define QMCPLUSPLUS_MCMILLANJ2_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
/** ModPade Jastrow functional
|
||||
*
|
||||
* \f[ u(r) = \frac{1}{2*A}\left[1-\exp(-A*r)\right], \f]
|
||||
*/
|
||||
template<class T>
|
||||
struct McMillanJ2Functor: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
|
||||
///coefficients
|
||||
real_type A;
|
||||
real_type B;
|
||||
|
||||
std::string ID_A;
|
||||
std::string ID_B;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
McMillanJ2Functor(real_type a=5.0, real_type b=4.9133):ID_A("McA"),ID_B("McB")
|
||||
{
|
||||
reset(a,b);
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new McMillanJ2Functor<T>(*this);
|
||||
}
|
||||
inline void reset()
|
||||
{
|
||||
reset(A,B);
|
||||
}
|
||||
|
||||
/** reset the internal variables.
|
||||
*
|
||||
* USE_resetParameters
|
||||
*/
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ia=myVars.where(0);
|
||||
if(ia>-1)
|
||||
A=active[ia];
|
||||
int ib=myVars.where(1);
|
||||
if(ib>-1)
|
||||
B=active[ib];
|
||||
reset(A,B);
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
|
||||
/** reset the internal variables.
|
||||
*@param a New Jastrow parameter a
|
||||
*/
|
||||
inline void reset(real_type a, real_type b)
|
||||
{
|
||||
A = a;
|
||||
B = b;
|
||||
}
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
return std::pow(B/r,A);
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
real_type wert = std::pow(B/r,A);
|
||||
dudr = -A*wert/r;
|
||||
d2udr2= (A+1.)*A*wert/(r*r);
|
||||
return wert;
|
||||
}
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
return -A*std::pow(B/r,A)/r;
|
||||
}
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if(cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
// std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if(aname == "a")
|
||||
{
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if(aname == "b")
|
||||
{
|
||||
ID_B = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(B,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_B,B);
|
||||
reset(A,B);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
// template<class T>
|
||||
// struct ModMcMillanJ2Functor: public OptimizableFunctorBase
|
||||
// {
|
||||
//
|
||||
// typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
//
|
||||
// ///coefficients
|
||||
// real_type A;
|
||||
// real_type B;
|
||||
// real_type RC;
|
||||
// real_type c0,c1,c2,c3,c4,c5,c6;
|
||||
//
|
||||
// std::string ID_A,ID_B,ID_RC;
|
||||
//
|
||||
// /** constructor
|
||||
// * @param a A coefficient
|
||||
// * @param samespin boolean to indicate if this function is for parallel spins
|
||||
// */
|
||||
// ModMcMillanJ2Functor(real_type a=4.9133, real_type b=5, real_type rc=1.0):ID_RC("Rcutoff"),ID_A("A"),ID_B("B")
|
||||
// {
|
||||
// A = a;
|
||||
// B = b;
|
||||
// RC = rc;
|
||||
// reset();
|
||||
// }
|
||||
//
|
||||
//
|
||||
// void checkInVariables(opt_variables_type& active)
|
||||
// {
|
||||
// active.insertFrom(myVars);
|
||||
// }
|
||||
//
|
||||
// void checkOutVariables(const opt_variables_type& active)
|
||||
// {
|
||||
// myVars.getIndex(active);
|
||||
// }
|
||||
//
|
||||
// OptimizableFunctorBase* makeClone() const
|
||||
// {
|
||||
// return new ModMcMillanJ2Functor<T>(*this);
|
||||
// }
|
||||
//
|
||||
// void resetParameters(const opt_variables_type& active)
|
||||
// {
|
||||
// int ia=myVars.where(0); if(ia>-1) A=active[ia];
|
||||
// int ib=myVars.where(1); if(ib>-1) B=active[ib];
|
||||
// int ic=myVars.where(2); if(ic>-1) RC=active[ic];
|
||||
// reset();
|
||||
// }
|
||||
//
|
||||
// inline void reset() {
|
||||
// c0 = std::pow(A,B);
|
||||
// c1 = -1.0*std::pow(A/RC,B);
|
||||
// c2 = B*c0*std::pow(RC,-B-1);
|
||||
// c3 = -0.5*(B+1)*B*c0*std::pow(RC,-B-2);
|
||||
//
|
||||
// c4 = -B*c0;
|
||||
// c5 = 2.0*c3;
|
||||
// c6 = B*(B+1.0)*c0;
|
||||
// }
|
||||
//
|
||||
//
|
||||
// /** evaluate the value at r
|
||||
// * @param r the distance
|
||||
// * @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
// */
|
||||
// inline real_type evaluate(real_type r) {
|
||||
// if (r<RC) {
|
||||
// real_type r1 = (r-RC);
|
||||
// return c0*std::pow(r,-B)+c1 + r1*(c2+r1*c3) ;
|
||||
// }
|
||||
// else {return 0.0;};
|
||||
// }
|
||||
//
|
||||
// /**@param r the distance
|
||||
// @param dudr first derivative
|
||||
// @param d2udr2 second derivative
|
||||
// @return the value
|
||||
// */
|
||||
// inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) {
|
||||
// if (r<RC){
|
||||
// real_type r1 = (r-RC);
|
||||
// real_type overr = 1.0/r;
|
||||
// real_type wert = std::pow(1.0/r,B);
|
||||
// dudr = c4*wert*overr + c2 + c5*r1;
|
||||
// d2udr2= c6*wert*overr*overr+ c5 ;
|
||||
// return c0*wert + c1 + r1*(c2+r1*c3) ;
|
||||
// } else {
|
||||
// dudr = 0.0;
|
||||
// d2udr2= 0.0;
|
||||
// return 0.0;
|
||||
// };
|
||||
// }
|
||||
//
|
||||
// /** return a value at r
|
||||
// */
|
||||
// real_type f(real_type r) {
|
||||
// return evaluate(r);
|
||||
// }
|
||||
//
|
||||
// /** return a derivative at r
|
||||
// */
|
||||
// real_type df(real_type r) {
|
||||
// if (r<RC) {return c4*std::pow(1.0/r,B+1) + c2 + c5*(r-RC) ;}
|
||||
// else {return 0.0;};
|
||||
// }
|
||||
//
|
||||
// /** Read in the parameter from the xml input file.
|
||||
// * @param cur current xmlNode from which the data members are reset
|
||||
// */
|
||||
// bool put(xmlNodePtr cur) {
|
||||
// RC = cutoff_radius;
|
||||
// xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
// while(tcur != NULL) {
|
||||
// //@todo Var -> <param(eter) role="opt"/>
|
||||
// std::string cname((const char*)(tcur->name));
|
||||
// if(cname == "parameter" || cname == "Var") {
|
||||
// std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
// std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
// if(aname == "A") {
|
||||
// ID_A = idname;
|
||||
// putContent(A,tcur);
|
||||
// } else if(aname == "B") {
|
||||
// ID_B = idname;
|
||||
// putContent(B,tcur);
|
||||
// }else if(aname == "RC")
|
||||
// {
|
||||
// ID_RC = idname;
|
||||
// putContent(RC,tcur);
|
||||
// }
|
||||
// }
|
||||
// tcur = tcur->next;
|
||||
// }
|
||||
// myVars.insert(ID_A,A);
|
||||
// myVars.insert(ID_B,B);
|
||||
// myVars.insert(ID_RC,RC);
|
||||
// reset();
|
||||
// LOGMSG(" Modified McMillan Parameters ")
|
||||
// LOGMSG(" A (" << ID_A << ") = " << A << " B (" << ID_B << ") = " << B<< " R_c (" << ID_RC << ") = " << RC)
|
||||
// return true;
|
||||
// }
|
||||
//
|
||||
// };
|
||||
|
||||
template<class T>
|
||||
struct ModMcMillanJ2Functor: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
///f(r) = A*r^-5 - A*r_c^-5 + A/(B*r_c^6) * tanh(B*(r-r_c))
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
|
||||
///coefficients
|
||||
real_type A;
|
||||
real_type B;
|
||||
real_type RC;
|
||||
real_type cA,c0,c1,c2,c3,c4,c5;
|
||||
|
||||
std::string ID_A,ID_B,ID_RC;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
ModMcMillanJ2Functor(real_type a=0.0, real_type b=0.0, real_type rc=0.0):ID_RC("Rcutoff"),ID_A("A"),ID_B("B")
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new ModMcMillanJ2Functor<T>(*this);
|
||||
}
|
||||
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ia=myVars.where(0);
|
||||
if(ia>-1)
|
||||
A=active[ia];
|
||||
int ib=myVars.where(1);
|
||||
if(ib>-1)
|
||||
B=active[ib];
|
||||
int ic=myVars.where(2);
|
||||
if(ic>-1)
|
||||
RC=active[ic];
|
||||
reset();
|
||||
}
|
||||
|
||||
inline void reset()
|
||||
{
|
||||
cA = std::pow(A,5);
|
||||
c0 = cA*std::pow(RC,-5);
|
||||
c1 = 5*cA*std::pow(RC,-6)/B;
|
||||
c2 = -5*cA;
|
||||
c3 = B*c1;
|
||||
c4 = 30*cA;
|
||||
c5 = -2*B*c3;
|
||||
}
|
||||
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
real_type r1 = (r-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
return cA*std::pow(r,-5)-c0 + c1*tanhZ;
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
real_type prt = cA*std::pow(r,-5);
|
||||
real_type r1 = (r-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
real_type sechZ = 2.0/(plZ+miZ);
|
||||
real_type sechZ2 = sechZ*sechZ;
|
||||
real_type ret = prt - c0 + c1*tanhZ;
|
||||
prt = prt/r;
|
||||
dudr = -5*prt + c3*sechZ2;
|
||||
prt = prt/r;
|
||||
d2udr2 = 30*prt + c5*sechZ2*tanhZ;
|
||||
return ret;
|
||||
}
|
||||
else
|
||||
{
|
||||
dudr = 0.0;
|
||||
d2udr2= 0.0;
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
if (r<RC)
|
||||
{
|
||||
real_type r1 = (r-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type sechZ = 2.0/(plZ+miZ);
|
||||
real_type sechZ2 = sechZ*sechZ;
|
||||
return c2*std::pow(r,-6) + c3*sechZ2;
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
RC = cutoff_radius;
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if(cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if(aname == "A")
|
||||
{
|
||||
ID_A = idname;
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if(aname == "B")
|
||||
{
|
||||
ID_B = idname;
|
||||
putContent(B,tcur);
|
||||
}
|
||||
else
|
||||
if(aname == "RC")
|
||||
{
|
||||
ID_RC = idname;
|
||||
putContent(RC,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_B,B);
|
||||
myVars.insert(ID_RC,RC);
|
||||
reset();
|
||||
LOGMSG(" Modified McMillan Parameters ")
|
||||
LOGMSG(" A (" << ID_A << ") = " << A << " B (" << ID_B << ") = " << B<< " R_c (" << ID_RC << ") = " << RC)
|
||||
return true;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
template<class T>
|
||||
struct comboMcMillanJ2Functor: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
///f(r) = A*r^-5 - A*r_c^-5 + A/(B*r_c^6) * tanh(B*(r-r_c))
|
||||
///made to be used with another jastrow within a cutoff distance R0 due to divergent cusp.
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
|
||||
///coefficients
|
||||
real_type A;
|
||||
real_type B;
|
||||
real_type RC;
|
||||
real_type R0;
|
||||
real_type cA,c0,c1,c2,c3,c4,c5;
|
||||
real_type b0,b1,b2,b3;
|
||||
|
||||
std::string ID_A,ID_B,ID_RC,ID_R0;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
comboMcMillanJ2Functor(real_type a=0.0, real_type b=0.0, real_type r0=0.0 , real_type rc=0.0):ID_RC("Rcutoff"),ID_R0("Rinner"),ID_A("A"),ID_B("B")
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new comboMcMillanJ2Functor<T>(*this);
|
||||
}
|
||||
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ia=myVars.where(0);
|
||||
if(ia>-1)
|
||||
A=active[ia];
|
||||
int ib=myVars.where(1);
|
||||
if(ib>-1)
|
||||
B=active[ib];
|
||||
int ic=myVars.where(2);
|
||||
if(ic>-1)
|
||||
RC=active[ic];
|
||||
int id=myVars.where(3);
|
||||
if(id>-1)
|
||||
R0=active[id];
|
||||
reset();
|
||||
}
|
||||
|
||||
inline void reset()
|
||||
{
|
||||
cA = std::pow(A,5);
|
||||
c0 = cA*std::pow(RC,-5);
|
||||
c1 = 5*cA*std::pow(RC,-6)/B;
|
||||
c2 = -5*cA;
|
||||
c3 = B*c1;
|
||||
c4 = 30*cA;
|
||||
c5 = -2*B*c3;
|
||||
real_type r1 = (R0-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
real_type sechZ = 2.0/(plZ+miZ);
|
||||
real_type sechZ2 = sechZ*sechZ;
|
||||
real_type btemp = cA/std::pow(R0,-7);
|
||||
b3 = 30.0*btemp + c5*sechZ2*tanhZ;
|
||||
b2 = 0.5*b3;
|
||||
b1 = -b3*R0 - 5.0*btemp*R0 + c3*sechZ2;
|
||||
b0 = -b2*R0*R0 -b1*R0 + btemp*R0*R0 - c0 +c1*tanhZ;
|
||||
}
|
||||
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
if ((r<RC)&&(r>R0))
|
||||
{
|
||||
real_type r1 = (r-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
return cA*std::pow(r,-5)-c0 + c1*tanhZ;
|
||||
}
|
||||
else
|
||||
if ((r<RC)&&(r<=R0))
|
||||
{
|
||||
// real_type r1 = (r-RC);
|
||||
// real_type plZ = std::exp(B*r1);
|
||||
// real_type miZ = 1.0/plZ;
|
||||
// real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
// real_type sechZ = 2.0/(plZ+miZ);
|
||||
// real_type sechZ2 = sechZ*sechZ;
|
||||
// real_type ret = b0 + b1*r + b2*r*r - c0 + c1*tanhZ;
|
||||
real_type ret = b0 + b1*r + b2*r*r;
|
||||
return ret;
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
if ((r<RC)&&(r>R0))
|
||||
{
|
||||
real_type prt = cA*std::pow(r,-5);
|
||||
real_type r1 = (r-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
real_type sechZ = 2.0/(plZ+miZ);
|
||||
real_type sechZ2 = sechZ*sechZ;
|
||||
real_type ret = prt - c0 + c1*tanhZ;
|
||||
prt = prt/r;
|
||||
dudr = -5*prt + c3*sechZ2;
|
||||
prt = prt/r;
|
||||
d2udr2 = 30*prt + c5*sechZ2*tanhZ;
|
||||
return ret;
|
||||
}
|
||||
else
|
||||
if ((r<RC)&&(r<=R0))
|
||||
{
|
||||
// real_type r1 = (r-RC);
|
||||
// real_type plZ = std::exp(B*r1);
|
||||
// real_type miZ = 1.0/plZ;
|
||||
// real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
// real_type sechZ = 2.0/(plZ+miZ);
|
||||
// real_type sechZ2 = sechZ*sechZ;
|
||||
// real_type ret = b0+b1*r+b2*r*r - c0 + c1*tanhZ;
|
||||
// dudr = b1 + b3*r + c3*sechZ2;
|
||||
// d2udr2 = b3 + c5*sechZ2*tanhZ;
|
||||
real_type ret = b0 + b1*r + b2*r*r ;
|
||||
dudr = b1 + b3*r ;
|
||||
d2udr2 = b3 ;
|
||||
return ret;
|
||||
}
|
||||
else
|
||||
{
|
||||
dudr = 0.0;
|
||||
d2udr2= 0.0;
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
if ((r<RC)&&(r>R0))
|
||||
{
|
||||
real_type r1 = (r-RC);
|
||||
real_type plZ = std::exp(B*r1);
|
||||
real_type miZ = 1.0/plZ;
|
||||
real_type sechZ = 2.0/(plZ+miZ);
|
||||
real_type sechZ2 = sechZ*sechZ;
|
||||
return c2*std::pow(r,-6) + c3*sechZ2;
|
||||
}
|
||||
else
|
||||
if ((r<RC)&&(r<R0))
|
||||
{
|
||||
// real_type r1 = (r-RC);
|
||||
// real_type plZ = std::exp(B*r1);
|
||||
// real_type miZ = 1.0/plZ;
|
||||
// real_type tanhZ = (plZ-miZ)/(plZ+miZ);
|
||||
// real_type sechZ = 2.0/(plZ+miZ);
|
||||
// real_type sechZ2 = sechZ*sechZ;
|
||||
// real_type ret = b1 + b3*r + c3*sechZ2;
|
||||
real_type ret = b1 + b3*r ;
|
||||
return ret;
|
||||
}
|
||||
else
|
||||
{
|
||||
return 0.0;
|
||||
};
|
||||
}
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
RC = cutoff_radius;
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if(cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if(aname == "A")
|
||||
{
|
||||
ID_A = idname;
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if(aname == "B")
|
||||
{
|
||||
ID_B = idname;
|
||||
putContent(B,tcur);
|
||||
}
|
||||
else
|
||||
if(aname == "RC")
|
||||
{
|
||||
ID_RC = idname;
|
||||
putContent(RC,tcur);
|
||||
if (RC==0.0)
|
||||
RC=cutoff_radius;
|
||||
}
|
||||
else
|
||||
if(aname == "R0")
|
||||
{
|
||||
ID_R0 = idname;
|
||||
putContent(R0,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_B,B);
|
||||
myVars.insert(ID_RC,RC);
|
||||
myVars.insert(ID_R0,R0);
|
||||
reset();
|
||||
LOGMSG(" Combo McMillan Parameters ")
|
||||
LOGMSG(" A (" << ID_A << ") = " << A << " B (" << ID_B << ") = " << B<< " R_c (" << ID_RC << ") = "<< RC <<" R_0 (" << ID_R0 << ") = " << R0)
|
||||
return true;
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,211 +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: D.C. Yang, 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: D.C. Yang, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_MCMILLANJ2G_H
|
||||
#define QMCPLUSPLUS_MCMILLANJ2G_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
/**
|
||||
* McMillan Jastrow functor with cusp condition that matches HFD-BHE
|
||||
* N = 64, density equilibrium
|
||||
* Will be tuned up for general N, rho values
|
||||
*/
|
||||
template<class T>
|
||||
struct McMillanJ2GFunctor: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
|
||||
///coefficients
|
||||
real_type A, B, R0, c1, c2;
|
||||
// setting R0 = 2.5 by default, in the constructor
|
||||
|
||||
std::string ID_A;
|
||||
std::string ID_B;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
McMillanJ2GFunctor(real_type a=5.0, real_type b=5.7448):ID_A("McA"),ID_B("McB"),R0(2.5)
|
||||
{
|
||||
reset(a,b);
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new McMillanJ2GFunctor<T>(*this);
|
||||
}
|
||||
inline void reset()
|
||||
{
|
||||
reset(A,B);
|
||||
real_type Y;
|
||||
c1 = evaluate(R0, Y, c2);
|
||||
c2 = -Y/(2.0*R0*c1);
|
||||
c1 *= std::exp(-c2*R0*R0);
|
||||
}
|
||||
|
||||
/** reset the internal variables.
|
||||
*
|
||||
* USE_resetParameters
|
||||
*/
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ia=myVars.where(0);
|
||||
if(ia>-1)
|
||||
A=active[ia];
|
||||
int ib=myVars.where(1);
|
||||
if(ib>-1)
|
||||
B=active[ib];
|
||||
reset(A,B);
|
||||
real_type Y;
|
||||
c1 = evaluate(R0, Y, c2);
|
||||
c2 = -Y/(2.0*R0*c1);
|
||||
c1 *= std::exp(-c2*R0*R0);
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
|
||||
/** reset the internal variables.
|
||||
*@param a New Jastrow parameter a
|
||||
*/
|
||||
inline void reset(real_type a, real_type b)
|
||||
{
|
||||
A = a;
|
||||
B = b;
|
||||
}
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
// 64 particles, equil density for now
|
||||
return ( r < 2.5 ? c1*std::exp(-c2*r*r) :
|
||||
( r < cutoff_radius ? std::pow(B,A)*(std::pow(r,-A) + std::pow(2.0*cutoff_radius-r,-A) - 2.0*std::pow(cutoff_radius,-A)) :
|
||||
0.0 ) );
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
real_type wert;
|
||||
if (r < 2.5)
|
||||
{
|
||||
wert = c1*std::exp(-c2*r*r);
|
||||
dudr = -2.0*c1*c2*r*std::exp(-c2*r*r);
|
||||
d2udr2 = -2.0*c1*c2*std::exp(-c2*r*r) + 4.0*c1*c2*c2*r*r*std::exp(-c2*r*r);
|
||||
/*
|
||||
wert = 17318.0912945689*std::exp(-1.00041061270143*r*r);
|
||||
dudr = -34650.4046456380*r*std::exp(-1.00041061270143*r*r);
|
||||
d2udr2 = -34650.4046456380*std::exp(-1.00041061270143*r*r) + 69329.2650837902*r*r*std::exp(-1.00041061270143*r*r);
|
||||
*/
|
||||
}
|
||||
else
|
||||
if (r < cutoff_radius)
|
||||
{
|
||||
real_type rd = 2.0*cutoff_radius - r;
|
||||
wert = std::pow(B,A)*(std::pow(r,-A) + std::pow(rd,-A) - 2.0*std::pow(cutoff_radius,-A));
|
||||
dudr = A*std::pow(B,A)*(std::pow(rd,-(A+1.0))-std::pow(r,-(A+1.0)));
|
||||
d2udr2 = (A+1.0)*A*std::pow(B,A)*(std::pow(rd,-(A+2.0))+std::pow(r,-(A+2.0)));
|
||||
}
|
||||
else
|
||||
{
|
||||
wert = 0.0;
|
||||
dudr = 0.0;
|
||||
d2udr2 = 0.0;
|
||||
}
|
||||
return wert;
|
||||
/*
|
||||
real_type rinv = 1.0/r, rdinv = 1.0/(2.0*cutoff_radius - r), wert1 = std::pow(B*rinv,A), wert2 = std::pow(B*rdinv,A);
|
||||
dudr = -A*(wert1*rinv - wert2*rdinv);
|
||||
d2udr2= (A+1.)*A*(wert1*rinv*rinv + wert2*rdinv*rdinv);
|
||||
return (wert1 + wert2 - 2.0*std::pow(B/cutoff_radius,A));
|
||||
*/
|
||||
}
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
return ( r < 2.5 ? -2.0*c1*c2*r*std::exp(-c2*r*r) :
|
||||
( r < cutoff_radius ? A*std::pow(B,A)*(std::pow(2.0*cutoff_radius-r,-(A+1.0))-std::pow(r,-(A+1.0))) :
|
||||
0.0 ) );
|
||||
}
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if(cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
// std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if(aname == "a")
|
||||
{
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if(aname == "b")
|
||||
{
|
||||
ID_B = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(B,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_B,B);
|
||||
reset(A,B);
|
||||
real_type Y;
|
||||
c1 = evaluate(R0, Y, c2);
|
||||
c2 = -Y/(2.0*R0*c1);
|
||||
c1 *= std::exp(c2*R0*R0);
|
||||
// std::cout << "ChangMo's test: " << R0 << ", " << cutoff_radius << ", " << c1 << ", " << c2 << std::endl << evaluate(R0-0.0001) << ", " << evaluate(R0+0.0001) << std::endl;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,181 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Ken Esler, kpesler@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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_MODIFIED_PADEFUNCTION_H
|
||||
#define QMCPLUSPLUS_MODIFIED_PADEFUNCTION_H
|
||||
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
/** ModPade Jastrow functional
|
||||
*
|
||||
* \f[ u(r) = \frac{1}{2*A}\left[1-\exp(-A*r)\right], \f]
|
||||
*/
|
||||
template<class T>
|
||||
struct ModPadeFunctor: public OptimizableFunctorBase
|
||||
{
|
||||
|
||||
///coefficients
|
||||
real_type A;
|
||||
real_type B;
|
||||
real_type Zeff;
|
||||
|
||||
real_type Coeff;
|
||||
real_type mAB;
|
||||
|
||||
std::string ID_B;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
ModPadeFunctor(real_type a=-0.5, real_type b=1): A(a),B(b),Zeff(1.0)
|
||||
{
|
||||
reset();
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new ModPadeFunctor<T>(*this);
|
||||
}
|
||||
|
||||
/** reset the internal variables.
|
||||
*@param a New Jastrow parameter a
|
||||
*/
|
||||
inline void reset()
|
||||
{
|
||||
Coeff=A/B;
|
||||
mAB = -A*B;
|
||||
}
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
return Coeff*(1.0-std::exp(-B*r));
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
real_type expmar=std::exp(-B*r);
|
||||
dudr=A*expmar;
|
||||
d2udr2=mAB*expmar;
|
||||
return Coeff*(1.0-expmar);
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr,
|
||||
real_type& d2udr2, real_type d3udr3)
|
||||
{
|
||||
std::cerr << "Third derivative not implemented for ModPadeFunctor.\n";
|
||||
real_type expmar=std::exp(-B*r);
|
||||
dudr=A*expmar;
|
||||
d2udr2=mAB*expmar;
|
||||
return Coeff*(1.0-expmar);
|
||||
}
|
||||
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
return A*std::exp(-B*r);
|
||||
}
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
Zeff=1.0;
|
||||
//jastrow[iab]->put(cur->xmlChildrenNode,wfs_ref.RealVars);
|
||||
cur = cur->xmlChildrenNode;
|
||||
while(cur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(cur->name));
|
||||
if(cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname("B"),idname("0");
|
||||
OhmmsAttributeSet p;
|
||||
p.add(aname,"name");
|
||||
p.add(idname,"id");
|
||||
p.put(cur);
|
||||
if(aname == "A")
|
||||
{
|
||||
putContent(A,cur);
|
||||
}
|
||||
else
|
||||
if(aname == "B")
|
||||
{
|
||||
ID_B = idname;
|
||||
putContent(B,cur);
|
||||
}
|
||||
else
|
||||
if(aname == "Z")
|
||||
{
|
||||
putContent(Zeff,cur);
|
||||
}
|
||||
}
|
||||
cur = cur->next;
|
||||
}
|
||||
reset();
|
||||
myVars.insert(ID_B,B);
|
||||
return true;
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
inline void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int loc=myVars.where(0);
|
||||
if(loc>-1)
|
||||
B = myVars[0] = active[loc];
|
||||
Coeff=A/B;
|
||||
mAB = -A*B;
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -1,181 +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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_GAUSSIANSLATER_H
|
||||
#define QMCPLUSPLUS_GAUSSIANSLATER_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
template<class T>
|
||||
struct OpenGaussianSlaterFunctor: public OptimizableFunctorBase
|
||||
{
|
||||
/// f(r) = A exp(-(r)^2/C^2)
|
||||
///for smooth truncation f(r)->f(r)+f(2rc-r)-2f(rc)
|
||||
typedef typename OptimizableFunctorBase::real_type real_type;
|
||||
|
||||
///coefficients
|
||||
real_type A;
|
||||
real_type B;
|
||||
std::string ID_A;
|
||||
std::string ID_B;
|
||||
|
||||
/** constructor
|
||||
* @param a A coefficient
|
||||
* @param samespin boolean to indicate if this function is for parallel spins
|
||||
*/
|
||||
OpenGaussianSlaterFunctor(real_type a=1.0, real_type b=1.0 ):ID_A("G_A"), ID_B("G_B")
|
||||
{
|
||||
A=a;
|
||||
B=b;
|
||||
reset();
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new OpenGaussianSlaterFunctor<T>(*this);
|
||||
}
|
||||
|
||||
inline void reset()
|
||||
{
|
||||
}
|
||||
|
||||
/** reset the internal variables.
|
||||
*
|
||||
* USE_resetParameters
|
||||
*/
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ia=myVars.where(0);
|
||||
if (ia>-1)
|
||||
A=myVars[0]=active[ia];
|
||||
int ib=myVars.where(1);
|
||||
if (ib>-1)
|
||||
B=myVars[1]=active[ib];
|
||||
reset();
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
}
|
||||
|
||||
/** evaluate the value at r
|
||||
* @param r the distance
|
||||
* @return \f$ u(r) = \frac{A}{r}\left[1-\exp(-\frac{r}{F})\right]\f$
|
||||
*/
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
return A*r*r/(1.0+B*r);
|
||||
}
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@param d3udr3 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr,
|
||||
real_type& d2udr2, real_type& d3udr3)
|
||||
{
|
||||
real_type part1 = 1.0/(1.0+B*r);
|
||||
dudr = A*r*(2.0+B*r)*part1*part1;
|
||||
d2udr2 = 2.0*A*part1*part1*part1;
|
||||
d3udr3 = -6.0*A*B*part1*part1*part1*part1;
|
||||
return A*r*r*part1;
|
||||
}
|
||||
|
||||
|
||||
/**@param r the distance
|
||||
@param dudr first derivative
|
||||
@param d2udr2 second derivative
|
||||
@return the value
|
||||
*/
|
||||
inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
real_type part1 = 1.0/(1.0+B*r);
|
||||
dudr = A*r*(2.0+B*r)*part1*part1;
|
||||
d2udr2 = 2.0*A*part1*part1*part1;
|
||||
return A*r*r*part1;
|
||||
}
|
||||
|
||||
/** return a value at r
|
||||
*/
|
||||
real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
/** return a derivative at r
|
||||
*/
|
||||
real_type df(real_type r)
|
||||
{
|
||||
real_type part1 = 1.0/(1.0+B*r);
|
||||
return A*r*(2.0+B*r)*part1*part1;
|
||||
}
|
||||
|
||||
inline bool
|
||||
evaluateDerivatives(real_type r, std::vector<TinyVector<real_type,3> >& derivs)
|
||||
{
|
||||
real_type u = 1.0/(1.0+B*r);
|
||||
derivs[0][0]= r*r*u;
|
||||
derivs[1][0]= -A*r*r*r*u*u;
|
||||
derivs[0][1]= r*(2.0+B*r)*u*u;
|
||||
derivs[1][1]= -A*r*r*(3.0+B*r)*u*u*u;
|
||||
derivs[0][2]= 2.0*u*u*u;
|
||||
derivs[1][2]= -6.0*A*r*u*u*u*u;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
/** Read in the parameter from the xml input file.
|
||||
* @param cur current xmlNode from which the data members are reset
|
||||
*/
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
xmlNodePtr tcur = cur->xmlChildrenNode;
|
||||
while (tcur != NULL)
|
||||
{
|
||||
//@todo Var -> <param(eter) role="opt"/>
|
||||
std::string cname((const char*)(tcur->name));
|
||||
if (cname == "parameter" || cname == "Var")
|
||||
{
|
||||
std::string aname((const char*)(xmlGetProp(tcur,(const xmlChar *)"name")));
|
||||
// std::string idname((const char*)(xmlGetProp(tcur,(const xmlChar *)"id")));
|
||||
if (aname == "A")
|
||||
{
|
||||
ID_A = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(A,tcur);
|
||||
}
|
||||
else
|
||||
if (aname == "B")
|
||||
{
|
||||
ID_B = (const char*)(xmlGetProp(tcur,(const xmlChar *)"id"));
|
||||
putContent(B,tcur);
|
||||
}
|
||||
}
|
||||
tcur = tcur->next;
|
||||
}
|
||||
myVars.insert(ID_A,A);
|
||||
myVars.insert(ID_B,B);
|
||||
reset();
|
||||
return true;
|
||||
}
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,516 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#include "OhmmsPETE/OhmmsMatrix.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ThreeBodyBlockSparse.h"
|
||||
#include "Particle/DistanceTable.h"
|
||||
#include "Particle/DistanceTableData.h"
|
||||
#include "Numerics/DeterminantOperators.h"
|
||||
#include "Numerics/OhmmsBlas.h"
|
||||
#include "Numerics/BlockMatrixFunctions.h"
|
||||
#include "Numerics/LibxmlNumericIO.h"
|
||||
#include "Utilities/IteratorUtility.h"
|
||||
#include "Numerics/MatrixOperators.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
ThreeBodyBlockSparse::ThreeBodyBlockSparse(const ParticleSet& ions, ParticleSet& els):
|
||||
CenterRef(ions), GeminalBasis(0), IndexOffSet(1), ID_Lambda("j3"), SameBlocksForGroup(true)
|
||||
{
|
||||
myTableIndex=els.addTable(ions,DT_AOS);
|
||||
//d_table = DistanceTable::add(ions,els);
|
||||
NumPtcls=els.getTotalNum();
|
||||
Optimizable=true;
|
||||
}
|
||||
|
||||
ThreeBodyBlockSparse::~ThreeBodyBlockSparse()
|
||||
{
|
||||
delete_iter(LambdaBlocks.begin(), LambdaBlocks.end());
|
||||
}
|
||||
|
||||
OrbitalBase::RealType
|
||||
ThreeBodyBlockSparse::evaluateLog(ParticleSet& P,
|
||||
ParticleSet::ParticleGradient_t& G,
|
||||
ParticleSet::ParticleLaplacian_t& L)
|
||||
{
|
||||
evaluateLogAndStore(P);
|
||||
//GeminalBasis->evaluateForWalkerMove(P);
|
||||
////this could be better but it is used only sparsely
|
||||
////MatrixOperators::product(GeminalBasis->Y, Lambda, V);
|
||||
//MatrixOperators::product(GeminalBasis->Y, Lambda, V, BlockOffset);
|
||||
//Uk=0.0;
|
||||
//LogValue=ValueType();
|
||||
//for(int i=0; i< NumPtcls-1; i++) {
|
||||
// const BasisSetType::RealType* restrict yptr=GeminalBasis->Y[i];
|
||||
// for(int j=i+1; j<NumPtcls; j++) {
|
||||
// RealType x= dot(V[j],yptr,BasisSize);
|
||||
// LogValue += x;
|
||||
// Uk[i]+= x;
|
||||
// Uk[j]+= x;
|
||||
// }
|
||||
//}
|
||||
//for(int i=0; i<NumPtcls; i++) {
|
||||
// const BasisSetType::GradType* restrict dptr=GeminalBasis->dY[i];
|
||||
// const BasisSetType::ValueType* restrict d2ptr=GeminalBasis->d2Y[i];
|
||||
// const BasisSetType::ValueType* restrict vptr=V[0];
|
||||
// BasisSetType::GradType grad(0.0);
|
||||
// BasisSetType::ValueType lap(0.0);
|
||||
// for(int j=0; j<NumPtcls; j++, vptr+=BasisSize) {
|
||||
// if(j!=i) {
|
||||
// grad += dot(vptr,dptr,BasisSize);
|
||||
// lap += dot(vptr,d2ptr,BasisSize);
|
||||
// }
|
||||
// }
|
||||
// G(i)+=grad;
|
||||
// L(i)+=lap;
|
||||
//}
|
||||
return LogValue;
|
||||
}
|
||||
|
||||
OrbitalBase::ValueType
|
||||
ThreeBodyBlockSparse::ratio(ParticleSet& P, int iat)
|
||||
{
|
||||
UpdateMode=ORB_PBYP_RATIO;
|
||||
GeminalBasis->evaluateForPtclMove(P,iat);
|
||||
const BasisSetType::ValueType* restrict y_ptr=GeminalBasis->Phi.data();
|
||||
//This is the only difference from ThreeBodyGeminal
|
||||
RealType* restrict cv_ptr = curV.data();
|
||||
for(int b=0; b<Blocks.size(); b++)
|
||||
{
|
||||
int nb=Blocks[b];
|
||||
if(nb)
|
||||
{
|
||||
//GEMV<RealType,0>::apply(LambdaBlocks[b]->data(),y_ptr,cv_ptr,nb,nb);
|
||||
GEMV<RealType,0>::apply(LambdaBlocks[BlockID[b]]->data(),y_ptr,cv_ptr,nb,nb);
|
||||
for(int ib=0,k=BlockOffset[b]; ib<nb; k++,ib++)
|
||||
delV[k] = (*cv_ptr++)-V[iat][k];
|
||||
y_ptr+=nb;
|
||||
}
|
||||
}
|
||||
diffVal=0.0;
|
||||
const RealType* restrict vptr=V[0];
|
||||
for(int j=0; j<NumPtcls; j++, vptr+=BasisSize)
|
||||
{
|
||||
if(j == iat)
|
||||
continue;
|
||||
diffVal+= (curVal[j]=simd::dot(delV.data(),Y[j],BasisSize));
|
||||
}
|
||||
curVal[iat]=diffVal;
|
||||
return std::exp(diffVal);
|
||||
}
|
||||
|
||||
void ThreeBodyBlockSparse::restore(int iat)
|
||||
{
|
||||
//nothing to do here
|
||||
}
|
||||
|
||||
void ThreeBodyBlockSparse::acceptMove(ParticleSet& P, int iat)
|
||||
{
|
||||
//add the differential
|
||||
LogValue += diffVal;
|
||||
Uk+=curVal; //accumulate the differences
|
||||
if(UpdateMode == ORB_PBYP_RATIO)
|
||||
{
|
||||
Y.replaceRow(GeminalBasis->Phi.data(),iat);
|
||||
V.replaceRow(curV.begin(),iat);
|
||||
}
|
||||
else
|
||||
{
|
||||
dUk.replaceRow(curGrad.begin(),iat);
|
||||
d2Uk.replaceRow(curLap.begin(),iat);
|
||||
dUk.add2Column(tGrad.begin(),iat);
|
||||
d2Uk.add2Column(tLap.begin(),iat);
|
||||
//Y.replaceRow(GeminalBasis->y(0),iat);
|
||||
//dY.replaceRow(GeminalBasis->dy(0),iat);
|
||||
//d2Y.replaceRow(GeminalBasis->d2y(0),iat);
|
||||
Y.replaceRow(GeminalBasis->Phi.data(),iat);
|
||||
dY.replaceRow(GeminalBasis->dPhi.data(),iat);
|
||||
d2Y.replaceRow(GeminalBasis->d2Phi.data(),iat);
|
||||
V.replaceRow(curV.begin(),iat);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
ThreeBodyBlockSparse::registerData(ParticleSet& P, WFBufferType& buf)
|
||||
{
|
||||
evaluateLogAndStore(P);
|
||||
FirstAddressOfdY=&(dY(0,0)[0]);
|
||||
LastAddressOfdY=FirstAddressOfdY+NumPtcls*BasisSize*DIM;
|
||||
FirstAddressOfgU=&(dUk(0,0)[0]);
|
||||
LastAddressOfgU = FirstAddressOfgU + NumPtcls*NumPtcls*DIM;
|
||||
buf.add(LogValue);
|
||||
buf.add(V.begin(), V.end());
|
||||
buf.add(Y.begin(), Y.end());
|
||||
buf.add(FirstAddressOfdY,LastAddressOfdY);
|
||||
buf.add(d2Y.begin(),d2Y.end());
|
||||
buf.add(Uk.begin(), Uk.end());
|
||||
buf.add(FirstAddressOfgU,LastAddressOfgU);
|
||||
buf.add(d2Uk.begin(), d2Uk.end());
|
||||
}
|
||||
|
||||
void
|
||||
ThreeBodyBlockSparse::evaluateLogAndStore(ParticleSet& P)
|
||||
{
|
||||
GeminalBasis->evaluateForWalkerMove(P);
|
||||
//this could be better but it is used only sparsely
|
||||
//MatrixOperators::product(GeminalBasis->Y, Lambda, V);
|
||||
MatrixOperators::product(GeminalBasis->Y, Lambda, V, BlockOffset);
|
||||
Y=GeminalBasis->Y;
|
||||
dY=GeminalBasis->dY;
|
||||
d2Y=GeminalBasis->d2Y;
|
||||
Uk=0.0;
|
||||
LogValue=RealType();
|
||||
for(int i=0; i< NumPtcls-1; i++)
|
||||
{
|
||||
const RealType* restrict yptr=GeminalBasis->Y[i];
|
||||
for(int j=i+1; j<NumPtcls; j++)
|
||||
{
|
||||
RealType x= simd::dot(V[j],yptr,BasisSize);
|
||||
LogValue += x;
|
||||
Uk[i]+= x;
|
||||
Uk[j]+= x;
|
||||
}
|
||||
}
|
||||
for(int i=0; i<NumPtcls; i++)
|
||||
{
|
||||
const BasisSetType::GradType* restrict dptr=GeminalBasis->dY[i];
|
||||
const BasisSetType::ValueType* restrict d2ptr=GeminalBasis->d2Y[i];
|
||||
const RealType* restrict vptr=V[0];
|
||||
BasisSetType::GradType grad(0.0);
|
||||
BasisSetType::ValueType lap(0.0);
|
||||
for(int j=0; j<NumPtcls; j++, vptr+=BasisSize)
|
||||
{
|
||||
if(j==i)
|
||||
{
|
||||
dUk(i,j) = 0.0;
|
||||
d2Uk(i,j)= 0.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
grad+= (dUk(i,j) = simd::dot(vptr,dptr,BasisSize));
|
||||
lap += (d2Uk(i,j)= simd::dot(vptr,d2ptr,BasisSize));
|
||||
}
|
||||
}
|
||||
P.G[i]+=grad;
|
||||
P.L[i]+=lap;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
ThreeBodyBlockSparse::copyFromBuffer(ParticleSet& P, WFBufferType& buf)
|
||||
{
|
||||
buf.get(LogValue);
|
||||
buf.get(V.begin(), V.end());
|
||||
buf.get(Y.begin(), Y.end());
|
||||
buf.get(FirstAddressOfdY,LastAddressOfdY);
|
||||
buf.get(d2Y.begin(),d2Y.end());
|
||||
buf.get(Uk.begin(), Uk.end());
|
||||
buf.get(FirstAddressOfgU,LastAddressOfgU);
|
||||
buf.get(d2Uk.begin(), d2Uk.end());
|
||||
}
|
||||
|
||||
OrbitalBase::RealType
|
||||
ThreeBodyBlockSparse::updateBuffer(ParticleSet& P, WFBufferType& buf,
|
||||
bool fromscratch)
|
||||
{
|
||||
evaluateLogAndStore(P);
|
||||
buf.put(LogValue);
|
||||
buf.put(V.begin(), V.end());
|
||||
buf.put(Y.begin(), Y.end());
|
||||
buf.put(FirstAddressOfdY,LastAddressOfdY);
|
||||
buf.put(d2Y.begin(),d2Y.end());
|
||||
buf.put(Uk.begin(), Uk.end());
|
||||
buf.put(FirstAddressOfgU,LastAddressOfgU);
|
||||
buf.put(d2Uk.begin(), d2Uk.end());
|
||||
return LogValue;
|
||||
}
|
||||
|
||||
void ThreeBodyBlockSparse::resetTargetParticleSet(ParticleSet& P)
|
||||
{
|
||||
//d_table = DistanceTable::add(CenterRef,P);
|
||||
GeminalBasis->resetTargetParticleSet(P);
|
||||
}
|
||||
|
||||
void ThreeBodyBlockSparse::checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
int ncur=active.size();
|
||||
GeminalBasis->checkInVariables(active);
|
||||
if(active.size()>ncur)
|
||||
GeminalBasis->checkInVariables(myVars);
|
||||
}
|
||||
|
||||
void ThreeBodyBlockSparse::checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
GeminalBasis->checkOutVariables(active);
|
||||
Optimizable=myVars.is_optimizable();
|
||||
app_log() << "<j3-variables>"<< std::endl;
|
||||
myVars.print(app_log());
|
||||
app_log() << "</j3-variables>"<< std::endl;
|
||||
}
|
||||
|
||||
///reset the value of all the Two-Body Jastrow functions
|
||||
void ThreeBodyBlockSparse::resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ii=0;
|
||||
for(int b=0; b<LambdaBlocks.size(); ++b)
|
||||
{
|
||||
Matrix<RealType>& m(*LambdaBlocks[b]);
|
||||
for(int ib=0; ib<m.rows(); ++ib)
|
||||
{
|
||||
int loc=myVars.where(ii++);
|
||||
if(loc>=0)
|
||||
m(ib,ib)=active[loc];
|
||||
}
|
||||
for(int ib=0; ib<m.rows(); ++ib)
|
||||
for(int jb=ib+1; jb<m.rows(); ++jb)
|
||||
{
|
||||
int loc=myVars.where(ii++);
|
||||
if(loc>=0)
|
||||
m(jb,ib)=m(ib,jb)=active[loc];
|
||||
}
|
||||
}
|
||||
//this is not necessary
|
||||
if(SameBlocksForGroup)
|
||||
checkLambda();
|
||||
GeminalBasis->resetParameters(active);
|
||||
for(int i=0; i<myVars.size(); ++i)
|
||||
if(myVars.where(i)>=0)
|
||||
myVars[i]=active[myVars.where(i)];
|
||||
//app_log() << "ThreeBodyBlockSparse::resetParameters " << std::endl << Lambda << std::endl;
|
||||
//for(int b=0; b<LambdaBlocks.size();++b)
|
||||
//{
|
||||
// app_log() << "3Blody Block " << b << std::endl << *LambdaBlocks[b] << std::endl;
|
||||
//}
|
||||
}
|
||||
|
||||
void ThreeBodyBlockSparse::reportStatus(std::ostream& os)
|
||||
{
|
||||
myVars.print(os);
|
||||
}
|
||||
|
||||
bool ThreeBodyBlockSparse::put(xmlNodePtr cur)
|
||||
{
|
||||
//BasisSize = GeminalBasis->TotalBasis;
|
||||
BasisSize = GeminalBasis->getBasisSetSize();
|
||||
app_log() << " The number of Geminal functions "
|
||||
<<"for Three-body Jastrow " << BasisSize << std::endl;
|
||||
app_log() << " The number of particles " << NumPtcls << std::endl;
|
||||
Lambda.resize(BasisSize,BasisSize);
|
||||
FreeLambda.resize(BasisSize,BasisSize);
|
||||
FreeLambda=false;
|
||||
//identity is the default
|
||||
for(int ib=0; ib<BasisSize; ib++)
|
||||
Lambda(ib,ib)=1.0;
|
||||
if(cur == NULL)
|
||||
{
|
||||
FreeLambda=true;
|
||||
}
|
||||
else
|
||||
{
|
||||
//read from an input nodes
|
||||
char coeffname[16];
|
||||
std::string aname("j3");
|
||||
std::string datatype("lambda");
|
||||
std::string sameblocks("yes");
|
||||
IndexOffSet=0;//defaults are c
|
||||
OhmmsAttributeSet attrib;
|
||||
attrib.add(aname,"id");
|
||||
attrib.add(aname,"name");
|
||||
attrib.add(datatype,"type");
|
||||
attrib.add(IndexOffSet,"offset");
|
||||
attrib.add(sameblocks,"sameBlocksForGroup");
|
||||
attrib.put(cur);
|
||||
SameBlocksForGroup = (sameblocks == "yes");
|
||||
ID_Lambda=aname;
|
||||
if(datatype.find("rray")<datatype.size())
|
||||
{
|
||||
putContent(Lambda,cur);
|
||||
FreeLambda=true;
|
||||
}
|
||||
else
|
||||
{
|
||||
xmlNodePtr tcur=cur->xmlChildrenNode;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
if(xmlStrEqual(tcur->name,(const xmlChar*)"lambda"))
|
||||
{
|
||||
int iIn=atoi((const char*)(xmlGetProp(tcur,(const xmlChar*)"i")));
|
||||
int jIn=atoi((const char*)(xmlGetProp(tcur,(const xmlChar*)"j")));
|
||||
int i=iIn-IndexOffSet;
|
||||
int j=jIn-IndexOffSet;
|
||||
double c=atof((const char*)(xmlGetProp(tcur,(const xmlChar*)"c")));
|
||||
Lambda(i,j)=c;
|
||||
FreeLambda(i,j)=true;
|
||||
if(i != j)
|
||||
Lambda(j,i)=c;
|
||||
//sprintf(coeffname,"%s_%d_%d",aname.c_str(),iIn,jIn);
|
||||
//varlist[coeffname]=c;
|
||||
}
|
||||
tcur=tcur->next;
|
||||
}
|
||||
}
|
||||
}
|
||||
V.resize(NumPtcls,BasisSize);
|
||||
Y.resize(NumPtcls,BasisSize);
|
||||
dY.resize(NumPtcls,BasisSize);
|
||||
d2Y.resize(NumPtcls,BasisSize);
|
||||
curGrad.resize(NumPtcls);
|
||||
curLap.resize(NumPtcls);
|
||||
curVal.resize(NumPtcls);
|
||||
tGrad.resize(NumPtcls);
|
||||
tLap.resize(NumPtcls);
|
||||
curV.resize(BasisSize);
|
||||
delV.resize(BasisSize);
|
||||
Uk.resize(NumPtcls);
|
||||
dUk.resize(NumPtcls,NumPtcls);
|
||||
d2Uk.resize(NumPtcls,NumPtcls);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
void ThreeBodyBlockSparse::setBlocks(const std::vector<int>& blockspergroup)
|
||||
{
|
||||
//copyc this
|
||||
BlocksPerGroup=blockspergroup;
|
||||
//test with three blocks
|
||||
Blocks.resize(CenterRef.getTotalNum());
|
||||
for(int i=0; i<CenterRef.getTotalNum(); i++)
|
||||
Blocks[i]=blockspergroup[CenterRef.GroupID[i]];
|
||||
BlockOffset.resize(Blocks.size()+1,0);
|
||||
for(int i=0; i<Blocks.size(); i++)
|
||||
BlockOffset[i+1]=BlockOffset[i]+Blocks[i];
|
||||
BlockID.resize(Blocks.size(),-1);
|
||||
if(LambdaBlocks.empty())
|
||||
{
|
||||
std::vector<int> firstK,lastK;
|
||||
if(SameBlocksForGroup)
|
||||
{
|
||||
std::vector<int> need2copy(CenterRef.getSpeciesSet().getTotalNum(),-1);
|
||||
firstK.resize(need2copy.size(),-1);
|
||||
lastK.resize(need2copy.size(),-1);
|
||||
LambdaBlocks.resize(need2copy.size(),0);
|
||||
for(int b=0; b<Blocks.size(); b++)
|
||||
{
|
||||
if(Blocks[b] == 0)
|
||||
continue;
|
||||
int gid=need2copy[CenterRef.GroupID[b]];
|
||||
if(gid<0)//assign the current block index
|
||||
{
|
||||
gid=CenterRef.GroupID[b];
|
||||
need2copy[gid] = gid;
|
||||
LambdaBlocks[gid]= new Matrix<RealType>(Blocks[b],Blocks[b]);
|
||||
//populate LambdaBlocks with input Lambda
|
||||
Matrix<RealType>& m(*LambdaBlocks[gid]);
|
||||
firstK[gid]=BlockOffset[b];
|
||||
lastK[gid]=BlockOffset[b+1];
|
||||
for(int k=BlockOffset[b],ib=0; k<BlockOffset[b+1]; ++k,++ib)
|
||||
for(int kp=BlockOffset[b],jb=0; kp<BlockOffset[b+1]; ++kp,++jb)
|
||||
m(ib,jb)=Lambda(k,kp);
|
||||
}
|
||||
BlockID[b]=gid;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
LambdaBlocks.resize(Blocks.size(),0);
|
||||
firstK.resize(Blocks.size(),-1);
|
||||
lastK.resize(Blocks.size(),-1);
|
||||
for(int b=0; b<Blocks.size(); b++)
|
||||
{
|
||||
if(Blocks[b] ==0)
|
||||
continue;
|
||||
BlockID[b]=b;
|
||||
LambdaBlocks[b]=new Matrix<RealType>(Blocks[b],Blocks[b]);
|
||||
Matrix<RealType>& m(*LambdaBlocks[b]);
|
||||
firstK[b]=BlockOffset[b];
|
||||
lastK[b]=BlockOffset[b+1];
|
||||
for(int k=BlockOffset[b],ib=0; k<BlockOffset[b+1]; ++k,++ib)
|
||||
for(int kp=BlockOffset[b],jb=0; kp<BlockOffset[b+1]; ++kp,++jb)
|
||||
m(ib,jb)=Lambda(k,kp);
|
||||
}
|
||||
}
|
||||
myVars.clear();
|
||||
for(int b=0; b<LambdaBlocks.size(); ++b)
|
||||
{
|
||||
int first=firstK[b];
|
||||
int last=lastK[b];
|
||||
Matrix<RealType>& m(*LambdaBlocks[b]);
|
||||
for(int k=first,ib=0; k<last; ++k,++ib)
|
||||
{
|
||||
std::ostringstream sstr;
|
||||
sstr << ID_Lambda<<"_"<<k<< "_"<<k;
|
||||
myVars.insert(sstr.str(),m(ib,ib),true);
|
||||
}
|
||||
for(int k=first,ib=0; k<last; ++k,++ib)
|
||||
for(int kp=k+1,jb=ib+1; kp<last; ++kp,++jb)
|
||||
{
|
||||
std::ostringstream sstr;
|
||||
sstr << ID_Lambda<<"_"<<k<< "_"<<kp;
|
||||
myVars.insert(sstr.str(),m(ib,jb),true);
|
||||
}
|
||||
}
|
||||
}
|
||||
if(SameBlocksForGroup)
|
||||
checkLambda();
|
||||
app_log() << " ThreeBodyBlockSparse handles optimizable variables " << std::endl;
|
||||
myVars.print(app_log());
|
||||
}
|
||||
|
||||
/** set dependent Lambda
|
||||
*/
|
||||
void ThreeBodyBlockSparse::checkLambda()
|
||||
{
|
||||
for(int b=0; b<Blocks.size(); b++)
|
||||
{
|
||||
if(Blocks[b] ==0)
|
||||
continue;
|
||||
Matrix<RealType>& m(*LambdaBlocks[BlockID[b]]);
|
||||
for(int k=BlockOffset[b],ib=0; k<BlockOffset[b+1]; k++,ib++)
|
||||
for(int kp=BlockOffset[b],jb=0; kp<BlockOffset[b+1]; kp++,jb++)
|
||||
Lambda(k,kp)=m(ib,jb);
|
||||
}
|
||||
}
|
||||
|
||||
OrbitalBasePtr ThreeBodyBlockSparse::makeClone(ParticleSet& tqp) const
|
||||
{
|
||||
ThreeBodyBlockSparse* myclone=new ThreeBodyBlockSparse(CenterRef,tqp);
|
||||
myclone->GeminalBasis=GeminalBasis->makeClone();
|
||||
myclone->GeminalBasis->resetTargetParticleSet(tqp);
|
||||
myclone->put(NULL);//to resize data
|
||||
myclone->Lambda=Lambda;
|
||||
myclone->FreeLambda=FreeLambda;
|
||||
myclone->ID_Lambda=ID_Lambda;
|
||||
myclone->setBlocks(BlocksPerGroup);
|
||||
return myclone;
|
||||
}
|
||||
|
||||
//void ThreeBodyBlockSparse::addOptimizables(OptimizableSetType& varlist)
|
||||
//{
|
||||
// char coeffname[16];
|
||||
// for(int ib=0; ib<BasisSize; ib++) {
|
||||
// sprintf(coeffname,"%s_%d_%d",ID_Lambda.c_str(),ib+IndexOffSet,ib+IndexOffSet);
|
||||
// varlist[coeffname]=Lambda(ib,ib);
|
||||
// for(int jb=ib+1; jb<BasisSize; jb++) {
|
||||
// sprintf(coeffname,"%s_%d_%d",ID_Lambda.c_str(),ib+IndexOffSet,jb+IndexOffSet);
|
||||
// Lambda(jb,ib) = Lambda(ib,jb);
|
||||
// varlist[coeffname]=Lambda(ib,jb);
|
||||
// }
|
||||
// }
|
||||
//}
|
||||
}
|
||||
|
||||
|
|
@ -1,151 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_THREEBODY_BLOCKSPARSE_H
|
||||
#define QMCPLUSPLUS_THREEBODY_BLOCKSPARSE_H
|
||||
#include "Configuration.h"
|
||||
#include "QMCWaveFunctions/OrbitalBase.h"
|
||||
#include "OhmmsPETE/OhmmsVector.h"
|
||||
#include "OhmmsPETE/OhmmsMatrix.h"
|
||||
#include "QMCWaveFunctions/BasisSetBase.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** @ingroup OrbitalComponent
|
||||
* @brief ThreeBodyBlockSparse functions
|
||||
*/
|
||||
class ThreeBodyBlockSparse: public OrbitalBase
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
typedef BasisSetBase<RealType> BasisSetType;
|
||||
|
||||
///constructor
|
||||
ThreeBodyBlockSparse(const ParticleSet& ions, ParticleSet& els);
|
||||
|
||||
~ThreeBodyBlockSparse();
|
||||
|
||||
//implement virtual functions for optimizations
|
||||
void checkInVariables(opt_variables_type& active);
|
||||
void checkOutVariables(const opt_variables_type& active);
|
||||
void resetParameters(const opt_variables_type& active);
|
||||
void reportStatus(std::ostream& os);
|
||||
//evaluate the distance table with els
|
||||
void resetTargetParticleSet(ParticleSet& P);
|
||||
|
||||
RealType evaluateLog(ParticleSet& P,
|
||||
ParticleSet::ParticleGradient_t& G,
|
||||
ParticleSet::ParticleLaplacian_t& L);
|
||||
|
||||
ValueType ratio(ParticleSet& P, int iat);
|
||||
|
||||
void restore(int iat);
|
||||
|
||||
void acceptMove(ParticleSet& P, int iat);
|
||||
|
||||
void registerData(ParticleSet& P, WFBufferType& buf);
|
||||
|
||||
RealType updateBuffer(ParticleSet& P, WFBufferType& buf, bool fromscratch=false);
|
||||
|
||||
void copyFromBuffer(ParticleSet& P, WFBufferType& buf);
|
||||
|
||||
OrbitalBasePtr makeClone(ParticleSet& tqp) const;
|
||||
|
||||
void setBasisSet(BasisSetType* abasis)
|
||||
{
|
||||
GeminalBasis=abasis;
|
||||
}
|
||||
|
||||
bool put(xmlNodePtr cur);
|
||||
|
||||
//set blocks
|
||||
void setBlocks(const std::vector<int>& blockspergroup);
|
||||
|
||||
///reference to the center
|
||||
const ParticleSet& CenterRef;
|
||||
///assign same blocks for the group
|
||||
bool SameBlocksForGroup;
|
||||
///index of the table for source-target
|
||||
int myTableIndex;
|
||||
///size of the localized basis set
|
||||
int BasisSize;
|
||||
///number of particles
|
||||
int NumPtcls;
|
||||
///offset of the index
|
||||
int IndexOffSet;
|
||||
/** temporary value for update */
|
||||
RealType diffVal;
|
||||
///root name for Lambda compoenents
|
||||
std::string ID_Lambda;
|
||||
/** Y(iat,ibasis) value of the iat-th ortbial, the basis index ibasis
|
||||
*/
|
||||
Matrix<RealType> Y;
|
||||
/** dY(iat,ibasis) value of the iat-th ortbial, the basis index ibasis
|
||||
*/
|
||||
Matrix<PosType> dY;
|
||||
/** d2Y(iat,ibasis) value of the iat-th ortbial, the basis index ibasis
|
||||
*/
|
||||
Matrix<RealType> d2Y;
|
||||
/** V(i,j) = Lambda(k,kk) U(i,kk)
|
||||
*/
|
||||
Matrix<RealType> V;
|
||||
|
||||
/** Symmetric matrix connecting Geminal Basis functions */
|
||||
Matrix<RealType> Lambda;
|
||||
/** boolean to enable/disable optmization of Lambda(i,j) component */
|
||||
Matrix<int> FreeLambda;
|
||||
std::vector<IndexType> BlocksPerGroup;
|
||||
std::vector<IndexType> Blocks;
|
||||
std::vector<IndexType> BlockOffset;
|
||||
std::vector<IndexType> BlockID;
|
||||
std::vector<Matrix<RealType>* > LambdaBlocks;
|
||||
|
||||
/** Uk[i] = \sum_j dot(U[i],V[j]) */
|
||||
Vector<RealType> Uk;
|
||||
|
||||
/** Gradient for update mode */
|
||||
Matrix<PosType> dUk;
|
||||
|
||||
/** Laplacian for update mode */
|
||||
Matrix<RealType> d2Uk;
|
||||
|
||||
/** temporary Laplacin for update */
|
||||
Vector<RealType> curLap, tLap;
|
||||
/** temporary Gradient for update */
|
||||
Vector<PosType> curGrad, tGrad;
|
||||
/** tempory Lambda*newY for update */
|
||||
Vector<RealType> curV;
|
||||
/** tempory Lambda*(newY-Y(iat)) for update */
|
||||
Vector<RealType> delV;
|
||||
/** tempory Lambda*(newY-Y(iat)) for update */
|
||||
Vector<RealType> curVal;
|
||||
|
||||
RealType *FirstAddressOfdY;
|
||||
RealType *LastAddressOfdY;
|
||||
RealType *FirstAddressOfgU;
|
||||
RealType *LastAddressOfgU;
|
||||
|
||||
/** Geminal basis function */
|
||||
BasisSetType *GeminalBasis;
|
||||
|
||||
/** evaluateLog and store data for particle-by-particle update */
|
||||
void evaluateLogAndStore(ParticleSet& P);
|
||||
|
||||
void checkLambda();
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -1,448 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#include "OhmmsPETE/OhmmsMatrix.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ThreeBodyGeminal.h"
|
||||
#include "Particle/DistanceTable.h"
|
||||
#include "Particle/DistanceTableData.h"
|
||||
#include "Numerics/DeterminantOperators.h"
|
||||
#include "Numerics/OhmmsBlas.h"
|
||||
#include "Numerics/LibxmlNumericIO.h"
|
||||
#include "Utilities/RandomGenerator.h"
|
||||
#include "Numerics/MatrixOperators.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
ThreeBodyGeminal::ThreeBodyGeminal(const ParticleSet& ions, ParticleSet& els):
|
||||
CenterRef(ions), GeminalBasis(0), IndexOffSet(1), ID_Lambda("j3")
|
||||
{
|
||||
d_table = DistanceTable::add(ions,els,DT_AOS);
|
||||
NumPtcls=els.getTotalNum();
|
||||
NormFac=1.0/static_cast<RealType>(NumPtcls*NumPtcls);
|
||||
Optimizable=true;
|
||||
}
|
||||
|
||||
ThreeBodyGeminal::~ThreeBodyGeminal()
|
||||
{
|
||||
//clean up
|
||||
}
|
||||
|
||||
void ThreeBodyGeminal::resetTargetParticleSet(ParticleSet& P)
|
||||
{
|
||||
d_table = DistanceTable::add(CenterRef,P,DT_AOS);
|
||||
GeminalBasis->resetTargetParticleSet(P);
|
||||
}
|
||||
|
||||
void ThreeBodyGeminal::checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
active.insertFrom(myVars);
|
||||
int ncur=active.size();
|
||||
GeminalBasis->checkInVariables(active);
|
||||
if(ncur!= active.size())
|
||||
GeminalBasis->checkInVariables(myVars);
|
||||
}
|
||||
|
||||
void ThreeBodyGeminal::checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
myVars.getIndex(active);
|
||||
GeminalBasis->checkOutVariables(active);
|
||||
app_log() << "<j3-variables>"<< std::endl;
|
||||
myVars.print(app_log());
|
||||
app_log() << "</j3-variables>"<< std::endl;
|
||||
}
|
||||
|
||||
///reset the value of all the Two-Body Jastrow functions
|
||||
void ThreeBodyGeminal::resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
int ii=0; //counter for (i,j) pairs for i<j<BasisSize
|
||||
int aii=0;//counter for the variables that are meant to be optimized
|
||||
for(int ib=0; ib<BasisSize; ib++)
|
||||
{
|
||||
if(FreeLambda[ii++])
|
||||
{
|
||||
int loc=myVars.where(aii++);
|
||||
if(loc>=0)
|
||||
Lambda(ib,ib)=active[loc];
|
||||
}
|
||||
for(int jb=ib+1; jb<BasisSize; jb++)
|
||||
{
|
||||
if(FreeLambda[ii++])
|
||||
{
|
||||
int loc=myVars.where(aii++);
|
||||
if(loc>=0)
|
||||
Lambda(jb,ib)=Lambda(ib,jb)=active[loc];
|
||||
}
|
||||
}
|
||||
}
|
||||
GeminalBasis->resetParameters(active);
|
||||
for(int i=0; i<myVars.size(); ++i)
|
||||
if(myVars.where(i)>=0)
|
||||
myVars[i]=active[myVars.where(i)];
|
||||
//app_log() << "ThreeBodyGeminal::resetParameters" << std::endl;
|
||||
//app_log() << Lambda << std::endl;
|
||||
}
|
||||
|
||||
void ThreeBodyGeminal::reportStatus(std::ostream& os)
|
||||
{
|
||||
myVars.print(os);
|
||||
}
|
||||
|
||||
OrbitalBase::RealType
|
||||
ThreeBodyGeminal::evaluateLog(ParticleSet& P,
|
||||
ParticleSet::ParticleGradient_t& G,
|
||||
ParticleSet::ParticleLaplacian_t& L)
|
||||
{
|
||||
//this is necessary to handle ratio for pseudopotentials
|
||||
evaluateLogAndStore(P);
|
||||
////GeminalBasis->evaluate(P);
|
||||
//GeminalBasis->evaluateForWalkerMove(P);
|
||||
//
|
||||
//MatrixOperators::product(GeminalBasis->Y, Lambda, V);
|
||||
////Rewrite with gemm
|
||||
////for(int i=0; i<NumPtcls; i++) {
|
||||
//// for(int k=0; k<BasisSize; k++) {
|
||||
//// V(i,k) = BLAS::dot(BasisSize,GeminalBasis->Y[i],Lambda[k]);
|
||||
//// }
|
||||
////}
|
||||
//Uk=0.0;
|
||||
//LogValue=ValueType();
|
||||
//for(int i=0; i< NumPtcls-1; i++) {
|
||||
// const RealType* restrict yptr=GeminalBasis->Y[i];
|
||||
// for(int j=i+1; j<NumPtcls; j++) {
|
||||
// RealType x= dot(V[j],yptr,BasisSize);
|
||||
// LogValue += x;
|
||||
// Uk[i]+= x;
|
||||
// Uk[j]+= x;
|
||||
// }
|
||||
//}
|
||||
//for(int i=0; i<NumPtcls; i++) {
|
||||
// const PosType* restrict dptr=GeminalBasis->dY[i];
|
||||
// const RealType* restrict d2ptr=GeminalBasis->d2Y[i];
|
||||
// const RealType* restrict vptr=V[0];
|
||||
// GradType grad(0.0);
|
||||
// ValueType lap(0.0);
|
||||
// for(int j=0; j<NumPtcls; j++, vptr+=BasisSize) {
|
||||
// if(j!=i) {
|
||||
// grad += dot(vptr,dptr,BasisSize);
|
||||
// lap += dot(vptr,d2ptr,BasisSize);
|
||||
// }
|
||||
// }
|
||||
// G(i)+=grad;
|
||||
// L(i)+=lap;
|
||||
//}
|
||||
return LogValue;
|
||||
}
|
||||
|
||||
OrbitalBase::ValueType
|
||||
ThreeBodyGeminal::ratio(ParticleSet& P, int iat)
|
||||
{
|
||||
UpdateMode=ORB_PBYP_RATIO;
|
||||
GeminalBasis->evaluateForPtclMove(P,iat);
|
||||
const BasisSetType::RealType* restrict y_ptr=GeminalBasis->Phi.data();
|
||||
for(int k=0; k<BasisSize; k++)
|
||||
{
|
||||
//curV[k] = BLAS::dot(BasisSize,y_ptr,Lambda[k]);
|
||||
curV[k] = simd::dot(y_ptr,Lambda[k],BasisSize);
|
||||
delV[k] = curV[k]-V[iat][k];
|
||||
}
|
||||
diffVal=0.0;
|
||||
const RealType* restrict vptr=V[0];
|
||||
for(int j=0; j<NumPtcls; j++, vptr+=BasisSize)
|
||||
{
|
||||
if(j==iat)
|
||||
continue;
|
||||
diffVal+= (curVal[j]=simd::dot(delV.data(),Y[j],BasisSize));
|
||||
}
|
||||
curVal[iat]=diffVal;
|
||||
return std::exp(diffVal);
|
||||
}
|
||||
|
||||
void ThreeBodyGeminal::restore(int iat)
|
||||
{
|
||||
//nothing to do here
|
||||
}
|
||||
|
||||
void ThreeBodyGeminal::acceptMove(ParticleSet& P, int iat)
|
||||
{
|
||||
//add the differential
|
||||
LogValue += diffVal;
|
||||
Uk+=curVal; //accumulate the differences
|
||||
if(UpdateMode == ORB_PBYP_RATIO)
|
||||
{
|
||||
Y.replaceRow(GeminalBasis->Phi.data(),iat);
|
||||
V.replaceRow(curV.begin(),iat);
|
||||
}
|
||||
else
|
||||
{
|
||||
dUk.replaceRow(curGrad.begin(),iat);
|
||||
d2Uk.replaceRow(curLap.begin(),iat);
|
||||
dUk.add2Column(tGrad.begin(),iat);
|
||||
d2Uk.add2Column(tLap.begin(),iat);
|
||||
//Y.replaceRow(GeminalBasis->y(0),iat);
|
||||
//dY.replaceRow(GeminalBasis->dy(0),iat);
|
||||
//d2Y.replaceRow(GeminalBasis->d2y(0),iat);
|
||||
Y.replaceRow(GeminalBasis->Phi.data(),iat);
|
||||
dY.replaceRow(GeminalBasis->dPhi.data(),iat);
|
||||
d2Y.replaceRow(GeminalBasis->d2Phi.data(),iat);
|
||||
V.replaceRow(curV.begin(),iat);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
ThreeBodyGeminal::registerData(ParticleSet& P, WFBufferType& buf)
|
||||
{
|
||||
evaluateLogAndStore(P);
|
||||
FirstAddressOfdY=&(dY(0,0)[0]);
|
||||
LastAddressOfdY=FirstAddressOfdY+NumPtcls*BasisSize*DIM;
|
||||
FirstAddressOfgU=&(dUk(0,0)[0]);
|
||||
LastAddressOfgU = FirstAddressOfgU + NumPtcls*NumPtcls*DIM;
|
||||
buf.add(LogValue);
|
||||
buf.add(V.begin(), V.end());
|
||||
buf.add(Y.begin(), Y.end());
|
||||
buf.add(FirstAddressOfdY,LastAddressOfdY);
|
||||
buf.add(d2Y.begin(),d2Y.end());
|
||||
buf.add(Uk.begin(), Uk.end());
|
||||
buf.add(FirstAddressOfgU,LastAddressOfgU);
|
||||
buf.add(d2Uk.begin(), d2Uk.end());
|
||||
}
|
||||
|
||||
void
|
||||
ThreeBodyGeminal::evaluateLogAndStore(ParticleSet& P)
|
||||
{
|
||||
GeminalBasis->evaluateForWalkerMove(P);
|
||||
MatrixOperators::product(GeminalBasis->Y, Lambda, V);
|
||||
Y=GeminalBasis->Y;
|
||||
dY=GeminalBasis->dY;
|
||||
d2Y=GeminalBasis->d2Y;
|
||||
Uk=0.0;
|
||||
LogValue=RealType();
|
||||
for(int i=0; i< NumPtcls-1; i++)
|
||||
{
|
||||
const RealType* restrict yptr=GeminalBasis->Y[i];
|
||||
for(int j=i+1; j<NumPtcls; j++)
|
||||
{
|
||||
RealType x= simd::dot(V[j],yptr,BasisSize);
|
||||
LogValue += x;
|
||||
Uk[i]+= x;
|
||||
Uk[j]+= x;
|
||||
}
|
||||
}
|
||||
for(int i=0; i<NumPtcls; i++)
|
||||
{
|
||||
const PosType* restrict dptr=GeminalBasis->dY[i];
|
||||
const RealType* restrict d2ptr=GeminalBasis->d2Y[i];
|
||||
const RealType* restrict vptr=V[0];
|
||||
BasisSetType::GradType grad(0.0);
|
||||
BasisSetType::ValueType lap(0.0);
|
||||
for(int j=0; j<NumPtcls; j++, vptr+=BasisSize)
|
||||
{
|
||||
if(j==i)
|
||||
{
|
||||
dUk(i,j) = 0.0;
|
||||
d2Uk(i,j)= 0.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
grad+= (dUk(i,j) = simd::dot(vptr,dptr,BasisSize));
|
||||
lap += (d2Uk(i,j)= simd::dot(vptr,d2ptr,BasisSize));
|
||||
}
|
||||
}
|
||||
P.G[i]+=grad;
|
||||
P.L[i]+=lap;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
ThreeBodyGeminal::copyFromBuffer(ParticleSet& P, WFBufferType& buf)
|
||||
{
|
||||
buf.get(LogValue);
|
||||
buf.get(V.begin(), V.end());
|
||||
buf.get(Y.begin(), Y.end());
|
||||
buf.get(FirstAddressOfdY,LastAddressOfdY);
|
||||
buf.get(d2Y.begin(),d2Y.end());
|
||||
buf.get(Uk.begin(), Uk.end());
|
||||
buf.get(FirstAddressOfgU,LastAddressOfgU);
|
||||
buf.get(d2Uk.begin(), d2Uk.end());
|
||||
}
|
||||
|
||||
OrbitalBase::RealType
|
||||
ThreeBodyGeminal::updateBuffer(ParticleSet& P, WFBufferType& buf,
|
||||
bool fromscratch)
|
||||
{
|
||||
evaluateLogAndStore(P);
|
||||
buf.put(LogValue);
|
||||
buf.put(V.begin(), V.end());
|
||||
buf.put(Y.begin(), Y.end());
|
||||
buf.put(FirstAddressOfdY,LastAddressOfdY);
|
||||
buf.put(d2Y.begin(),d2Y.end());
|
||||
buf.put(Uk.begin(), Uk.end());
|
||||
buf.put(FirstAddressOfgU,LastAddressOfgU);
|
||||
buf.put(d2Uk.begin(), d2Uk.end());
|
||||
return LogValue;
|
||||
}
|
||||
|
||||
bool ThreeBodyGeminal::put(xmlNodePtr cur)
|
||||
{
|
||||
//BasisSize = GeminalBasis->TotalBasis;
|
||||
BasisSize = GeminalBasis->getBasisSetSize();
|
||||
app_log() << " The number of Geminal functions "
|
||||
<<"for Three-body Jastrow " << BasisSize << std::endl;
|
||||
app_log() << " The number of particles " << NumPtcls << std::endl;
|
||||
Lambda.resize(BasisSize,BasisSize);
|
||||
//disable lambda's so that element-by-element input can be handled
|
||||
FreeLambda.resize(BasisSize*(BasisSize+1)/2);
|
||||
FreeLambda=false;
|
||||
//zero is default
|
||||
Lambda=0.0;
|
||||
//for(int ib=0; ib<BasisSize; ib++) Lambda(ib,ib)=NormFac;
|
||||
//for(int ib=0; ib<BasisSize; ib++)
|
||||
// for(int jb=ib; jb<BasisSize; ++jb)
|
||||
// {
|
||||
// Lambda(ib,jb)=Random();
|
||||
// if(jb!=ib) Lambda(jb,ib)=Lambda(ib,jb);
|
||||
// }
|
||||
if(cur == NULL)
|
||||
{
|
||||
FreeLambda=true;
|
||||
}
|
||||
else
|
||||
{
|
||||
//read from an input nodes
|
||||
std::string aname("j3");
|
||||
std::string datatype("no");
|
||||
int sizeIn(0);
|
||||
IndexOffSet=1;
|
||||
OhmmsAttributeSet attrib;
|
||||
attrib.add(aname,"id");
|
||||
attrib.add(sizeIn,"size");
|
||||
attrib.add(aname,"name");
|
||||
attrib.add(datatype,"type");
|
||||
attrib.add(IndexOffSet,"offset");
|
||||
attrib.put(cur);
|
||||
ID_Lambda=aname;
|
||||
if(datatype.find("rray")<datatype.size())
|
||||
{
|
||||
if (sizeIn==Lambda.rows())
|
||||
{
|
||||
putContent(Lambda,cur);
|
||||
}
|
||||
FreeLambda=true;
|
||||
//addOptimizables(varlist);
|
||||
//symmetrize it
|
||||
//for(int ib=0; ib<BasisSize; ib++) {
|
||||
// sprintf(coeffname,"%s_%d_%d",aname.c_str(),ib+IndexOffSet,ib+IndexOffSet);
|
||||
// varlist[coeffname]=Lambda(ib,ib);
|
||||
// for(int jb=ib+1; jb<BasisSize; jb++) {
|
||||
// sprintf(coeffname,"%s_%d_%d",aname.c_str(),ib+IndexOffSet,jb+IndexOffSet);
|
||||
// Lambda(jb,ib) = Lambda(ib,jb);
|
||||
// varlist[coeffname]=Lambda(ib,jb);
|
||||
// }
|
||||
//}
|
||||
}
|
||||
else
|
||||
{
|
||||
xmlNodePtr tcur=cur->xmlChildrenNode;
|
||||
while(tcur != NULL)
|
||||
{
|
||||
if(xmlStrEqual(tcur->name,(const xmlChar*)"lambda"))
|
||||
{
|
||||
int iIn=atoi((const char*)(xmlGetProp(tcur,(const xmlChar*)"i")));
|
||||
int jIn=atoi((const char*)(xmlGetProp(tcur,(const xmlChar*)"j")));
|
||||
int i=iIn-IndexOffSet;
|
||||
int j=jIn-IndexOffSet;
|
||||
double c=atof((const char*)(xmlGetProp(tcur,(const xmlChar*)"c")));
|
||||
Lambda(i,j)=c;
|
||||
FreeLambda[i*BasisSize+j]=true;
|
||||
if(i != j)
|
||||
Lambda(j,i)=c;
|
||||
//sprintf(coeffname,"%s_%d_%d",aname.c_str(),iIn,jIn);
|
||||
//varlist[coeffname]=c;
|
||||
}
|
||||
tcur=tcur->next;
|
||||
}
|
||||
}
|
||||
}
|
||||
//myVars are set
|
||||
myVars.clear();
|
||||
std::string coeffname;
|
||||
int ii=0;
|
||||
for(int ib=0; ib<BasisSize; ib++)
|
||||
{
|
||||
if(FreeLambda[ii++])
|
||||
{
|
||||
coeffname=ID_Lambda+"_"+std::to_string(ib)+"_"+std::to_string(ib);
|
||||
myVars.insert(coeffname,Lambda(ib,ib));
|
||||
}
|
||||
for(int jb=ib+1; jb<BasisSize; jb++)
|
||||
{
|
||||
if(FreeLambda[ii++])
|
||||
{
|
||||
coeffname=ID_Lambda+"_"+std::to_string(ib)+"_"+std::to_string(jb);
|
||||
myVars.insert(coeffname,Lambda(ib,jb));
|
||||
}
|
||||
}
|
||||
}
|
||||
//app_log() << " Lambda Variables " << std::endl;
|
||||
//myVars.print(app_log());
|
||||
//app_log() << std::endl;
|
||||
V.resize(NumPtcls,BasisSize);
|
||||
Y.resize(NumPtcls,BasisSize);
|
||||
dY.resize(NumPtcls,BasisSize);
|
||||
d2Y.resize(NumPtcls,BasisSize);
|
||||
curGrad.resize(NumPtcls);
|
||||
curLap.resize(NumPtcls);
|
||||
curVal.resize(NumPtcls);
|
||||
tGrad.resize(NumPtcls);
|
||||
tLap.resize(NumPtcls);
|
||||
curV.resize(BasisSize);
|
||||
delV.resize(BasisSize);
|
||||
Uk.resize(NumPtcls);
|
||||
dUk.resize(NumPtcls,NumPtcls);
|
||||
d2Uk.resize(NumPtcls,NumPtcls);
|
||||
//app_log() << " Three-body Geminal coefficients " << std::endl;
|
||||
//app_log() << Lambda << std::endl;
|
||||
//GeminalBasis->resize(NumPtcls);
|
||||
return true;
|
||||
}
|
||||
|
||||
OrbitalBasePtr ThreeBodyGeminal::makeClone(ParticleSet& tqp) const
|
||||
{
|
||||
ThreeBodyGeminal* myclone=new ThreeBodyGeminal(CenterRef,tqp);
|
||||
myclone->GeminalBasis=GeminalBasis->makeClone();
|
||||
myclone->GeminalBasis->resetTargetParticleSet(tqp);
|
||||
myclone->put(NULL);
|
||||
myclone->myVars=myVars;
|
||||
myclone->Lambda=Lambda;
|
||||
myclone->ID_Lambda=ID_Lambda;
|
||||
myclone->FreeLambda=FreeLambda;
|
||||
return myclone;
|
||||
}
|
||||
|
||||
//void ThreeBodyGeminal::addOptimizables(OptimizableSetType& varlist)
|
||||
//{
|
||||
// char coeffname[16];
|
||||
// for(int ib=0; ib<BasisSize; ib++) {
|
||||
// sprintf(coeffname,"%s_%d_%d",ID_Lambda.c_str(),ib+IndexOffSet,ib+IndexOffSet);
|
||||
// varlist[coeffname]=Lambda(ib,ib);
|
||||
// for(int jb=ib+1; jb<BasisSize; jb++) {
|
||||
// sprintf(coeffname,"%s_%d_%d",ID_Lambda.c_str(),ib+IndexOffSet,jb+IndexOffSet);
|
||||
// Lambda(jb,ib) = Lambda(ib,jb);
|
||||
// varlist[coeffname]=Lambda(ib,jb);
|
||||
// }
|
||||
// }
|
||||
//}
|
||||
}
|
||||
|
||||
|
|
@ -1,140 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_THREEBODY_GEMINAL_H
|
||||
#define QMCPLUSPLUS_THREEBODY_GEMINAL_H
|
||||
#include "Configuration.h"
|
||||
#include "QMCWaveFunctions/OrbitalBase.h"
|
||||
#include "OhmmsPETE/OhmmsVector.h"
|
||||
#include "OhmmsPETE/OhmmsMatrix.h"
|
||||
#include "QMCWaveFunctions/BasisSetBase.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** @ingroup OrbitalComponent
|
||||
* @brief ThreeBodyGeminal functions
|
||||
*/
|
||||
class ThreeBodyGeminal: public OrbitalBase
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
typedef BasisSetBase<RealType> BasisSetType;
|
||||
|
||||
///constructor
|
||||
ThreeBodyGeminal(const ParticleSet& ions, ParticleSet& els);
|
||||
|
||||
~ThreeBodyGeminal();
|
||||
|
||||
//implement virtual functions for optimizations
|
||||
void checkInVariables(opt_variables_type& active);
|
||||
void checkOutVariables(const opt_variables_type& active);
|
||||
void resetParameters(const opt_variables_type& active);
|
||||
void reportStatus(std::ostream& os);
|
||||
void resetTargetParticleSet(ParticleSet& P);
|
||||
|
||||
RealType evaluateLog(ParticleSet& P,
|
||||
ParticleSet::ParticleGradient_t& G, ParticleSet::ParticleLaplacian_t& L);
|
||||
|
||||
ValueType ratio(ParticleSet& P, int iat);
|
||||
|
||||
void restore(int iat);
|
||||
|
||||
void acceptMove(ParticleSet& P, int iat);
|
||||
|
||||
void registerData(ParticleSet& P, WFBufferType& buf);
|
||||
|
||||
RealType updateBuffer(ParticleSet& P, WFBufferType& buf, bool fromscratch=false);
|
||||
|
||||
void copyFromBuffer(ParticleSet& P, WFBufferType& buf);
|
||||
|
||||
void setBasisSet(BasisSetType* abasis)
|
||||
{
|
||||
GeminalBasis=abasis;
|
||||
}
|
||||
|
||||
bool put(xmlNodePtr cur);
|
||||
|
||||
OrbitalBasePtr makeClone(ParticleSet& tqp) const;
|
||||
|
||||
private:
|
||||
|
||||
///reference to the center
|
||||
const ParticleSet& CenterRef;
|
||||
///distance table
|
||||
const DistanceTableData* d_table;
|
||||
///size of the localized basis set
|
||||
int BasisSize;
|
||||
///number of particles
|
||||
int NumPtcls;
|
||||
///offset of the index
|
||||
int IndexOffSet;
|
||||
///normalization factor 1.0/N*N
|
||||
RealType NormFac;
|
||||
/** temporary value for update */
|
||||
RealType diffVal;
|
||||
///root name for Lambda compoenents
|
||||
std::string ID_Lambda;
|
||||
/** Y(iat,ibasis) value of the iat-th ortbial, the basis index ibasis
|
||||
*/
|
||||
Matrix<RealType> Y;
|
||||
/** dY(iat,ibasis) value of the iat-th ortbial, the basis index ibasis
|
||||
*/
|
||||
Matrix<PosType> dY;
|
||||
/** d2Y(iat,ibasis) value of the iat-th ortbial, the basis index ibasis
|
||||
*/
|
||||
Matrix<RealType> d2Y;
|
||||
/** V(i,j) = Lambda(k,kk) U(i,kk)
|
||||
*/
|
||||
Matrix<RealType> V;
|
||||
|
||||
/** Symmetric matrix connecting Geminal Basis functions */
|
||||
Matrix<RealType> Lambda;
|
||||
/** boolean to enable/disable optmization of Lambda(i,j) component */
|
||||
Vector<int> FreeLambda;
|
||||
/** Uk[i] = \sum_j dot(U[i],V[j]) */
|
||||
Vector<RealType> Uk;
|
||||
|
||||
/** Gradient for update mode */
|
||||
Matrix<PosType> dUk;
|
||||
|
||||
/** Laplacian for update mode */
|
||||
Matrix<RealType> d2Uk;
|
||||
|
||||
/** temporary Laplacin for update */
|
||||
Vector<RealType> curLap, tLap;
|
||||
/** temporary Gradient for update */
|
||||
Vector<PosType> curGrad, tGrad;
|
||||
/** tempory Lambda*newY for update */
|
||||
Vector<RealType> curV;
|
||||
/** tempory Lambda*(newY-Y(iat)) for update */
|
||||
Vector<RealType> delV;
|
||||
/** tempory Lambda*(newY-Y(iat)) for update */
|
||||
Vector<RealType> curVal;
|
||||
|
||||
RealType *FirstAddressOfdY;
|
||||
RealType *LastAddressOfdY;
|
||||
RealType *FirstAddressOfgU;
|
||||
RealType *LastAddressOfgU;
|
||||
|
||||
/** Geminal basis function */
|
||||
BasisSetType *GeminalBasis;
|
||||
|
||||
/** evaluateLog and store data for particle-by-particle update */
|
||||
void evaluateLogAndStore(ParticleSet& P);
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -1,69 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
/**@file ThreeBodyGeminalBuilder.cpp
|
||||
*@brief definition of three-body jastrow of Geminal functions
|
||||
*/
|
||||
#include "QMCWaveFunctions/Jastrow/ThreeBodyGeminalBuilder.h"
|
||||
#include "QMCWaveFunctions/Jastrow/JastrowBasisBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/NGOBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/GTOBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/STOBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/MolecularBasisBuilder.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/GTOMolecularOrbitals.h"
|
||||
#include "QMCWaveFunctions/Jastrow/ThreeBodyGeminal.h"
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
ThreeBodyGeminalBuilder::ThreeBodyGeminalBuilder(ParticleSet& els,
|
||||
TrialWaveFunction& wfs,
|
||||
ParticleSet& ions):
|
||||
OrbitalBuilderBase(els,wfs), sourcePtcl(ions), basisBuilder(0)
|
||||
{
|
||||
}
|
||||
|
||||
bool ThreeBodyGeminalBuilder::put(xmlNodePtr cur)
|
||||
{
|
||||
ThreeBodyGeminal::BasisSetType *basisSet=0;
|
||||
bool foundBasisSet=false;
|
||||
xmlNodePtr basisPtr=NULL;
|
||||
xmlNodePtr coeffPtr=NULL;
|
||||
cur = cur->xmlChildrenNode;
|
||||
while(cur != NULL)
|
||||
{
|
||||
std::string cname((const char*)(cur->name));
|
||||
if(cname == basisset_tag)
|
||||
{
|
||||
basisPtr=cur;
|
||||
}
|
||||
else
|
||||
if(cname == "coefficient" || cname == "coefficients")
|
||||
{
|
||||
coeffPtr=cur;
|
||||
}
|
||||
cur=cur->next;
|
||||
}
|
||||
if(basisPtr != NULL)
|
||||
{
|
||||
ThreeBodyGeminal* J3 = new ThreeBodyGeminal(sourcePtcl, targetPtcl);
|
||||
basisBuilder = new JastrowBasisBuilder(targetPtcl,sourcePtcl,"gto",false);
|
||||
basisBuilder->put(basisPtr);
|
||||
J3->setBasisSet(basisBuilder->myBasisSet);
|
||||
J3->put(coeffPtr,targetPsi.VarList);
|
||||
targetPsi.addOrbital(J3);
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
|
@ -1,46 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
/**@file ThreeBodyGeminalBuilder.h
|
||||
*@brief declaration of three-body jastrow of Geminal functions
|
||||
*/
|
||||
#ifndef QMCPLUSPLUS_THREEBODY_GEMINALBUILDER_H
|
||||
#define QMCPLUSPLUS_THREEBODY_GEMINALBUILDER_H
|
||||
#include "QMCWaveFunctions/OrbitalBuilderBase.h"
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
class ThreeBodyGeminal;
|
||||
class BasisSetBuilder;
|
||||
|
||||
/**@ingroup WFSBuilder
|
||||
* @brief An abstract class for wave function builders
|
||||
*/
|
||||
class ThreeBodyGeminalBuilder: public OrbitalBuilderBase
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
ThreeBodyGeminalBuilder(ParticleSet& els, TrialWaveFunction& wfs,
|
||||
ParticleSet& ions);
|
||||
|
||||
/// process a xml node at cur
|
||||
bool put(xmlNodePtr cur);
|
||||
|
||||
protected:
|
||||
|
||||
ParticleSet& sourcePtcl;
|
||||
BasisSetBuilder* basisBuilder;
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,301 +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: 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_WMFUNCTOR_WITHPADECORRECTION_H
|
||||
#define QMCPLUSPLUS_WMFUNCTOR_WITHPADECORRECTION_H
|
||||
#include "Numerics/OptimizableFunctorBase.h"
|
||||
#include "OhmmsPETE/TinyVector.h"
|
||||
#include "OhmmsData/AttributeSet.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** Implements a screened Function \f$u[r]=(1-z(r/rcut))/(1+B*z(r/rcut)\f$
|
||||
*
|
||||
* Short-range functor introduced by Wagner and Mitas, cond-mat/0610088
|
||||
*/
|
||||
template<class T>
|
||||
struct WMFunctor: public OptimizableFunctorBase
|
||||
{
|
||||
///input B
|
||||
real_type B0;
|
||||
///input Rcut
|
||||
real_type Rcut;
|
||||
///1/Rcut
|
||||
real_type OneOverRc;
|
||||
///12/Rcut
|
||||
real_type DxDr;
|
||||
///12/Rcut/Rcut
|
||||
real_type D2xDr2;
|
||||
///id
|
||||
std::string ID_B;
|
||||
///name of B-attribute
|
||||
std::string attribName;
|
||||
///constructor
|
||||
WMFunctor(real_type b, real_type rc=7.5, const std::string& bname="exponent"):
|
||||
B0(b),Rcut(rc),attribName(bname)
|
||||
{
|
||||
reset();
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new WMFunctor<T>(*this);
|
||||
}
|
||||
|
||||
void reset()
|
||||
{
|
||||
//B0=b; Rcut=rc;
|
||||
OneOverRc=1.0/Rcut;
|
||||
DxDr=12*OneOverRc;
|
||||
D2xDr2=OneOverRc*DxDr;
|
||||
}
|
||||
|
||||
inline real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
|
||||
inline real_type df(real_type r)
|
||||
{
|
||||
if(r>Rcut)
|
||||
return 0.0;
|
||||
real_type x=r*OneOverRc;
|
||||
real_type z=x*x*(6.0-8*x+3.0*x*x);
|
||||
return -(1+B0)/(1+B0*z)/(1+B0*z)*OneOverRc*12*x*(1-2.0*x+x*x);
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
if(r>Rcut)
|
||||
return 0.0;
|
||||
real_type x=r*OneOverRc;
|
||||
real_type z=x*x*(6.0-8*x+3.0*x*x);
|
||||
return (1-z)/(1+B0*z);
|
||||
}
|
||||
|
||||
inline real_type
|
||||
evaluate(real_type r, real_type& dudr, real_type& d2udr2)
|
||||
{
|
||||
if(r>Rcut)
|
||||
{
|
||||
dudr=0.0;
|
||||
d2udr2=0.0;
|
||||
return 0.0;
|
||||
}
|
||||
real_type x=r*OneOverRc;
|
||||
real_type xx=x*x;
|
||||
real_type z=xx*(6.0-8*x+3.0*xx);
|
||||
real_type u=1.0/(1+B0*z);
|
||||
real_type du=(1.0+B0)*u*u;
|
||||
real_type dzdr=DxDr*x*(1.0-2.0*x+xx);
|
||||
dudr=-du*dzdr;
|
||||
d2udr2=du*(2.0*B0*u*dzdr*dzdr-D2xDr2*(1.0-4.0*x+3.0*xx));
|
||||
return (1.0-z)*u;
|
||||
}
|
||||
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
OhmmsAttributeSet rAttrib;
|
||||
rAttrib.add(ID_B,"id");
|
||||
rAttrib.add(ID_B,"ref");
|
||||
rAttrib.add(B0,"exponent");
|
||||
rAttrib.put(cur);
|
||||
ID_B.append("_E");
|
||||
myVars.insert(ID_B,B0,true);
|
||||
return true;
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
//disable optimization of E
|
||||
//active.insertFrom(myVars);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
//disable optimization of E
|
||||
//myVars.getIndex(active);
|
||||
}
|
||||
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
//disable optimization of E
|
||||
//int loc=myVars.where(0);
|
||||
//if(loc>=0) {myVars[0]=B0=active[loc];}
|
||||
}
|
||||
};
|
||||
|
||||
template<class T>
|
||||
struct WMFunctorSet: public OptimizableFunctorBase
|
||||
{
|
||||
//contraction and exponent pair
|
||||
typedef TinyVector<real_type,2> variable_type;
|
||||
///input Rcut
|
||||
real_type Rcut;
|
||||
///1/Rcut
|
||||
real_type OneOverRc;
|
||||
///12/Rcut
|
||||
real_type DxDr;
|
||||
///12/Rcut/Rcut
|
||||
real_type D2xDr2;
|
||||
///input Param[i][0]=B; Param[i][1]=C for i-th component
|
||||
std::vector<variable_type> Params;
|
||||
|
||||
WMFunctorSet(real_type rc=7.5):Rcut(rc)
|
||||
{
|
||||
OneOverRc=1.0/Rcut;
|
||||
DxDr=12*OneOverRc;
|
||||
D2xDr2=OneOverRc*DxDr;
|
||||
Params.reserve(8);//limit to 8
|
||||
}
|
||||
|
||||
OptimizableFunctorBase* makeClone() const
|
||||
{
|
||||
return new WMFunctorSet<T>(*this);
|
||||
}
|
||||
|
||||
inline void reset() {}
|
||||
|
||||
inline void addComponent(real_type c, real_type b, const std::string& aname)
|
||||
{
|
||||
int i=Params.size();
|
||||
Params.push_back(variable_type(c,b));
|
||||
std::string cname=aname+"_C";
|
||||
std::string ename=aname+"_E";
|
||||
myVars.insert(cname,Params[i][0]);
|
||||
myVars.insert(ename,Params[i][1]);
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r)
|
||||
{
|
||||
if(r>Rcut)
|
||||
return 0.0;
|
||||
real_type x=r*OneOverRc;
|
||||
real_type xx=x*x;
|
||||
real_type z=xx*(6.0-8*x+3.0*xx);
|
||||
real_type value=0.0;
|
||||
for(int i=0; i<Params.size(); ++i)
|
||||
value += Params[i][0]*(1.0-z)/(1+Params[i][1]*z);
|
||||
return value;
|
||||
}
|
||||
|
||||
inline real_type evaluate(real_type r, real_type& dv, real_type& d2v)
|
||||
{
|
||||
dv=0.0;
|
||||
d2v=0.0;
|
||||
if(r>Rcut)
|
||||
return 0.0;
|
||||
real_type x=r*OneOverRc;
|
||||
real_type xx=x*x;
|
||||
real_type z=xx*(6.0-8*x+3.0*xx);
|
||||
real_type dzdr=DxDr*x*(1.0-2.0*x+xx);
|
||||
real_type dzdrsq=dzdr*dzdr;
|
||||
real_type d2zdr2=D2xDr2*(1.0-4.0*x+3.0*xx);
|
||||
real_type v=0.0;
|
||||
for(int i=0; i<Params.size(); ++i)
|
||||
{
|
||||
real_type c=Params[i][0];
|
||||
real_type b=Params[i][1];
|
||||
real_type u=1.0/(1+b*z);
|
||||
real_type du=c*(1.0+b)*u*u;
|
||||
v += c*(1.0-z)*u;
|
||||
dv -= du*dzdr;
|
||||
d2v+= du*(2.0*b*u*dzdrsq-d2zdr2);
|
||||
}
|
||||
return v;
|
||||
}
|
||||
|
||||
inline bool
|
||||
evaluateDerivatives(real_type r, std::vector<TinyVector<real_type,3> >& derivs)
|
||||
{
|
||||
if(r>=Rcut)
|
||||
return false;
|
||||
real_type x=r*OneOverRc;
|
||||
real_type xx=x*x;
|
||||
real_type z=xx*(6.0-8*x+3.0*xx);
|
||||
real_type dzdr=DxDr*x*(1.0-2.0*x+xx);
|
||||
real_type dzdrsq=dzdr*dzdr;
|
||||
real_type d2zdr2=D2xDr2*(1.0-4.0*x+3.0*xx);
|
||||
typedef TinyVector<real_type,3> e_type;
|
||||
typename std::vector<e_type>::iterator dit(derivs.begin());
|
||||
for(int i=0,ii=0; i<Params.size(); ++i)
|
||||
{
|
||||
real_type c=Params[i][0];
|
||||
real_type b=Params[i][1];
|
||||
real_type u=1.0/(1+b*z);
|
||||
real_type uu=u*u;
|
||||
real_type du=(1.0+b)*uu;
|
||||
(*dit++)=e_type((1.0-z)*u,-du*dzdr,du*(2.0*b*u*dzdrsq-d2zdr2));
|
||||
//Disable derivatives with respect to exponents: incorrect and unstable
|
||||
//real_type dudz=(z*(b+2)-1.0)*uu*u;
|
||||
//(*dit++)=e_type(c*z*(z-1)*uu,c*dudz*dzdr,c*uu*uu*((b+2)*(1.0-2*b*z)+3*b)*dzdrsq+dudz*d2zdr2);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
inline real_type f(real_type r)
|
||||
{
|
||||
return evaluate(r);
|
||||
}
|
||||
inline real_type df(real_type r)
|
||||
{
|
||||
if(r>Rcut)
|
||||
return 0.0;
|
||||
real_type x=r*OneOverRc;
|
||||
real_type xx=x*x;
|
||||
real_type z=xx*(6.0-8*x+3.0*xx);
|
||||
real_type dzdr=DxDr*x*(1.0-2.0*x+xx);
|
||||
real_type dv=0.0;
|
||||
for(int i=0; i<Params.size(); ++i)
|
||||
{
|
||||
real_type c=Params[i][0];
|
||||
real_type b=Params[i][1];
|
||||
real_type u=1.0/(1+b*z);
|
||||
dv -= c*(1.0+b)*u*u*dzdr;
|
||||
}
|
||||
return dv;
|
||||
}
|
||||
|
||||
bool put(xmlNodePtr cur)
|
||||
{
|
||||
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 ii=0;
|
||||
for(int i=0; i<Params.size(); ++i)
|
||||
{
|
||||
int loc_c=myVars.where(ii++);
|
||||
int loc_e=myVars.where(ii++);
|
||||
if(loc_c>=0)
|
||||
Params[i][0]=active[loc_c];
|
||||
if(loc_e>=0)
|
||||
Params[i][1]=active[loc_e];
|
||||
}
|
||||
}
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,47 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/BsplineAOBuilder.h"
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
BsplineAOBuilder::BsplineAOBuilder(xmlNodePtr cur): m_orbitals(0)
|
||||
{
|
||||
if(cur != NULL)
|
||||
putCommon(cur);
|
||||
}
|
||||
|
||||
bool BsplineAOBuilder::putCommon(xmlNodePtr cur)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
bool
|
||||
BsplineAOBuilder::addRadialOrbital(xmlNodePtr cur, const QuantumNumberType& nlms, bool useSphericalHarmonicsNormalization)
|
||||
{
|
||||
if(!m_orbitals)
|
||||
{
|
||||
ERRORMSG("m_orbitals, SphericalOrbitals<ROT,GT>*, is not initialized")
|
||||
return false;
|
||||
}
|
||||
std::cout << "#### BsplineAOBuilder::addRadialOrbital " << std::endl;
|
||||
RadialOrbitalType* radorb= new RadialOrbitalType(0.0);
|
||||
radorb->put(cur);
|
||||
m_orbitals->Rnl.push_back(radorb);
|
||||
m_orbitals->RnlID.push_back(nlms);
|
||||
return true;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -1,58 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
#ifndef QMCPLUSPLUS_BSPLINE_AO_BUILDER_H
|
||||
#define QMCPLUSPLUS_BSPLINE_AO_BUILDER_H
|
||||
|
||||
#include "Configuration.h"
|
||||
#include "QMCWaveFunctions/MolecularOrbitals/SphericalBasisSet.h"
|
||||
#include "QMCWaveFunctions/Jastrow/BsplineFunctor.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
class BsplineAOBuilder: public QMCTraits
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
typedef BsplineFunctor<RealType> RadialOrbitalType;
|
||||
typedef SphericalBasisSet<RadialOrbitalType> CenteredOrbitalType;
|
||||
|
||||
///the radial orbitals
|
||||
CenteredOrbitalType* m_orbitals;
|
||||
///the species
|
||||
std::string m_species;
|
||||
///constructor
|
||||
BsplineAOBuilder(xmlNodePtr cur=NULL);
|
||||
|
||||
///assign a CenteredOrbitalType to work on
|
||||
void setOrbitalSet(CenteredOrbitalType* oset, const std::string& acenter)
|
||||
{
|
||||
m_orbitals = oset;
|
||||
m_species = acenter;
|
||||
}
|
||||
|
||||
bool addGrid(xmlNodePtr cur)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
bool addRadialOrbital(xmlNodePtr cur, const QuantumNumberType& nlms, bool useSphericalHarmonicsNormalization=true);
|
||||
|
||||
bool putCommon(xmlNodePtr cur);
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,298 +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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
|
||||
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
|
||||
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
|
||||
//
|
||||
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
|
||||
//////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
/** @file SparseSparseLocalizsedBasisSet.h
|
||||
* @author Jeongnim Kim
|
||||
* @brief A derived class from BasisSetBase which implements a sparse vector method to handle locazed basis sets.
|
||||
*/
|
||||
#ifndef QMCPLUSPLUS_SPARSELOCALIZEDBASISSET_H
|
||||
#define QMCPLUSPLUS_SPARSELOCALIZEDBASISSET_H
|
||||
|
||||
#include "QMCWaveFunctions/BasisSetBase.h"
|
||||
#include "Particle/DistanceTable.h"
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** A localized basis set derived from BasisSetBase<typename COT::value_type>
|
||||
*
|
||||
* This class performs the evaluation of the basis functions and their
|
||||
* derivatives for each of the N-particles in a configuration.
|
||||
* The template parameter COT denotes Centered-Orbital-Type which provides
|
||||
* a set of localized orbitals associated with a center. Unlike LocalizedBasisSet<COT>
|
||||
* this class does not assume that all the centers have basis functions, as the name
|
||||
* SparseLocalizedBasisSet implies.
|
||||
*/
|
||||
template<class COT>
|
||||
struct SparseLocalizedBasisSet: public BasisSetBase<typename COT::value_type>
|
||||
{
|
||||
typedef BasisSetBase<typename COT::value_type> BasisSetType;
|
||||
typedef typename BasisSetType::RealType RealType;
|
||||
typedef typename BasisSetType::ValueType ValueType;
|
||||
typedef typename BasisSetType::IndexType IndexType;
|
||||
typedef typename BasisSetType::IndexVector_t IndexVector_t;
|
||||
typedef typename BasisSetType::ValueVector_t ValueVector_t;
|
||||
typedef typename BasisSetType::ValueMatrix_t ValueMatrix_t;
|
||||
typedef typename BasisSetType::GradVector_t GradVector_t;
|
||||
typedef typename BasisSetType::GradMatrix_t GradMatrix_t;
|
||||
|
||||
using BasisSetType::BasisSetSize;
|
||||
using BasisSetType::Phi;
|
||||
using BasisSetType::dPhi;
|
||||
using BasisSetType::d2Phi;
|
||||
using BasisSetType::grad_grad_Phi;
|
||||
using BasisSetType::Y;
|
||||
using BasisSetType::dY;
|
||||
using BasisSetType::d2Y;
|
||||
///Reference to the center
|
||||
const ParticleSet& CenterRef;
|
||||
///number of centers, e.g., ions
|
||||
int NumCenters;
|
||||
///number of quantum particles
|
||||
int NumTargets;
|
||||
///number of groups
|
||||
int NumGroups;
|
||||
///BasisOffset
|
||||
std::vector<int> BasisOffset;
|
||||
|
||||
/** class to handle a COT and particle indices which share the COT
|
||||
*/
|
||||
struct LOForCenter
|
||||
{
|
||||
///Centered Orbital Type for this center
|
||||
COT* myO;
|
||||
///IDs for the particle which belong to the same center
|
||||
std::vector<int> myP;
|
||||
inline LOForCenter(COT* o=0):myO(o) {}
|
||||
LOForCenter(const LOForCenter& old):myP(old.myP)
|
||||
{
|
||||
myO=old.myO->makeClone();
|
||||
}
|
||||
};
|
||||
|
||||
typedef std::vector<LOForCenter*> ContainerType;
|
||||
ContainerType LOBasisSet;
|
||||
|
||||
/** distance table, e.g., ion-electron
|
||||
*
|
||||
* Localized basis sets require a pair relationship between CenterRef
|
||||
* and the quantum particle set.
|
||||
*/
|
||||
const DistanceTableData* myTable;
|
||||
|
||||
/** constructor
|
||||
* @param ions ionic system
|
||||
* @param els electronic system
|
||||
*/
|
||||
SparseLocalizedBasisSet(const ParticleSet& ions, ParticleSet& els):
|
||||
CenterRef(ions), myTable(0)
|
||||
{
|
||||
myTable = DistanceTable::add(ions,els,DT_AOS);
|
||||
NumCenters=CenterRef.getTotalNum();
|
||||
NumTargets=els.getTotalNum();
|
||||
NumGroups=CenterRef.getSpeciesSet().getTotalNum();
|
||||
BasisOffset.resize(NumCenters+1);
|
||||
LOBasisSet.resize(NumGroups,0);
|
||||
}
|
||||
|
||||
BasisSetBase<typename COT::value_type>* makeClone() const
|
||||
{
|
||||
SparseLocalizedBasisSet<COT>* myclone=new SparseLocalizedBasisSet<COT>(*this);
|
||||
for(int i=0; i<LOBasisSet.size(); ++i)
|
||||
myclone->LOBasisSet[i]=new LOForCenter(*LOBasisSet[i]);
|
||||
return myclone;
|
||||
}
|
||||
|
||||
void checkInVariables(opt_variables_type& active)
|
||||
{
|
||||
for(int i=0; i<LOBasisSet.size(); ++i)
|
||||
LOBasisSet[i]->myO->checkInVariables(active);
|
||||
}
|
||||
|
||||
void checkOutVariables(const opt_variables_type& active)
|
||||
{
|
||||
for(int i=0; i<LOBasisSet.size(); ++i)
|
||||
LOBasisSet[i]->myO->checkOutVariables(active);
|
||||
}
|
||||
|
||||
void resetParameters(const opt_variables_type& active)
|
||||
{
|
||||
for(int i=0; i<LOBasisSet.size(); ++i)
|
||||
LOBasisSet[i]->myO->resetParameters(active);
|
||||
}
|
||||
|
||||
/**
|
||||
@param atable the distance table (ion-electron)
|
||||
@brief Assign the distance table (ion-electron) and
|
||||
*determine the total number of basis states.
|
||||
*/
|
||||
void setBasisSetSize(int nbs)
|
||||
{
|
||||
if(nbs == BasisSetSize)
|
||||
return;
|
||||
if(myTable ==0)
|
||||
{
|
||||
APP_ABORT("SparseLocalizsedBasisSet cannot function without a distance table. Abort");
|
||||
}
|
||||
for(int ig=0; ig<NumGroups; ig++)
|
||||
if(LOBasisSet[ig])
|
||||
LOBasisSet[ig]->myO->setTable(myTable);
|
||||
BasisOffset[0]=0;
|
||||
for(int c=0; c<NumCenters; c++)
|
||||
{
|
||||
int ig=CenterRef.GroupID[c];
|
||||
if(LOBasisSet[ig])
|
||||
{
|
||||
LOBasisSet[ig]->myP.push_back(c);//add this center to the basis set group
|
||||
BasisOffset[c+1] = BasisOffset[c]+LOBasisSet[ig]->myO->getBasisSetSize();
|
||||
}
|
||||
else
|
||||
{
|
||||
BasisOffset[c+1] = BasisOffset[c];
|
||||
}
|
||||
}
|
||||
BasisSetSize = BasisOffset[NumCenters];
|
||||
//for(int ig=0; ig<NumGroups; ig++)
|
||||
//{
|
||||
// if(LOBasisSet[ig])
|
||||
// {
|
||||
// std::cout << "### SparseLOBasis " << ig << " : ";
|
||||
// const std::vector<int>& plist(LOBasisSet[ig]->myP);
|
||||
// for(int i=0; i<plist.size(); i++)
|
||||
// {
|
||||
// std::cout << plist[i] << "[" << BasisOffset[plist[i]] << "] ";
|
||||
// }
|
||||
// std::cout << std::endl;
|
||||
// }
|
||||
//}
|
||||
this->resize(NumTargets);
|
||||
}
|
||||
|
||||
/** reset the distance table with a new target P
|
||||
*/
|
||||
void resetTargetParticleSet(ParticleSet& P)
|
||||
{
|
||||
LOGMSG("SparseLocalizsedBasisSet::resetTargetParticleSet")
|
||||
myTable = DistanceTable::add(CenterRef,P,DT_AOS);
|
||||
for(int ig=0; ig<NumGroups; ig++)
|
||||
if(LOBasisSet[ig])
|
||||
LOBasisSet[ig]->myO->setTable(myTable);
|
||||
}
|
||||
|
||||
inline void
|
||||
evaluateWithHessian(const ParticleSet& P, int iat)
|
||||
{
|
||||
APP_ABORT("SparseLocalizsedBasisSet::evaluateWithHessian has not been implemented. \n");
|
||||
}
|
||||
|
||||
void evaluateWithThirdDeriv(const ParticleSet& P, int iat)
|
||||
{
|
||||
APP_ABORT("SparseLocalizsedBasisSet::evaluateWithThirdDeriv has not been implemented. \n");
|
||||
}
|
||||
|
||||
void evaluateThirdDerivOnly(const ParticleSet& P, int iat)
|
||||
{
|
||||
APP_ABORT("SparseLocalizsedBasisSet::evaluateThirdDerivOnly has not been implemented. \n");
|
||||
}
|
||||
|
||||
inline void
|
||||
evaluateForWalkerMove(const ParticleSet& P)
|
||||
{
|
||||
for(int ig=0; ig<NumGroups; ig++)
|
||||
{
|
||||
if(LOBasisSet[ig])
|
||||
{
|
||||
const std::vector<int>& plist(LOBasisSet[ig]->myP);
|
||||
COT& o(*(LOBasisSet[ig]->myO));
|
||||
for(int i=0; i<plist.size(); i++)
|
||||
{
|
||||
o.evaluateForWalkerMove(plist[i],0,NumTargets,BasisOffset[plist[i]],Y,dY,d2Y);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
inline void
|
||||
evaluateForWalkerMove(const ParticleSet& P, int iat)
|
||||
{
|
||||
for(int ig=0; ig<NumGroups; ig++)
|
||||
{
|
||||
if(LOBasisSet[ig])
|
||||
{
|
||||
const std::vector<int>& plist(LOBasisSet[ig]->myP);
|
||||
COT& o(*(LOBasisSet[ig]->myO));
|
||||
for(int i=0; i<plist.size(); i++)
|
||||
{
|
||||
o.evaluateForWalkerMove(plist[i],iat,BasisOffset[plist[i]],Phi,dPhi,d2Phi);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
inline void
|
||||
evaluateForPtclMove(const ParticleSet& P, int iat)
|
||||
{
|
||||
for(int ig=0; ig<NumGroups; ig++)
|
||||
{
|
||||
if(LOBasisSet[ig])
|
||||
{
|
||||
const std::vector<int>& plist(LOBasisSet[ig]->myP);
|
||||
COT& o(*(LOBasisSet[ig]->myO));
|
||||
for(int i=0; i<plist.size(); i++)
|
||||
{
|
||||
o.evaluateForPtclMove(plist[i],iat,BasisOffset[plist[i]],Phi);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
inline void
|
||||
evaluateAllForPtclMove(const ParticleSet& P, int iat)
|
||||
{
|
||||
for(int ig=0; ig<NumGroups; ig++)
|
||||
{
|
||||
if(LOBasisSet[ig])
|
||||
{
|
||||
const std::vector<int>& plist(LOBasisSet[ig]->myP);
|
||||
COT& o(*(LOBasisSet[ig]->myO));
|
||||
for(int i=0; i<plist.size(); i++)
|
||||
{
|
||||
o.evaluateAllForPtclMove(plist[i],iat,BasisOffset[plist[i]],Phi,dPhi,d2Phi);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
inline void
|
||||
evaluateForPtclMoveWithHessian(const ParticleSet& P, int iat)
|
||||
{
|
||||
APP_ABORT("SparseLocalizsedBasisSet::evaluateForPtclMoveWithHessian has not been implemented. \n");
|
||||
}
|
||||
|
||||
/** add a new set of Centered Atomic Orbitals
|
||||
* @param icenter the index of the center
|
||||
* @param aos a set of Centered Atomic Orbitals
|
||||
*/
|
||||
void add(int icenter, COT* aos)
|
||||
{
|
||||
aos->setTable(myTable);
|
||||
LOBasisSet[icenter]=new LOForCenter(aos);
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
Loading…
Reference in New Issue