scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Remove variable arrays and cleanup multi quintic.

This commit is contained in:
jnkim 2017-09-27 16:00:52 -07:00
parent e77437d09a
commit 097b16a1a2
3 changed files with 258 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
src/QMCWaveFunctions/lcao/MultiQuinticSpline1D.h
View File

@ -25,6 +25,40 @@
namespace qmcplusplus
{
template<typename T>
struct LogGridLight
{
T lower_bound;
T upper_bound;
T LogDelta;
T OneOverLogDelta;
template<typename T1>
inline void set(T1 ri, T1 rf, int n)
{
lower_bound=static_cast<T>(ri);
upper_bound=static_cast<T>(rf);
//num_points=n;
double ratio = rf/ri;
double log_ratio = std::log(ratio);
double dlog_ratio = log_ratio/static_cast<double>(n-1);
LogDelta = dlog_ratio;
OneOverLogDelta = 1.0/dlog_ratio;
}
inline int locate(T r) const
{
return static_cast<int>(std::log(r/lower_bound)*OneOverLogDelta);
}
inline T getCLForQuintic(T r, int& loc) const
{
constexpr T cone(1);
loc=static_cast<int>(std::log(r/lower_bound)*OneOverLogDelta);
return r-lower_bound*std::exp(loc*LogDelta);
}
};
/** multivalue implementation for OneDimQuintic
*/
template<typename T>
@ -32,6 +66,7 @@ struct MultiQuinticSpline1D
{
typedef T value_type;
typedef OneDimGridBase<T> grid_type;
typedef Matrix<value_type, aligned_allocator<value_type> > coeff_type;
///number of splines
size_t num_splines;
@ -40,11 +75,15 @@ struct MultiQuinticSpline1D
///default is false
bool own_spline;
/** shared_ptr
* Coeffs[6*spline_points][num_splines]
*/
///will be the real grid
LogGridLight<T> myGrid;
grid_type* m_grid;
Matrix<value_type, aligned_allocator<value_type> > *Coeffs;
/** Need to use shared_ptr
*
* Coeffs[6*spline_points][num_splines+padding]
*/
coeff_type* Coeffs;
aligned_vector<value_type> first_deriv;
MultiQuinticSpline1D():own_spline(false),m_grid(nullptr),Coeffs(nullptr){}
@ -54,7 +93,9 @@ struct MultiQuinticSpline1D
MultiQuinticSpline1D<T>* makeClone() const
{
MultiQuinticSpline1D<T>* myclone= new MultiQuinticSpline1D<T>(*this);
#if !defined(USE_MULTIQUINTIC)
myclone->m_grid=m_grid->makeClone(); //wasting X, to remove this
#endif
myclone->own_spline=false;
return myclone;
}
@ -66,6 +107,20 @@ struct MultiQuinticSpline1D
inline void evaluate(T r, T* restrict u)
{
#if defined(USE_MULTIQUINTIC)
if(r<myGrid.lower_bound)
{
const value_type dr=r-myGrid.lower_bound;
const value_type* restrict a=(*Coeffs)[0];
for(size_t i=0; i<num_splines; ++i)
u[i]=a[i]+first_deriv[i]*dr;
}
else
{
int loc;
const auto cL=myGrid.getCLForQuintic(r,loc);
const size_t offset=loc*6;
#else
if(r<m_grid->rmin())
{
const value_type dr=r-m_grid->rmin();
@ -75,16 +130,16 @@ struct MultiQuinticSpline1D
}
else
{
//m_grid->locate(r);
m_grid->updateForQuintic(r,false);
const size_t offset=m_grid->Loc*6;
int loc;
const auto cL=myGrid.getCLForQuintic(r,loc);
const size_t offset=loc*6;
#endif
const value_type* restrict a=(*Coeffs)[offset+0];
const value_type* restrict b=(*Coeffs)[offset+1];
const value_type* restrict c=(*Coeffs)[offset+2];
const value_type* restrict d=(*Coeffs)[offset+3];
const value_type* restrict e=(*Coeffs)[offset+4];
const value_type* restrict f=(*Coeffs)[offset+5];
const auto cL=m_grid->cL;
for(size_t i=0; i<num_splines; ++i)
u[i]=a[i]+cL*(b[i]+cL*(c[i]+cL*(d[i]+cL*(e[i]+cL*f[i]))));
}
@ -93,6 +148,24 @@ struct MultiQuinticSpline1D
inline void evaluate(T r, T* restrict u, T* restrict du, T* restrict d2u)
{
constexpr value_type czero(0);
#if defined(USE_MULTIQUINTIC)
if(r<myGrid.lower_bound)
{
const value_type dr=r-myGrid.lower_bound;
const value_type* restrict a=(*Coeffs)[0];
for(size_t i=0; i<num_splines; ++i)
{
u[i]=a[i]+first_deriv[i]*dr;
du[i]=first_deriv[i]*dr;
d2u[i]=czero;
}
}
else
{
int loc;
const auto cL=myGrid.getCLForQuintic(r,loc);
const size_t offset=loc*6;
#else
if(r<m_grid->rmin())
{
const T dr=r-m_grid->rmin();
@ -107,6 +180,9 @@ struct MultiQuinticSpline1D
else
{
m_grid->updateForQuintic(r,true);
const size_t offset=m_grid->Loc*6;
const auto cL=m_grid->cL;
#endif
constexpr value_type ctwo(2);
constexpr value_type cthree(3);
constexpr value_type cfour(4);
@ -115,15 +191,13 @@ struct MultiQuinticSpline1D
constexpr value_type c12(5);
constexpr value_type c20(6);
const size_t offset=m_grid->Loc*6;
const value_type* restrict a=(*Coeffs)[offset+0];
const value_type* restrict b=(*Coeffs)[offset+1];
const value_type* restrict c=(*Coeffs)[offset+2];
const value_type* restrict d=(*Coeffs)[offset+3];
const value_type* restrict e=(*Coeffs)[offset+4];
const value_type* restrict f=(*Coeffs)[offset+5];
const auto cR=m_grid->cR;
const auto cL=m_grid->cL;
for(size_t i=0; i<num_splines; ++i)
{
u[i] = a[i]+cL*(b[i]+cL*(c[i]+cL*(d[i]+cL*(e[i]+cL*f[i]))));
@ -133,6 +207,13 @@ struct MultiQuinticSpline1D
}
}
/** initialize grid */
template<typename T1>
void setGrid(T1 rmin, T1 rmax, int npts)
{
myGrid.set(rmin,rmax,npts);
}
/** initialize grid and container
* @param ri minimum grid point
* @param rf maximum grid point
@ -146,9 +227,9 @@ struct MultiQuinticSpline1D
//this is very shaky
if(Coeffs==nullptr && !own_spline)
{
spline_order=order;
spline_order=order+1;
num_splines=norbs;
Coeffs=new Matrix<value_type, aligned_allocator<value_type> >((order+1)*m_grid->size(),getAlignedSize<T>(norbs));
Coeffs=new coeff_type((order+1)*m_grid->size(),getAlignedSize<T>(norbs));
first_deriv.resize(num_splines);
own_spline=true;
}
@ -157,28 +238,28 @@ struct MultiQuinticSpline1D
template<typename T1>
void add_spline(int ispline, OneDimQuinticSpline<T1>& in)
{
const T1* restrict A=in.m_Y.data();
const T1* restrict B=in.B.data();
const T1* restrict C=in.m_Y2.data();
const T1* restrict D=in.D.data();
const T1* restrict E=in.E.data();
const T1* restrict F=in.F.data();
first_deriv[ispline]=in.first_deriv;
value_type* restrict out=Coeffs->data();
const size_t ncols=Coeffs->cols();
const size_t num_points=in.size();
for(size_t i=0; i<num_points; ++i)
//if(spline_order==QUINTIC)
{
out[(i*6+0)*ncols+ispline]=A[i];
out[(i*6+1)*ncols+ispline]=B[i];
out[(i*6+2)*ncols+ispline]=C[i];
out[(i*6+3)*ncols+ispline]=D[i];
out[(i*6+4)*ncols+ispline]=E[i];
out[(i*6+5)*ncols+ispline]=F[i];
const T1* restrict A=in.m_Y.data();
const T1* restrict B=in.B.data();
const T1* restrict C=in.m_Y2.data();
const T1* restrict D=in.D.data();
const T1* restrict E=in.E.data();
const T1* restrict F=in.F.data();
value_type* restrict out=Coeffs->data();
const size_t ncols=Coeffs->cols();
const size_t num_points=in.size();
for(size_t i=0; i<num_points; ++i)
{
out[(i*6+0)*ncols+ispline]=static_cast<T>(A[i]);
out[(i*6+1)*ncols+ispline]=static_cast<T>(B[i]);
out[(i*6+2)*ncols+ispline]=static_cast<T>(C[i]);
out[(i*6+3)*ncols+ispline]=static_cast<T>(D[i]);
out[(i*6+4)*ncols+ispline]=static_cast<T>(E[i]);
out[(i*6+5)*ncols+ispline]=static_cast<T>(F[i]);
}
}
}
#if 0
template<typename InType, typename VV>

159
src/QMCWaveFunctions/lcao/RadialOrbitalSetBuilder.h
View File

@ -31,26 +31,73 @@
namespace qmcplusplus
{
/** helper classes and function to handle transformations of the radial functions
*/
/** temporary function to compute the cutoff without constructing NGFunctor */
template<typename Fin, typename T>
inline double find_cutoff(Fin& in, T rmax)
{
//myInFunc=in;
LogGrid<T> agrid;
const T delta=0.1;
const T eps=1e-5;
bool too_small=true;
agrid.set(eps,rmax,1001);
int i=1000;
T r=rmax;
while(too_small && i>0)
template<typename Fin, typename T>
inline double find_cutoff(Fin& in, T rmax)
{
//myInFunc=in;
LogGrid<T> agrid;
const T delta=0.1;
const T eps=1e-0;
bool too_small=true;
agrid.set(eps,rmax,1001);
int i=1000;
T r=rmax;
while(too_small && i>0)
{
r=agrid[i--];
T x=in.f(r);
too_small=(std::abs(x)<eps);
}
return static_cast<double>(r);
}
/** abstract class to defer generation of spline functors
* @tparam T precision of the final result
*/
template<typename T>
struct TransformerBase
{
///temporary grid in double precision
typedef OneDimGridBase<double> grid_type;
///the multiple set
typedef MultiQuinticSpline1D<T> FnOut;
/** convert input 1D functor to the multi set
* @param agrid original grid
* @param multiset the object that should be populated
* @param ispline index of the this analytic function
* @param int order quintic (or cubic) only quintic is used
*/
virtual void convert(grid_type* agrid, FnOut* multiset, int ispline, int order)=0;
virtual ~TransformerBase(){}
};
template<typename T, typename FnIn>
struct A2NTransformer: TransformerBase<T>
{
r=agrid[i--];
T x=in.f(r);
too_small=(std::abs(x)<eps);
}
return static_cast<double>(r);
}
typedef typename TransformerBase<T>::grid_type grid_type;
typedef typename TransformerBase<T>::FnOut FnOut;
FnIn* m_ref; //candidate for unique_ptr
A2NTransformer(FnIn* in): m_ref(in){}
~A2NTransformer() { delete m_ref; }
inline double find_cutoff(T rmax) const
{
return find_cutoff(*m_ref,rmax);
}
void convert(grid_type* agrid, FnOut* multiset, int ispline, int order)
{
typedef OneDimQuinticSpline<double> spline_type;
spline_type radorb(agrid);
Transform2GridFunctor<FnIn,spline_type> transform(*m_ref,radorb);
transform.generate(agrid->rmin(),agrid->rmax(),agrid->size());
multiset->add_spline(ispline,radorb);
}
};
/** Build a set of radial orbitals at the origin
@ -109,6 +156,8 @@ public:
void finalize();
private:
OneDimGridBase<double>* myGrid;
//only check the cutoff
void addGaussian(xmlNodePtr cur);
void addSlater(xmlNodePtr cur);
@ -118,16 +167,22 @@ private:
///safe common cutoff radius
double m_rcut_safe;
/** radial functors to be finalized
*/
std::vector<TransformerBase<RealType>*> radTemp;
//std::vector<std::pair<int,OptimizableFunctorBase*> > radFuncTemp;
///store the temporary analytic data
std::vector<GaussianCombo<RealType>*> gtoTemp;
std::vector<SlaterCombo<RealType>*> stoTemp;
std::tuple<int,double,double> grid_param_in;
};
template<typename COT>
RadialOrbitalSetBuilder<COT>::RadialOrbitalSetBuilder(xmlNodePtr cur)
: Normalized(true),m_orbitals(0),input_grid(nullptr),m_rcut(-1.0),
m_infunctype("Gaussian"), m_fileid(-1)
m_infunctype("Gaussian"), m_fileid(-1), myGrid(nullptr)
{
if(cur != NULL)
putCommon(cur);
@ -194,7 +249,16 @@ private:
{
GridType *agrid = OneDimGridFactory::createGrid(cur);
m_orbitals->Grids.push_back(agrid);
#if 0
if(agrid->GridTag == LINEAR_1DGRID)
myGrid=new LinearGrid;
else
myGrid=new LogGrid;
myGrid->set(static_cast<double>(agrid->rmin()),
static_cast<double>(agrid->rmax()),agrid->size());
#endif
}
#if !defined(ENABLE_SOA)
else
{
app_log() << " Grid is created by the input paremters in h5" << std::endl;
@ -243,6 +307,7 @@ private:
//}
//using 0.01 for the time being
}
#endif
//set zero to use std::max
m_rcut_safe=0;
@ -305,12 +370,14 @@ private:
{
int L= m_nlms[1];
#if defined(USE_MULTIQUINTIC)
GaussianCombo<RealType>* gset=new GaussianCombo<RealType>(L,Normalized);
using gto_type=GaussianCombo<double>;
gto_type* gset=new gto_type(L,Normalized);
gset->putBasisGroup(cur);
auto r0=find_cutoff(*gset,m_orbitals->Grids[0]->rmax());
m_rcut_safe=std::max(m_rcut_safe,r0);
radTemp.push_back(new A2NTransformer<RealType,gto_type>(gset));
//add this basisGroup
gtoTemp.push_back(gset);
//gtoTemp.push_back(gset);
m_orbitals->Rnl.push_back(nullptr); //CLEANUP
m_orbitals->RnlID.push_back(m_nlms);
#else
@ -337,19 +404,32 @@ private:
template<typename COT>
void RadialOrbitalSetBuilder<COT>::finalize()
{
std::cout << "Going to finalize " << m_rcut_safe << std::endl;
//std::cout << "Going to finalize " << m_rcut_safe << std::endl;
#if defined(USE_MULTIQUINTIC)
GridType* agrid = m_orbitals->Grids[0];
agrid->locate(static_cast<RealType>(m_rcut_safe));
agrid->set(agrid->rmin(),m_rcut_safe,agrid->Loc);
MultiQuinticSpline1D<RealType>* multiset=new MultiQuinticSpline1D<RealType>;
m_orbitals->setRmax(m_rcut_safe);
int norbs=gtoTemp.size();
MultiQuinticSpline1D<RealType>* multiset=new MultiQuinticSpline1D<RealType>;;
int npts=agrid->size();
if(myGrid==nullptr)
{
myGrid=new LogGrid<double>;
myGrid->set(1.e-6,1.e2,1001);
}
int norbs=radTemp.size();
multiset->initialize(agrid,norbs);
multiset->setGrid(myGrid->rmin(),myGrid->rmax(),myGrid->size());
for(int ib=0; ib<norbs; ++ib)
radTemp[ib]->convert(myGrid,multiset,ib,5);
m_orbitals->MultiRnl=multiset;
m_orbitals->setRmax(static_cast<RealType>(myGrid->rmax()));
#if 0
int npts=agrid->size();
int norbs=gtoTemp.size();
multiset->initialize(agrid,norbs);
if(gtoTemp.size())
{
RadialOrbitalType radorb(agrid);
@ -363,6 +443,7 @@ private:
}
}
m_orbitals->MultiRnl=multiset;
#endif
#else
m_orbitals->setRmax(m_rcut);
#endif
@ -372,6 +453,17 @@ private:
template<typename COT>
void RadialOrbitalSetBuilder<COT>::addSlater(xmlNodePtr cur)
{
#if defined(USE_MULTIQUINTIC)
using sto_type=SlaterCombo<RealType>;
sto_type* gset=new sto_type(m_nlms[1],Normalized);
gset->putBasisGroup(cur);
auto r0=find_cutoff(*gset,m_orbitals->Grids[0]->rmax());
m_rcut_safe=std::max(m_rcut_safe,r0);
radTemp.push_back(new A2NTransformer<RealType,sto_type>(gset));
m_orbitals->RnlID.push_back(m_nlms);
#else
////pointer to the grid
GridType* agrid = m_orbitals->Grids[0];
RadialOrbitalType *radorb = new RadialOrbitalType(agrid);
@ -384,11 +476,15 @@ private:
//add the radial orbital to the list
m_orbitals->Rnl.push_back(radorb);
m_orbitals->RnlID.push_back(m_nlms);
#endif
}
template<typename COT>
void RadialOrbitalSetBuilder<COT>::addNumerical(xmlNodePtr cur, const std::string& dsname)
{
#if defined(USE_MULTIQUINTIC)
APP_ABORT("RadialOrbitalSetBuilder<COT>::addNumerical is not working with multiquintic");
#else
int imin = 0;
OhmmsAttributeSet aAttrib;
aAttrib.add(imin,"imin");
@ -447,6 +543,7 @@ private:
//}
//cout << " Power " << rinv_p << " size=" << rad_orb.size() << std::endl;
//APP_ABORT("NGOBuilder::addNumerical");
#endif
}
}
#endif

26
src/QMCWaveFunctions/lcao/SoaAtomicBasisSet.h
View File

@ -49,6 +49,9 @@ namespace qmcplusplus
///Replace Rnl, Grids
MultiQuinticSpline1D<value_type> *MultiRnl;
#endif
///temporary storage
VectorSoaContainer<value_type,4> tempS;
///set of grids : keep this until completion
std::vector<grid_type*> Grids;
///the constructor
@ -110,6 +113,7 @@ namespace qmcplusplus
inline void setBasisSetSize(int n)
{
BasisSetSize=LM.size();
tempS.resize(std::max(Ylm.size(),RnlID.size()));
}
/** Set Rmax */
@ -151,14 +155,18 @@ namespace qmcplusplus
const T x=-dr[0], y=-dr[1], z=-dr[2];
Ylm.evaluateVGL(x,y,z);
const size_t nl_max=RnlID.size();
T phi[nl_max];
T dphi[nl_max];
T d2phi[nl_max];
//T phi[nl_max];
//T dphi[nl_max];
//T d2phi[nl_max];
//one can assert the alignment
value_type* restrict phi=tempS.data(0);
value_type* restrict dphi=tempS.data(1);
value_type* restrict d2phi=tempS.data(2);
#if defined(USE_MULTIQUINTIC)
MultiRnl->evaluate(r,phi,dphi,d2phi);
#else
const size_t nl_max=RnlID.size();
for(size_t nl=0; nl<nl_max; ++nl)
{
phi[nl]=Rnl[nl]->evaluate(r,dphi[nl],d2phi[nl]);
@ -195,7 +203,7 @@ namespace qmcplusplus
template<typename T, typename PosType>
inline void
evaluateV(const T r, const PosType& dr, T* restrict psi) const
evaluateV(const T r, const PosType& dr, T* restrict psi)
{
if(r>Rmax)
{
@ -203,14 +211,14 @@ namespace qmcplusplus
return;
}
T ylm_v[Ylm.size()];
Ylm.evaluateV(-dr[0],-dr[1],-dr[2],ylm_v);
value_type* restrict ylm_v=tempS.data(0);
value_type* restrict phi_r=tempS.data(1);
const int nl_max=RnlID.size();
T phi_r[nl_max];
Ylm.evaluateV(-dr[0],-dr[1],-dr[2],ylm_v);
#if defined(USE_MULTIQUINTIC)
MultiRnl->evaluate(r,phi_r);
#else
const int nl_max=RnlID.size();
for(size_t nl=0; nl<nl_max; ++nl)
phi_r[nl]=Rnl[nl]->evaluate(r);
#endif