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

Merge branch 'develop' into Manual_LCAO_features

This commit is contained in:
Paul R. C. Kent 2018-08-01 11:53:07 -04:00 committed by GitHub
parent ad6d6a0ef2 86b3d08719
commit b18d0f2424
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
19 changed files with 76 additions and 236 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/Particle/AsymmetricDistanceTableData.h
View File

@ -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

1
src/Particle/DistanceTableAA.cpp
View File

@ -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)

1
src/Particle/DistanceTableAB.cpp
View File

@ -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)

73
src/Particle/DistanceTableData.h
View File

@ -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

7
src/Particle/ParticleSet.cpp
View File

@ -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);

13
src/Particle/ParticleSet.h
View File

@ -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();

2
src/Particle/SoaDistanceTableAA.h
View File

@ -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);

36
src/Particle/SoaDistanceTableBA.h
View File

@ -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

4
src/Particle/SymmetricDistanceTableData.h
View File

@ -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

8
src/QMCDrivers/tests/test_dmc.cpp
View File

@ -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();

4
src/QMCDrivers/tests/test_vmc.cpp
View File

@ -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();

1
src/QMCHamiltonians/CoulombPBCAB.cpp
View File

@ -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;

17
src/QMCHamiltonians/ECPotentialBuilder.cpp
View File

@ -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;
}

35
src/QMCHamiltonians/NonLocalECPotential.cpp
View File

@ -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 jel=0; jel<P.getTotalNum(); jel++)
{
const auto &dist = myTable->Distances[jel];
const auto &displ = myTable->Displacements[jel];
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)
Value += PP[iat]->evaluateOne(P,iat,Psi,J[nj],dist[nj],displ[nj],Tmove,Txy);
}
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
{

5
src/Sandbox/FakeWaveFunction.h
View File

@ -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;

24
src/Sandbox/diff_distancetables.cpp
View File

@ -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;

24
src/Sandbox/diff_j1.cpp
View File

@ -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,30 +276,27 @@ int main(int argc, char** argv)
r_ratio=0.0;
constexpr int nknots=12;
int nsphere=0;
for(int jel=0; jel<nels; ++jel)
{
const auto &dist = d_ie->Distances[jel];
for(int iat=0; iat<nions; ++iat)
if(dist[iat]<Rmax)
{
for(int nj=0, jmax=d_ie->nadj(iat); nj<jmax; ++nj)
{
const auto r=d_ie->distance(iat,nj);
if(r<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;
}
} //end of omp parallel

24
src/Sandbox/diff_wfs.cpp
View File

@ -253,31 +253,27 @@ int main(int argc, char** argv)
r_ratio=0.0;
constexpr int nknots=12;
int nsphere=0;
for(int jel=0; jel<nels; ++jel)
{
const auto &dist = d_ie->Distances[jel];
for(int iat=0; iat<nions; ++iat)
if(dist[iat]<Rmax)
{
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 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;
}
} //end of omp parallel

27
src/Sandbox/miniqmc.cpp
View File

@ -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,24 +285,20 @@ int main(int argc, char** argv)
const DistanceTableData* d_ie=Jastrow->d_ie;
clock_mc.restart();
for(int jel=0; jel<nels; ++jel)
{
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 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 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);
}
}
}