Add a SPO set scanner which allows to print the values of SPOs into files on the selected path.

Usage:
add the following block at the end of determinantset block.
<spo_scanner source="i">
  <path name="a" from_atom="0" to_atom="1" nknots="41"/>
  <path name="b" from_pos="0.0 0.0 0.0" to_pos="6.74632230 6.74632230 6.74632230" nknots="61"/>
</spo_scanner>



git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@7390 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Ye Luo 2017-01-04 02:01:40 +00:00
parent 28be5d29d8
commit 03375a745f
8 changed files with 226 additions and 32 deletions

View File

@ -17,6 +17,7 @@
#include "QMCWaveFunctions/BasisSetFactory.h"
#include "QMCWaveFunctions/SPOSetScanner.h"
#include "QMCWaveFunctions/Fermion/SlaterDetBuilder.h"
#include "Utilities/ProgressReportEngine.h"
#include "OhmmsData/AttributeSet.h"
@ -97,7 +98,7 @@ bool SlaterDetBuilder::put(xmlNodePtr cur)
if (myBasisSetFactory == 0)
{//always create one, using singleton and just to access the member functions
myBasisSetFactory = new BasisSetFactory(targetPtcl,targetPsi, ptclPool);
myBasisSetFactory = new BasisSetFactory(targetPtcl, targetPsi, ptclPool);
myBasisSetFactory->setReportLevel(ReportLevel);
}
@ -224,7 +225,15 @@ bool SlaterDetBuilder::put(xmlNodePtr cur)
while (cur != NULL)//check the basis set
{
getNodeName(cname,cur);
if (cname == sd_tag)
if (cname == sposcanner_tag)
{
if(myComm->rank()==0)
{
SPOSetScanner ascanner(spomap, targetPtcl, ptclPool);
ascanner.put(cur);
}
}
else if (cname == sd_tag)
{
multiDet=false;
if(slaterdet_0)

View File

@ -59,6 +59,8 @@ std::string OrbitalBuilderBase::backflow_tag="backflow";
std::string OrbitalBuilderBase::multisd_tag="multideterminant";
std::string OrbitalBuilderBase::sposcanner_tag="spo_scanner";
OrbitalBuilderBase::OrbitalBuilderBase(ParticleSet& p, TrialWaveFunction& psi):
MPIObjectBase(psi.getCommunicator()),
targetPtcl(p), targetPsi(psi), myNode(NULL)

View File

@ -82,6 +82,8 @@ public:
static std::string backflow_tag;
/// the element name for a multi slater determinant wavefunction
static std::string multisd_tag;
/// the element name for a SPO scanner
static std::string sposcanner_tag;
//@}

View File

@ -21,6 +21,7 @@
#include "QMCWaveFunctions/PlaneWave/PWParameterSet.h"
#include "QMCWaveFunctions/Fermion/SlaterDet.h"
#include "QMCWaveFunctions/DummyBasisSet.h"
#include "QMCWaveFunctions/SPOSetScanner.h"
#include "OhmmsData/ParameterSet.h"
#include "OhmmsData/AttributeSet.h"
#include "Numerics/HDFSTLAttrib.h"
@ -30,8 +31,8 @@
namespace qmcplusplus
{
PWOrbitalBuilder::PWOrbitalBuilder(ParticleSet& els, TrialWaveFunction& psi)
: OrbitalBuilderBase(els,psi),hfileID(-1), rootNode(NULL)
PWOrbitalBuilder::PWOrbitalBuilder(ParticleSet& els, TrialWaveFunction& psi, PtclPoolType& psets)
: OrbitalBuilderBase(els,psi), ptclPool(psets), hfileID(-1), rootNode(NULL)
#if !defined(EANBLE_SMARTPOINTER)
,myBasisSet(0)
#endif
@ -73,27 +74,30 @@ bool PWOrbitalBuilder::put(xmlNodePtr cur)
if(aptr != NULL)
myParam->Ecut=atof((const char*)aptr);
}
else
if(cname == "coefficients")
else if(cname == "coefficients")
{
//close
if(hfileID>0)
H5Fclose(hfileID);
hfileID=getH5(cur,"hdata");
}
else if(cname == sd_tag)
{
if(hfileID<0)
hfileID=getH5(cur,"href");
if(hfileID<0)
{
//close
if(hfileID>0)
H5Fclose(hfileID);
hfileID=getH5(cur,"hdata");
app_error() << " Cannot create a SlaterDet due to missing h5 file" << std::endl;
OHMMS::Controller->abort();
}
else
if(cname == "slaterdeterminant")
{
if(hfileID<0)
hfileID=getH5(cur,"href");
if(hfileID<0)
{
app_error() << " Cannot create a SlaterDet due to missing h5 file" << std::endl;
OHMMS::Controller->abort();
}
success = createPWBasis(cur);
success = putSlaterDet(cur);
}
success = createPWBasis(cur);
success = putSlaterDet(cur);
}
else if (cname == sposcanner_tag)
{
SPOSetScanner ascanner(spomap, targetPtcl, ptclPool);
ascanner.put(cur);
}
cur=cur->next;
}
H5Fclose(hfileID);
@ -108,7 +112,6 @@ bool PWOrbitalBuilder::putSlaterDet(xmlNodePtr cur)
typedef SlaterDet SlaterDeterminant_t;
typedef DiracDeterminantBase Det_t;
SlaterDeterminant_t* sdet(new SlaterDeterminant_t(targetPtcl));
std::map<std::string,SPOSetBasePtr>& spo_ref(sdet->mySPOSet);
int spin_group=0;
cur=cur->children;
while(cur != NULL)
@ -126,20 +129,21 @@ bool PWOrbitalBuilder::putSlaterDet(xmlNodePtr cur)
if(ref == "0")
ref=id;
int firstIndex=targetPtcl.first(spin_group);
std::map<std::string,SPOSetBasePtr>::iterator lit(spo_ref.find(ref));
std::map<std::string,SPOSetBasePtr>::iterator lit(spomap.find(ref));
Det_t* adet=0;
//int spin_group=0;
if(lit == spo_ref.end())
if(lit == spomap.end())
{
app_log() << " Create a PWOrbitalSet" << std::endl;;
SPOSetBasePtr psi(createPW(cur,spin_group));
sdet->add(psi,ref);
adet= new Det_t(psi,firstIndex);
spomap[ref] = psi;
adet = new Det_t(psi,firstIndex);
}
else
{
app_log() << " Reuse a PWOrbitalSet" << std::endl;
adet= new Det_t((*lit).second,firstIndex);
adet = new Det_t((*lit).second,firstIndex);
}
app_log()<< " spin=" << spin_group << " id=" << id << " ref=" << ref << std::endl;
if(adet)

View File

@ -44,6 +44,9 @@ private:
typedef PWRealOrbitalSet::PWBasisPtr PWBasisPtr;
#endif
std::map<std::string,SPOSetBasePtr> spomap;
PtclPoolType& ptclPool;
///Read routine for HDF wavefunction file version 0.10
void ReadHDFWavefunction(hid_t hfile);
@ -62,7 +65,7 @@ private:
public:
///constructor
PWOrbitalBuilder(ParticleSet& els, TrialWaveFunction& wfs);
PWOrbitalBuilder(ParticleSet& els, TrialWaveFunction& wfs, PtclPoolType& psets);
~PWOrbitalBuilder();
///implement vritual function

View File

@ -0,0 +1,174 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
//
// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_SPOSET_SCANNER_H
#define QMCPLUSPLUS_SPOSET_SCANNER_H
#include "Particle/ParticleSet.h"
#include "QMCWaveFunctions/OrbitalSetTraits.h"
#include "QMCWaveFunctions/SPOSetBase.h"
#include "OhmmsData/AttributeSet.h"
namespace qmcplusplus
{
/** a scanner for all the SPO sets.
*/
class SPOSetScanner
{
public:
typedef std::map<std::string,ParticleSet*> PtclPoolType;
typedef std::map<std::string,SPOSetBasePtr> SPOMapType;
typedef QMCTraits::ValueType ValueType;
typedef OrbitalSetTraits<ValueType>::ValueVector_t ValueVector_t;
typedef OrbitalSetTraits<ValueType>::GradVector_t GradVector_t;
typedef OrbitalSetTraits<ValueType>::HessVector_t HessVector_t;
SPOMapType& SPOMap;
PtclPoolType& PtclPool;
ParticleSet* ions;
ParticleSet& target;
// construction/destruction
SPOSetScanner(SPOMapType& spomap, ParticleSet& targetPtcl, PtclPoolType& psets): SPOMap(spomap), target(targetPtcl), PtclPool(psets), ions(0){};
//~SPOSetScanner(){};
// processing scanning
void put(xmlNodePtr cur)
{
app_log() << "Entering the SPO set scanner!" << std::endl;
// check in the source particle set and search for it in the pool.
std::string sourcePtcl("ion0");
OhmmsAttributeSet aAttrib;
aAttrib.add(sourcePtcl,"source");
aAttrib.put(cur);
PtclPoolType::iterator pit(PtclPool.find(sourcePtcl));
if(pit == PtclPool.end())
app_log() << "Source particle set not found. Can not be used as reference point." << std::endl;
else
ions=(*pit).second;
// scanning the SPO sets
xmlNodePtr cur_save=cur;
SPOSetBasePtr mySPOSet;
for (SPOMapType::iterator spo_iter=SPOMap.begin(); spo_iter!=SPOMap.end(); spo_iter++)
{
app_log() << " Processing SPO " << spo_iter->first << std::endl;
mySPOSet = spo_iter->second;
// scanning the paths
cur = cur_save->children;
while (cur != NULL)
{
std::string trace_name("no name");
OhmmsAttributeSet aAttrib;
aAttrib.add(trace_name,"name");
aAttrib.put(cur);
std::string cname;
getNodeName(cname,cur);
std::string prefix(spo_iter->first+"_"+cname+"_"+trace_name);
if (cname == "path")
{
app_log() << " Scanning a " << cname << " called " << trace_name << " and writing to " << prefix+"_v/g/l.dat" << std::endl;
scan_path(cur, mySPOSet, prefix);
}
else
{
if(cname !="text" && cname !="comment") app_log() << " Unknown type of scanning " << cname << std::endl;
}
cur = cur->next;
}
}
app_log() << "Exiting the SPO set scanner!" << std::endl << std::endl;
}
// scanning a path
void scan_path(xmlNodePtr cur, SPOSetBasePtr mySPOSet, std::string prefix)
{
std::string file_name;
file_name=prefix+"_v.dat";
std::ofstream output_v(file_name.c_str());
file_name=prefix+"_g.dat";
std::ofstream output_g(file_name.c_str());
file_name=prefix+"_l.dat";
std::ofstream output_l(file_name.c_str());
int nknots(2);
int from_atom(-1);
int to_atom(-1);
TinyVector<double,OHMMS_DIM> from_pos(0.0, 0.0, 0.0);
TinyVector<double,OHMMS_DIM> to_pos(0.0, 0.0, 0.0);
OhmmsAttributeSet aAttrib;
aAttrib.add(nknots,"nknots");
aAttrib.add(from_atom,"from_atom");
aAttrib.add(to_atom,"to_atom");
aAttrib.add(from_pos,"from_pos");
aAttrib.add(to_pos,"to_pos");
aAttrib.put(cur);
// sanity check
if ( nknots < 2 ) nknots=2;
// check out the reference atom coordinates
if (ions)
{
if ( from_atom >= 0 && from_atom < ions->R.size() )
from_pos = ions->R[from_atom];
if ( to_atom >= 0 && to_atom < ions->R.size() )
to_pos = ions->R[to_atom];
}
// prepare a fake particle set
ValueVector_t SPO_v, SPO_l;
GradVector_t SPO_g;
int OrbitalSize(mySPOSet->size());
SPO_v.resize(OrbitalSize);
SPO_g.resize(OrbitalSize);
SPO_l.resize(OrbitalSize);
double Delta= 1.0/(nknots-1);
int elec_count=target.R.size();
ParticleSet::SingleParticlePos_t zero_pos(0.0,0.0,0.0);
for(int icount=0, ind=0; icount<nknots; icount++, ind++)
{
if ( ind == elec_count ) ind = 0 ;
target.R[ind][0] = (to_pos[0]-from_pos[0]) * Delta * icount + from_pos[0];
target.R[ind][1] = (to_pos[1]-from_pos[1]) * Delta * icount + from_pos[1];
target.R[ind][2] = (to_pos[2]-from_pos[2]) * Delta * icount + from_pos[2];
target.makeMoveAndCheck(ind, zero_pos);
mySPOSet->evaluate(target, ind, SPO_v, SPO_g, SPO_l);
std::ostringstream o;
o << "x,y,z " << std::fixed << std::setprecision(7) << target.R[ind][0] << " " << target.R[ind][1] << " " << target.R[ind][2] ;
output_v << o.str() << " : " << std::scientific << std::setprecision(12);
output_g << o.str() << " : " << std::scientific << std::setprecision(12);
output_l << o.str() << " : " << std::scientific << std::setprecision(12);
for(int iorb=0; iorb<OrbitalSize; iorb++)
{
output_v << SPO_v[iorb] << " ";
output_g << SPO_g[iorb][0] << " " << SPO_g[iorb][1] << " " << SPO_g[iorb][2] << " ";
output_l << SPO_l[iorb] << " ";
}
output_v << std::endl;
output_g << std::endl;
output_l << std::endl;
}
output_v.close();
output_g.close();
output_l.close();
}
};
}
#endif
/***************************************************************************
* $RCSfile$ $Author: abenali $
* $Revision: 7138 $ $Date: 2016-09-27 18:45:29 -0500 (Tue, 27 Sep 2016) $
* $Id: OrbitalBuilderBase.h 7138 2016-09-27 23:45:29Z abenali $
***************************************************************************/

View File

@ -186,7 +186,7 @@ bool WaveFunctionFactory::addFermionTerm(xmlNodePtr cur)
//#if OHMMS_DIM == 3 && QMC_BUILD_LEVEL>1
else if(orbtype == "PWBasis" || orbtype == "PW" || orbtype == "pw")
{
detbuilder = new PWOrbitalBuilder(*targetPtcl,*targetPsi);
detbuilder = new PWOrbitalBuilder(*targetPtcl,*targetPsi,ptclPool);
}
//#endif /* QMC_BUILD_LEVEL>1 */
else

View File

@ -127,7 +127,7 @@ const char *particles =
xmlNodePtr pw1 = xmlFirstElementChild(root);
PWOrbitalBuilder pw_builder(elec, psi);
PWOrbitalBuilder pw_builder(elec, psi, ptcl.getPool());
pw_builder.put(pw1);
REQUIRE(psi.getOrbitals().size() == 1);
@ -286,7 +286,7 @@ const char *particles =
xmlNodePtr pw1 = xmlFirstElementChild(root);
PWOrbitalBuilder pw_builder(elec, psi);
PWOrbitalBuilder pw_builder(elec, psi, ptcl.getPool());
pw_builder.put(pw1);
REQUIRE(psi.getOrbitals().size() == 1);