mirror of https://github.com/QMCPACK/qmcpack.git
Merge pull request #939 from ye-luo/interchange-ion-electron-in-NLPP
Remove neighbourlist in distance tables
This commit is contained in:
commit
86b3d08719
|
@ -172,7 +172,7 @@ struct AsymmetricDTD
|
|||
}
|
||||
|
||||
///not so useful inline but who knows
|
||||
inline void evaluate(ParticleSet& P, bool update_neighbor_list)
|
||||
inline void evaluate(ParticleSet& P)
|
||||
{
|
||||
const int ns=N[SourceIndex];
|
||||
const int nt=N[VisitorIndex];
|
||||
|
@ -186,7 +186,7 @@ struct AsymmetricDTD
|
|||
{
|
||||
APP_ABORT(" No need to call AsymmetricDTD::evaluate(ParticleSet& P, int jat)");
|
||||
//based on full evaluation. Only compute it if jat==0
|
||||
if(jat==0) evaluate(P,true);
|
||||
if(jat==0) evaluate(P);
|
||||
}
|
||||
|
||||
///evaluate the temporary pair relations
|
||||
|
|
|
@ -86,7 +86,6 @@ DistanceTableData* createDistanceTable(ParticleSet& s, int dt_type)
|
|||
dt = new SymmetricDTD<RealType,DIM,PPPS>(s,s);
|
||||
}
|
||||
}
|
||||
o << " Setting Rmax = " << s.Lattice.SimulationCellRadius << std::endl;
|
||||
}
|
||||
}
|
||||
else if(sc == SUPERCELL_SLAB)
|
||||
|
|
|
@ -88,7 +88,6 @@ DistanceTableData* createDistanceTable(const ParticleSet& s, ParticleSet& t, int
|
|||
dt = new AsymmetricDTD<RealType,DIM,PPPS>(s,t);
|
||||
}
|
||||
}
|
||||
o << " Setting Rmax = " << s.Lattice.SimulationCellRadius << std::endl;
|
||||
}
|
||||
}
|
||||
else if(sc == SUPERCELL_SLAB)
|
||||
|
|
|
@ -31,43 +31,6 @@
|
|||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** container for the pair data
|
||||
*/
|
||||
template<class T, unsigned D>
|
||||
struct PairDataType
|
||||
{
|
||||
///distance-related data
|
||||
T r, rr, rinv;
|
||||
///displacement vector
|
||||
TinyVector<T,D> dr;
|
||||
///default constructor
|
||||
inline PairDataType() {}
|
||||
///copy constructor
|
||||
inline PairDataType(const PairDataType<T,D>& p):r(p.r),rr(p.rr),rinv(p.rinv),dr(p.dr) {}
|
||||
///copy operator
|
||||
inline PairDataType<T,D>& operator=(const PairDataType<T,D>& p)
|
||||
{
|
||||
r=p.r;
|
||||
rr=p.rr;
|
||||
rinv=p.rinv;
|
||||
dr=p.dr;
|
||||
return *this;
|
||||
}
|
||||
///set the values
|
||||
inline void set(const TinyVector<T,D>& displ, T sep2)
|
||||
{
|
||||
r=sqrt(sep2);
|
||||
rr=sep2;
|
||||
rinv=1.0/r;
|
||||
dr=displ;
|
||||
}
|
||||
///clear the value
|
||||
inline void reset()
|
||||
{
|
||||
r=0.0;
|
||||
}
|
||||
};
|
||||
|
||||
/** @defgroup nnlist Distance-table group
|
||||
* @brief class to manage a set of data for distance relations between ParticleSet objects.
|
||||
*/
|
||||
|
@ -142,12 +105,6 @@ struct DistanceTableData
|
|||
int DTType;
|
||||
///size of indicies
|
||||
TinyVector<IndexType,4> N;
|
||||
///true, if ratio computations need displacement, e.g. LCAO type
|
||||
bool NeedDisplacement;
|
||||
///Maximum radius set by a Hamiltonian
|
||||
RealType Rmax;
|
||||
///** Maximum square */
|
||||
//RealType Rmax2;
|
||||
|
||||
/** @brief M.size() = N[SourceIndex]+1
|
||||
*
|
||||
|
@ -206,11 +163,8 @@ struct DistanceTableData
|
|||
std::string Name;
|
||||
///constructor using source and target ParticleSet
|
||||
DistanceTableData(const ParticleSet& source, const ParticleSet& target)
|
||||
: Origin(&source), N(0), NeedDisplacement(false), Need_full_table_loadWalker(false)
|
||||
{
|
||||
Rmax=0; //set 0
|
||||
//Rmax=source.Lattice.WignerSeitzRadius;
|
||||
}
|
||||
: Origin(&source), N(0), Need_full_table_loadWalker(false)
|
||||
{ }
|
||||
|
||||
///virutal destructor
|
||||
virtual ~DistanceTableData() { }
|
||||
|
@ -225,8 +179,6 @@ struct DistanceTableData
|
|||
{
|
||||
Name = tname;
|
||||
}
|
||||
///set the maximum radius
|
||||
inline void setRmax(RealType rc) { Rmax=std::max(Rmax,rc);}
|
||||
|
||||
///returns the reference the origin particleset
|
||||
const ParticleSet& origin() const
|
||||
|
@ -285,24 +237,24 @@ struct DistanceTableData
|
|||
/// return the distance |R[iadj(i,nj)]-R[i]|
|
||||
inline RealType distance(int i, int nj) const
|
||||
{
|
||||
return (DTType)? r_m2(i,nj): r_m[M[i]+nj];
|
||||
return r_m[M[i]+nj];
|
||||
}
|
||||
|
||||
/// return the displacement R[iadj(i,nj)]-R[i]
|
||||
inline PosType displacement(int i, int nj) const
|
||||
{
|
||||
return (DTType)? dr_m2(i,nj): dr_m[M[i]+nj];
|
||||
return dr_m[M[i]+nj];
|
||||
}
|
||||
|
||||
//!< Returns a number of neighbors of the i-th ptcl.
|
||||
inline IndexType nadj(int i) const
|
||||
{
|
||||
return (DTType)? M[i]:M[i+1]-M[i];
|
||||
return M[i+1]-M[i];
|
||||
}
|
||||
//!< Returns the id of nj-th neighbor for i-th ptcl
|
||||
inline IndexType iadj(int i, int nj) const
|
||||
{
|
||||
return (DTType)? J2(i,nj):J[M[i] +nj];
|
||||
return J[M[i] +nj];
|
||||
}
|
||||
//!< Returns the id of j-th neighbor for i-th ptcl
|
||||
inline IndexType loc(int i, int j) const
|
||||
|
@ -347,11 +299,8 @@ struct DistanceTableData
|
|||
return -1;
|
||||
}
|
||||
|
||||
///update internal data after completing particle-by-particle moves
|
||||
virtual void donePbyP() { }
|
||||
|
||||
///evaluate the full Distance Table and neighbor_list if requested
|
||||
virtual void evaluate(ParticleSet& P, bool update_neighbor_list=true) = 0;
|
||||
///evaluate the full Distance Table
|
||||
virtual void evaluate(ParticleSet& P) = 0;
|
||||
|
||||
/// evaluate the Distance Table for a given electron
|
||||
virtual void evaluate(ParticleSet& P, int jat)=0;
|
||||
|
@ -454,12 +403,6 @@ struct DistanceTableData
|
|||
std::vector<RealType> rinv_m;
|
||||
/** displacement vectors \f$dr(i,j) = R(j)-R(i)\f$ */
|
||||
std::vector<PosType> dr_m;
|
||||
/** full distance AB or AA return r2_m(iat,jat) */
|
||||
Matrix<RealType,aligned_allocator<RealType> > r_m2;
|
||||
/** full displacement AB or AA */
|
||||
Matrix<PosType> dr_m2;
|
||||
/** J2 for compact neighbors */
|
||||
Matrix<int,aligned_allocator<int> > J2;
|
||||
/*@}*/
|
||||
|
||||
/**resize the storage
|
||||
|
|
|
@ -87,10 +87,7 @@ ParticleSet::ParticleSet(const ParticleSet& p)
|
|||
addTable(p.DistTables[i]->origin(),p.DistTables[i]->DTType);
|
||||
}
|
||||
for(int i=0; i<p.DistTables.size(); ++i)
|
||||
{
|
||||
DistTables[i]->Need_full_table_loadWalker = p.DistTables[i]->Need_full_table_loadWalker;
|
||||
DistTables[i]->Rmax = p.DistTables[i]->Rmax;
|
||||
}
|
||||
if(p.SK)
|
||||
{
|
||||
LRBox=p.LRBox; //copy LRBox
|
||||
|
@ -754,8 +751,6 @@ void ParticleSet::rejectMove(Index_t iat)
|
|||
void ParticleSet::donePbyP()
|
||||
{
|
||||
myTimers[2]->start();
|
||||
for (size_t i=0; i<DistTables.size(); i++)
|
||||
DistTables[i]->donePbyP();
|
||||
if (SK && !SK->DoUpdate)
|
||||
SK->UpdateAllPart(*this);
|
||||
activePtcl=-1;
|
||||
|
@ -785,7 +780,7 @@ void ParticleSet::loadWalker(Walker_t& awalker, bool pbyp)
|
|||
// in certain cases, full tables must be ready
|
||||
for (int i=0; i< DistTables.size(); i++)
|
||||
if(DistTables[i]->DTType==DT_AOS||DistTables[i]->Need_full_table_loadWalker)
|
||||
DistTables[i]->evaluate(*this,false);
|
||||
DistTables[i]->evaluate(*this);
|
||||
//computed so that other objects can use them, e.g., kSpaceJastrow
|
||||
if(SK && SK->DoUpdate)
|
||||
SK->UpdateAllPart(*this);
|
||||
|
|
|
@ -278,7 +278,7 @@ public:
|
|||
*
|
||||
* Ensure that the distance for this-this is always created first.
|
||||
*/
|
||||
int addTable(const ParticleSet& psrc, int dt_type);
|
||||
int addTable(const ParticleSet& psrc, int dt_type);
|
||||
|
||||
/** returns index of a distance table, -1 if not present
|
||||
* @param psrc source particle set
|
||||
|
@ -453,15 +453,12 @@ public:
|
|||
*/
|
||||
void saveWalker(Walker_t& awalker);
|
||||
|
||||
/** update neighbor list, structure factor and unmark activePtcl
|
||||
/** update structure factor and unmark activePtcl
|
||||
*
|
||||
* Currently the trial wave function depends only on distances and
|
||||
* doesn't use any neightbor lists from ParticleSet. However, the
|
||||
* evaluation of non-local pseudopotential relies on the neighbor
|
||||
* list of electron-ion and the Coulomb interaction needs the
|
||||
* structure factor. For these reason, donePbyP after the loop of
|
||||
* single electron moves before evaluating the Hamiltonian.
|
||||
* unmark activePtcl is more of a safety measure probably not needed.
|
||||
* The Coulomb interaction evaluation needs the structure factor.
|
||||
* For these reason, call donePbyP after the loop of single
|
||||
* electron moves before evaluating the Hamiltonian. Unmark
|
||||
* activePtcl is more of a safety measure probably not needed.
|
||||
*/
|
||||
void donePbyP();
|
||||
|
||||
|
|
|
@ -61,7 +61,7 @@ struct SoaDistanceTableAA: public DTD_BConds<T,D,SC>, public DistanceTableData
|
|||
Temp_dr.resize(Ntargets);
|
||||
}
|
||||
|
||||
inline void evaluate(ParticleSet& P, bool update_neighbor_list)
|
||||
inline void evaluate(ParticleSet& P)
|
||||
{
|
||||
CONSTEXPR T BigR= std::numeric_limits<T>::max();
|
||||
//P.RSoA.copyIn(P.R);
|
||||
|
|
|
@ -50,12 +50,6 @@ struct SoaDistanceTableBA: public DTD_BConds<T,D,SC>, public DistanceTableData
|
|||
|
||||
Temp_r.resize(Nsources);
|
||||
Temp_dr.resize(Nsources);
|
||||
|
||||
//this is used to build a compact list
|
||||
M.resize(Nsources);
|
||||
r_m2.resize(Nsources,Ntargets_padded);
|
||||
dr_m2.resize(Nsources,Ntargets_padded);
|
||||
J2.resize(Nsources,Ntargets_padded);
|
||||
}
|
||||
|
||||
#if (__cplusplus >= 201103L)
|
||||
|
@ -65,12 +59,11 @@ struct SoaDistanceTableBA: public DTD_BConds<T,D,SC>, public DistanceTableData
|
|||
~SoaDistanceTableBA() {}
|
||||
|
||||
/** evaluate the full table */
|
||||
inline void evaluate(ParticleSet& P, bool update_neighbor_list)
|
||||
inline void evaluate(ParticleSet& P)
|
||||
{
|
||||
//be aware of the sign of Displacement
|
||||
for(int iat=0; iat<Ntargets; ++iat)
|
||||
DTD_BConds<T,D,SC>::computeDistances(P.R[iat],Origin->RSoA, Distances[iat], Displacements[iat], 0, Nsources);
|
||||
if(update_neighbor_list) donePbyP();
|
||||
}
|
||||
|
||||
/** evaluate the iat-row with the current position
|
||||
|
@ -169,33 +162,6 @@ struct SoaDistanceTableBA: public DTD_BConds<T,D,SC>, public DistanceTableData
|
|||
return nn;
|
||||
}
|
||||
|
||||
inline void donePbyP()
|
||||
{
|
||||
//Rmax is zero: no need to transpose the table.
|
||||
if(Rmax<std::numeric_limits<T>::epsilon()) return;
|
||||
|
||||
CONSTEXPR T cminus(-1);
|
||||
for(int iat=0; iat<Nsources; ++iat)
|
||||
{
|
||||
int nn=0;
|
||||
int* restrict jptr=J2[iat];
|
||||
RealType* restrict rptr=r_m2[iat];
|
||||
PosType* restrict dptr=dr_m2[iat];
|
||||
for(int jat=0; jat<Ntargets; ++jat)
|
||||
{
|
||||
const RealType rij=Distances[jat][iat];
|
||||
if(rij<Rmax)
|
||||
{//make the compact list
|
||||
rptr[nn]=rij;
|
||||
dptr[nn]=cminus*Displacements[jat][iat];
|
||||
jptr[nn]=jat;
|
||||
nn++;
|
||||
}
|
||||
}
|
||||
M[iat]=nn;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
}
|
||||
#endif
|
||||
|
|
|
@ -124,7 +124,7 @@ struct SymmetricDTD
|
|||
// }
|
||||
// }
|
||||
|
||||
inline void evaluate(ParticleSet& P, bool update_neighbor_list)
|
||||
inline void evaluate(ParticleSet& P)
|
||||
{
|
||||
const int n = N[SourceIndex];
|
||||
for(int i=0,ij=0; i<n; i++)
|
||||
|
@ -139,7 +139,7 @@ struct SymmetricDTD
|
|||
{
|
||||
APP_ABORT(" No need to call SymmetricDTD::evaluate(ParticleSet& P, int jat)");
|
||||
//based on full evaluation. Only compute it if jat==0
|
||||
if(jat==0) evaluate(P,true);
|
||||
if(jat==0) evaluate(P);
|
||||
}
|
||||
|
||||
///evaluate the temporary pair relations
|
||||
|
|
|
@ -80,7 +80,11 @@ TEST_CASE("DMC Particle-by-Particle advanceWalkers ConstantOrbital", "[drivers][
|
|||
tspecies(massIdx, upIdx) = 1.0;
|
||||
tspecies(massIdx, downIdx) = 1.0;
|
||||
|
||||
#ifdef ENABLE_SOA
|
||||
elec.addTable(ions,DT_SOA);
|
||||
#else
|
||||
elec.addTable(ions,DT_AOS);
|
||||
#endif
|
||||
elec.update();
|
||||
|
||||
|
||||
|
@ -181,7 +185,11 @@ TEST_CASE("DMC Particle-by-Particle advanceWalkers LinearOrbital", "[drivers][dm
|
|||
tspecies(massIdx, upIdx) = 1.0;
|
||||
tspecies(massIdx, downIdx) = 1.0;
|
||||
|
||||
#ifdef ENABLE_SOA
|
||||
elec.addTable(ions,DT_SOA);
|
||||
#else
|
||||
elec.addTable(ions,DT_AOS);
|
||||
#endif
|
||||
elec.update();
|
||||
|
||||
|
||||
|
|
|
@ -80,7 +80,11 @@ TEST_CASE("VMC Particle-by-Particle advanceWalkers", "[drivers][vmc]")
|
|||
tspecies(massIdx, upIdx) = 1.0;
|
||||
tspecies(massIdx, downIdx) = 1.0;
|
||||
|
||||
#ifdef ENABLE_SOA
|
||||
elec.addTable(ions,DT_SOA);
|
||||
#else
|
||||
elec.addTable(ions,DT_AOS);
|
||||
#endif
|
||||
elec.update();
|
||||
|
||||
|
||||
|
|
|
@ -37,7 +37,6 @@ CoulombPBCAB::CoulombPBCAB(ParticleSet& ions, ParticleSet& elns,
|
|||
//AB = new LRHandlerType(ions);
|
||||
myTableIndex=elns.addTable(ions,DT_SOA_PREFERRED);
|
||||
initBreakup(elns);
|
||||
elns.DistTables[myTableIndex]->setRmax(myRcut);
|
||||
prefix="Flocal";
|
||||
app_log() << " Rcut " << myRcut << std::endl;
|
||||
app_log() << " Maximum K shell " << AB->MaxKshell << std::endl;
|
||||
|
|
|
@ -73,8 +73,6 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
|
|||
useSimpleTableFormat();
|
||||
}
|
||||
|
||||
RealType rcore_max=0;
|
||||
|
||||
///create LocalECPotential
|
||||
bool usePBC =
|
||||
!(IonConfig.Lattice.SuperCellEnum == SUPERCELL_OPEN || pbc =="no");
|
||||
|
@ -91,13 +89,7 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
|
|||
LocalECPotential* apot = new LocalECPotential(IonConfig,targetPtcl);
|
||||
#endif
|
||||
for(int i=0; i<localPot.size(); i++)
|
||||
{
|
||||
if(localPot[i])
|
||||
{
|
||||
apot->add(i,localPot[i],localZeff[i]);
|
||||
rcore_max=std::max(rcore_max,localPot[i]->r_max);
|
||||
}
|
||||
}
|
||||
if(localPot[i]) apot->add(i,localPot[i],localZeff[i]);
|
||||
targetH.addOperator(apot,"LocalECP");
|
||||
}
|
||||
else
|
||||
|
@ -145,18 +137,11 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
|
|||
<< nknot_max << std::endl;
|
||||
if(NLPP_algo=="batched") app_log() << " Using batched ratio computing in NonLocalECP" << std::endl;
|
||||
|
||||
rcore_max=std::max(rc2,rcore_max);
|
||||
targetPtcl.checkBoundBox(2*rc2);
|
||||
|
||||
targetH.addOperator(apot,"NonLocalECP");
|
||||
}
|
||||
|
||||
//app_log() << "Checking THIS " << std::endl;
|
||||
//DEV::OPTIMIZE_SOA
|
||||
int tid=targetPtcl.addTable(IonConfig,DT_SOA_PREFERRED);
|
||||
targetPtcl.DistTables[tid]->setRmax(rcore_max); //For the optimization only
|
||||
app_log() << " ECPotential::Rmax " << targetPtcl.DistTables[tid]->Rmax << std::endl;
|
||||
|
||||
app_log().flush();
|
||||
return true;
|
||||
}
|
||||
|
|
|
@ -153,20 +153,16 @@ NonLocalECPotential::evaluate(ParticleSet& P, bool Tmove)
|
|||
}
|
||||
else
|
||||
{
|
||||
const DistanceTableData* myTable = P.DistTables[myTableIndex];
|
||||
const auto myTable = P.DistTables[myTableIndex];
|
||||
if(myTable->DTType == DT_SOA)
|
||||
{
|
||||
for(int iat=0; iat<NumIons; iat++)
|
||||
for(int jel=0; jel<P.getTotalNum(); jel++)
|
||||
{
|
||||
if(PP[iat]==nullptr) continue;
|
||||
const int* restrict J=myTable->J2[iat];
|
||||
const RealType* restrict dist=myTable->r_m2[iat];
|
||||
const PosType* restrict displ=myTable->dr_m2[iat];
|
||||
for(size_t nj=0; nj<myTable->M[iat]; ++nj)
|
||||
{
|
||||
if(dist[nj]<PP[iat]->Rmax)
|
||||
Value += PP[iat]->evaluateOne(P,iat,Psi,J[nj],dist[nj],displ[nj],Tmove,Txy);
|
||||
}
|
||||
const auto &dist = myTable->Distances[jel];
|
||||
const auto &displ = myTable->Displacements[jel];
|
||||
for(int iat=0; iat<NumIons; iat++)
|
||||
if(PP[iat]!=nullptr && dist[iat]<PP[iat]->Rmax)
|
||||
Value += PP[iat]->evaluateOne(P,iat,Psi,jel,dist[iat],RealType(-1)*displ[iat],Tmove,Txy);
|
||||
}
|
||||
}
|
||||
else
|
||||
|
@ -212,21 +208,14 @@ void
|
|||
NonLocalECPotential::computeOneElectronTxy(ParticleSet& P, const int ref_elec)
|
||||
{
|
||||
std::vector<NonLocalData>& Txy(nonLocalOps.Txy);
|
||||
const DistanceTableData* myTable = P.DistTables[myTableIndex];
|
||||
const auto myTable = P.DistTables[myTableIndex];
|
||||
if(myTable->DTType == DT_SOA)
|
||||
{
|
||||
const auto &dist = myTable->Distances[ref_elec];
|
||||
const auto &displ = myTable->Displacements[ref_elec];
|
||||
for(int iat=0; iat<NumIons; iat++)
|
||||
{
|
||||
if(PP[iat]==nullptr) continue;
|
||||
const int* restrict J=myTable->J2[iat];
|
||||
const RealType* restrict dist=myTable->r_m2[iat];
|
||||
const PosType* restrict displ=myTable->dr_m2[iat];
|
||||
for(size_t nj=0; nj<myTable->M[iat]; ++nj)
|
||||
{
|
||||
if(dist[nj]<PP[iat]->Rmax && J[nj]==ref_elec)
|
||||
PP[iat]->evaluateOne(P,iat,Psi,J[nj],dist[nj],displ[nj],true,Txy);
|
||||
}
|
||||
}
|
||||
if(PP[iat]!=nullptr && dist[iat]<PP[iat]->Rmax)
|
||||
PP[iat]->evaluateOne(P,iat,Psi,ref_elec,dist[iat],RealType(-1)*displ[iat],true,Txy);
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
|
@ -29,11 +29,6 @@ namespace qmcplusplus
|
|||
DistanceTableData* d_ee;
|
||||
DistanceTableData* d_ie;
|
||||
|
||||
inline void setRmax(valT x)
|
||||
{
|
||||
d_ie->setRmax(x);
|
||||
}
|
||||
|
||||
virtual ~FakeWaveFunctionBase(){}
|
||||
virtual void evaluateLog(ParticleSet& P)=0;
|
||||
virtual posT evalGrad(ParticleSet& P, int iat)=0;
|
||||
|
|
|
@ -120,9 +120,7 @@ int main(int argc, char** argv)
|
|||
DistanceTableData* d_ee=DistanceTable::add(els,DT_SOA);
|
||||
DistanceTableData* d_ie=DistanceTable::add(ions,els,DT_SOA);
|
||||
|
||||
d_ie->setRmax(els.Lattice.WignerSeitzRadius);
|
||||
RealType Rsim=els.Lattice.WignerSeitzRadius;
|
||||
std::cout << "Setting 1-body cutoff " << d_ie->Rmax << std::endl;
|
||||
|
||||
DistanceTableData* d_ee_aos=DistanceTable::add(els_aos,DT_AOS);
|
||||
DistanceTableData* d_ie_aos=DistanceTable::add(ions_aos,els_aos,DT_AOS);
|
||||
|
@ -191,7 +189,6 @@ int main(int argc, char** argv)
|
|||
els.donePbyP();
|
||||
els_aos.donePbyP();
|
||||
|
||||
std::cout << "WignerSeitzRadius " << d_ie->Rmax << std::endl;
|
||||
{
|
||||
double r_err=0.0;
|
||||
for(int iat=0; iat<nels; ++iat)
|
||||
|
@ -275,27 +272,6 @@ int main(int argc, char** argv)
|
|||
cout << "AB SoA-AoS in displacement = " << disp_err/nn << endl;
|
||||
}
|
||||
|
||||
|
||||
{
|
||||
double dist_err=0.0;
|
||||
double displ_err=0.0;
|
||||
int nn=0;
|
||||
for(int i=0; i<nions; ++i)
|
||||
for(int nj=0, jmax=d_ie->M[i]; nj<jmax; ++nj)
|
||||
{
|
||||
int j=d_ie->J2(i,nj);
|
||||
int ij=i*nels+j;
|
||||
PosType diff_dr=d_ie->dr_m2(i,nj)-d_ie_aos->dr(ij);
|
||||
RealType diff_r=d_ie->r_m2(i,nj)- d_ie_aos->r(ij);
|
||||
dist_err += abs(diff_r);
|
||||
displ_err += sqrt(dot(diff_dr,diff_dr));
|
||||
nn++;
|
||||
}
|
||||
cout << "---------------------------------" << endl;
|
||||
cout << "AB SoA-AoS compact = " << dist_err/nn << " " << displ_err/nn;
|
||||
cout << " Total number of nn = " << nn << " / " << nels*nions << endl;
|
||||
}
|
||||
|
||||
OHMMS::Controller->finalize();
|
||||
|
||||
return 0;
|
||||
|
|
|
@ -141,7 +141,6 @@ int main(int argc, char** argv)
|
|||
OneBodyJastrowOrbital<BsplineFunctor<RealType> > J_aos(ions,els_aos);
|
||||
|
||||
DistanceTableData* d_ie=DistanceTable::add(ions,els,DT_SOA);
|
||||
d_ie->setRmax(Rmax);
|
||||
|
||||
RealType r1_cut=std::min(RealType(6.4),els.Lattice.WignerSeitzRadius);
|
||||
|
||||
|
@ -277,29 +276,26 @@ int main(int argc, char** argv)
|
|||
r_ratio=0.0;
|
||||
constexpr int nknots=12;
|
||||
int nsphere=0;
|
||||
for(int iat=0; iat<nions; ++iat)
|
||||
for(int jel=0; jel<nels; ++jel)
|
||||
{
|
||||
for(int nj=0, jmax=d_ie->nadj(iat); nj<jmax; ++nj)
|
||||
{
|
||||
const auto r=d_ie->distance(iat,nj);
|
||||
if(r<Rmax)
|
||||
const auto &dist = d_ie->Distances[jel];
|
||||
for(int iat=0; iat<nions; ++iat)
|
||||
if(dist[iat]<Rmax)
|
||||
{
|
||||
const int iel=d_ie->iadj(iat,nj);
|
||||
nsphere++;
|
||||
random_th.generate_uniform(&delta[0][0],nknots*3);
|
||||
for(int k=0; k<nknots;++k)
|
||||
{
|
||||
els.makeMoveOnSphere(iel,delta[k]);
|
||||
RealType r_soa=J.ratio(els,iel);
|
||||
els.rejectMove(iel);
|
||||
els.makeMoveOnSphere(jel,delta[k]);
|
||||
RealType r_soa=J.ratio(els,jel);
|
||||
els.rejectMove(jel);
|
||||
|
||||
els_aos.makeMoveOnSphere(iel,delta[k]);
|
||||
RealType r_aos=J_aos.ratio(els_aos,iel);
|
||||
els_aos.rejectMove(iel);
|
||||
els_aos.makeMoveOnSphere(jel,delta[k]);
|
||||
RealType r_aos=J_aos.ratio(els_aos,jel);
|
||||
els_aos.rejectMove(jel);
|
||||
r_ratio += abs(r_soa/r_aos-1);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
cout << "ratio with SphereMove Error = " << r_ratio/nsphere << " # of moves =" << nsphere << endl;
|
||||
}
|
||||
|
|
|
@ -253,30 +253,26 @@ int main(int argc, char** argv)
|
|||
r_ratio=0.0;
|
||||
constexpr int nknots=12;
|
||||
int nsphere=0;
|
||||
for(int iat=0; iat<nions; ++iat)
|
||||
for(int jel=0; jel<nels; ++jel)
|
||||
{
|
||||
const auto centerP=ions.R[iat];
|
||||
for(int nj=0, jmax=d_ie->nadj(iat); nj<jmax; ++nj)
|
||||
{
|
||||
const auto r=d_ie->distance(iat,nj);
|
||||
if(r<Rmax)
|
||||
const auto &dist = d_ie->Distances[jel];
|
||||
for(int iat=0; iat<nions; ++iat)
|
||||
if(dist[iat]<Rmax)
|
||||
{
|
||||
const int iel=d_ie->iadj(iat,nj);
|
||||
nsphere++;
|
||||
random_th.generate_uniform(&delta[0][0],nknots*3);
|
||||
for(int k=0; k<nknots;++k)
|
||||
{
|
||||
els.makeMoveOnSphere(iel,delta[k]);
|
||||
RealType r_soa=J.ratio(els,iel);
|
||||
els.rejectMove(iel);
|
||||
els.makeMoveOnSphere(jel,delta[k]);
|
||||
RealType r_soa=J.ratio(els,jel);
|
||||
els.rejectMove(jel);
|
||||
|
||||
els_aos.makeMoveOnSphere(iel,delta[k]);
|
||||
RealType r_aos=J_aos.ratio(els_aos,iel);
|
||||
els_aos.rejectMove(iel);
|
||||
els_aos.makeMoveOnSphere(jel,delta[k]);
|
||||
RealType r_aos=J_aos.ratio(els_aos,jel);
|
||||
els_aos.rejectMove(jel);
|
||||
r_ratio += abs(r_soa/r_aos-1);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
cout << "ratio with SphereMove Error = " << r_ratio/nsphere << " # of moves =" << nsphere << endl;
|
||||
}
|
||||
|
|
|
@ -207,9 +207,6 @@ int main(int argc, char** argv)
|
|||
else
|
||||
Jastrow=new AoSWaveFunction(ions,els);
|
||||
|
||||
//set Rmax for ion-el distance table for PP
|
||||
Jastrow->setRmax(Rmax);
|
||||
|
||||
//create pseudopp
|
||||
NonLocalPP<OHMMS_PRECISION> ecp(random_th);
|
||||
|
||||
|
@ -288,26 +285,22 @@ int main(int argc, char** argv)
|
|||
const DistanceTableData* d_ie=Jastrow->d_ie;
|
||||
|
||||
clock_mc.restart();
|
||||
for(int iat=0; iat<nions; ++iat)
|
||||
for(int jel=0; jel<nels; ++jel)
|
||||
{
|
||||
const auto centerP=ions.R[iat];
|
||||
for(int nj=0, jmax=d_ie->nadj(iat); nj<jmax; ++nj)
|
||||
{
|
||||
const auto r=d_ie->distance(iat,nj);
|
||||
if(r<Rmax)
|
||||
const auto &dist = d_ie->Distances[jel];
|
||||
const auto &displ = d_ie->Displacements[jel];
|
||||
for(int iat=0; iat<nions; ++iat)
|
||||
if(dist[iat]<Rmax)
|
||||
{
|
||||
const int iel=d_ie->iadj(iat,nj);
|
||||
const auto dr=d_ie->displacement(iat,nj);
|
||||
for (int k=0; k < nknots ; k++)
|
||||
{
|
||||
PosType deltar(r*rOnSphere[k]-dr);
|
||||
els.makeMoveOnSphere(iel,deltar);
|
||||
spo.evaluate_v(els.R[iel]);
|
||||
Jastrow->ratio(els,iel);
|
||||
els.rejectMove(iel);
|
||||
PosType deltar(dist[iat]*rOnSphere[k]-displ[iat]);
|
||||
els.makeMoveOnSphere(jel,deltar);
|
||||
spo.evaluate_v(els.R[jel]);
|
||||
Jastrow->ratio(els,jel);
|
||||
els.rejectMove(jel);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
t_pseudo_loc+=clock_mc.elapsed();
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue