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

Complete Pade for SoA. Some tests still fail.

This commit is contained in:
Ye Luo 2017-10-03 09:37:14 -05:00
parent d7f48f6b9b
commit 97356367dc
6 changed files with 41 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/QMCWaveFunctions/Jastrow/BsplineFunctor.h
View File

@ -121,6 +121,7 @@ struct BsplineFunctor: public OptimizableFunctorBase
}
/** compute value, gradient and laplacian for [iStart, iEnd) pairs
* @param iat dummy
* @param iStart starting particle index
* @param iEnd ending particle index
* @param _distArray distance arrUay
@ -130,7 +131,7 @@ struct BsplineFunctor: public OptimizableFunctorBase
* @param distArrayCompressed temp storage to filter r_j < cutoff_radius
* @param distIndices temp storage for the compressed index
*/
void evaluateVGL(const int iStart, const int iEnd,
void evaluateVGL(const int iat, const int iStart, const int iEnd,
const T* _distArray,
T* restrict _valArray,
T* restrict _gradArray,
@ -138,13 +139,14 @@ struct BsplineFunctor: public OptimizableFunctorBase
T* restrict distArrayCompressed, int* restrict distIndices ) const;
/** evaluate sum of the pair potentials for [iStart,iEnd)
* @param iat dummy
* @param iStart starting particle index
* @param iEnd ending particle index
* @param _distArray distance arrUay
* @param distArrayCompressed temp storage to filter r_j < cutoff_radius
* @return \f$\sum u(r_j)\f$ for r_j < cutoff_radius
*/
T evaluateV(const int iStart, const int iEnd,
T evaluateV(const int iat, const int iStart, const int iEnd,
const T* restrict _distArray,
T* restrict distArrayCompressed) const;
@ -644,7 +646,7 @@ struct BsplineFunctor: public OptimizableFunctorBase
template<typename T>
inline T
BsplineFunctor<T>::evaluateV(const int iStart, const int iEnd,
BsplineFunctor<T>::evaluateV(const int iat, const int iStart, const int iEnd,
const T* restrict _distArray, T* restrict distArrayCompressed ) const
{
const real_type* restrict distArray = _distArray + iStart;
@ -681,7 +683,7 @@ BsplineFunctor<T>::evaluateV(const int iStart, const int iEnd,
}
template<typename T>
inline void BsplineFunctor<T>::evaluateVGL(const int iStart, const int iEnd,
inline void BsplineFunctor<T>::evaluateVGL(const int iat, const int iStart, const int iEnd,
const T* _distArray, T* restrict _valArray,
T* restrict _gradArray, T* restrict _laplArray,
T* restrict distArrayCompressed, int* restrict distIndices ) const

4
src/QMCWaveFunctions/Jastrow/J1OrbitalSoA.h
View File

@ -141,7 +141,7 @@ struct J1OrbitalSoA : public OrbitalBase
for(int jg=0; jg<NumGroups; ++jg)
{
if(F[jg]!=nullptr)
curAt += F[jg]->evaluateV(Ions.first(jg), Ions.last(jg), dist, DistCompressed.data() );
curAt += F[jg]->evaluateV(-1, Ions.first(jg), Ions.last(jg), dist, DistCompressed.data());
}
}
else
@ -209,7 +209,7 @@ struct J1OrbitalSoA : public OrbitalBase
for(int jg=0; jg<NumGroups; ++jg)
{
if(F[jg]==nullptr) continue;
F[jg]->evaluateVGL( Ions.first(jg), Ions.last(jg), dist,
F[jg]->evaluateVGL(-1, Ions.first(jg), Ions.last(jg), dist,
U.data(), dU.data(), d2U.data(), DistCompressed.data(), DistIndice.data());
}
}

4
src/QMCWaveFunctions/Jastrow/J2OrbitalSoA.h
View File

@ -414,7 +414,7 @@ J2OrbitalSoA<FT>::computeU3(ParticleSet& P, int iat, const RealType* restrict di
const FuncType& f2(*F[igt+jg]);
int iStart = P.first(jg);
int iEnd = P.last(jg);
f2.evaluateVGL(iStart, iEnd, dist, u, du, d2u, DistCompressed.data(), DistIndice.data());
f2.evaluateVGL(iat, iStart, iEnd, dist, u, du, d2u, DistCompressed.data(), DistIndice.data());
}
//u[iat]=czero;
//du[iat]=czero;
@ -437,7 +437,7 @@ J2OrbitalSoA<FT>::ratio(ParticleSet& P, int iat)
const FuncType& f2(*F[igt+jg]);
int iStart = P.first(jg);
int iEnd = P.last(jg);
cur_Uat += f2.evaluateV( iStart, iEnd, dist, DistCompressed.data() );
cur_Uat += f2.evaluateV(iat, iStart, iEnd, dist, DistCompressed.data());
}
return std::exp(Uat[iat]-cur_Uat);

25
src/QMCWaveFunctions/Jastrow/PadeFunctors.h
View File

@ -97,13 +97,13 @@ struct PadeFunctor:public OptimizableFunctorBase
AoverB=A/B;
}
inline real_type evaluate(real_type r)
inline real_type evaluate(real_type r) const
{
return A*r/(1.0+B*r)-AoverB;
}
inline real_type
evaluate(real_type r, real_type& dudr, real_type& d2udr2)
evaluate(real_type r, real_type& dudr, real_type& d2udr2) const
{
real_type u = 1.0/(1.0+B*r);
dudr = A*u*u;
@ -112,7 +112,7 @@ struct PadeFunctor:public OptimizableFunctorBase
}
inline real_type
evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3)
evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3) const
{
real_type u = 1.0/(1.0+B*r);
dudr = A*u*u;
@ -121,6 +121,25 @@ struct PadeFunctor:public OptimizableFunctorBase
return A*u*r-AoverB;
}
inline real_type evaluateV(const int iat, const int iStart, const int iEnd,
const T* restrict _distArray, T* restrict distArrayCompressed ) const
{
real_type sum(0);
for(int idx=iStart; idx<iEnd; idx++)
if (idx!=iat) sum += evaluate(_distArray[idx]);
return sum;
}
inline void evaluateVGL(const int iat, const int iStart, const int iEnd,
const T* distArray, T* restrict valArray,
T* restrict gradArray, T* restrict laplArray,
T* restrict distArrayCompressed, int* restrict distIndices ) const
{
for(int idx=iStart; idx<iEnd; idx++)
valArray[idx] = evaluate(distArray[idx], gradArray[idx], laplArray[idx]);
if ( iat>=iStart && iat<iEnd )
valArray[iat] = gradArray[iat] = laplArray[iat] = T(0);
}
inline real_type f(real_type r)
{

5
src/QMCWaveFunctions/Jastrow/PadeJastrowBuilder.cpp
View File

@ -14,6 +14,7 @@
#include "QMCWaveFunctions/Jastrow/PadeJastrowBuilder.h"
#include "QMCWaveFunctions/Jastrow/TwoBodyJastrowOrbital.h"
#include "QMCWaveFunctions/Jastrow/J2OrbitalSoA.h"
#include "QMCWaveFunctions/Jastrow/OneBodyJastrowOrbital.h"
#include "QMCWaveFunctions/DiffOrbitalBase.h"
#include "Utilities/IteratorUtility.h"
@ -53,7 +54,11 @@ bool PadeJastrowBuilder::put(xmlNodePtr cur)
if(sourceOpt == targetPtcl.getName())
{
//two-body
#if defined(ENABLE_SOA)
typedef J2OrbitalSoA<RadFuncType> J2Type;
#else
typedef TwoBodyJastrowOrbital<RadFuncType> J2Type;
#endif
typedef DiffTwoBodyJastrowOrbital<RadFuncType> dJ2Type;
int taskid=(targetPsi.is_manager())?targetPsi.getGroupID():-1;
J2Type *J2 = new J2Type(targetPtcl,taskid);

6
src/QMCWaveFunctions/tests/test_wf.cpp
View File

@ -47,7 +47,6 @@ TEST_CASE("Pade functor", "[wavefunction]")
REQUIRE(u == Approx(2.232142857142857));
}
#ifndef ENABLE_SOA
TEST_CASE("Pade Jastrow", "[wavefunction]")
{
@ -82,7 +81,11 @@ TEST_CASE("Pade Jastrow", "[wavefunction]")
tspecies(chargeIdx, upIdx) = -1;
tspecies(chargeIdx, downIdx) = -1;
#ifdef ENABLE_SOA
elec_.addTable(ions_,DT_SOA);
#else
elec_.addTable(ions_,DT_AOS);
#endif
elec_.update();
@ -118,6 +121,5 @@ const char *particles = \
REQUIRE(logpsi == Approx(-1.862821769493147));
}
#endif
}