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

Remove unused QMCTools

This commit is contained in:
Jaron Krogel 2018-05-03 09:41:05 -04:00
parent fa53be33ae
commit fb5b9267ee
32 changed files with 1 additions and 13899 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
src/QMCTools/Any2Slater.h
View File

@ -1,145 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_NONLINEAR_FITTING_H
#define QMCPLUSPLUS_NONLINEAR_FITTING_H
#include "Numerics/SlaterTypeOrbital.h"
#include "Numerics/OneDimIntegration.h"
#include "Optimize/Minimize.h"
/** A class to optimize a radial function by a Slater-type orbital
*
* scalar is typedef (double) inherited from MinimizeFunction
* Implements the virtual functions of MinizeFunction for optimization
* - int NumParams() ; returns the number of parameters to be optimized
* - scalar& Params(int i); assign the i-th value of the optimizable parameters
* - scalar Params(int i); returns the i-th value of the optimizable parameters
* - scalar Cost(); returns the cost function
* - void WriteStuff
*/
class Any2Slater : public MinimizeFunction
{
public:
///typedef of the source function
typedef OneDimGridFunctor<scalar> SourceType;
///typedef of the grid, using LogGrid
typedef LogGrid<scalar> GridType;
Any2Slater(SourceType& in): Source(in),
cg_tolerance(1e-6),cg_stepsize(0.001),cg_epsilon(1e-6)
{
OptParams.resize(1);
psi_sq.resize(Source.size());
scalar r2psi2=-10000.0;
int igMax=-1;
for(int ig=1; ig<Source.size(); ig++)
{
psi_sq[ig]=Source(ig)*Source(ig);
scalar r=Source.r(ig);
scalar t=r*r*psi_sq[ig];
if(t>r2psi2)
{
r2psi2=t;
igMax=ig;
}
}
//integrate_RK2_forward(psi_sq,psi_norm);
igMax+=2;
Target.N=1;
Target.Power=0;
Target.Z=1/Source.r(igMax);
Target.reset();
Source.grid().locate(0.2*Source.r(igMax));
minIndex=Source.grid().Loc;
std::cout << "Initial exponent " << Target.Z << " at " << Source.r(igMax) << std::endl;
std::cout << "The index of the cutoff " << minIndex << " at " << Source.r(minIndex) << std::endl;
OptParams[0]=Target.Z;
}
bool put(xmlNodePtr cur)
{
mPtr=cur;
return true;
}
///return the number of optimizable parameters
int NumParams()
{
return OptParams.size();
}
///assign optimization parameter i
scalar& Params(int i)
{
return OptParams[i];
}
///return optimization parameter i
scalar Params(int i) const
{
return OptParams[i];
}
///use the OptParams modified the optimization library and evaluate the cost function
scalar Cost()
{
Target.Z=OptParams[0];
Target.reset();
scalar del=0.0;
for(int ig=minIndex; ig<Source.size(); ig++)
{
//scalar y= Target.f((*myGrid)(ig))-Source[ig];
scalar r=Source.r(ig);
scalar y= Target.f(r);
scalar t=r*r*(y*y*-psi_sq[ig]);
del += t*t;
}
return del;
}
void WriteStuff()
{
std::cout << "Slater Z = " << Target.Z << " Norm = " << Target.Norm << std::endl;
}
/** main optimization function using ConjugateGradient method
*/
bool optimize()
{
ConjugateGradient CG;
CG.Tolerance = cg_tolerance;
CG.StepSize = cg_stepsize;
CG.epsilon = cg_epsilon;
CG.Minimize(*this);
return true;
}
private:
SourceType& Source;
scalar cg_tolerance;
scalar cg_stepsize;
scalar cg_epsilon;
int minIndex;
xmlNodePtr mPtr;
GenericSTO<scalar> Target;
std::vector<scalar> OptParams;
Vector<scalar> psi_sq, psi_norm;
};
#endif

275
src/QMCTools/Assign.cpp
View File

@ -1,275 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jordan E. Vincent, University of Illinois at Urbana-Champaign
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jordan E. Vincent, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Configuration.h"
#include "ParticleBase/RandomSeqGenerator.h"
#include "ParticleIO/XMLParticleIO.h"
#include "Message/Communicate.h"
#include "Utilities/OhmmsInfo.h"
#include "Particle/MCWalkerConfiguration.h"
#include "Particle/HDFWalkerIO.h"
#include "Particle/DistanceTable.h"
#include "Particle/DistanceTableData.h"
#include "Numerics/OneDimGridBase.h"
#include "OhmmsPETE/OhmmsMatrix.h"
#include "OhmmsApp/ProjectData.h"
#include "OhmmsApp/RandomNumberControl.h"
#include "QMC/QMCUtilities.h"
#include <fstream>
int main(int argc, char **argv)
{
using namespace qmcplusplus;
xmlDocPtr m_doc;
xmlNodePtr m_root;
xmlXPathContextPtr m_context;
enum {SourceIndex = DistanceTableData::SourceIndex};
OHMMS::Controller->initialize(argc,argv);
OhmmsInfo welcome(argc,argv,OHMMS::Controller->mycontext());
///project description
OHMMS::ProjectData myProject;
///random number controller
OHMMS::RandomNumberControl myRandomControl;
if(argc>1)
{
// build an XML tree from a the file;
LOGMSG("Opening file " << argv[1])
m_doc = xmlParseFile(argv[1]);
if (m_doc == NULL)
{
ERRORMSG("File " << argv[1] << " is invalid")
xmlFreeDoc(m_doc);
return 1;
}
// Check the document is of the right kind
m_root = xmlDocGetRootElement(m_doc);
if (m_root == NULL)
{
ERRORMSG("Empty document");
xmlFreeDoc(m_doc);
return 1;
}
}
else
{
WARNMSG("No argument is given. Assume that does not need an input file")
}
m_context = xmlXPathNewContext(m_doc);
xmlXPathObjectPtr result
= xmlXPathEvalExpression((const xmlChar*)"//project",m_context);
if(xmlXPathNodeSetIsEmpty(result->nodesetval))
{
WARNMSG("Project is not defined")
myProject.reset();
}
else
{
myProject.put(result->nodesetval->nodeTab[0]);
}
xmlXPathFreeObject(result);
//initialize the random number generator
xmlNodePtr rptr = myRandomControl.initialize(m_context);
if(rptr)
{
xmlAddChild(m_root,rptr);
}
///the ions
ParticleSet ion;
MCWalkerConfiguration el;
el.setName("e");
int iu = el.Species.addSpecies("u");
int id = el.Species.addSpecies("d");
int icharge = el.Species.addAttribute("charge");
el.Species(icharge,iu) = -1;
el.Species(icharge,id) = -1;
bool init_els = determineNumOfElectrons(el,m_context);
result
= xmlXPathEvalExpression((const xmlChar*)"//particleset",m_context);
xmlNodePtr el_ptr=NULL, ion_ptr=NULL;
for(int i=0; i<result->nodesetval->nodeNr; i++)
{
xmlNodePtr cur=result->nodesetval->nodeTab[i];
xmlChar* aname= xmlGetProp(cur,(const xmlChar*)"name");
if(aname)
{
char fc = aname[0];
if(fc == 'e')
{
el_ptr=cur;
}
else
if(fc == 'i')
{
ion_ptr=cur;
}
}
}
bool donotresize = false;
if(init_els)
{
el.setName("e");
XMLReport("The configuration for electrons is already determined by the wave function")
donotresize = true;
}
if(el_ptr)
{
XMLParticleParser pread(el,donotresize);
pread.put(el_ptr);
}
if(ion_ptr)
{
XMLParticleParser pread(ion);
pread.put(ion_ptr);
}
xmlXPathFreeObject(result);
if(!ion.getTotalNum())
{
ion.setName("i");
ion.create(1);
ion.R[0] = 0.0;
}
//The ion-ion distance-table
DistanceTableData* d_ii = DistanceTable::getTable(DistanceTable::add(ion));
d_ii->create(1);
d_ii->evaluate(ion);
std::vector<double> Cut, Core;
int Centers = ion.getTotalNum();
//attribute id for cut
int icut = ion.Species.addAttribute("cut");
//store the max distance from atom
Cut.resize(Centers);
for(int iat=0; iat<Centers; iat++)
{
int id = ion.GroupID[iat];
Cut[iat] = ion.Species(icut,id);
}
int icore = ion.Species.addAttribute("core");
//store the max distance from atom
Core.resize(Centers);
for(int iat=0; iat<Centers; iat++)
{
Core[iat]=ion.Species(icore,ion.GroupID[iat]);
}
//3N-dimensional Gaussian
ParticleSet::ParticlePos_t chi(el.getTotalNum());
makeGaussRandom(chi);
//determine if odd or even number of particles
int irem = el.getTotalNum()%2;
int ihalf = el.getTotalNum()/2;
//assign the core
int ncore(0);
for(int iat=0; iat<Centers; iat++)
{
double sep=0.8*Cut[iat];
for(int iel=0; iel<Core[iat]/2; iel++,ncore++)
{
el.R[ncore]=ion.R[iat]+sep*chi[ncore];
el.R[ncore+ihalf]=ion.R[iat]+sep*chi[ncore+ihalf];
}
}
int ipart = ncore;
int isave_iat=0;
for(int iat=0; iat<Centers; iat++)
{
for(int nn=d_ii->M[iat]; nn<d_ii->M[iat+1]; nn++)
{
double bondlength = d_ii->r(nn);
int jat = d_ii->J[nn];
//only assign if the half bond-length < cutoff
if(bondlength < Cut[iat]+Cut[jat])
{
if(ipart < ihalf)
{
XMLReport("Assigning particles = " << ipart << " and " << ipart+ihalf)
/*place 2 electrons (an up and a down) at half
the bond-length plus a random number multiplied
by 10% of the bond-length*/
el.R[ipart] = ion.R[iat]+0.5*d_ii->dr(nn)+0.1*bondlength*chi[ipart];
el.R[ipart+ihalf] = ion.R[iat]+0.5*d_ii->dr(nn)+0.1*bondlength*chi[ipart+ihalf];
ipart++;
isave_iat = iat;
}
}
}
}
//assign the last particle (if odd number of particles)
int flag = 1;
ipart = el.getTotalNum()-1;
if(irem)
{
XMLReport("Assigning last particle.")
for(int iat = isave_iat+1; iat<Centers; iat++)
{
for(int nn=d_ii->M[iat]; nn<d_ii->M[iat+1]; nn++)
{
double bondlength = d_ii->r(nn);
if((0.5*bondlength < Cut[iat]) && flag)
{
XMLReport("Assigning particle = " << ipart)
el.R[ipart] = ion.R[iat]+0.5*d_ii->dr(nn)+0.1*bondlength*chi[ipart];
flag = 0;
}
}
}
}
std::cout << "Ionic configuration : " << ion.getName() << std::endl;
ion.get(std::cout);
std::cout << "Electronic configuration : " << el.getName() << std::endl;
el.get(std::cout);
std::string newxml(myProject.CurrentRoot());
newxml.append(".ptcl.xml");
std::ofstream ptcl_out(newxml.c_str());
/*
std::ofstream molmol("assign.xyz");
molmol << Centers+el.getTotalNum() << std::endl;
molmol << std::endl;
for(int iat=0; iat<Centers; iat++)
molmol << ion.Species.speciesName[ion.GroupID[iat]] << 0.5292*ion.R[iat] << std::endl;
for(int ipart=0; ipart<el.getTotalNum(); ipart++)
molmol << "He" << 0.5292*el.R[ipart] << std::endl;
molmol.close();
*/
xmlXPathFreeContext(m_context);
xmlFreeDoc(m_doc);
int nup = el.last(0);
int ndown = el.last(1)-el.last(0);
ptcl_out << "<?xml version=\"1.0\"?>" << std::endl;
ptcl_out << "<particleset name=\"e\">" << std::endl;
ptcl_out << "<group name=\"u\" size=\"" << nup << "\">" << std::endl;
ptcl_out << "<parameter name=\"charge\">-1</parameter>" << std::endl;
ptcl_out << "<attrib name=\"position\" datatype=\"posArray\">" << std::endl;
for (int ipart=0; ipart<nup; ++ipart)
ptcl_out << el.R[ipart] << std::endl;
ptcl_out << "</attrib>" << std::endl;
ptcl_out << "</group>" << std::endl;
ptcl_out << "<group name=\"d\" size=\"" << ndown << "\">" << std::endl;
ptcl_out << "<parameter name=\"charge\">-1</parameter>" << std::endl;
ptcl_out << "<attrib name=\"position\" datatype=\"posArray\">" << std::endl;
for (int ipart=nup; ipart<el.getTotalNum(); ++ipart)
ptcl_out << el.R[ipart] << std::endl;
ptcl_out << "</attrib>" << std::endl;
ptcl_out << "</group>" << std::endl;
ptcl_out << "</particleset>" << std::endl;
OHMMS::Controller->finalize();
return 0;
}

2
src/QMCTools/CMakeLists.txt
View File

@ -49,7 +49,7 @@ SET(MOSRCS
# create libmocommon # create libmocommon
ADD_LIBRARY(mocommon ${MOSRCS}) ADD_LIBRARY(mocommon ${MOSRCS})
set(QTOOLS convert4qmc MSDgenerator extract-eshdf-kvectors trace_density) set(QTOOLS convert4qmc MSDgenerator extract-eshdf-kvectors)
ADD_EXECUTABLE(getSupercell getSupercell.cpp) ADD_EXECUTABLE(getSupercell getSupercell.cpp)
FOREACH(p ${QTOOLS}) FOREACH(p ${QTOOLS})

1
src/QMCTools/GTO2GridBuilder.cpp
View File

@ -20,7 +20,6 @@
#include "Numerics/Transform2GridFunctor.h" #include "Numerics/Transform2GridFunctor.h"
#include "Numerics/OneDimCubicSpline.h" #include "Numerics/OneDimCubicSpline.h"
#include "Numerics/GaussianBasisSet.h" #include "Numerics/GaussianBasisSet.h"
#include "QMCTools/GridMolecularOrbitals.h"
#include "QMCTools/GTO2GridBuilder.h" #include "QMCTools/GTO2GridBuilder.h"
#include "QMCFactory/OneDimGridFactory.h" #include "QMCFactory/OneDimGridFactory.h"

331
src/QMCTools/GridMolecularOrbitals.cpp
View File

@ -1,331 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Utilities/OhmmsInfo.h"
#include "Particle/DistanceTableData.h"
#include "Particle/DistanceTable.h"
#include "OhmmsData/AttributeSet.h"
#include "QMCWaveFunctions/DetSetBuilderWithBasisSet.h"
#include "QMCTools/GridMolecularOrbitals.h"
#include "QMCTools/RGFBuilderBase.h"
#include "QMCTools/STO2GridBuilder.h"
#include "QMCTools/GTO2GridBuilder.h"
#include "QMCTools/Any2GridBuilder.h"
#include "QMCTools/NumericalRGFBuilder.h"
namespace qmcplusplus
{
GridMolecularOrbitals::GridMolecularOrbitals(ParticleSet& els, TrialWaveFunction& psi,
ParticleSet& ions):
OrbitalBuilderBase(els,psi), IonSys(ions), BasisSet(0), d_table(0), rbuilder(0)
{
//int d_ie = DistanceTable::add(ions,els);
d_table = DistanceTable::add(ions,els);
nlms_id["n"] = q_n;
nlms_id["l"] = q_l;
nlms_id["m"] = q_m;
nlms_id["s"] = q_s;
}
bool GridMolecularOrbitals::put(xmlNodePtr cur)
{
LOGMSG("GridMolecularOrbitals::put")
DetSetBuilderWithBasisSet<GridMolecularOrbitals> spobuilder(targetPtcl,targetPsi,*this);
if(spobuilder.put(cur))
{
BasisSet->resize(spobuilder.NumPtcl);
return true;
}
else
{
return false;
}
}
GridMolecularOrbitals::BasisSetType*
GridMolecularOrbitals::addBasisSet(xmlNodePtr cur)
{
if(!BasisSet)
BasisSet = new BasisSetType(IonSys.getSpeciesSet().getTotalNum());
QuantumNumberType nlms;
std::string rnl;
//current number of centers
int ncenters = CenterID.size();
int activeCenter;
int gridmode = -1;
bool addsignforM = false;
std::string sph("default"), Morder("gaussian");
//go thru the tree
cur = cur->xmlChildrenNode;
std::map<std::string,RGFBuilderBase*> rbuilderlist;
while(cur!=NULL)
{
std::string cname((const char*)(cur->name));
if(cname == basis_tag || cname == "atomicBasisSet")
{
int expandlm = GAUSSIAN_EXPAND;
std::string abasis("invalid"), btype("Numerical");
//Register valid attributes attributes
OhmmsAttributeSet aAttrib;
aAttrib.add(abasis,"elementType");
aAttrib.add(abasis,"species");
aAttrib.add(btype,"type");
aAttrib.add(sph,"angular");
aAttrib.add(addsignforM,"expM");
aAttrib.add(Morder,"expandYlm");
aAttrib.put(cur);
if(abasis == "invalid")
continue;
if(sph == "spherical")
addsignforM=1; //include (-1)^m
if(Morder == "gaussian")
{
expandlm = GAUSSIAN_EXPAND;
}
else
if(Morder == "natural")
{
expandlm = NATURAL_EXPAND;
}
else
if(Morder == "no")
{
expandlm = DONOT_EXPAND;
}
if(addsignforM)
LOGMSG("Spherical Harmonics contain (-1)^m factor")
else
LOGMSG("Spherical Harmonics DO NOT contain (-1)^m factor")
//search the species name
std::map<std::string,int>::iterator it = CenterID.find(abasis);
if(it == CenterID.end())
//add the name to the map CenterID
{
if(btype == "Numerical" || btype == "NG" || btype == "HFNG")
{
rbuilder = new NumericalRGFBuilder(cur);
}
else
{
rbuilder = new Any2GridBuilder(cur);
}
//CenterID[abasis] = activeCenter = ncenters++;
CenterID[abasis]=activeCenter=IonSys.getSpeciesSet().findSpecies(abasis);
int Lmax(0); //maxmimum angular momentum of this center
int num(0);//the number of localized basis functions of this center
//process the basic property: maximun angular momentum, the number of basis functions to be added
std::vector<xmlNodePtr> radGroup;
xmlNodePtr cur1 = cur->xmlChildrenNode;
xmlNodePtr gptr=0;
while(cur1 != NULL)
{
std::string cname1((const char*)(cur1->name));
if(cname1 == basisfunc_tag || cname1 == "basisGroup")
{
radGroup.push_back(cur1);
int l=atoi((const char*)(xmlGetProp(cur1, (const xmlChar *)"l")));
Lmax = std::max(Lmax,l);
//expect that only Rnl is given
if(expandlm)
num += 2*l+1;
else
num++;
}
else
if(cname1 == "grid")
{
gptr = cur1;
}
cur1 = cur1->next;
}
XMLReport("Adding a center " << abasis << " centerid "<< CenterID[abasis])
XMLReport("Maximum angular momentum = " << Lmax)
XMLReport("Number of centered orbitals = " << num)
//create a new set of atomic orbitals sharing a center with (Lmax, num)
//if(addsignforM) the basis function has (-1)^m sqrt(2)Re(Ylm)
CenteredOrbitalType* aos = new CenteredOrbitalType(Lmax,addsignforM);
aos->LM.resize(num);
aos->NL.resize(num);
//Now, add distinct Radial Orbitals and (l,m) channels
num=0;
rbuilder->setOrbitalSet(aos,abasis); //assign radial orbitals for the new center
rbuilder->addGrid(gptr); //assign a radial grid for the new center
std::vector<xmlNodePtr>::iterator it(radGroup.begin());
std::vector<xmlNodePtr>::iterator it_end(radGroup.end());
while(it != it_end)
{
cur1 = (*it);
xmlAttrPtr att = cur1->properties;
while(att != NULL)
{
std::string aname((const char*)(att->name));
if(aname == "rid" || aname == "id")
//accept id/rid
{
rnl = (const char*)(att->children->content);
}
else
{
std::map<std::string,int>::iterator iit = nlms_id.find(aname);
if(iit != nlms_id.end())
//valid for n,l,m,s
{
nlms[(*iit).second] = atoi((const char*)(att->children->content));
}
}
att = att->next;
}
XMLReport("\n(n,l,m,s) " << nlms[0] << " " << nlms[1] << " " << nlms[2] << " " << nlms[3])
//add Ylm channels
num = expandYlm(rnl,nlms,num,aos,cur1,expandlm);
++it;
}
//add the new atomic basis to the basis set
BasisSet->add(aos,activeCenter);
#if !defined(HAVE_MPI)
rbuilder->print(abasis,1);
#endif
if(rbuilder)
{
delete rbuilder;
rbuilder=0;
}
}
else
{
WARNMSG("Species " << abasis << " is already initialized. Ignore the input.")
}
}
cur = cur->next;
}
if(BasisSet)
{
BasisSet->setTable(d_table);
LOGMSG("The total number of basis functions " << BasisSet->TotalBasis)
return BasisSet;
}
else
{
ERRORMSG("BasisSet is not initialized.")
return NULL;
}
}
int
GridMolecularOrbitals::expandYlm(const std::string& rnl, const QuantumNumberType& nlms,
int num, CenteredOrbitalType* aos, xmlNodePtr cur1,
int expandlm)
{
if(expandlm == GAUSSIAN_EXPAND)
{
XMLReport("Expanding Ylm according to Gaussian98")
std::map<std::string,int>::iterator rnl_it = RnlID.find(rnl);
if(rnl_it == RnlID.end())
{
int nl = aos->Rnl.size();
if(rbuilder->addRadialOrbital(cur1,nlms))
{
RnlID[rnl] = nl;
int l = nlms[q_l];
XMLReport("Adding " << 2*l+1 << " spherical orbitals for l= " << l)
switch (l)
{
case(0):
aos->LM[num] = aos->Ylm.index(0,0);
aos->NL[num] = nl;
num++;
break;
case(1)://px(1),py(-1),pz(0)
aos->LM[num] = aos->Ylm.index(1,1);
aos->NL[num] = nl;
num++;
aos->LM[num] = aos->Ylm.index(1,-1);
aos->NL[num] = nl;
num++;
aos->LM[num] = aos->Ylm.index(1,0);
aos->NL[num] = nl;
num++;
break;
default://0,1,-1,2,-2,...,l,-l
aos->LM[num] = aos->Ylm.index(l,0);
aos->NL[num] = nl;
num++;
for(int tm=1; tm<=l; tm++)
{
aos->LM[num] = aos->Ylm.index(l,tm);
aos->NL[num] = nl;
num++;
aos->LM[num] = aos->Ylm.index(l,-tm);
aos->NL[num] = nl;
num++;
}
break;
}
}
}
}
else
if(expandlm == NATURAL_EXPAND)
{
XMLReport("Expanding Ylm as -l,-l+1,...,l-1,l")
std::map<std::string,int>::iterator rnl_it = RnlID.find(rnl);
if(rnl_it == RnlID.end())
{
int nl = aos->Rnl.size();
if(rbuilder->addRadialOrbital(cur1,nlms))
{
RnlID[rnl] = nl;
int l = nlms[q_l];
XMLReport("Adding " << 2*l+1 << " spherical orbitals")
for(int tm=-l; tm<=l; tm++,num++)
{
aos->LM[num] = aos->Ylm.index(l,tm);
aos->NL[num] = nl;
}
}
}
}
else
{
XMLReport("Ylm is NOT expanded.")
//assign the index for real Spherical Harmonic with (l,m)
aos->LM[num] = aos->Ylm.index(nlms[q_l],nlms[q_m]);
//radial orbitals: add only distinct orbitals
std::map<std::string,int>::iterator rnl_it = RnlID.find(rnl);
if(rnl_it == RnlID.end())
{
int nl = aos->Rnl.size();
if(rbuilder->addRadialOrbital(cur1,nlms))
//assign the index for radial orbital with (n,l)
{
aos->NL[num] = nl;
RnlID[rnl] = nl;
}
}
else
{
//assign the index for radial orbital with (n,l) if repeated
XMLReport("Already added radial function for id: " << rnl)
aos->NL[num] = (*rnl_it).second;
}
//increment number of basis functions
num++;
}
return num;
}
}

99
src/QMCTools/GridMolecularOrbitals.h
View File

@ -1,99 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_CONVERT2_RADIALGRID_H
#define QMCPLUSPLUS_CONVERT2_RADIALGRID_H
#include "QMCWaveFunctions/OrbitalBuilderBase.h"
#include "QMCTools/MolecularOrbitalBasis.h"
#include "QMCTools/RGFBuilderBase.h"
namespace qmcplusplus
{
/** derived class from OrbitalBuilderBase
*
* Create a basis set of molecular orbital types as defined in MolecularOrbitalBasis
* with radial wave functions on the radial grids.
*/
class GridMolecularOrbitals: public OrbitalBuilderBase
{
public:
//@typedef radial grid type
typedef OneDimGridBase<RealType> GridType;
//@typedef \f$R_{nl}\f$ radial functor type defined on a grid
typedef OneDimGridFunctor<RealType> RadialOrbitalType;
//@typedef centered orbital type which represents a set of \f$R_{nl}Y_{lm}\f$ for a center
typedef SphericalBasisSet<RadialOrbitalType,GridType> CenteredOrbitalType;
//@typedef molecuar orbital basis composed of multiple CenteredOrbitalType s
typedef MolecularOrbitalBasis<CenteredOrbitalType> BasisSetType;
/** constructor
* \param els reference to the electrons
* \param psi reference to the wavefunction
* \param ions reference to the ions
*/
GridMolecularOrbitals(ParticleSet& els, TrialWaveFunction& psi, ParticleSet& ions);
/** initialize the Antisymmetric wave function for electrons
*@param cur the current xml node
*
*/
bool put(xmlNodePtr cur);
/** process basis element to build the basis set
*@param cur the current xml node
*@return a pointer to the BasisSet
*
*This member function is necessary in order to use SDSetBuilderWithBasisSet
*to initialize fermion wavefunctions with slater determinants.
*/
BasisSetType* addBasisSet(xmlNodePtr cur);
private:
enum {DONOT_EXPAND=0, GAUSSIAN_EXPAND=1, NATURAL_EXPAND};
///reference to the ionic system with which the basis set are associated
ParticleSet& IonSys;
///pointer to the BasisSet built by GridMolecularOrbitals
BasisSetType* BasisSet;
///distance table pointer
DistanceTableData* d_table;
///Current RadiaaGridFunctorBuilder
RGFBuilderBase* rbuilder;
///map for the radial orbitals
std::map<std::string,int> RnlID;
///map for the centers
std::map<std::string,int> CenterID;
///map for (n,l,m,s) to its quantum number index
std::map<std::string,int> nlms_id;
///append Ylm channels
int expandYlm(const std::string& rnl, const QuantumNumberType& nlms, int num,
CenteredOrbitalType* aos, xmlNodePtr cur1,
int expandlm=DONOT_EXPAND);
};
}
#endif

208
src/QMCTools/HDFWalkerMerger.cpp
View File

@ -1,208 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <iostream>
#include "QMCTools/HDFWalkerMerger.h"
#include "Utilities/RandomGeneratorIO.h"
using namespace qmcplusplus;
HDFWalkerMerger::HDFWalkerMerger(const std::string& aroot, int ncopy)
: NumCopy(ncopy), FileRoot(aroot)
{
}
HDFWalkerMerger::~HDFWalkerMerger() { }
void HDFWalkerMerger::init()
{
char fname[128];
char GrpName[128];
hsize_t dimin[3], summary_size=3;
std::vector<double> summaryIn;
//determine the minimum configuration and walkers for each config
int min_config=10000;
MaxNumWalkers=0;
for(int ip=0; ip<NumCopy; ip++)
{
int nconf=1;
numWalkersIn.push_back(new std::vector<hsize_t>);
sprintf(fname,"%s.p%03d.config.h5",FileRoot.c_str(),ip);
hid_t hfile = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
hid_t h_config = H5Gopen(hfile,"config_collection");
hid_t h1=H5Dopen(h_config,"NumOfConfigurations");
if(h1>-1)
{
H5Dread(h1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&(nconf));
H5Dclose(h1);
}
for(int iconf=0; iconf<nconf; iconf++)
{
sprintf(GrpName,"config%04d/coord",iconf);
h1 = H5Dopen(h_config,GrpName);
hid_t dataspace = H5Dget_space(h1);
int rank = H5Sget_simple_extent_ndims(dataspace);
int status_n = H5Sget_simple_extent_dims(dataspace, dimin, NULL);
H5Sclose(dataspace);
H5Dclose(h1);
numWalkersIn[ip]->push_back(dimin[0]);
MaxNumWalkers = std::max(MaxNumWalkers,dimin[0]);
}
h1 = H5Dopen(h_config,"Summary");
hid_t summary_space = H5Dget_space(h1);
H5Sget_simple_extent_ndims(summary_space);
H5Sget_simple_extent_dims(summary_space,&summary_size, NULL);
if(Summary.size() != summary_size)
{
Summary.resize(summary_size,0.0);
summaryIn.resize(summary_size,0.0);
}
H5Dread(h1, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &summaryIn[0]);
for(hsize_t k=0; k<summary_size; k++)
{
Summary[k] += summaryIn[k];
}
H5Sclose(summary_space);
H5Dclose(h1);
H5Fclose(hfile);
min_config = std::min(min_config,nconf);
}
NumPtcl = dimin[1];
Dimension = dimin[2];
NumConfig = min_config;
OffSet.resize(NumCopy+1,NumConfig);
OffSet=0;
for(int ip=0; ip<NumCopy; ip++)
{
for(int iconf=0; iconf<NumConfig; iconf++)
{
OffSet(ip+1,iconf) = OffSet(ip,iconf)+(*numWalkersIn[ip])[iconf];
}
}
}
void HDFWalkerMerger::writeHeader(hid_t gid)
{
int cur_version=1;
hsize_t dim=1;
hid_t header_space= H5Screate_simple(1, &dim, NULL);
hid_t header_set= H5Dcreate(gid, "version", H5T_NATIVE_INT, header_space, H5P_DEFAULT);
H5Dwrite(header_set, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&cur_version);
H5Dclose(header_set);
header_set= H5Dcreate(gid, "ncontexts", H5T_NATIVE_INT, header_space, H5P_DEFAULT);
H5Dwrite(header_set, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&NumCopy);
H5Dclose(header_set);
H5Sclose(header_space);
}
void HDFWalkerMerger::merge()
{
init();
char fname[128];
char GrpName[128];
char ConfName[128];
#if H5_VERS_RELEASE < 4
hssize_t offset[]= {0,0,0};
#else
hsize_t offset[]= {0,0,0};
#endif
typedef std::vector<OHMMS_PRECISION> Container_t;
Container_t vecIn(MaxNumWalkers*NumPtcl*Dimension);
sprintf(fname,"%s.config.h5",FileRoot.c_str());
hid_t masterfile= H5Fcreate(fname,H5F_ACC_TRUNC, H5P_DEFAULT,H5P_DEFAULT);
writeHeader(masterfile);
hid_t mastercf = H5Gcreate(masterfile,"config_collection",0);
hsize_t onedim=1;
hid_t dataspace= H5Screate_simple(1, &onedim, NULL);
hid_t dataset= H5Dcreate(mastercf, "NumOfConfigurations", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&NumConfig);
H5Dclose(dataset);
H5Sclose(dataspace);
hsize_t dimout[3], dimin[3];
dimin[1]=NumPtcl;
dimin[2]=Dimension;
dimout[1]=NumPtcl;
dimout[2]=Dimension;
//create new dataspace for each config
for(int iconf=0; iconf<NumConfig; iconf++)
{
dimout[0]=OffSet(NumCopy,iconf);
sprintf(ConfName,"config%04d",iconf);
hid_t cf_id = H5Gcreate(mastercf,ConfName,0);
dataspace = H5Screate_simple(3,dimout,NULL);
dataset = H5Dcreate(cf_id,"coord",H5T_NATIVE_DOUBLE, dataspace,H5P_DEFAULT);
H5Dclose(dataset);
H5Sclose(dataspace);
H5Gclose(cf_id);
}
//use random_states opposed to random_states
for(int ip=0; ip<NumCopy; ip++)
{
sprintf(fname,"%s.p%03d.config.h5",FileRoot.c_str(),ip);
hid_t h0 = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
hid_t cfcollect = H5Gopen(h0,"config_collection");
const std::vector<hsize_t>& nw(*(numWalkersIn[ip]));
for(int iconf=0; iconf<NumConfig; iconf++)
{
//get the number of walkers for (ip,iconf)
dimin[0] = nw[iconf];
sprintf(ConfName,"config%04d",iconf);
hid_t cf_id = H5Gopen(cfcollect,ConfName);
HDFAttribIO<Container_t> in(vecIn, 3, dimin);
in.read(cf_id,"coord");
H5Gclose(cf_id);
//change the offset
offset[0]=OffSet(ip,iconf);
//create simple data to write
hid_t memspace = H5Screate_simple(3, dimin, NULL);
cf_id = H5Gopen(mastercf,ConfName);
hid_t dataset = H5Dopen(cf_id,"coord");
hid_t filespace = H5Dget_space(dataset);
herr_t status = H5Sselect_hyperslab(filespace,H5S_SELECT_SET, offset,NULL,dimin,NULL);
status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace, filespace, H5P_DEFAULT, &vecIn[0]);
H5Dclose(dataset);
H5Sclose(filespace);
H5Sclose(memspace);
H5Gclose(cf_id);
}
H5Gclose(cfcollect);
sprintf(GrpName,"context%04d",ip);
hid_t mycontext = H5Gcreate(masterfile,GrpName,0);
hid_t ranIn = H5Gopen(h0,"random_state");
HDFAttribIO<RandomGenerator_t> r(Random);
r.read(ranIn,"dummy");
H5Gclose(ranIn);
hid_t ranOut=H5Gcreate(mycontext, "random_state",0);
r.write(ranOut,"dummy");
H5Gclose(ranOut);
H5Gclose(mycontext);
H5Fclose(h0);
}
//write the summaryset
double d=1.0/static_cast<double>(NumCopy);
for(int k=0; k<Summary.size(); k++)
Summary[k]*=d;
hsize_t dim=Summary.size();
hid_t summary_space = H5Screate_simple(1, &dim, NULL);
hid_t summary_set = H5Dcreate(mastercf, "Summary", H5T_NATIVE_DOUBLE, summary_space, H5P_DEFAULT);
ret = H5Dwrite(summary_set, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,&Summary[0]);
H5Dclose(summary_set);
H5Sclose(summary_space);
H5Gclose(mastercf);
H5Fclose(masterfile);
}

81
src/QMCTools/HDFWalkerMerger.h
View File

@ -1,81 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_WALKER_MERGE_IO_H
#define QMCPLUSPLUS_WALKER_MERGE_IO_H
#include <string>
#include "OhmmsPETE/OhmmsMatrix.h"
#include "Numerics/HDFSTLAttrib.h"
class Communicate;
namespace qmcplusplus
{
/** Writes a set of walker configurations to an HDF5 file. */
class HDFWalkerMerger
{
///the physical dimension
int Dimension;
///the number of separate files
int NumCopy;
///the number of configuration
int NumConfig;
///the number of particles
int NumPtcl;
///communicator
Communicate* myComm;
///maxmimum number of walkers for any config
hsize_t MaxNumWalkers;
///file root
std::string FileRoot;
///numWalkerIn[node]->operator[](iconfig) returns the number of walkers
std::vector<std::vector<hsize_t>*> numWalkersIn;
/** offset of the walkers
*
* OffSet(NumCopy,*) is the total number of walkers for the merged
* configuration.
*/
Matrix<hsize_t> OffSet;
/** summary is accumulated */
std::vector<double> Summary;
/** write header information */
void writeHeader(hid_t git);
/** initialize sizes */
void init();
public:
HDFWalkerMerger(const std::string& froot, int ncopy);
~HDFWalkerMerger();
/** set communicator */
void setCommunicator(Communicate* c);
void merge();
};
}
#endif

273
src/QMCTools/PC.cpp
View File

@ -1,273 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jordan E. Vincent, University of Illinois at Urbana-Champaign
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jordan E. Vincent, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Configuration.h"
#include "ParticleBase/RandomSeqGenerator.h"
#include "ParticleIO/XMLParticleIO.h"
#include "Message/Communicate.h"
#include "Utilities/OhmmsInfo.h"
#include "Particle/MCWalkerConfiguration.h"
#include "Particle/HDFWalkerIO.h"
#include "Particle/DistanceTable.h"
#include "Particle/DistanceTableData.h"
#include "Numerics/OneDimGridBase.h"
#include "OhmmsPETE/OhmmsMatrix.h"
#include "OhmmsApp/ProjectData.h"
#include "OhmmsApp/RandomNumberControl.h"
#include "QMC/QMCUtilities.h"
#include "OhmmsData/ParameterSet.h"
#include <fstream>
int main(int argc, char **argv)
{
using namespace qmcplusplus;
xmlDocPtr m_doc;
xmlNodePtr m_root;
xmlXPathContextPtr m_context;
enum {SourceIndex = DistanceTableData::SourceIndex};
OHMMS::Controller->initialize(argc,argv);
OhmmsInfo welcome(argc,argv,OHMMS::Controller->mycontext());
///project description
OHMMS::ProjectData myProject;
///random number controller
OHMMS::RandomNumberControl myRandomControl;
if(argc>1)
{
// build an XML tree from a the file;
LOGMSG("Opening file " << argv[1])
m_doc = xmlParseFile(argv[1]);
if (m_doc == NULL)
{
ERRORMSG("File " << argv[1] << " is invalid")
xmlFreeDoc(m_doc);
return 1;
}
// Check the document is of the right kind
m_root = xmlDocGetRootElement(m_doc);
if (m_root == NULL)
{
ERRORMSG("Empty document");
xmlFreeDoc(m_doc);
return 1;
}
}
else
{
WARNMSG("No argument is given. Assume that does not need an input file")
}
m_context = xmlXPathNewContext(m_doc);
xmlXPathObjectPtr result
= xmlXPathEvalExpression((const xmlChar*)"//project",m_context);
if(xmlXPathNodeSetIsEmpty(result->nodesetval))
{
WARNMSG("Project is not defined")
myProject.reset();
}
else
{
myProject.put(result->nodesetval->nodeTab[0]);
}
xmlXPathFreeObject(result);
//initialize the random number generator
xmlNodePtr rptr = myRandomControl.initialize(m_context);
if(rptr)
{
xmlAddChild(m_root,rptr);
}
///the ions
ParticleSet ion;
MCWalkerConfiguration el;
el.setName("e");
int iu = el.Species.addSpecies("u");
int id = el.Species.addSpecies("d");
int icharge = el.Species.addAttribute("charge");
el.Species(icharge,iu) = -1;
el.Species(icharge,id) = -1;
bool init_els = determineNumOfElectrons(el,m_context);
result
= xmlXPathEvalExpression((const xmlChar*)"//particleset",m_context);
xmlNodePtr el_ptr=NULL, ion_ptr=NULL;
for(int i=0; i<result->nodesetval->nodeNr; i++)
{
xmlNodePtr cur=result->nodesetval->nodeTab[i];
xmlChar* aname= xmlGetProp(cur,(const xmlChar*)"name");
if(aname)
{
char fc = aname[0];
if(fc == 'e')
{
el_ptr=cur;
}
else
if(fc == 'i')
{
ion_ptr=cur;
}
}
}
bool donotresize = false;
if(init_els)
{
el.setName("e");
XMLReport("The configuration for electrons is already determined by the wave function")
donotresize = true;
}
if(el_ptr)
{
XMLParticleParser pread(el,donotresize);
pread.put(el_ptr);
}
if(ion_ptr)
{
XMLParticleParser pread(ion);
pread.put(ion_ptr);
}
xmlXPathFreeObject(result);
if(!ion.getTotalNum())
{
ion.setName("i");
ion.create(1);
ion.R[0] = 0.0;
}
double ri = 0.0;
double rf = 20.0;
int npts = 100;
ParameterSet m_param;
m_param.add(ri,"ri","AU");
m_param.add(rf,"rf","AU");
m_param.add(npts,"npts","int");
std::string type = "elel";
result
= xmlXPathEvalExpression((const xmlChar*)"//paircorrelation",m_context);
xmlNodePtr cur=result->nodesetval->nodeTab[0];
type = ((const char*)xmlGetProp(cur,(const xmlChar*)"type"));
XMLReport("Type = " << type)
xmlNodePtr tcur = cur->children;
m_param.put(cur);
std::vector<xmlNodePtr> wset;
int pid=OHMMS::Controller->mycontext();
while(tcur != NULL)
{
std::string cname((const char*)(tcur->name));
if(cname == "mcwalkerset")
{
int pid_target=pid;
xmlChar* anode=xmlGetProp(tcur,(const xmlChar*)"node");
if(anode)
{
pid_target = atoi((const char*)anode);
}
if(pid_target == pid)
wset.push_back(tcur);
}
tcur=tcur->next;
}
std::vector<std::string> ConfigFile;
if(wset.size())
{
for(int i=0; i<wset.size(); i++)
{
xmlChar* att= xmlGetProp(wset[i],(const xmlChar*)"file");
if(att)
ConfigFile.push_back((const char*)att);
}
}
xmlXPathFreeObject(result);
int nup = el.last(0);
int ndown = el.last(1)-el.last(0);
int factor = 1;
if(nup)
factor*=nup;
if(ndown)
factor*=ndown;
DistanceTableData* d_table;
//create a grid
LinearGrid<double> grid;
grid.set(ri,rf,npts);
double rcut = grid.rmax();
XMLReport("Grid Values: ri = " << ri << " rf = " << rf << " npts = " << npts)
int nsp = el.groups();
int ng = nsp*nsp;
Matrix<int> histogram;
histogram.resize(grid.size(),ng);
d_table = DistanceTable::getTable(DistanceTable::add(el));
d_table->create(1);
int count = 0;
XMLReport("List of the configuration files:")
for(int i=0; i<ConfigFile.size(); i++)
XMLReport(ConfigFile[i]);
for(int i=0; i<ConfigFile.size(); i++)
{
HDFWalkerInput wReader(ConfigFile[i]);
while(wReader.put(el))
{
for(MCWalkerConfiguration::iterator it = el.begin();
it != el.end(); ++it)
{
el.R = (*it)->R;
DistanceTable::update(el);
for(int i=0; i<d_table->size(SourceIndex); i++)
{
for(int nn=d_table->M[i]; nn<d_table->M[i+1]; nn++)
{
double r = d_table->r(nn);
if(r < rcut)
{
count++;
int index = grid.index(r);
int id = d_table->PairID[nn];
histogram(index,id)++;
}
}
}
}
}
}
std::string froot(myProject.CurrentRoot());
froot.append(".pc");
std::ofstream pc_out(froot.c_str());
double vol = 4.0*4.0*atan(1.0)/3.0*pow(rcut,3);
double norm = static_cast<double>(count*factor)/(vol);
for(int ng=0; ng<histogram.rows()-1; ng++)
{
double r1 = grid[ng];
//double r2 = (i<(grid.size()-1)) ? grid[i+1]:(2.0*grid[i]-grid[i-1]);
double r2 = grid[ng+1];
// double r = 0.5*(r1+r2);
double r = 0.75*(pow(r2,4)-pow(r1,4))/(pow(r2,3)-pow(r1,3));
pc_out << std::setw(20) << r;
for(int j=0; j<histogram.cols(); j++)
{
//volume element dr = 4/3\pi[(r+dr)^3-r^3]
//double binVol = 4.0*4.0*atan(1.0)/3.0*(pow(r2,3)-pow(r1,3));
double binVol = 1.0;
pc_out << std::setw(20) << static_cast<double>(histogram(ng,j))/(binVol*norm);
}
pc_out << std::endl;
}
std::cout << "Ionic configuration : " << ion.getName() << std::endl;
ion.get(std::cout);
std::cout << "Electronic configuration : " << el.getName() << std::endl;
el.get(std::cout);
xmlXPathFreeContext(m_context);
xmlFreeDoc(m_doc);
OHMMS::Controller->finalize();
return 0;
}

212
src/QMCTools/PHDFWalkerMerger.cpp
View File

@ -1,212 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <iostream>
#include "QMCTools/HDFWalkerMerger.h"
#include "Utilities/RandomGeneratorIO.h"
#include "Utilities/UtilityFunctions.h"
#include "Message/Communicate.h"
#include "Message/CommOperators.h"
using namespace qmcplusplus;
HDFWalkerMerger::HDFWalkerMerger(const std::string& aroot, int ncopy)
: NumCopy(ncopy), FileRoot(aroot)
{
}
HDFWalkerMerger::~HDFWalkerMerger() { }
void HDFWalkerMerger::setCommunicator(Communicate* c)
{
myComm=c;
}
void HDFWalkerMerger::init()
{
char fname[128];
char GrpName[128];
int nprocs=myComm->ncontexts();
int mynode=myComm->mycontext();
hsize_t dimin[3], summary_size=3;
std::vector<int> nw_pernode(nprocs,0);
std::vector<int> nc_offset(nprocs+1,0);
FairDivideLow(NumCopy,nprocs,nc_offset);
std::cout << "Number of copies " << NumCopy << std::endl;
std::cout << "Number of proces " << nprocs << std::endl;
for(int i=0; i<=nprocs; i++)
std::cout << nc_offset[i] << std::endl;
numWalkersIn.push_back(new std::vector<hsize_t>);
//determine the minimum configuration and walkers for each config
int min_config=1;
for(int ip=nc_offset[mynode]; ip<nc_offset[mynode+1]; ip++)
{
int nconf=1;
sprintf(fname,"%s.p%03d.config.h5",FileRoot.c_str(),ip);
hid_t hfile = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
hid_t h_config = H5Gopen(hfile,"config_collection");
hid_t h1=H5Dopen(h_config,"NumOfConfigurations");
if(h1>-1)
{
H5Dread(h1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&(nconf));
H5Dclose(h1);
}
//iconf is the last one
int iconf=nconf-1;
sprintf(GrpName,"config%04d/coord",iconf);
h1 = H5Dopen(h_config,GrpName);
hid_t dataspace = H5Dget_space(h1);
int rank = H5Sget_simple_extent_ndims(dataspace);
int status_n = H5Sget_simple_extent_dims(dataspace, dimin, NULL);
H5Sclose(dataspace);
H5Dclose(h1);
nw_pernode[mynode]+=dimin[0];
numWalkersIn[0]->push_back(dimin[0]);
H5Fclose(hfile);
}
myComm->allreduce(nw_pernode);
NumPtcl = dimin[1];
Dimension = dimin[2];
OffSet.resize(2,nprocs+1);
OffSet=0;
int ip_in=0;
for(int ip=0; ip<nprocs; ip++)
{
OffSet(0,ip+1)=OffSet(0,ip)+nw_pernode[ip];
if(mynode ==0)
{
std::cout << ip << " has " << nw_pernode[ip] << std::endl;
}
}
for(int ip=0; ip<=nprocs; ip++)
OffSet(1,ip)=nc_offset[ip];
OffSet(nprocs,1)=nc_offset[nprocs];
MaxNumWalkers=nw_pernode[mynode];
NumConfig=1;
}
void HDFWalkerMerger::writeHeader(hid_t gid)
{
int cur_version=1;
hsize_t dim=1;
hid_t header_space= H5Screate_simple(1, &dim, NULL);
hid_t header_set= H5Dcreate(gid, "version", H5T_NATIVE_INT, header_space, H5P_DEFAULT);
H5Dwrite(header_set, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&cur_version);
H5Dclose(header_set);
header_set= H5Dcreate(gid, "ncontexts", H5T_NATIVE_INT, header_space, H5P_DEFAULT);
H5Dwrite(header_set, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&NumCopy);
H5Dclose(header_set);
H5Sclose(header_space);
}
void HDFWalkerMerger::merge()
{
init();
char fname[128];
char GrpName[128];
char ConfName[128];
#if H5_VERS_RELEASE < 4
hssize_t offset[]= {0,0,0};
#else
hsize_t offset[]= {0,0,0};
#endif
typedef std::vector<OHMMS_PRECISION> Container_t;
Container_t vecLoc(MaxNumWalkers*NumPtcl*Dimension);
hsize_t dimout[3], dimin[3];
dimin[1]=NumPtcl;
dimin[2]=Dimension;
dimout[1]=NumPtcl;
dimout[2]=Dimension;
////create new dataspace for each config
//for(int iconf=0; iconf<NumConfig; iconf++) {
// dimout[0]=OffSet(NumCopy,iconf);
// sprintf(ConfName,"config%04d",iconf);
// hid_t cf_id = H5Gcreate(mastercf,ConfName,0);
// dataspace = H5Screate_simple(3,dimout,NULL);
// dataset = H5Dcreate(cf_id,"coord",H5T_NATIVE_DOUBLE, dataspace,H5P_DEFAULT);
// H5Dclose(dataset);
// H5Sclose(dataspace);
// H5Gclose(cf_id);
//}
int nprocs=myComm->ncontexts();
int mynode=myComm->mycontext();
int ndacc=0; //local data count
for(int ip=OffSet(1,mynode),ipL=0; ip<OffSet(1,mynode+1); ip++,ipL++)
{
sprintf(fname,"%s.p%03d.config.h5",FileRoot.c_str(),ip);
hid_t h0 = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
hid_t cfcollect = H5Gopen(h0,"config_collection");
//get the number of walkers for (ip,iconf)
dimin[0] = (*numWalkersIn[0])[ipL];
sprintf(ConfName,"config%04d",0);
std::cout << "Node " << mynode << " reads " << dimin[0] << " walkers from " << fname << std::endl;
int nd=dimin[0]*NumPtcl*Dimension;
hid_t cf_id = H5Gopen(cfcollect,ConfName);
Container_t vecIn(nd);
HDFAttribIO<Container_t> in(vecIn, 3, dimin);
in.read(cf_id,"coord");
H5Gclose(cf_id);
copy(vecIn.begin(),vecIn.begin()+nd,vecLoc.begin()+ndacc);
ndacc+=nd;
H5Fclose(h0);
}
MPI_Info info=MPI_INFO_NULL;
hid_t acc_tpl1=H5Pcreate(H5P_FILE_ACCESS);
herr_t ret=H5Pset_fapl_mpio(acc_tpl1,myComm->getMPI(),info);
sprintf(fname,"%s.config.h5",FileRoot.c_str());
hid_t masterfile= H5Fcreate(fname,H5F_ACC_TRUNC, H5P_DEFAULT,acc_tpl1);
ret=H5Pclose(acc_tpl1);
dimout[0]=OffSet(0,nprocs);
hid_t mastercf = H5Gcreate(masterfile,"config_collection",0);
hid_t c1 = H5Gcreate(mastercf,"config0000",0);
//create global dimensionality
hid_t sid1=H5Screate_simple(3,dimout,NULL);
//create dataset collectively
hid_t dataset1=H5Dcreate(c1,"coord",H5T_NATIVE_DOUBLE,sid1,H5P_DEFAULT);
hsize_t start[3],count[3],stride[3]= {1,1,1};
start[0]=OffSet(0,mynode);
start[1]=0;
start[2]=0;
count[0]=OffSet(0,mynode+1)-OffSet(0,mynode);
count[1]=NumPtcl;
count[2]=Dimension;
//get the database
hid_t file_dataspace=H5Dget_space(dataset1);
ret=H5Sselect_hyperslab(file_dataspace,H5S_SELECT_SET,start,stride,count,NULL);
//create data independently
hid_t mem_dataspace=H5Screate_simple(3,count,NULL);
ret=H5Dwrite(dataset1,H5T_NATIVE_DOUBLE,mem_dataspace,file_dataspace,H5P_DEFAULT, &vecLoc[0]);
H5Sclose(file_dataspace);
ret = H5Dclose(dataset1);
H5Sclose(sid1);
H5Gclose(c1);
H5Gclose(mastercf);
H5Fclose(masterfile);
////write the summaryset
//double d=1.0/static_cast<double>(NumCopy);
//for(int k=0; k<Summary.size(); k++) Summary[k]*=d;
//hsize_t dim=Summary.size();
//hid_t summary_space = H5Screate_simple(1, &dim, NULL);
//hid_t summary_set = H5Dcreate(mastercf, "Summary", H5T_NATIVE_DOUBLE, summary_space, H5P_DEFAULT);
//ret = H5Dwrite(summary_set, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,&Summary[0]);
//H5Dclose(summary_set);
//H5Sclose(summary_space);
//// H5Gclose(mastercf);
//H5Fclose(masterfile);
}

137
src/QMCTools/QMCUtilities.h
View File

@ -1,137 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_SIMPLE_UTILITIES_H
#define QMCPLUSPLUS_SIMPLE_UTILITIES_H
#include "QMCWaveFunctions/OrbitalBuilderBase.h"
namespace qmcplusplus
{
/** A convenient function to evaluate the electron configuration
*@param el the electron configuration
*@param acontext xpath context
*/
inline
bool
determineNumOfElectrons(ParticleSet& el, xmlXPathContextPtr acontext)
{
//initialize el with the wave function information
//This is to determine the number of up and down electrons
std::vector<int> N;
std::string sdname("//");
sdname.append(OrbitalBuilderBase::sd_tag);
xmlXPathObjectPtr result
= xmlXPathEvalExpression((const xmlChar*)(sdname.c_str()),acontext);
bool found_el=false;
int nsd= result->nodesetval->nodeNr;
XMLReport("Found " << nsd << " SlaterDeterminant.")
if(nsd)
{
std::vector<xmlNodePtr> dset;
xmlNodePtr cur=result->nodesetval->nodeTab[0]->children;
while(cur != NULL)
{
std::string cname((const char*)(cur->name));
if(cname == OrbitalBuilderBase::det_tag)
dset.push_back(cur);
cur=cur->next;
}
if(dset.size())
{
XMLReport("Found " << dset.size() << OrbitalBuilderBase::det_tag)
for(int i=0; i<dset.size(); i++)
{
int norb = 0;
xmlChar* orb=xmlGetProp(dset[i],(const xmlChar*)"orbitals");
if(orb)
{
norb = atoi((const char*)orb);
XMLReport("Using attribute orbitals " << norb)
}
else
{
cur = dset[i]->children;
while(cur != NULL)
{
std::string cname((const char*)(cur->name));
if(cname == OrbitalBuilderBase::spo_tag)
norb++;
cur=cur->next;
}
XMLReport("Counting the number or ortbials " << norb)
}
N.push_back(norb);
}
el.create(N);
found_el=true;
}
}
xmlXPathFreeObject(result);
if(N.empty())
{
result = xmlXPathEvalExpression((const xmlChar*)"//particleset",acontext);
int n= result->nodesetval->nodeNr;
int pset=0;
while(!found_el && pset<n)
{
xmlChar* s= xmlGetProp(result->nodesetval->nodeTab[pset],(const xmlChar*)"name");
if(s)
{
if(s[0] == 'e')
{
xmlNodePtr cur = result->nodesetval->nodeTab[pset]->children;
found_el=true;
while(cur != NULL)
{
std::string cname((const char*)(cur->name));
if(cname == "group")
{
xmlChar* g= xmlGetProp(cur,(const xmlChar*)"size");
if(s)
{
N.push_back(atoi((const char*)g));
}
}
cur = cur->next;
}
}
}
pset++;
}
xmlXPathFreeObject(result);
}
if(N.empty())
{
ERRORMSG("Could not determine the number of electrons. Assuming He(1,1)")
N.resize(2);
N[0]=1;
N[1]=1;
}
else
{
XMLReport("The electron configured by //slaterdeterminant or //particle/group/@name=\'e\'")
}
el.create(N);
return true;
}
}
#endif

1
src/QMCTools/RGFBuilderBase.cpp
View File

@ -20,7 +20,6 @@
#if defined(HAVE_LIBHDF5) #if defined(HAVE_LIBHDF5)
#include "Numerics/HDFNumericAttrib.h" #include "Numerics/HDFNumericAttrib.h"
#endif #endif
#include "QMCTools/GridMolecularOrbitals.h"
#include "QMCTools/RGFBuilderBase.h" #include "QMCTools/RGFBuilderBase.h"
#include "QMCFactory/OneDimGridFactory.h" #include "QMCFactory/OneDimGridFactory.h"
namespace qmcplusplus namespace qmcplusplus

83
src/QMCTools/STO2GridBuilder.cpp
View File

@ -1,83 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Utilities/OhmmsInfo.h"
#include "Numerics/LibxmlNumericIO.h"
#include "Numerics/SlaterTypeOrbital.h"
#include "Numerics/Transform2GridFunctor.h"
#include "Numerics/OneDimCubicSpline.h"
#include "QMCWaveFunctions/MolecularOrbitals/GridMolecularOrbitals.h"
#include "QMCWaveFunctions/MolecularOrbitals/STO2GridBuilder.h"
namespace qmcplusplus
{
bool
STO2GridBuilder::putCommon(xmlNodePtr cur)
{
return true;
}
/** Add a new Slater Type Orbital with quantum numbers \f$(n,l,m,s)\f$
* \param cur the current xmlNode to be processed
* \param nlms a vector containing the quantum numbers \f$(n,l,m,s)\f$
* \return true is succeeds
*
This function puts the STO on a logarithmic grid and calculates the boundary
conditions for the 1D Cubic Spline. The derivates at the endpoint
are assumed to be all zero. Note: for the radial orbital we use
\f[ f(r) = \frac{R(r)}{r^l}, \f] where \f$ R(r) \f$ is the usual
radial orbital and \f$ l \f$ is the angular momentum.
*/
bool
STO2GridBuilder::addRadialOrbital(xmlNodePtr cur,
const QuantumNumberType& nlms)
{
if(!m_orbitals)
{
ERRORMSG("m_orbitals, SphericalOrbitals<ROT,GT>*, is not initialized")
return false;
}
RadialOrbitalType *radorb = NULL;
int n=nlms[0];
int l=nlms[1];
RealType zeta = 1.0;
xmlNodePtr s = cur->xmlChildrenNode;
while(s != NULL)
{
std::string cname((const char*)(s->name));
if(cname == "parameter" || cname =="Var")
{
putContent(zeta,s);
}
s=s->next;
}
XMLReport("Zeta = " << zeta)
STONorm<RealType> anorm(n);
GenericSTO<RealType> sto(n-l-1,zeta,anorm(n-1,zeta));
XMLReport("Calculating 1D-Cubic spline.")
//pointer to the grid
GridType* agrid = m_orbitals->Grids[0];
radorb = new OneDimCubicSpline<RealType>(agrid);
//spline the slater type orbital
Transform2GridFunctor<GenericSTO<RealType>,RadialOrbitalType> transform(sto, *radorb);
transform.generate(agrid->rmin(), agrid->rmax(),agrid->size());
//add the radial orbital to the list
m_orbitals->Rnl.push_back(radorb);
m_orbitals->RnlID.push_back(nlms);
return true;
}
}

41
src/QMCTools/STO2GridBuilder.h
View File

@ -1,41 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#ifndef QMCPLUSPLUS_STO2RADIALGRIDFUNCTOR_H
#define QMCPLUSPLUS_STO2RADIALGRIDFUNCTOR_H
#include "QMCWaveFunctions/MolecularOrbitals/RGFBuilderBase.h"
namespace qmcplusplus
{
/**Class to convert SlaterTypeOrbital to a radial orbital on a log grid.
*
* For a center,
* - only one grid is used
* - any number of radial orbitals
*/
struct STO2GridBuilder: public RGFBuilderBase
{
///constructor
STO2GridBuilder() {}
bool addRadialOrbital(xmlNodePtr cur, const QuantumNumberType& nlms);
bool putCommon(xmlNodePtr cur);
};
}
#endif

179
src/QMCTools/densityplot.cpp
View File

@ -1,179 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <sstream>
#include "OhmmsData/OhmmsElementBase.h"
#include "OhmmsData/HDFAttribIO.h"
#include "OhmmsPETE/TinyVector.h"
#include "OhmmsPETE/OhmmsMatrix.h"
using namespace qmcplusplus;
int print_help()
{
std::cout << "Usage: density [--fileroot std::string | --diff root1 root2 ] --np int " << std::endl;
return 1;
}
struct DensityObserver
{
typedef TinyVector<double,3> PosType;
int npts;
int nsamples;
double rmin;
double rmax;
double delta;
double deltainv;
Matrix<double> value;
DensityObserver():rmin(-7.92),rmax(7.92),delta(0.16),nsamples(0)
{
deltainv=1.0/delta;
npts=static_cast<int>((rmax-rmin)/delta)+1;
value.resize(npts,npts);
value=0.0;
std::cout << " Density grid " << npts << std::endl;
}
inline void accumulate(const std::vector<PosType>& pos, int nat)
{
nsamples+=pos.size();
for(int i=0; i<nat; i++)
{
double x=pos[i][0];
double y=pos[i][1];
if(x>rmin && x<rmax && y >rmin && y<rmax)
{
value(static_cast<int>((x-rmin)*deltainv),static_cast<int>((y-rmin)*deltainv))+=1.0;
}
}
}
void getData(const char* fname, int nproc);
void print(const char* fname)
{
std::ofstream fout(fname);
fout.setf(std::ios::scientific, std::ios::floatfield);
fout.setf(std::ios::left,std::ios::adjustfield);
fout.precision(6);
double norm=1.0/static_cast<double>(nsamples);
double x=rmin;
for(int ix=0; ix<npts; ix++,x+=delta)
{
double y=rmin;
for(int jx=0; jx<npts; jx++,y+=delta)
fout << std::setw(20) << y << std::setw(20) << x << std::setw(20) << norm*value(ix,jx) << std::endl;
fout << std::endl;
}
}
};
int main(int argc, char** argv)
{
if(argc<2)
return print_help();
int iargc=0;
int nproc=1;
std::vector<std::string> fnamein(2);
std::string h5fileroot("0");
bool getdiff=false;
while(iargc<argc)
{
std::string c(argv[iargc]);
if(c == "--fileroot")
{
h5fileroot=argv[++iargc];
}
else
if(c == "--diff")
{
fnamein[0]=argv[++iargc];
fnamein[1]=argv[++iargc];
getdiff=true;
}
else
if(c == "--np")
{
nproc=atoi(argv[++iargc]);
}
else
if(c == "--help" || c == "-h")
{
return print_help();
}
++iargc;
}
if(getdiff)
{
std::cout << "Going to compare two density " << fnamein[0] << " " << fnamein[1] << std::endl;
DensityObserver recorder1;
recorder1.getData(fnamein[0].c_str(),nproc);
DensityObserver recorder2;
recorder2.getData(fnamein[1].c_str(),nproc);
recorder1.value -= recorder2.value;
recorder1.print("rho_diff.dat");
}
else
{
DensityObserver recorder;
recorder.getData(h5fileroot.c_str(),nproc);
}
return 0;
}
void DensityObserver::getData(const char* froot, int nproc)
{
typedef TinyVector<double,3> PosType;
hsize_t dimIn[3],dimTot[3];
char fname[256],coordname[128];
for(int ip=0; ip<nproc; ip++)
{
if(nproc>1)
sprintf(fname,"%s.p%03d.config.h5",froot,ip);
else
sprintf(fname,"%s.config.h5",froot);
hid_t fileID = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
int nconf=1;
hid_t mastercf = H5Gopen(fileID,"config_collection");
hid_t h1=H5Dopen(mastercf,"NumOfConfigurations");
H5Dread(h1, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&(nconf));
H5Dclose(h1);
std::cout << "Adding data in " << fname << " number of configurations = " << nconf << std::endl;
for(int iconf=0; iconf<nconf; iconf++)
{
sprintf(coordname,"config%04d/coord",iconf);
hid_t dataset = H5Dopen(mastercf,coordname);
hid_t dataspace = H5Dget_space(dataset);
int rank = H5Sget_simple_extent_ndims(dataspace);
int status_n = H5Sget_simple_extent_dims(dataspace, dimTot, NULL);
std::vector<PosType> pos(dimTot[0]*dimTot[1]);
hid_t ret = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(pos[0][0]));
H5Dclose(dataset);
H5Sclose(dataspace);
std::accumulate(pos,dimTot[0]*dimTot[1]);
}
}
std::string outfname(froot);
outfname.append(".den.dat");
print(outfname.c_str());
}

174
src/QMCTools/ecp4qmc.cpp
View File

@ -1,174 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <fstream>
#include "Message/Communicate.h"
#include "Utilities/OhmmsInfo.h"
#include "Numerics/GaussianTimesRN.h"
#include "Numerics/OneDimGridBase.h"
#include "Numerics/OneDimGridFunctor.h"
#include "Numerics/OneDimCubicSpline.h"
#include "QMCFactory/OneDimGridFactory.h"
using namespace qmcplusplus;
struct ECPTest
{
typedef OneDimGridFactory::RealType RealType;
typedef OneDimGridFactory::GridType GridType;
typedef GaussianTimesRN<RealType> InFuncType;
typedef OneDimCubicSpline<RealType> OutFuncType;
std::string Name;
RealType Zeff;
ECPTest(const std::string& fname);
void buildSemiLocal(xmlNodePtr cur);
void buildLocal(xmlNodePtr cur);
};
int main(int argc, char** argv)
{
OHMMS::Controller->initialize(argc,argv);
OhmmsInfo welcome(argc,argv,OHMMS::Controller->mycontext());
int iargc=0;
ECPTest ecp(argv[1]);
return 0;
}
ECPTest::ECPTest(const std::string& fname)
{
// build an XML tree from a the file;
xmlDocPtr m_doc = xmlParseFile(fname.c_str());
if (m_doc == NULL)
{
ERRORMSG("File " << fname << " is invalid")
xmlFreeDoc(m_doc);
}
// Check the document is of the right kind
xmlNodePtr cur = xmlDocGetRootElement(m_doc);
if (cur == NULL)
{
ERRORMSG("Empty document");
xmlFreeDoc(m_doc);
}
cur=cur->children;
while(cur != NULL)
{
std::string cname((const char*)cur->name);
if(cname == "header")
{
Zeff = atoi((const char*)xmlGetProp(cur,(const xmlChar*)"zval"));
Name = (const char*)xmlGetProp(cur,(const xmlChar*)"symbol");
}
else
if(cname == "semilocal")
{
buildSemiLocal(cur);
}
else
if(cname == "local")
{
buildLocal(cur);
}
cur=cur->next;
}
xmlFreeDoc(m_doc);
}
void ECPTest::buildSemiLocal(xmlNodePtr cur)
{
GridType* grid_semilocal=0;
std::vector<InFuncType*> semilocal;
cur=cur->children;
while(cur != NULL)
{
std::string cname((const char*)cur->name);
if(cname == "grid")
{
grid_semilocal=OneDimGridFactory::createGrid(cur);
}
else
if(cname == "vps")
{
xmlNodePtr cur1=cur->children;
while(cur1 != NULL)
{
if(xmlStrEqual(cur1->name,(const xmlChar*)"basisGroup"))
{
InFuncType* a=new InFuncType;
a->putBasisGroup(cur1);
semilocal.push_back(a);
}
cur1=cur1->next;
}
}
cur=cur->next;
}
std::string fname(Name);
fname.append(".semilocal.dat");
std::ofstream fout(fname.c_str());
fout.setf(std::ios::scientific, std::ios::floatfield);
fout.precision(12);
int ig=0;
while(ig<grid_semilocal->size())
{
double r=(*grid_semilocal)[ig++];
fout << std::setw(20) << setprecision(12) << r;
for(int i=0; i<semilocal.size(); i++)
{
fout << std::setw(20) << semilocal[i]->f(r);
}
fout << std::endl;
}
}
void ECPTest::buildLocal(xmlNodePtr cur)
{
GridType* grid_semilocal=0;
InFuncType* localpp;
cur=cur->children;
while(cur != NULL)
{
std::string cname((const char*)cur->name);
if(cname == "grid")
{
grid_semilocal=OneDimGridFactory::createGrid(cur);
}
else
if(cname == "basisGroup")
{
localpp=new InFuncType;
localpp->putBasisGroup(cur);
}
cur=cur->next;
}
std::cout << " Effective Z = " << Zeff << std::endl;
std::string fname(Name);
fname.append(".local.dat");
std::ofstream fout(fname.c_str());
fout.setf(std::ios::scientific, std::ios::floatfield);
fout.precision(12);
int ig=0;
while(ig<grid_semilocal->size())
{
double r=(*grid_semilocal)[ig++];
fout << std::setw(20) << setprecision(12) << r << std::setw(20) << localpp->f(r)-Zeff/r<< std::endl;
}
delete localpp;
}

75
src/QMCTools/gemtest.cpp
View File

@ -1,75 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <iomanip>
#include "Numerics/MatrixOperators.h"
#include "Utilities/RandomGenerator.h"
#include "Utilities/Timer.h"
using namespace qmcplusplus;
int main(int argc, char** argv)
{
Random.init(0,1,-1);
typedef double scalar_t;
int N=100;
int M=2000;
int L=5;
Matrix<scalar_t> A(N,M), B(M,L), C(N,L);
Vector<scalar_t> y(M),z(N);
double gemm_t=0.0;
double gemv_t=0.0;
Timer myTimer;
int niter=100;
for(int iter=0; iter<niter; iter++)
{
for(int i=0; i<N; i++)
for(int j=0; j<M; j++)
A(i,j)=Random();
//A(i,j)=scalar_t(Random(),Random());
for(int i=0; i<M; i++)
for(int j=0; j<L; j++)
B(i,j)=Random();
//B(i,j)=scalar_t(Random(),Random());
for(int i=0; i<M; i++)
y(i)=Random();
//y(i)=scalar_t(Random(),Random());
myTimer.restart();
MatrixOperators::product(A,B,C);
gemm_t+=myTimer.elapsed();
myTimer.restart();
MatrixOperators::product(A,y,z.data());
gemv_t+=myTimer.elapsed();
// std::cout << "<<<<<< GEMM TEST " << std::endl;
// for(int i=0; i<N; i++) {
// for(int j=0; j<L; j++) {
// scalar_t v=0.0;
// for(int k=0; k<M; k++) v+=A(i,k)*B(k,j);
// std::cout << i << "," << j << " " << v-C(i,j) << " " << v << std::endl;
// }
// }
// std::cout << "<<<<<< GEMV TEST " << std::endl;
// for(int i=0; i<N; i++) {
// scalar_t v=0.0;
// for(int j=0; j<M; j++) v+=A(i,j)*y(j);
// std::cout << i << " " << v-z[i] << " " << v << std::endl;
// }
}
std::cout << "|"<<N << "|"<<M << "|" << L << "|" << gemm_t/niter << "|" << gemv_t/niter<< "|" << std::endl;
}

178
src/QMCTools/gofr.cpp
View File

@ -1,178 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <sstream>
#include "Configuration.h"
#include "Numerics/HDFNumericAttrib.h"
#include "Utilities/IteratorUtility.h"
using namespace qmcplusplus;
int print_help()
{
std::cout << "Usage: gofr --fileroot std::string --np int " << std::endl;
return 1;
}
struct GofRObserver
{
int NumSamples;
int NumNodes;
std::string DataSetName;
Vector<double> gofr;
Vector<double> gofr2;
Matrix<double> gofr_dat;
GofRObserver(const std::string& aname):DataSetName(aname),NumSamples(0)
{
}
inline void accumulate(int count)
{
int nbin=gofr_dat.cols();
for(int i=0; i<count; i++)
{
double* restrict dptr=gofr_dat[i];
for(int k=0; k<nbin; k++,dptr++)
{
gofr[k] += *dptr;
gofr2[k] += (*dptr)*(*dptr);
}
}
NumSamples = count;
}
inline void resize(int ns, int nbin)
{
NumSamples=ns;
gofr_dat.resize(ns,nbin);
gofr_dat=0.0;
gofr.resize(nbin);
gofr2.resize(nbin);
gofr=0.0;
gofr2=0.0;
}
void getData(const char* fname, int nproc);
void print(const char* fname)
{
std::ofstream fout(fname);
fout.setf(std::ios::scientific, std::ios::floatfield);
fout.setf(std::ios::left,std::ios::adjustfield);
fout.precision(6);
//double norm=static_cast<double>(NumNodes)/static_cast<double>(NumSamples);
double norm=1.0/static_cast<double>(NumSamples);
double sqrtnorm=sqrt(norm);
for(int i=0; i<gofr.size(); i++)
{
double avg=gofr[i]*norm;
double var= sqrt(gofr2[i]*norm-avg*avg);
fout << std::setw(3) << i << std::setw(20) << avg << std::setw(20) << var*sqrtnorm << std::setw(20) << var << std::endl;
}
}
};
int main(int argc, char** argv)
{
if(argc<2)
return print_help();
int iargc=0;
int nproc=1;
std::vector<std::string> gofr_name;
std::string h5fileroot("0");
while(iargc<argc)
{
std::string c(argv[iargc]);
if(c == "--fileroot")
{
h5fileroot=argv[++iargc];
}
else
if(c == "--np" || c == "-np")
{
nproc=atoi(argv[++iargc]);
}
else
if(c == "--help" || c == "-h")
{
return print_help();
}
else
if(c == "--dataset" || c == "-d")
{
gofr_name.push_back(argv[++iargc]);
}
++iargc;
}
for(int i=0; i<gofr_name.size(); i++)
{
GofRObserver recorder(gofr_name[i]);
recorder.getData(h5fileroot.c_str(),nproc);
}
return 0;
}
void GofRObserver::getData(const char* froot, int nproc)
{
NumNodes=nproc;
char fname[128];
int count_max=1000000;
for(int ip=0; ip<nproc; ip++)
{
if(nproc>1)
sprintf(fname,"%s.p%03d.config.h5",froot,ip);
else
sprintf(fname,"%s.config.h5",froot);
std::cout << "Getting data from " << fname << std::endl;
hid_t fileid = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
hid_t obsid = H5Gopen(fileid,"observables");
int count=0;
HDFAttribIO<int> i(count);
i.read(obsid,"count");
count_max = std::min(count_max,count);
std::string gname(DataSetName);
gname.append("/value");
//check dimension
hid_t h1 = H5Dopen(obsid, gname.c_str());
hid_t dataspace = H5Dget_space(h1);
hsize_t dims_out[2];
int rank = H5Sget_simple_extent_ndims(dataspace);
int status_n = H5Sget_simple_extent_dims(dataspace, dims_out, NULL);
H5Dclose(h1);
if(gofr_dat.size() ==0)
{
resize(dims_out[0],dims_out[1]);
}
Matrix<double> tmp(dims_out[0],dims_out[1]);
HDFAttribIO<Matrix<double> > o(tmp);
o.read(obsid,gname.c_str());
gofr_dat += tmp;
H5Gclose(obsid);
H5Fclose(fileid);
}
std::cout << "Number of blocks " << count_max << std::endl;
std::accumulate(count_max);
std::string outfname(froot);
outfname += "."+ DataSetName+".dat";
print(outfname.c_str());
}

110
src/QMCTools/gto2grid.cpp
View File

@ -1,110 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Message/Communicate.h"
#include "Utilities/OhmmsInfo.h"
#include "QMCWaveFunctions/MolecularOrbitals/GTO2GridBuilder.h"
using namespace qmcplusplus;
void buildBasisSet(xmlNodePtr cur);
int main(int argc, char** argv)
{
OHMMS::Controller->initialize(argc,argv);
OhmmsInfo welcome(argc,argv,OHMMS::Controller->mycontext());
// build an XML tree from a the file;
xmlDocPtr m_doc = xmlParseFile(argv[1]);
if (m_doc == NULL)
{
ERRORMSG("File " << argv[1] << " is invalid")
xmlFreeDoc(m_doc);
return 1;
}
// Check the document is of the right kind
xmlNodePtr cur = xmlDocGetRootElement(m_doc);
if (cur == NULL)
{
ERRORMSG("Empty document");
xmlFreeDoc(m_doc);
return 1;
}
xmlXPathContextPtr m_context = xmlXPathNewContext(m_doc);
xmlXPathObjectPtr result
= xmlXPathEvalExpression((const xmlChar*)"//atomicBasisSet",m_context);
if(!xmlXPathNodeSetIsEmpty(result->nodesetval))
{
for(int ic=0; ic<result->nodesetval->nodeNr; ic++)
{
buildBasisSet(result->nodesetval->nodeTab[ic]);
}
}
xmlXPathFreeObject(result);
return 0;
}
void buildBasisSet(xmlNodePtr cur)
{
xmlNodePtr anchor = cur;
xmlNodePtr grid_ptr = 0;
std::vector<xmlNodePtr> phi_ptr;
std::vector<QuantumNumberType> nlms;
int Lmax = 0;
int current = 0;
std::string acenter("none");
const xmlChar* aptr = xmlGetProp(cur,(const xmlChar*)"species");
if(aptr)
acenter = (const char*)aptr;
std::cout << "Building basis set for " << acenter << std::endl;
cur = anchor->children;
while(cur != NULL)
{
std::string cname((const char*)(cur->name));
if(cname == "grid")
grid_ptr = cur;
else
if(cname == "basisSet")
{
phi_ptr.push_back(cur);
nlms.push_back(QuantumNumberType());
int n=1,l=0,m=0;
const xmlChar* lptr = xmlGetProp(cur,(const xmlChar*)"l");
if(lptr)
l = atoi((const char*)lptr);
const xmlChar* nptr = xmlGetProp(cur,(const xmlChar*)"n");
if(nptr)
n = atoi((const char*)nptr);
const xmlChar* mptr = xmlGetProp(cur,(const xmlChar*)"m");
if(mptr)
m = atoi((const char*)mptr);
Lmax = std::max(l,Lmax);
nlms[current][0]=n;
nlms[current][1]=l;
nlms[current][2]=m;
current++;
}
cur = cur->next;
}
RGFBuilderBase::CenteredOrbitalType aos(Lmax);
RGFBuilderBase* rbuilder = new GTO2GridBuilder;
rbuilder->setOrbitalSet(&aos,acenter);
rbuilder->addGrid(grid_ptr);
for(int i=0; i<nlms.size(); i++)
{
rbuilder->addRadialOrbital(phi_ptr[i],nlms[i]);
}
rbuilder->print(acenter,1);
}

172
src/QMCTools/gto2slater.cpp
View File

@ -1,172 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Message/Communicate.h"
#include "Utilities/OhmmsInfo.h"
#include "Numerics/OneDimGridBase.h"
#include "Numerics/OneDimCubicSpline.h"
#include "Numerics/Transform2GridFunctor.h"
#include "Numerics/GaussianBasisSet.h"
#include "QMCTools/Any2Slater.h"
struct GTO2Slater
{
typedef GaussianCombo<double> GTOType;
typedef LogGrid<double> GridType;
bool Normalized;
GridType myGrid;
xmlNodePtr gridPtr;
std::map<std::string,xmlNodePtr> sPtr;
GTO2Slater():Normalized(false),gridPtr(0) {}
int parse(const char* fname);
bool put(xmlNodePtr);
void optimize();
};
int main(int argc, char** argv)
{
OHMMS::Controller->initialize(argc,argv);
OhmmsInfo welcome(argc,argv,OHMMS::Controller->mycontext());
GTO2Slater transformer;
transformer.parse(argv[1]);
return 0;
}
int GTO2Slater::parse(const char* fname)
{
// build an XML tree from a the file;
xmlDocPtr m_doc = xmlParseFile(fname);
if (m_doc == NULL)
{
ERRORMSG("File " << fname << " is invalid")
xmlFreeDoc(m_doc);
return 1;
}
// Check the document is of the right kind
xmlNodePtr cur = xmlDocGetRootElement(m_doc);
if (cur == NULL)
{
ERRORMSG("Empty document");
xmlFreeDoc(m_doc);
return 1;
}
xmlXPathContextPtr m_context = xmlXPathNewContext(m_doc);
xmlXPathObjectPtr result
= xmlXPathEvalExpression((const xmlChar*)"//atomicBasisSet",m_context);
if(!xmlXPathNodeSetIsEmpty(result->nodesetval))
{
for(int ic=0; ic<result->nodesetval->nodeNr; ic++)
{
std::cout << "Going to optimize" << std::endl;
if(put(result->nodesetval->nodeTab[ic]))
{
optimize();
}
}
}
xmlXPathFreeObject(result);
return 1;
}
bool GTO2Slater::put(xmlNodePtr cur)
{
cur = cur->children;
while(cur != NULL)
{
std::string cname((const char*)(cur->name));
if(cname == "grid")
gridPtr = cur;
else
if(cname == "basisGroup")
{
std::string rid("invalid");
std::string rtype("Gaussian");
std::string norm("no");
int l=0;
OhmmsAttributeSet inAttrib;
inAttrib.add(rid,"rid");
inAttrib.add(l,"l");
inAttrib.add(rtype,"type");
inAttrib.add(norm,"normalized");
inAttrib.put(cur);
if(rtype == "Gaussian" && l == 0)
//pick only S
{
//if Ngto==1, don't do it
if(norm == "yes")
Normalized=true;
else
Normalized=false;
std::map<std::string,xmlNodePtr>::iterator it(sPtr.find(rid));
if(it == sPtr.end())
{
sPtr[rid]=cur;
}
}
}
cur=cur->next;
}
if(sPtr.empty())
return false;
return
true;
}
/** main functio to optimize multiple contracted S orbitals
*/
void GTO2Slater::optimize()
{
//construct one-dim grid
double ri = 1e-5;
double rf = 10.0;
int npts = 101;
std::string gridType("log");
if(gridPtr)
{
OhmmsAttributeSet radAttrib;
radAttrib.add(gridType,"type");
radAttrib.add(npts,"npts");
radAttrib.add(ri,"ri");
radAttrib.add(rf,"rf");
radAttrib.put(gridPtr);
}
myGrid.set(ri,rf,npts);
//create a numerical grid funtor
typedef OneDimCubicSpline<double> RadialOrbitalType;
RadialOrbitalType radorb(&myGrid);
int L= 0;
//Loop over all the constracted S orbitals
std::map<std::string,xmlNodePtr>::iterator it(sPtr.begin()),it_end(sPtr.end());
while(it != it_end)
{
//create contracted gaussian
GTOType gset(L,Normalized);
//read the radfunc's of basisGroup
gset.putBasisGroup((*it).second);
//convert to a radial functor
Transform2GridFunctor<GTOType,RadialOrbitalType> transform(gset, radorb);
transform.generate(myGrid.rmin(),myGrid.rmax(),myGrid.size());
//optimize it with the radial functor
Any2Slater gto2slater(radorb);
gto2slater.optimize();
++it;
}
sPtr.clear();
}

55
src/QMCTools/h5merge.cpp
View File

@ -1,55 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include "Utilities/RandomGenerator.h"
#include <vector>
#include "QMCTools/HDFWalkerMerger.h"
#include "Utilities/OhmmsInfo.h"
int main(int argc, char **argv)
{
OHMMS::Controller->initialize(argc,argv);
OhmmsInfo welcome(argc,argv,OHMMS::Controller->rank());
qmcplusplus::Random.init(0,1,-1);
if(argc<2)
{
std::cerr << " Usage: h5merge <rootname> -n <number-of-processor> -o[utfile] <outfile> " << std::endl;
return 1;
}
int ic=0;
int np=1;
std::string ofile(argv[1]);
while(ic<argc)
{
std::string w(argv[ic]);
if(w.find("-n")<w.size())
{
np=atoi(argv[++ic]);
}
else
if(w.find("-o")<w.size())
{
ofile=argv[++ic];
}
++ic;
}
std::cout << "Number of processors = " << np << std::endl;
//qmcplusplus::HDFWalkerMerger merger(argv[1],atoi(argv[2]));
//qmcplusplus::HDFWalkerMerger merger(argv[1],np);
//merger.merge();
return 0;
}

7317
src/QMCTools/makefun.f90
View File

File diff suppressed because it is too large Load Diff

61
src/QMCTools/numberOfOrbitalsForCI.pl
View File

@ -1,61 +0,0 @@
#//////////////////////////////////////////////////////////////////////////////////////
#// 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: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
#// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
#//
#//
#// File created by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
#//////////////////////////////////////////////////////////////////////////////////////
#! /usr/bin/perl
# only works for csf's right now
if($#ARGV != 1)
{
print "Usage: script filename cutoff \n";
exit(1)
}
$filename = $ARGV[0];
$dbl1 = $ARGV[1];
# check file exists
open(INN,"< $filename");
@alldata = <INN>;
close(INN);
$dcnt = 0;
$int1 = 0;
$ndet = 0;
while($dcnt < $#alldata+1)
{
if($alldata[$dcnt] =~ /<csf/)
{
$alldata[$dcnt] =~ /qchem_coeff.\"(\S*)\"/;
$somevar = $1;
#print "my num $somevar\n";
if(abs($somevar) > abs($dbl1))
{
$ndet++;
$alldata[$dcnt] =~ /occ.\"(.*)[12](0*)"/;
$mylength = length($1)+1;
#print "dol1 $1 dol2 $2 \n";
#print "my length is $mylength\n";
if($mylength > $int1)
{
$int1 = $mylength;
}
}
}
$dcnt++;
}
# must add number of core orbitals to this
print "Number of determinants: $ndet \n";
print "Required number of orbitals: $int1 \n";

164
src/QMCTools/observable.cpp
View File

@ -1,164 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <sstream>
#include <OhmmsData/OhmmsElementBase.h>
#include <Numerics/HDFSTLAttrib.h>
using namespace qmcplusplus;
int print_help()
{
std::cerr << "Usage: observable [--fileroot std::string ] list-of-obervables " << std::endl;
std::cerr << " Example: observable --fileroot myhdf --first 10 density " << std::endl;
return 1;
}
struct ObservableFile
{
int first_block;
std::string h5name;
std::vector<std::string> olist;
ObservableFile(int argc, char** argv);
int dump_averages();
template<typename IT>
inline void accumulate_n(IT first, IT last, IT result)
{
for(; first<last;)
*result++ += *first++;
}
template<typename IT, typename T>
inline void scale_n(IT first, IT last, T fac)
{
for(; first<last;)
*first++ *=fac;
}
};
int main(int argc, char** argv)
{
if(argc<2)
return print_help();
ObservableFile in(argc,argv);
return in.dump_averages();
}
ObservableFile::ObservableFile(int argc, char** argv)
:first_block(0)
{
int iargc=1;
while(iargc<argc)
{
std::string c(argv[iargc]);
if(c.find("fileroot")<c.size())
{
h5name=argv[++iargc];
}
else
if(c.find("first")<c.size())
{
first_block=atoi(argv[++iargc]);
}
else
if(c.find("help")<c.size())
{
print_help();
}
else
{
olist.push_back(argv[iargc]);
}
++iargc;
}
}
int ObservableFile::dump_averages()
{
if(h5name.empty())
{
std::cerr << " Invalid file. Provide --fileroot root-of-h5-file " << std::endl;
return 1;
}
herr_t status = H5Eset_auto(NULL, NULL);
char fname[256];
if(h5name.find("h5")<h5name.size())
sprintf(fname,"%s",h5name.c_str());
else
sprintf(fname,"%s.h5",h5name.c_str());
hid_t fileID = H5Fopen(fname,H5F_ACC_RDWR,H5P_DEFAULT);
if(fileID<0)
{
std::cerr << " Cannot find " << fname << std::endl;
return 1;
}
if(olist.empty())
{
std::cerr << "No observables given." << std::endl;
print_help();
return 1;
}
for(int i=0; i<olist.size(); ++i)
{
hid_t gid=H5Gopen(fileID,olist[i].c_str());
if(gid<0)
{
std::cerr << "Cannot find " << olist[i] << " group. Abort" << std::endl;
return 1;
}
hid_t dataset = H5Dopen(gid,"value");
hid_t dataspace = H5Dget_space(dataset);
//check the dimension of the active observable
int rank = H5Sget_simple_extent_ndims(dataspace);
std::vector<hsize_t> in_dims(rank);
int status_n = H5Sget_simple_extent_dims(dataspace, &in_dims[0], NULL);
std::vector<int> dims(rank-1);
int num_blocks=static_cast<int>(in_dims[0]);
int nb=num_blocks-first_block;
//skip it
if(nb<0)
continue;
std::cout << "<<<< Input dimension " << num_blocks;
int tot_size=1;
for(int dim=0; dim<rank-1; ++dim)
{
tot_size *= dims[dim]=static_cast<int>(in_dims[dim+1]);
std::cout << " "<< dims[dim] ;
}
std::cout << "\n Total size = " << tot_size << std::endl;
std::vector<double> odata_in;
odata_in.resize(tot_size*num_blocks);
hid_t ret = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(odata_in[0]));
H5Dclose(dataset);
H5Sclose(dataspace);
std::vector<double> odata_tot(tot_size,0.0), odata2_tot(tot_size,0.0);
std::vector<double>::iterator first=odata_in.begin()+first_block*tot_size;
for(int b=first_block; b<num_blocks; ++b,first+=tot_size)
accumulate_n(first,first+tot_size,odata_tot.begin());
double fac=1.0/static_cast<double>(nb);
scale_n(odata_tot.begin(),odata_tot.end(),fac);
HDFAttribIO<std::vector<double> > dummy(odata_tot,dims);
dummy.write(gid,"average");
}
return 0;
}

593
src/QMCTools/pw2qmcpack.qe3.0.f90
View File

@ -1,593 +0,0 @@
!//////////////////////////////////////////////////////////////////////////////////////
!// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
!// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
!//
!// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
!//////////////////////////////////////////////////////////////////////////////////////
!
! Copyright (C) 2004 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
PROGRAM pw2qmcpack
!-----------------------------------------------------------------------
! This subroutine writes the file "prefix".pwscf.xml and "prefix".pwscf.h5
! containing the plane wave coefficients and other stuff needed by the
! QMC code QMCPACK.
#include "f_defs.h"
USE io_files, ONLY : nd_nmbr, prefix, outdir, tmp_dir, trimcheck
USE io_global, ONLY : ionode, ionode_id
USE mp, ONLY : mp_bcast
!
IMPLICIT NONE
INTEGER :: ios
NAMELIST / inputpp / prefix, outdir
CALL start_postproc(nd_nmbr)
!
! set default values for variables in namelist
!
prefix = 'pwscf'
CALL get_env( 'ESPRESSO_TMPDIR', outdir )
IF ( TRIM( outdir ) == ' ' ) outdir = './'
IF ( ionode ) THEN
!
READ (5, inputpp, err=200, iostat=ios)
200 CALL errore('pw2qmcpack', 'reading inputpp namelist', ABS(ios))
tmp_dir = trimcheck (outdir)
!
END IF
!
! ... Broadcast variables
!
CALL mp_bcast( prefix, ionode_id )
CALL mp_bcast(tmp_dir, ionode_id )
!
CALL read_file
CALL openfil_pp
!
CALL compute_qmcpack
!
CALL stop_pp
STOP
END PROGRAM pw2qmcpack
SUBROUTINE compute_qmcpack
USE kinds, ONLY: DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, zv, atm
USE cell_base, ONLY: omega, alat, tpiba2, at, bg
USE char, ONLY: title
USE constants, ONLY: tpi
USE ener, ONLY: ewld, ehart, etxc, vtxc, etot, etxcc
USE gvect, ONLY: ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, g, gg, ecutwfc, gcutm, nl, igtongl
USE klist , ONLY: nks, nelec, nelup, neldw, xk
USE lsda_mod, ONLY: lsda, nspin
USE scf, ONLY: rho, rho_core, rhog, rhog_core
USE vlocal, ONLY: vloc, vnew, strf
USE wvfct, ONLY: npw, npwx, nbnd, gamma_only, igk, g2kin, wg, et
USE uspp, ONLY: nkb, vkb, dvan
USE uspp_param, ONLY: nh
USE becmod, ONLY: becp
USE io_global, ONLY: stdout
USE io_files, ONLY: nd_nmbr, nwordwfc, iunwfc, iun => iunsat, tmp_dir, prefix
USE wavefunctions_module, ONLY : evc, psic
use gsmooth, only: nls, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s
use iotk_module
use iotk_xtox_interf
IMPLICIT NONE
INTEGER :: ig, ibnd, ik, io, na, j, ispin, nbndup, nbnddown, &
nk, ngtot, ig7, ikk, nt, ijkb0, ikb, ih, jh, jkb, at_num, &
nelec_tot, nelec_up, nelec_down, ii, igx, igy, igz
INTEGER, ALLOCATABLE :: INDEX(:), igtog(:)
LOGICAL :: exst, found
REAL(DP) :: ek, eloc, enl, charge, etotefield
REAL(DP) :: bg_qmc(3,3), tau_r(3), g_red(3)
COMPLEX(DP), ALLOCATABLE :: phase(:),aux(:), hpsi(:,:), eigpacked(:)
COMPLEX(DP), ALLOCATABLE :: psitr(:)
REAL(DP), ALLOCATABLE :: g_cart(:,:),psireal(:)
INTEGER :: ios, ierr, h5len,ig_c,save_complex, nup,ndown
INTEGER, EXTERNAL :: atomic_number, is_complex
INTEGER, ALLOCATABLE:: g_qmc(:,:)
REAL (DP), EXTERNAL :: ewald
CHARACTER(256) :: tmp,h5name,eigname
CHARACTER(iotk_attlenx) :: attr
CALL init_us_1
CALL newd
!####io = 77
!####WRITE (6,'(/,5x,''Writing file pwfn.data for program CASINO'')')
!####CALL seqopn( 77, 'pwfn.data', 'formatted',exst)
!ALLOCATE (hpsi(npwx, nbnd))
ALLOCATE (aux(nrxx))
ALLOCATE (becp (nkb,nbnd))
! four times npwx should be enough
ALLOCATE (INDEX (4*npwx) )
ALLOCATE (igtog (4*npwx) )
!hpsi (:,:) = (0.d0, 0.d0)
INDEX(:) = 0
igtog(:) = 0
IF( lsda ) THEN
nbndup = nbnd
nbnddown = nbnd
nk = nks/2
! nspin = 2
ELSE
nbndup = nbnd
nbnddown = 0
nk = nks
! nspin = 1
ENDIF
! if(nks > 1) rewind(iunigk)
! do ik=1,nks
! if(nks > 1) read(iunigk) npw, igk
!
! if(nks > 1) rewind(iunigk)
ek = 0.d0
eloc= 0.d0
enl = 0.d0
!
DO ispin = 1, nspin
!
! calculate the local contribution to the total energy
!
! bring rho to G-space
!
aux(:) = CMPLX ( rho(:,ispin), 0.d0)
CALL cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
!
DO nt=1,ntyp
DO ig = 1, ngm
eloc = eloc + vloc(igtongl(ig),nt) * strf(ig,nt) &
* CONJG(aux(nl(ig)))
ENDDO
ENDDO
DO ik = 1, nk
ikk = ik + nk*(ispin-1)
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio (evc, nwordwfc, iunwfc, ikk, - 1)
CALL init_us_2 (npw, igk, xk (1, ikk), vkb)
CALL ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
DO ig =1, npw
IF( igk(ig) > 4*npwx ) &
CALL errore ('pw2qmcpack','increase allocation of index', ig)
INDEX( igk(ig) ) = 1
ENDDO
ENDDO
ENDDO
#ifdef __PARA
CALL reduce(1,eloc)
CALL reduce(1,ek)
CALL poolreduce(1,ek)
CALL poolreduce(1,enl)
#endif
eloc = eloc * omega
ek = ek * tpiba2
ngtot = 0
DO ig = 1, 4*npwx
IF( INDEX(ig) == 1 ) THEN
ngtot = ngtot + 1
igtog(ngtot) = ig
ENDIF
ENDDO
ALLOCATE (g_qmc(3,ngtot))
ALLOCATE (g_cart(3,ngtot))
! get the number of electrons
nelec_tot= NINT(nelec)
nup=NINT(nelup)
ndown=NINT(neldw)
if(nup .eq. 0) then
ndown=nelec_tot/2
nup=nelec_tot-ndown
endif
h5name = TRIM( prefix ) // '.pwscf.h5'
eigname = "eigenstates_"//trim(iotk_itoa(nr1s))//'_'//trim(iotk_itoa(nr2s))//'_'//trim(iotk_itoa(nr3s))
bg_qmc(:,:)=bg(:,:)/alat
!! create a file for particle set
tmp = TRIM( tmp_dir ) // TRIM( prefix )// '.ptcl.xml'
CALL iotk_open_write(iun, FILE=TRIM(tmp), ROOT="qmcsystem", IERR=ierr )
CALL iotk_write_attr (attr,"name","global",first=.true.)
CALL iotk_write_begin(iun, "simulationcell",ATTR=attr)
CALL iotk_write_attr (attr,"name","lattice",first=.true.)
CALL iotk_write_attr (attr,"units","bohr")
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,100) alat*at(1,1), alat*at(2,1), alat*at(3,1)
WRITE(iun,100) alat*at(1,2), alat*at(2,2), alat*at(3,2)
WRITE(iun,100) alat*at(1,3), alat*at(2,3), alat*at(3,3)
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_attr (attr,"name","reciprocal",first=.true.)
CALL iotk_write_attr (attr,"units","2pi/bohr")
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,100) bg_qmc(1,1), bg_qmc(2,1), bg_qmc(3,1)
WRITE(iun,100) bg_qmc(1,2), bg_qmc(2,2), bg_qmc(3,2)
WRITE(iun,100) bg_qmc(1,3), bg_qmc(2,3), bg_qmc(3,3)
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_attr (attr,"name","bconds",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,'(a)') 'p p p'
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_attr (attr,"name","LR_dim_cutoff",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,'(a)') '10'
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "simulationcell")
! <particleset name="ions">
CALL iotk_write_attr (attr,"name","ion0",first=.true.)
CALL iotk_write_attr (attr,"size",nat)
CALL iotk_write_begin(iun, "particleset",ATTR=attr)
! ionic species --> group
DO na=1,ntyp
CALL iotk_write_attr (attr,"name",TRIM(atm(na)),first=.true.)
CALL iotk_write_begin(iun, "group",ATTR=attr)
CALL iotk_write_attr (attr,"name","charge",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
write(iun,*) zv(na)
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "group")
ENDDO
! <attrib name="ionid"/>
CALL iotk_write_attr (attr,"name","ionid",first=.true.)
CALL iotk_write_attr (attr,"datatype","stringArray")
CALL iotk_write_begin(iun, "attrib",ATTR=attr)
write(iun,'(a)') (TRIM(atm(ityp(na))),na=1,nat)
CALL iotk_write_end(iun, "attrib")
! <attrib name="position"/>
CALL iotk_write_attr (attr,"name","position",first=.true.)
CALL iotk_write_attr (attr,"datatype","posArray")
CALL iotk_write_attr (attr,"condition","0")
CALL iotk_write_begin(iun, "attrib",ATTR=attr)
! write in cartesian coordinates in bohr
! problem with xyz ordering inrelation to real-space grid
DO na = 1, nat
!tau_r(1)=alat*(tau(1,na)*bg_qmc(1,1)+tau(2,na)*bg_qmc(1,2)+tau(3,na)*bg_qmc(1,3))
!tau_r(2)=alat*(tau(1,na)*bg_qmc(2,1)+tau(2,na)*bg_qmc(2,2)+tau(3,na)*bg_qmc(2,3))
!tau_r(3)=alat*(tau(1,na)*bg_qmc(3,1)+tau(2,na)*bg_qmc(3,2)+tau(3,na)*bg_qmc(3,3))
tau_r(1)=alat*tau(1,na)
tau_r(2)=alat*tau(2,na)
tau_r(3)=alat*tau(3,na)
WRITE(iun,100) (tau_r(j),j=1,3)
ENDDO
!write(iun,100) tau
CALL iotk_write_end(iun, "attrib")
CALL iotk_write_end(iun, "particleset")
! </particleset>
! <particleset name="e">
CALL iotk_write_attr (attr,"name","e",first=.true.)
CALL iotk_write_attr (attr,"random","yes")
CALL iotk_write_begin(iun, "particleset",ATTR=attr)
! <group name="u" size="" >
CALL iotk_write_attr (attr,"name","u",first=.true.)
CALL iotk_write_attr (attr,"size",nup)
CALL iotk_write_begin(iun, "group",ATTR=attr)
CALL iotk_write_attr (attr,"name","charge",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
write(iun,*) -1
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "group")
! <group name="d" size="" >
CALL iotk_write_attr (attr,"name","d",first=.true.)
CALL iotk_write_attr (attr,"size",ndown)
CALL iotk_write_begin(iun, "group",ATTR=attr)
CALL iotk_write_attr (attr,"name","charge",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
write(iun,*) -1
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "group")
CALL iotk_write_end(iun, "particleset")
CALL iotk_close_write(iun)
!! close the file
DO ik = 0, nk-1
! create a xml input file for each k-point
IF(nk .gt. 1) THEN
tmp = TRIM( tmp_dir ) // TRIM( prefix ) //TRIM(iotk_index(ik))// '.wfs.xml'
ELSE
tmp = TRIM( tmp_dir ) // TRIM( prefix )// '.wfs.xml'
ENDIF
CALL iotk_open_write(iun, FILE=TRIM(tmp), ROOT="qmcsystem", IERR=ierr )
! <wavefunction name="psi0">
CALL iotk_write_attr (attr,"name","psi0",first=.true.)
CALL iotk_write_attr (attr,"target","e")
CALL iotk_write_begin(iun, "wavefunction",ATTR=attr)
write(iun,'(a)') '<!-- Uncomment this out to use plane-wave basis functions'
CALL iotk_write_attr (attr,"type","PW",first=.true.)
CALL iotk_write_attr (attr,"href",TRIM(h5name))
CALL iotk_write_attr (attr,"version","1.10")
CALL iotk_write_begin(iun, "determinantset",ATTR=attr)
write(iun,'(a)') '--> '
CALL iotk_write_attr (attr,"type","bspline",first=.true.)
CALL iotk_write_attr (attr,"href",TRIM(h5name))
CALL iotk_write_attr (attr,"version","0.10")
CALL iotk_write_begin(iun, "determinantset",ATTR=attr)
CALL iotk_write_attr (attr,"ecut",ecutwfc/2,first=.true.)
! basisset to overwrite cutoff to a smaller value
CALL iotk_write_begin(iun, "basisset",ATTR=attr)
! add grid to use spline on FFT grid
CALL iotk_write_attr (attr,"dir","0",first=.true.)
CALL iotk_write_attr (attr,"npts",nr1s)
CALL iotk_write_attr (attr,"closed","no")
CALL iotk_write_empty(iun, "grid",ATTR=attr)
CALL iotk_write_attr (attr,"dir","1",first=.true.)
CALL iotk_write_attr (attr,"npts",nr2s)
CALL iotk_write_attr (attr,"closed","no")
CALL iotk_write_empty(iun, "grid",ATTR=attr)
CALL iotk_write_attr (attr,"dir","2",first=.true.)
CALL iotk_write_attr (attr,"npts",nr3s)
CALL iotk_write_attr (attr,"closed","no")
CALL iotk_write_empty(iun, "grid",ATTR=attr)
CALL iotk_write_end(iun, "basisset")
!CALL iotk_write_attr (attr,"href",TRIM(h5name),first=.true.)
!CALL iotk_write_empty(iun, "coefficients",ATTR=attr)
! write the index of the twist angle
CALL iotk_write_attr (attr,"name","twistIndex",first=.true.)
CALL iotk_write_begin(iun, "h5tag",ATTR=attr)
write(iun,*) ik
CALL iotk_write_end(iun, "h5tag")
CALL iotk_write_attr (attr,"name","twistAngle",first=.true.)
CALL iotk_write_begin(iun, "h5tag",ATTR=attr)
g_red(1)=at(1,1)*xk(1,ik+1)+at(2,1)*xk(2,ik+1)+at(3,1)*xk(3,ik+1)
g_red(2)=at(1,2)*xk(1,ik+1)+at(2,2)*xk(2,ik+1)+at(3,2)*xk(3,ik+1)
g_red(3)=at(1,3)*xk(1,ik+1)+at(2,3)*xk(2,ik+1)+at(3,3)*xk(3,ik+1)
!write(iun,100) xk(1,ik+1),xk(2,ik+1),xk(3,ik+1)
write(iun,100) g_red(1),g_red(2),g_red(3)
CALL iotk_write_end(iun, "h5tag")
!write(iun,'(a)') '<!-- Uncomment this out for bspline wavefunctions '
CALL iotk_write_attr (attr,"name","eigenstates",first=.true.)
CALL iotk_write_begin(iun, "h5tag",ATTR=attr)
write(iun,'(a)') TRIM(eigname)
CALL iotk_write_end(iun, "h5tag")
!write(iun,'(a)') '--> '
CALL iotk_write_begin(iun, "slaterdeterminant")
! build determinant for up electrons
CALL iotk_write_attr (attr,"id","updet",first=.true.)
CALL iotk_write_attr (attr,"size",nup)
CALL iotk_write_begin(iun, "determinant",ATTR=attr)
CALL iotk_write_attr (attr,"mode","ground",first=.true.)
CALL iotk_write_attr (attr,"spindataset",0)
CALL iotk_write_begin(iun, "occupation",ATTR=attr)
CALL iotk_write_end(iun, "occupation")
CALL iotk_write_end(iun, "determinant")
! build determinant for down electrons
CALL iotk_write_attr (attr,"id","downdet",first=.true.)
CALL iotk_write_attr (attr,"size",ndown)
IF( lsda ) CALL iotk_write_attr (attr,"ref","updet")
CALL iotk_write_begin(iun, "determinant",ATTR=attr)
CALL iotk_write_attr (attr,"mode","ground",first=.true.)
IF( lsda ) THEN
CALL iotk_write_attr (attr,"spindataset",1)
ELSE
CALL iotk_write_attr (attr,"spindataset",0)
ENDIF
CALL iotk_write_begin(iun, "occupation",ATTR=attr)
CALL iotk_write_end(iun, "occupation")
CALL iotk_write_end(iun, "determinant")
CALL iotk_write_end(iun, "slaterdeterminant")
CALL iotk_write_end(iun, "determinantset")
CALL iotk_write_end(iun, "wavefunction")
CALL iotk_close_write(iun)
ENDDO
tmp = TRIM( tmp_dir )//TRIM( prefix ) //'.pwscf.h5'
h5len = LEN_TRIM(tmp)
DO ig=1, ngtot
ig_c =igtog(ig)
g_cart(1,ig)=tpi/alat*g(1,ig_c)
g_cart(2,ig)=tpi/alat*g(2,ig_c)
g_cart(3,ig)=tpi/alat*g(3,ig_c)
g_qmc(1,ig)=at(1,1)*g(1,ig_c)+at(2,1)*g(2,ig_c)+at(3,1)*g(3,ig_c)
g_qmc(2,ig)=at(1,2)*g(1,ig_c)+at(2,2)*g(2,ig_c)+at(3,2)*g(3,ig_c)
g_qmc(3,ig)=at(1,3)*g(1,ig_c)+at(2,3)*g(2,ig_c)+at(3,3)*g(3,ig_c)
ENDDO
! generate hdf5 file containing all the parameters
! print out 2*occupied bands
CALL pwhdf_open_file(tmp,h5len)
!CALL pwhdf_write_basis(g,igtog,ngtot)
CALL pwhdf_write_basis(g_qmc,g_cart,ngtot)
CALL pwhdf_write_parameters(nelec_tot,nspin,nbnd,nk,ecutwfc/2,alat,at)
!
100 FORMAT (3(1x,f20.15))
ALLOCATE (eigpacked(ngtot))
! start a main section to save eigen values and vector
CALL pwhdf_open_eigg
DO ik = 1, nk
g_red(1)=at(1,1)*xk(1,ik)+at(2,1)*xk(2,ik)+at(3,1)*xk(3,ik)
g_red(2)=at(1,2)*xk(1,ik)+at(2,2)*xk(2,ik)+at(3,2)*xk(3,ik)
g_red(3)=at(1,3)*xk(1,ik)+at(2,3)*xk(2,ik)+at(3,3)*xk(3,ik)
CALL pwhdf_open_twist(ik,g_red(1),nbnd,nspin)
!CALL pwhdf_open_twist(ik,xk(1,ik),2*nbnd,nspin)
DO ispin = 1, nspin
ikk = ik + nk*(ispin-1)
IF( nks > 1 ) THEN
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio(evc,nwordwfc,iunwfc,ikk,-1)
ENDIF
DO ibnd = 1, nbnd
DO ig=1, ngtot
! now for all G vectors find the PW coefficient for this k-point
found = .FALSE.
DO ig7 = 1, npw
IF( igk(ig7) == igtog(ig) )THEN
eigpacked(ig)=evc(ig7,ibnd)
found = .TRUE.
GOTO 17
ENDIF
ENDDO
! if can't find the coefficient this is zero
17 IF( .NOT. found ) eigpacked(ig)=(0.d0,0.d0)
ENDDO
CALL pwhdf_write_band(ibnd,ispin,et(ibnd,ikk)/2,eigpacked,ngtot)
ENDDO
ENDDO
CALL pwhdf_close_twist
ENDDO
CALL pwhdf_close_eigg
!ALLOCATE (phase(nrxxs))
ALLOCATE(psireal(nrx1s*nrx2s*nrx3s))
ALLOCATE(psitr(nrx1s*nrx2s*nrx3s))
! open real-space wavefunction on FFT grid
CALL pwhdf_open_eigr(nr1s,nr2s,nr3s)
DO ik = 1, nk
!! evaluate the phase
!phase(:) = (0.d0,0.d0)
!if ( ig_(ik,ib)>0) phase( nls(ig_(ik,ib)) ) = (1.d0,0.d0)
!call cft3s (phase, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2)
g_red(1)=at(1,1)*xk(1,ik)+at(2,1)*xk(2,ik)+at(3,1)*xk(3,ik)
g_red(2)=at(1,2)*xk(1,ik)+at(2,2)*xk(2,ik)+at(3,2)*xk(3,ik)
g_red(3)=at(1,3)*xk(1,ik)+at(2,3)*xk(2,ik)+at(3,3)*xk(3,ik)
IF(g_red(1)*g_red(1)+g_red(2)*g_red(2)+g_red(3)*g_red(3)<1e-12) THEN
save_complex=0
ELSE
save_complex=1
ENDIF
! open a twsit
CALL pwhdf_open_twist(ik,g_red(1),nbnd,nspin)
DO ispin = 1, nspin
ikk = ik + nk*(ispin-1)
IF( nks > 1 ) THEN
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio(evc,nwordwfc,iunwfc,ikk,-1)
ENDIF
DO ibnd = 1, nbnd !!transform G to R
psic(:)=(0.d0,0.d0)
psic(nls(igk(1:npw)))=evc(1:npw,ibnd)
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
!! THIS IS ORIGINAL: changes are made to convert real and reorder for C arrays
!!CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psic,1)
IF(save_complex .eq. 1) THEN
!psic(:)=psic(:)/omega
!CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psic,save_complex)
ii=1
DO igx=1,nr1s
DO igy=0,nr2s-1
DO igz=0,nr3s-1
psitr(ii)=psic(igx+nr1s*(igy+igz*nr2s))/omega
ii=ii+1
ENDDO
ENDDO
ENDDO
CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psitr,save_complex)
ELSE
!DO ig=1,nr1s*nr2s*nr3s
! psireal(ig)=real(psic(ig))
!ENDDO
!psireal(1:nr1s*nr2s*nr3s)=real(psic(1:nr1s*nr2s*nr3s))/omega
ii=1
DO igx=1,nr1s
DO igy=0,nr2s-1
DO igz=0,nr3s-1
psireal(ii)=real(psic(igx+nr1s*(igy+igz*nr2s)))/omega
ii=ii+1
ENDDO
ENDDO
ENDDO
CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psireal,save_complex)
ENDIF
!! conversion and output complete for each band
ENDDO
ENDDO
CALL pwhdf_close_twist
ENDDO
CALL pwhdf_close_eigr
! write charge density
! ignore spin for the time being
!CALL pwhdf_write_rho(rho,rhog(1,1),ngm)
ii=1
DO igx=1,nr1s
DO igy=0,nr2s-1
DO igz=0,nr3s-1
psireal(ii)=rho(igx+nr1s*(igy+igz*nr2s),1)
ii=ii+1
ENDDO
ENDDO
ENDDO
CALL pwhdf_write_rho(psireal,rhog(1,1),ngm)
! close the file
CALL pwhdf_close_file
DEALLOCATE (igtog)
DEALLOCATE (index)
DEALLOCATE (becp)
DEALLOCATE (aux)
!DEALLOCATE (hpsi)
DEALLOCATE (eigpacked)
DEALLOCATE (g_qmc)
DEALLOCATE (g_cart)
DEALLOCATE (psireal)
DEALLOCATE (psitr)
!DEALLOCATE (phase)
END SUBROUTINE compute_qmcpack

614
src/QMCTools/pw2qmcpack.qe4.0.f90
View File

@ -1,614 +0,0 @@
!//////////////////////////////////////////////////////////////////////////////////////
!// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
!// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
!//
!// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
!//////////////////////////////////////////////////////////////////////////////////////
!
! Copyright (C) 2004 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
PROGRAM pw2qmcpack
!-----------------------------------------------------------------------
! This subroutine writes the file "prefix".pwscf.xml and "prefix".pwscf.h5
! containing the plane wave coefficients and other stuff needed by QMCPACK.
#include "f_defs.h"
USE io_files, ONLY : nd_nmbr, prefix, outdir, tmp_dir, trimcheck
USE io_global, ONLY : ionode, ionode_id
USE mp, ONLY : mp_bcast
!
IMPLICIT NONE
INTEGER :: ios
NAMELIST / inputpp / prefix, outdir
CALL start_postproc(nd_nmbr)
!
! set default values for variables in namelist
!
prefix = 'pwscf'
CALL get_env( 'ESPRESSO_TMPDIR', outdir )
IF ( TRIM( outdir ) == ' ' ) outdir = './'
ios = 0
IF ( ionode ) THEN
!
!READ (5, inputpp, err=200, iostat=ios)
READ (5, inputpp, iostat=ios)
tmp_dir = trimcheck (outdir)
!
END IF
CALL mp_bcast( ios, ionode_id )
IF ( ios/=0 ) CALL errore('pw2qmcpack', 'reading inputpp namelist', ABS(ios))
!
! ... Broadcast variables
!
CALL mp_bcast( prefix, ionode_id )
CALL mp_bcast(tmp_dir, ionode_id )
!
CALL read_file
CALL openfil_pp
!
CALL compute_qmcpack
!
CALL stop_pp
STOP
END PROGRAM pw2qmcpack
SUBROUTINE compute_qmcpack
USE kinds, ONLY: DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, zv, atm
USE cell_base, ONLY: omega, alat, tpiba2, at, bg
USE char, ONLY: title
USE constants, ONLY: tpi
USE ener, ONLY: ewld, ehart, etxc, vtxc, etot, etxcc
USE gvect, ONLY: ngm, gstart, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, g, gg, ecutwfc, gcutm, nl, igtongl
USE klist , ONLY: nks, nelec, nelup, neldw, xk
USE lsda_mod, ONLY: lsda, nspin
!!V3.0
!!USE scf, ONLY: rho, rho_core, rhog, rhog_core
!!USE vlocal, ONLY: vloc, vnew, strf
!!USE wvfct, ONLY: npw, npwx, nbnd, gamma_only, igk, g2kin, wg, et
USE scf, ONLY: rho, rho_core, rhog_core, vnew
USE vlocal, ONLY: vloc, strf
USE wvfct, ONLY: npw, npwx, nbnd, igk, g2kin, wg, et
USE control_flags, ONLY: gamma_only
USE uspp, ONLY: nkb, vkb, dvan
USE uspp_param, ONLY: nh
USE becmod, ONLY: becp, calbec
USE io_global, ONLY: stdout
USE io_files, ONLY: nd_nmbr, nwordwfc, iunwfc, iun => iunsat, tmp_dir, prefix
USE wavefunctions_module, ONLY : evc, psic
use gsmooth, only: nls, nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s
use iotk_module
use iotk_xtox_interf
USE mp_global, ONLY: inter_pool_comm, intra_pool_comm
USE mp, ONLY: mp_sum
IMPLICIT NONE
INTEGER :: ig, ibnd, ik, io, na, j, ispin, nbndup, nbnddown, &
nk, ngtot, ig7, ikk, nt, ijkb0, ikb, ih, jh, jkb, at_num, &
nelec_tot, nelec_up, nelec_down, ii, igx, igy, igz
INTEGER, ALLOCATABLE :: INDEX(:), igtog(:)
LOGICAL :: exst, found
REAL(DP) :: ek, eloc, enl, charge, etotefield
REAL(DP) :: bg_qmc(3,3), tau_r(3), g_red(3)
COMPLEX(DP), ALLOCATABLE :: phase(:),aux(:), hpsi(:,:), eigpacked(:)
COMPLEX(DP), ALLOCATABLE :: psitr(:)
REAL(DP), ALLOCATABLE :: g_cart(:,:),psireal(:)
INTEGER :: ios, ierr, h5len,ig_c,save_complex, nup,ndown
INTEGER, EXTERNAL :: atomic_number, is_complex
INTEGER, ALLOCATABLE:: g_qmc(:,:)
REAL (DP), EXTERNAL :: ewald
CHARACTER(256) :: tmp,h5name,eigname
CHARACTER(iotk_attlenx) :: attr
CALL init_us_1
CALL newd
!####io = 77
!####WRITE (6,'(/,5x,''Writing file pwfn.data for program CASINO'')')
!####CALL seqopn( 77, 'pwfn.data', 'formatted',exst)
!ALLOCATE (hpsi(npwx, nbnd))
ALLOCATE (aux(nrxx))
ALLOCATE (becp (nkb,nbnd))
! four times npwx should be enough
ALLOCATE (INDEX (4*npwx) )
ALLOCATE (igtog (4*npwx) )
!hpsi (:,:) = (0.d0, 0.d0)
INDEX(:) = 0
igtog(:) = 0
IF( lsda ) THEN
nbndup = nbnd
nbnddown = nbnd
nk = nks/2
! nspin = 2
ELSE
nbndup = nbnd
nbnddown = 0
nk = nks
! nspin = 1
ENDIF
! if(nks > 1) rewind(iunigk)
! do ik=1,nks
! if(nks > 1) read(iunigk) npw, igk
!
! if(nks > 1) rewind(iunigk)
ek = 0.d0
eloc= 0.d0
enl = 0.d0
!
DO ispin = 1, nspin
!
! calculate the local contribution to the total energy
!
! bring rho to G-space
!
!aux(:) = CMPLX ( rho(:,ispin), 0.d0)
aux(:) = CMPLX ( rho%of_r(:,ispin), 0.d0)
CALL cft3(aux,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
!
DO nt=1,ntyp
DO ig = 1, ngm
eloc = eloc + vloc(igtongl(ig),nt) * strf(ig,nt) &
* CONJG(aux(nl(ig)))
ENDDO
ENDDO
DO ik = 1, nk
ikk = ik + nk*(ispin-1)
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio (evc, nwordwfc, iunwfc, ikk, - 1)
CALL init_us_2 (npw, igk, xk (1, ikk), vkb)
!CALL ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
CALL calbec ( npw, vkb, evc, becp )
DO ig =1, npw
IF( igk(ig) > 4*npwx ) &
CALL errore ('pw2qmcpack','increase allocation of index', ig)
INDEX( igk(ig) ) = 1
ENDDO
ENDDO
ENDDO
!#ifdef __PARA
! CALL reduce(1,eloc)
! CALL reduce(1,ek)
! CALL poolreduce(1,ek)
! CALL poolreduce(1,enl)
!#endif
#ifdef __PARA
call mp_sum( eloc, intra_pool_comm )
call mp_sum( ek, intra_pool_comm )
call mp_sum( ek, inter_pool_comm )
call mp_sum( enl, inter_pool_comm )
!call mp_sum( demet, inter_pool_comm )
#endif
eloc = eloc * omega
ek = ek * tpiba2
ngtot = 0
DO ig = 1, 4*npwx
IF( INDEX(ig) == 1 ) THEN
ngtot = ngtot + 1
igtog(ngtot) = ig
ENDIF
ENDDO
ALLOCATE (g_qmc(3,ngtot))
ALLOCATE (g_cart(3,ngtot))
! get the number of electrons
nelec_tot= NINT(nelec)
nup=NINT(nelup)
ndown=NINT(neldw)
if(nup .eq. 0) then
ndown=nelec_tot/2
nup=nelec_tot-ndown
endif
h5name = TRIM( prefix ) // '.pwscf.h5'
eigname = "eigenstates_"//trim(iotk_itoa(nr1s))//'_'//trim(iotk_itoa(nr2s))//'_'//trim(iotk_itoa(nr3s))
bg_qmc(:,:)=bg(:,:)/alat
!! create a file for particle set
tmp = TRIM( tmp_dir ) // TRIM( prefix )// '.ptcl.xml'
CALL iotk_open_write(iun, FILE=TRIM(tmp), ROOT="qmcsystem", IERR=ierr )
CALL iotk_write_attr (attr,"name","global",first=.true.)
CALL iotk_write_begin(iun, "simulationcell",ATTR=attr)
CALL iotk_write_attr (attr,"name","lattice",first=.true.)
CALL iotk_write_attr (attr,"units","bohr")
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,100) alat*at(1,1), alat*at(2,1), alat*at(3,1)
WRITE(iun,100) alat*at(1,2), alat*at(2,2), alat*at(3,2)
WRITE(iun,100) alat*at(1,3), alat*at(2,3), alat*at(3,3)
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_attr (attr,"name","reciprocal",first=.true.)
CALL iotk_write_attr (attr,"units","2pi/bohr")
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,100) bg_qmc(1,1), bg_qmc(2,1), bg_qmc(3,1)
WRITE(iun,100) bg_qmc(1,2), bg_qmc(2,2), bg_qmc(3,2)
WRITE(iun,100) bg_qmc(1,3), bg_qmc(2,3), bg_qmc(3,3)
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_attr (attr,"name","bconds",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,'(a)') 'p p p'
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_attr (attr,"name","LR_dim_cutoff",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
WRITE(iun,'(a)') '10'
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "simulationcell")
! <particleset name="ions">
CALL iotk_write_attr (attr,"name","ion0",first=.true.)
CALL iotk_write_attr (attr,"size",nat)
CALL iotk_write_begin(iun, "particleset",ATTR=attr)
! ionic species --> group
DO na=1,ntyp
CALL iotk_write_attr (attr,"name",TRIM(atm(na)),first=.true.)
CALL iotk_write_begin(iun, "group",ATTR=attr)
CALL iotk_write_attr (attr,"name","charge",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
write(iun,*) zv(na)
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "group")
ENDDO
! <attrib name="ionid"/>
CALL iotk_write_attr (attr,"name","ionid",first=.true.)
CALL iotk_write_attr (attr,"datatype","stringArray")
CALL iotk_write_begin(iun, "attrib",ATTR=attr)
write(iun,'(a)') (TRIM(atm(ityp(na))),na=1,nat)
CALL iotk_write_end(iun, "attrib")
! <attrib name="position"/>
CALL iotk_write_attr (attr,"name","position",first=.true.)
CALL iotk_write_attr (attr,"datatype","posArray")
CALL iotk_write_attr (attr,"condition","0")
CALL iotk_write_begin(iun, "attrib",ATTR=attr)
! write in cartesian coordinates in bohr
! problem with xyz ordering inrelation to real-space grid
DO na = 1, nat
!tau_r(1)=alat*(tau(1,na)*bg_qmc(1,1)+tau(2,na)*bg_qmc(1,2)+tau(3,na)*bg_qmc(1,3))
!tau_r(2)=alat*(tau(1,na)*bg_qmc(2,1)+tau(2,na)*bg_qmc(2,2)+tau(3,na)*bg_qmc(2,3))
!tau_r(3)=alat*(tau(1,na)*bg_qmc(3,1)+tau(2,na)*bg_qmc(3,2)+tau(3,na)*bg_qmc(3,3))
tau_r(1)=alat*tau(1,na)
tau_r(2)=alat*tau(2,na)
tau_r(3)=alat*tau(3,na)
WRITE(iun,100) (tau_r(j),j=1,3)
ENDDO
!write(iun,100) tau
CALL iotk_write_end(iun, "attrib")
CALL iotk_write_end(iun, "particleset")
! </particleset>
! <particleset name="e">
CALL iotk_write_attr (attr,"name","e",first=.true.)
CALL iotk_write_attr (attr,"random","yes")
CALL iotk_write_begin(iun, "particleset",ATTR=attr)
! <group name="u" size="" >
CALL iotk_write_attr (attr,"name","u",first=.true.)
CALL iotk_write_attr (attr,"size",nup)
CALL iotk_write_begin(iun, "group",ATTR=attr)
CALL iotk_write_attr (attr,"name","charge",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
write(iun,*) -1
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "group")
! <group name="d" size="" >
CALL iotk_write_attr (attr,"name","d",first=.true.)
CALL iotk_write_attr (attr,"size",ndown)
CALL iotk_write_begin(iun, "group",ATTR=attr)
CALL iotk_write_attr (attr,"name","charge",first=.true.)
CALL iotk_write_begin(iun, "parameter",ATTR=attr)
write(iun,*) -1
CALL iotk_write_end(iun, "parameter")
CALL iotk_write_end(iun, "group")
CALL iotk_write_end(iun, "particleset")
CALL iotk_close_write(iun)
!! close the file
DO ik = 0, nk-1
! create a xml input file for each k-point
IF(nk .gt. 1) THEN
tmp = TRIM( tmp_dir ) // TRIM( prefix ) //TRIM(iotk_index(ik))// '.wfs.xml'
ELSE
tmp = TRIM( tmp_dir ) // TRIM( prefix )// '.wfs.xml'
ENDIF
CALL iotk_open_write(iun, FILE=TRIM(tmp), ROOT="qmcsystem", IERR=ierr )
! <wavefunction name="psi0">
CALL iotk_write_attr (attr,"name","psi0",first=.true.)
CALL iotk_write_attr (attr,"target","e")
CALL iotk_write_begin(iun, "wavefunction",ATTR=attr)
write(iun,'(a)') '<!-- Uncomment this out to use plane-wave basis functions'
CALL iotk_write_attr (attr,"type","PW",first=.true.)
CALL iotk_write_attr (attr,"href",TRIM(h5name))
CALL iotk_write_attr (attr,"version","1.10")
CALL iotk_write_begin(iun, "determinantset",ATTR=attr)
write(iun,'(a)') '--> '
CALL iotk_write_attr (attr,"type","bspline",first=.true.)
CALL iotk_write_attr (attr,"href",TRIM(h5name))
CALL iotk_write_attr (attr,"version","0.10")
CALL iotk_write_begin(iun, "determinantset",ATTR=attr)
CALL iotk_write_attr (attr,"ecut",ecutwfc/2,first=.true.)
! basisset to overwrite cutoff to a smaller value
CALL iotk_write_begin(iun, "basisset",ATTR=attr)
! add grid to use spline on FFT grid
CALL iotk_write_attr (attr,"dir","0",first=.true.)
CALL iotk_write_attr (attr,"npts",nr1s)
CALL iotk_write_attr (attr,"closed","no")
CALL iotk_write_empty(iun, "grid",ATTR=attr)
CALL iotk_write_attr (attr,"dir","1",first=.true.)
CALL iotk_write_attr (attr,"npts",nr2s)
CALL iotk_write_attr (attr,"closed","no")
CALL iotk_write_empty(iun, "grid",ATTR=attr)
CALL iotk_write_attr (attr,"dir","2",first=.true.)
CALL iotk_write_attr (attr,"npts",nr3s)
CALL iotk_write_attr (attr,"closed","no")
CALL iotk_write_empty(iun, "grid",ATTR=attr)
CALL iotk_write_end(iun, "basisset")
!CALL iotk_write_attr (attr,"href",TRIM(h5name),first=.true.)
!CALL iotk_write_empty(iun, "coefficients",ATTR=attr)
! write the index of the twist angle
CALL iotk_write_attr (attr,"name","twistIndex",first=.true.)
CALL iotk_write_begin(iun, "h5tag",ATTR=attr)
write(iun,*) ik
CALL iotk_write_end(iun, "h5tag")
CALL iotk_write_attr (attr,"name","twistAngle",first=.true.)
CALL iotk_write_begin(iun, "h5tag",ATTR=attr)
g_red(1)=at(1,1)*xk(1,ik+1)+at(2,1)*xk(2,ik+1)+at(3,1)*xk(3,ik+1)
g_red(2)=at(1,2)*xk(1,ik+1)+at(2,2)*xk(2,ik+1)+at(3,2)*xk(3,ik+1)
g_red(3)=at(1,3)*xk(1,ik+1)+at(2,3)*xk(2,ik+1)+at(3,3)*xk(3,ik+1)
!write(iun,100) xk(1,ik+1),xk(2,ik+1),xk(3,ik+1)
write(iun,100) g_red(1),g_red(2),g_red(3)
CALL iotk_write_end(iun, "h5tag")
!write(iun,'(a)') '<!-- Uncomment this out for bspline wavefunctions '
CALL iotk_write_attr (attr,"name","eigenstates",first=.true.)
CALL iotk_write_begin(iun, "h5tag",ATTR=attr)
write(iun,'(a)') TRIM(eigname)
CALL iotk_write_end(iun, "h5tag")
!write(iun,'(a)') '--> '
CALL iotk_write_begin(iun, "slaterdeterminant")
! build determinant for up electrons
CALL iotk_write_attr (attr,"id","updet",first=.true.)
CALL iotk_write_attr (attr,"size",nup)
CALL iotk_write_begin(iun, "determinant",ATTR=attr)
CALL iotk_write_attr (attr,"mode","ground",first=.true.)
CALL iotk_write_attr (attr,"spindataset",0)
CALL iotk_write_begin(iun, "occupation",ATTR=attr)
CALL iotk_write_end(iun, "occupation")
CALL iotk_write_end(iun, "determinant")
! build determinant for down electrons
CALL iotk_write_attr (attr,"id","downdet",first=.true.)
CALL iotk_write_attr (attr,"size",ndown)
IF( lsda ) CALL iotk_write_attr (attr,"ref","updet")
CALL iotk_write_begin(iun, "determinant",ATTR=attr)
CALL iotk_write_attr (attr,"mode","ground",first=.true.)
IF( lsda ) THEN
CALL iotk_write_attr (attr,"spindataset",1)
ELSE
CALL iotk_write_attr (attr,"spindataset",0)
ENDIF
CALL iotk_write_begin(iun, "occupation",ATTR=attr)
CALL iotk_write_end(iun, "occupation")
CALL iotk_write_end(iun, "determinant")
CALL iotk_write_end(iun, "slaterdeterminant")
CALL iotk_write_end(iun, "determinantset")
CALL iotk_write_end(iun, "wavefunction")
CALL iotk_close_write(iun)
ENDDO
tmp = TRIM( tmp_dir )//TRIM( prefix ) //'.pwscf.h5'
h5len = LEN_TRIM(tmp)
DO ig=1, ngtot
ig_c =igtog(ig)
g_cart(1,ig)=tpi/alat*g(1,ig_c)
g_cart(2,ig)=tpi/alat*g(2,ig_c)
g_cart(3,ig)=tpi/alat*g(3,ig_c)
g_qmc(1,ig)=at(1,1)*g(1,ig_c)+at(2,1)*g(2,ig_c)+at(3,1)*g(3,ig_c)
g_qmc(2,ig)=at(1,2)*g(1,ig_c)+at(2,2)*g(2,ig_c)+at(3,2)*g(3,ig_c)
g_qmc(3,ig)=at(1,3)*g(1,ig_c)+at(2,3)*g(2,ig_c)+at(3,3)*g(3,ig_c)
ENDDO
! generate hdf5 file containing all the parameters
! print out 2*occupied bands
CALL pwhdf_open_file(tmp,h5len)
!CALL pwhdf_write_basis(g,igtog,ngtot)
CALL pwhdf_write_basis(g_qmc,g_cart,ngtot)
CALL pwhdf_write_parameters(nelec_tot,nspin,nbnd,nk,ecutwfc/2,alat,at)
!
100 FORMAT (3(1x,f20.15))
ALLOCATE (eigpacked(ngtot))
! start a main section to save eigen values and vector
CALL pwhdf_open_eigg
DO ik = 1, nk
g_red(1)=at(1,1)*xk(1,ik)+at(2,1)*xk(2,ik)+at(3,1)*xk(3,ik)
g_red(2)=at(1,2)*xk(1,ik)+at(2,2)*xk(2,ik)+at(3,2)*xk(3,ik)
g_red(3)=at(1,3)*xk(1,ik)+at(2,3)*xk(2,ik)+at(3,3)*xk(3,ik)
CALL pwhdf_open_twist(ik,g_red(1),nbnd,nspin)
!CALL pwhdf_open_twist(ik,xk(1,ik),2*nbnd,nspin)
DO ispin = 1, nspin
ikk = ik + nk*(ispin-1)
IF( nks > 1 ) THEN
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio(evc,nwordwfc,iunwfc,ikk,-1)
ENDIF
DO ibnd = 1, nbnd
DO ig=1, ngtot
! now for all G vectors find the PW coefficient for this k-point
found = .FALSE.
DO ig7 = 1, npw
IF( igk(ig7) == igtog(ig) )THEN
eigpacked(ig)=evc(ig7,ibnd)
found = .TRUE.
GOTO 17
ENDIF
ENDDO
! if can't find the coefficient this is zero
17 IF( .NOT. found ) eigpacked(ig)=(0.d0,0.d0)
ENDDO
CALL pwhdf_write_band(ibnd,ispin,et(ibnd,ikk)/2,eigpacked,ngtot)
ENDDO
ENDDO
CALL pwhdf_close_twist
ENDDO
CALL pwhdf_close_eigg
!ALLOCATE (phase(nrxxs))
ALLOCATE(psireal(nrx1s*nrx2s*nrx3s))
ALLOCATE(psitr(nrx1s*nrx2s*nrx3s))
! open real-space wavefunction on FFT grid
CALL pwhdf_open_eigr(nr1s,nr2s,nr3s)
DO ik = 1, nk
!! evaluate the phase
!phase(:) = (0.d0,0.d0)
!if ( ig_(ik,ib)>0) phase( nls(ig_(ik,ib)) ) = (1.d0,0.d0)
!call cft3s (phase, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, +2)
g_red(1)=at(1,1)*xk(1,ik)+at(2,1)*xk(2,ik)+at(3,1)*xk(3,ik)
g_red(2)=at(1,2)*xk(1,ik)+at(2,2)*xk(2,ik)+at(3,2)*xk(3,ik)
g_red(3)=at(1,3)*xk(1,ik)+at(2,3)*xk(2,ik)+at(3,3)*xk(3,ik)
IF(g_red(1)*g_red(1)+g_red(2)*g_red(2)+g_red(3)*g_red(3)<1e-12) THEN
save_complex=0
ELSE
save_complex=1
ENDIF
! open a twsit
CALL pwhdf_open_twist(ik,g_red(1),nbnd,nspin)
DO ispin = 1, nspin
ikk = ik + nk*(ispin-1)
IF( nks > 1 ) THEN
CALL gk_sort (xk (1, ikk), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
CALL davcio(evc,nwordwfc,iunwfc,ikk,-1)
ENDIF
DO ibnd = 1, nbnd !!transform G to R
psic(:)=(0.d0,0.d0)
psic(nls(igk(1:npw)))=evc(1:npw,ibnd)
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
!! THIS IS ORIGINAL: changes are made to convert real and reorder for C arrays
!!CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psic,1)
IF(save_complex .eq. 1) THEN
!psic(:)=psic(:)/omega
!CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psic,save_complex)
ii=1
DO igx=1,nr1s
DO igy=0,nr2s-1
DO igz=0,nr3s-1
psitr(ii)=psic(igx+nr1s*(igy+igz*nr2s))/omega
ii=ii+1
ENDDO
ENDDO
ENDDO
CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psitr,save_complex)
ELSE
!DO ig=1,nr1s*nr2s*nr3s
! psireal(ig)=real(psic(ig))
!ENDDO
!psireal(1:nr1s*nr2s*nr3s)=real(psic(1:nr1s*nr2s*nr3s))/omega
ii=1
DO igx=1,nr1s
DO igy=0,nr2s-1
DO igz=0,nr3s-1
psireal(ii)=real(psic(igx+nr1s*(igy+igz*nr2s)))/omega
ii=ii+1
ENDDO
ENDDO
ENDDO
CALL pwhdf_write_wfr(ibnd,ispin,et(ibnd,ikk)/2,psireal,save_complex)
ENDIF
!! conversion and output complete for each band
ENDDO
ENDDO
CALL pwhdf_close_twist
ENDDO
CALL pwhdf_close_eigr
! write charge density
! ignore spin for the time being
!CALL pwhdf_write_rho(rho,rhog(1,1),ngm)
ii=1
DO igx=1,nr1s
DO igy=0,nr2s-1
DO igz=0,nr3s-1
!psireal(ii)=rho(igx+nr1s*(igy+igz*nr2s),1)
psireal(ii)=rho%of_r(igx+nr1s*(igy+igz*nr2s),1)
ii=ii+1
ENDDO
ENDDO
ENDDO
!CALL pwhdf_write_rho(psireal,rhog(1,1),ngm)
CALL pwhdf_write_rho(psireal,rho%of_g(1,1),ngm)
! close the file
CALL pwhdf_close_file
DEALLOCATE (igtog)
DEALLOCATE (index)
DEALLOCATE (becp)
DEALLOCATE (aux)
!DEALLOCATE (hpsi)
DEALLOCATE (eigpacked)
DEALLOCATE (g_qmc)
DEALLOCATE (g_cart)
DEALLOCATE (psireal)
DEALLOCATE (psitr)
!DEALLOCATE (phase)
END SUBROUTINE compute_qmcpack

394
src/QMCTools/pw_hdf5.c
View File

@ -1,394 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
/*
* Copyright (C) 2004 PWSCF group
* Copyright (C) 2007 QMCPACK developers
*
* @author Jeongnim Kim http://www.mcc.uiuc.edu/qmcpack/
* @brief Implements generic hdf5 interfaces for plane wave codes and qmcpack
*
* - pwhdf_open_file: open hdf5 file
* - pwhdf_close_file : close hdf5 file
* - pwhdf_open_eigg : open eigenstates
* - pwhdf_close_eigg : close eigenstates
* - pwhdf_open_eigr : open eigenstates_nx_ny_nz
* - pwhdf_close_eigr : close eigenstates_nx_ny_nz
*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <sys/types.h>
#include "c_defs.h"
#include "hdf5.h"
/* file handler */
static hid_t h_file=-1;
/* main handler */
static hid_t h_main=-1;
/* twist angle handler */
static hid_t h_twist=-1;
/* number of fft grid */
static int h_ngrid[3];
/* number of real-space grids */
static int h_ngridtot=0;
/* check for gamma */
static int is_gamma=0;
void F77_FUNC_(pwhdf_open_file,PWHDF_OPEN_FILE)(const char* fname, const int* length)
{
char * hfname = ( char * ) malloc( (*length) + 1 ) ;
memcpy( hfname , fname , *length ) ;
hfname[*length] = '\0' ;
if(h_file>=0) H5Fclose(h_file);
h_file = H5Fcreate(hfname,H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
/* impelements version 1.00 hdf5 format */
int version[]={1,10};
hsize_t dim=2;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_file, "version", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,version);
H5Sclose(dataspace);
H5Dclose(dataset);
free(hfname);
}
void F77_FUNC_(pwhdf_close_file,PWHDF_CLOSE_FILE)()
{
if(h_file>=0) H5Fclose(h_file);
h_file=-1;
}
/** write basisset: number of plane waves, plane wave coefficients
*/
void F77_FUNC_(pwhdf_write_basis,PWHDF_WRITE_BASIS)(const int* ig,
const double* gcart, const int* ngtot)
{
int ng=*ngtot;
hid_t h_basis = H5Gcreate(h_file,"basis",0);
hsize_t dim=1;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_basis, "num_planewaves", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,ngtot);
H5Sclose(dataspace);
H5Dclose(dataset);
hsize_t dims[2];
dims[0] = ng;
dims[1] = 3;
dataspace = H5Screate_simple(2, dims, NULL);
dataset = H5Dcreate(h_basis, "multipliers", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,ig);
H5Dclose(dataset);
H5Sclose(dataspace);
dataspace = H5Screate_simple(2, dims, NULL);
dataset = H5Dcreate(h_basis, "planewaves", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,gcart);
H5Sclose(dataspace);
H5Dclose(dataset);
H5Gclose(h_basis);
}
/** write basisset: number of plane waves, plane wave coefficients
void F77_FUNC_(pwhdf_write_basis,PWHDF_WRITE_BASIS)(const double* g, const int* igtog, const int* ngtot)
{
int ng=*ngtot;
int *ig=(int*)malloc(3*ng*sizeof(int));
for(int i=0,i3=0; i<ng; i++)
{
int cur=3*(igtog[i]-1);
ig[i3++]=(int)g[cur++];
ig[i3++]=(int)g[cur++];
ig[i3++]=(int)g[cur++];
}
hid_t h_basis = H5Gcreate(h_file,"basis",0);
hsize_t dim=1;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_basis, "num_planewaves", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,ngtot);
H5Sclose(dataspace);
H5Dclose(dataset);
hsize_t dims[2];
dims[0] = ng;
dims[1] = 3;
dataspace = H5Screate_simple(2, dims, NULL);
dataset = H5Dcreate(h_basis, "planewaves", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,ig);
H5Sclose(dataspace);
H5Dclose(dataset);
H5Gclose(h_basis);
free(ig);
}
*/
void F77_FUNC_(pwhdf_write_parameters,PWHDF_WRITE_PARAMETERS)(
const int* nelec, const int* nspin, const int* nband, const int* nk,
const double* ecut, const double* alat, const double* at)
{
hid_t h_param = H5Gcreate(h_file,"parameters",0);
hsize_t dim=1;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_param, "num_spins", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,nspin);
H5Sclose(dataspace);
H5Dclose(dataset);
dataspace= H5Screate_simple(1, &dim, NULL);
dataset= H5Dcreate(h_param, "num_electrons", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,nelec);
H5Sclose(dataspace);
H5Dclose(dataset);
dataspace= H5Screate_simple(1, &dim, NULL);
dataset= H5Dcreate(h_param, "num_bands", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,nband);
H5Sclose(dataspace);
H5Dclose(dataset);
dataspace= H5Screate_simple(1, &dim, NULL);
dataset= H5Dcreate(h_param, "num_twists", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,nk);
H5Sclose(dataspace);
H5Dclose(dataset);
int iscomplex=1;
dataspace= H5Screate_simple(1, &dim, NULL);
dataset= H5Dcreate(h_param, "complex_coefficients", H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,&iscomplex);
H5Sclose(dataspace);
H5Dclose(dataset);
dataspace= H5Screate_simple(1, &dim, NULL);
dataset= H5Dcreate(h_param, "maximum_ecut", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,ecut);
H5Sclose(dataspace);
H5Dclose(dataset);
double lattice[9];
for(int i=0; i<9; i++) lattice[i]=(*alat)*at[i];
hsize_t dims[]={3,3};
dataspace = H5Screate_simple(2, dims, NULL);
dataset = H5Dcreate(h_param, "lattice", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,lattice);
H5Sclose(dataspace);
H5Dclose(dataset);
H5Gclose(h_param);
}
/* open mainbody:eigenstates */
void F77_FUNC_(pwhdf_open_eigg,PWHDF_OPEN_EIGG)()
{
if(h_main>=0) H5Gclose(h_main);
h_main = H5Gcreate(h_file,"eigenstates",0);
}
/* close eigenstates */
void F77_FUNC_(pwhdf_close_eigg,PWHDF_CLOSE_EIGG)()
{
if(h_main>=0) H5Gclose(h_main);
h_main=-1;
}
/* open twist# */
void F77_FUNC_(pwhdf_open_twist,PWHDF_OPEN_TWIST)(const int* ik, const double *xk,
const int* nband, const int* nspin)
{
char twistname[16];
sprintf(twistname,"twist%i",(*ik)-1);
if(h_twist>=0) H5Gclose(h_twist);
h_twist = H5Gcreate(h_main,twistname,0);
/* write twist_angle */
hsize_t dim=3;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_twist, "twist_angle", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,xk);
H5Sclose(dataspace);
H5Dclose(dataset);
if((xk[0]*xk[0]+xk[1]*xk[1]+xk[2]*xk[2]) < 1e-12)
is_gamma=1;
else
is_gamma=0;
/* create band#/spin# groups so that H5Gopen can be used */
for(int ib=0; ib<*nband; ib++)
{
char bandname[16];
sprintf(bandname,"band%i",ib);
hid_t h_band = H5Gcreate(h_twist,bandname,0);
for(int ispin=0; ispin<*nspin; ispin++)
{
char spinname[16];
sprintf(spinname,"spin%i",ispin);
hid_t h_spin = H5Gcreate(h_band,spinname,0);
H5Gclose(h_spin);
}
H5Gclose(h_band);
}
}
/* close twist# */
void F77_FUNC_(pwhdf_close_twist,PWHDF_CLOSE_TWIST)()
{
if(h_twist>=0) H5Gclose(h_twist);
h_twist=-1;
}
/* write eigen value and eigen vector for (ibnd, ispin) */
void F77_FUNC_(pwhdf_write_band,PWHDF_WRITE_BAND)(const int* ibnd,
const int* ispin, const double* e,
const double* eigv, const int* ngtot)
{
char spinname[16];
sprintf(spinname,"band%i/spin%i",(*ibnd)-1,(*ispin)-1);
hid_t h_spin = H5Gopen(h_twist,spinname);
hsize_t dim=1;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_spin, "eigenvalue", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,e);
H5Sclose(dataspace);
H5Dclose(dataset);
hsize_t dims[2];
dims[0] = *ngtot;
dims[1] = 2;
dataspace = H5Screate_simple(2, dims, NULL);
dataset = H5Dcreate(h_spin, "eigenvector", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,eigv);
H5Sclose(dataspace);
H5Dclose(dataset);
H5Gclose(h_spin);
}
/* open mainbody:eigenstates */
void F77_FUNC_(pwhdf_open_eigr,PWHDF_OPEN_EIGR)(const int* nr1, const int* nr2, const int*nr3)
{
/* swap the order : CHECK THIS! */
h_ngrid[0]=*nr1;
h_ngrid[1]=*nr2;
h_ngrid[2]=*nr3;
h_ngridtot=(*nr1)*(*nr2)*(*nr3);
char wfrname[16];
sprintf(wfrname,"eigenstates_%i_%i_%i",h_ngrid[0],h_ngrid[1],h_ngrid[2]);
if(h_main>=0) H5Gclose(h_main);
h_main = H5Gcreate(h_file,wfrname,0);
}
/* close twist# */
void F77_FUNC_(pwhdf_close_eigr,PWHDF_CLOSE_EIGR)()
{
if(h_main>=0) H5Gclose(h_main);
h_main=-1;
}
/* write eigen value and eigen vector for (ibnd, ispin) */
void F77_FUNC_(pwhdf_write_wfr,PWHDF_WRITE_WFR)(const int* ibnd,
const int* ispin, const double* e,
const double* eigr, const int* use_complex)
{
char spinname[16];
sprintf(spinname,"band%i/spin%i",(*ibnd)-1,(*ispin)-1);
/*sprintf(spinname,"band_%i",(*ibnd)-1);*/
hid_t h_spin = H5Gopen(h_twist,spinname);
/* write eigenvalue */
hsize_t dim=1;
hid_t dataspace= H5Screate_simple(1, &dim, NULL);
hid_t dataset= H5Dcreate(h_spin, "eigenvalue", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,e);
H5Sclose(dataspace);
H5Dclose(dataset);
hsize_t dims_out=(*use_complex)?4:3;
hsize_t dims[4];
dims[0] = h_ngrid[0];
dims[1] = h_ngrid[1];
dims[2] = h_ngrid[2];
dims[3] = 2;
dataspace = H5Screate_simple(dims_out, dims, NULL);
dataset = H5Dcreate(h_spin, "eigenvector", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,eigr);
H5Sclose(dataspace);
H5Dclose(dataset);
H5Gclose(h_spin);
}
/* write eigen value and eigen vector for (ibnd, ispin) */
void F77_FUNC_(pwhdf_write_rho,PWHDF_WRITE_RHO)(const double* rho, const double* rhog, int ngm)
{
/* write eigenvector */
hsize_t dims[3];
dims[0] = h_ngrid[0];
dims[1] = h_ngrid[1];
dims[2] = h_ngrid[2];
hid_t dataspace = H5Screate_simple(3, dims, NULL);
hid_t dataset = H5Dcreate(h_file, "chargedensity_r", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
hid_t ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,rho);
H5Sclose(dataspace);
H5Dclose(dataset);
/*
hsize_t gdims[2];
gdims[0]=ngm;
gdims[1]=2;
dataspace = H5Screate_simple(2, gdims, NULL);
dataset = H5Dcreate(h_file, "chargedensity_g", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT);
ret = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,rhog);
H5Sclose(dataspace);
H5Dclose(dataset);
*/
/* testing with paraview/vtk
if(is_gamma)
{
char vtkname[32];
sprintf(vtkname,"band%i.vtk",(*ibnd)-1);
FILE *vtk=fopen(vtkname,"w");
fprintf(vtk,"# vtk DataFile Version 3.0\n");
fprintf(vtk,"vtk output\n");
fprintf(vtk,"ASCII\n");
fprintf(vtk,"DATASET STRUCTURED_POINTS\n");
fprintf(vtk,"DIMENSIONS %i %i %i\n",h_ngrid[0],h_ngrid[1],h_ngrid[2]);
fprintf(vtk,"ORIGIN 0 0 0\n");
fprintf(vtk,"SPACING 1 1 1\n");
fprintf(vtk,"\nPOINT_DATA %i\n",h_ngridtot);
fprintf(vtk,"SCALARS scalars float\n");
fprintf(vtk,"LOOKUP_TABLE default\n");
for(int i=0,i2=0; i<h_ngridtot;i+=10)
{
for(int j=0; j<10; j++,i2+=2) fprintf(vtk,"%12.6e ",eigr[i2]*eigr[i2]);
fprintf(vtk,"\n");
}
fprintf(vtk,"\n");
fclose(vtk);
}
*/
}

92
src/QMCTools/spline.cpp
View File

@ -1,92 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include "OhmmsPETE/OhmmsVector.h"
#include "Numerics/OneDimGridBase.h"
#include "Numerics/OneDimCubicSpline.h"
#define FUNCTION(r) cos(r)
#define DFUNCTION(r) sin(r)
//#define USE_PBC
int main(int argc, char** argv)
{
double ri = 0.0;
double rf = 1.0;
int npts = 101;
const double nk=2;
const double twopi = 8.0*atan(1.0)*nk;
LinearGrid<double> agrid;
agrid.set(ri,rf,npts);
#if defined(USE_PBC)
typedef OneDimCubicSplinePBC<double> OrbitalType;
#else
typedef OneDimCubicSplineFirst<double> OrbitalType;
#endif
OrbitalType aorb(&agrid);
aorb.resize(agrid.size());
for(int i=0; i<agrid.size(); i++)
{
aorb(i)=FUNCTION(twopi*agrid(i));
}
aorb.spline(0,twopi*DFUNCTION(twopi*agrid.rmin()),agrid.size()-1, twopi*DFUNCTION(twopi*agrid.rmax()));
std::string fname("testpbc.dat");
std::ofstream dfile(fname.c_str());
dfile.setf(std::ios::scientific, std::ios::floatfield);
dfile.setf(std::ios::left,std::ios::adjustfield);
dfile.precision(15);
const double c1=1.0/twopi;
const double c2=c1*c1;
double du,d2u,_r,_rinv,y;
#if defined(USE_PBC)
for(int ig=agrid.size()/2; ig<agrid.size() ; ig++)
{
_r = agrid(ig)+0.5*agrid.dr(ig)-agrid.rmax();
_rinv = 1.0/_r;
//aorb.setgrid(_r);
y=aorb.evaluate(_r,_rinv,du,d2u);
dfile << std::setw(30) << _r << std::setw(30) << FUNCTION(twopi*_r) << std::setw(30) << y << std::setw(30) << du*c1 << std::setw(3) << d2u*c2 << std::endl;
}
#endif
for(int ig=0; ig<agrid.size()-1; ig++)
{
_r = agrid(ig)+0.5*agrid.dr(ig);
_rinv = 1.0/_r;
//aorb.setgrid(_r);
y=aorb.evaluate(_r,_rinv,du,d2u);
dfile << std::setw(30) << _r << std::setw(30) << FUNCTION(twopi*_r) << std::setw(30) << y
<< std::setw(30) << du*c1 << std::setw(3) << d2u*c2 << std::endl;
}
#if defined(USE_PBC)
for(int ig=0; ig<agrid.size()/2; ig++)
{
_r = agrid(ig)+0.5*agrid.dr(ig)+agrid.rmax();
_rinv = 1.0/_r;
//aorb.setgrid(_r);
y=aorb.evaluate(_r,_rinv,du,d2u);
dfile << std::setw(30) << _r << std::setw(30) << FUNCTION(twopi*_r) << std::setw(30) << y << std::setw(30) << du*c1 << std::setw(3) << d2u*c2 << std::endl;
}
#endif
dfile << std::endl;
return 0;
}

237
src/QMCTools/spline3.cpp
View File

@ -1,237 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <sstream>
#include "OhmmsData/OhmmsElementBase.h"
#include "Numerics/HDFSTLAttrib.h"
#include "Numerics/TriCubicSplineT.h"
#include "Numerics/HDFTriCubicSpline.h"
#include "Utilities/Clock.h"
struct TestFunc
{
double k0,k1,k2;
double d2factor;
TestFunc(int nk0=1, int nk1=1, int nk2=1)
{
const double twopi = 8.0*atan(1.0);
k0=twopi*static_cast<double>(nk0);
k1=twopi*static_cast<double>(nk1);
k2=twopi*static_cast<double>(nk2);
d2factor = -(k0*k0+k1*k1+k2*k2);
}
//inline double operator()(const TinyVector<double,3>& pos) {
// return sin(k0*pos[0])*sin(k1*pos[1])*sin(k2*pos[2]);
//}
//inline double operator()(double x, double y, double z) {
// return sin(k0*x)*sin(k1*y)*sin(k2*z);
//}
//
inline double f(const TinyVector<double,3>& pos)
{
return sin(k0*pos[0])*sin(k1*pos[1])*sin(k2*pos[2]);
}
inline double f(double x, double y, double z)
{
return sin(k0*x)*sin(k1*y)*sin(k2*z);
}
inline double d2f(const TinyVector<double,3>& pos)
{
return d2factor*f(pos);
}
inline double d2f(double x, double y, double z)
{
return d2factor*f(x,y,z);
}
};
struct ComboFunc
{
std::vector<double> C;
std::vector<TestFunc*> F;
ComboFunc() {}
void push_back(double c, TestFunc* fn)
{
C.push_back(c);
F.push_back(fn);
}
inline double f(double x, double y, double z)
{
double res=0;
for(int i=0; i<C.size(); i++)
res += C[i]*F[i]->f(x,y,z);
return res;
}
inline double d2f(double x, double y, double z)
{
double res=0;
for(int i=0; i<C.size(); i++)
res += C[i]*F[i]->d2f(x,y,z);
return res;
}
inline double f(const TinyVector<double,3>& pos)
{
return f(pos[0],pos[1],pos[2]);
}
inline double d2f(const TinyVector<double,3>& pos)
{
return d2f(pos[0],pos[1],pos[2]);
}
};
int main(int argc, char** argv)
{
double ri = 0.0;
double rf = 1.0;
std::vector<int> npts(3);
npts[0]=51;
npts[1]=51;
npts[2]=51;
double xcut=0.23;
double ycut=0.67;
const int nk0=1;
const int nk1=1;
const int nk2=1;
//Create one-dimensional grids for three orthogonal directions
typedef LinearGrid<double> GridType;
GridType gridX, gridY, gridZ;
gridX.set(ri,rf,npts[0]);
gridY.set(ri,rf,npts[1]);
gridZ.set(ri,rf,npts[2]);
//Create an analytic function for assignment
ComboFunc infunc;
infunc.push_back(0.5,new TestFunc(1,1,1));
//infunc.push_back(0.3,new TestFunc(1,1,2));
//infunc.push_back(0.1,new TestFunc(1,2,1));
//infunc.push_back(0.01,new TestFunc(2,1,1));
//infunc.push_back(0.01,new TestFunc(2,2,1));
//infunc.push_back(0.001,new TestFunc(2,1,2));
//infunc.push_back(0.001,new TestFunc(2,2,2));
//infunc.push_back(0.001,new TestFunc(5,5,5));
//infunc.push_back(-0.3,new TestFunc(7,2,3));
//infunc.push_back(0.01,new TestFunc(7,7,7));
//infunc.push_back(0.001,new TestFunc(5,5,5));
//Write to an array
std::vector<double> inData(npts[0]*npts[1]*npts[2]);
std::vector<double>::iterator it(inData.begin());
Pooma::Clock timer;
timer.start();
//Assign the values
for(int ix=0; ix<npts[0]; ix++)
{
double x(gridX(ix));
for(int iy=0; iy<npts[1]; iy++)
{
double y(gridY(iy));
for(int iz=0; iz<npts[2]; iz++)
{
(*it)=infunc.f(x,y,gridZ(iz));
++it;
}
}
}
timer.stop();
std::cout << "Time to evaluate " << timer.cpu_time() << std::endl;
//Test TriCubicSplineT function
//Create XYZCubicGrid
XYZCubicGrid<double> grid3(&gridX,&gridY,&gridZ);
//Create a TriCubicSpline with PBC: have to think more about fixed-boundary conditions
TriCubicSplineT<double> aorb(&grid3);
//Reset the coefficients
aorb.reset(inData.begin(), inData.end());
double lap,val;
TinyVector<double,3> grad;
//aorb.reset();
//Write for vtk ImageData
std::string fname("spline3d.vti");
std::ofstream dfile(fname.c_str());
dfile.setf(std::ios::scientific, std::ios::floatfield);
dfile.setf(std::ios::left,std::ios::adjustfield);
dfile.precision(10);
dfile << "<?xml version=\"1.0\"?>" << std::endl;
dfile << "<VTKFile type=\"ImageData\" version=\"0.1\">" << std::endl;
dfile << " <ImageData WholeExtent=\"0 " << npts[0]-2 << " 0 " << npts[1]-2 << " 0 " << npts[2]-2
<< "\" Origin=\"0 0 0\" Spacing=\"1 1 1\">"<< std::endl;
dfile << " <Piece Extent=\"0 " << npts[0]-2 << " 0 " << npts[1]-2 << " 0 " << npts[2]-2 << "\">" << std::endl;
dfile << " <PointData Scalars=\"wfs\">" << std::endl;
dfile << " <DataArray type=\"Float32\" Name=\"wfs\">" << std::endl;
timer.start();
int ng=0;
for(int ix=0; ix<npts[0]-1; ix++)
{
double x(gridX(ix));
for(int iy=0; iy<npts[1]-1; iy++)
{
double y(gridY(iy));
for(int iz=0; iz<npts[2]-1; iz++, ng++)
{
TinyVector<double,3> p(x,y,gridZ(iz));
//aorb.setgrid(p);
//Timing with the std::ofstream is not correct.
//Uncomment the line below and comment out the next two line.
//double t=aorb.evaluate(p,grad,lap);
dfile << std::setw(20) << aorb.evaluate(p,grad,lap);
if(ng%5 == 4)
dfile << std::endl;
}
}
}
timer.stop();
std::cout << "Time to evaluate with spline " << timer.cpu_time() << std::endl;
dfile << " </DataArray>" << std::endl;
dfile << " </PointData>" << std::endl;
dfile << " </Piece>" << std::endl;
dfile << " </ImageData>" << std::endl;
dfile << "</VTKFile>" << std::endl;
hid_t h_file = H5Fcreate("spline3d.h5",H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
HDFAttribIO<std::vector<double> > dump(inData,npts);
dump.write(h_file,"orb0000");
HDFAttribIO<TriCubicSplineT<double> > dump1(aorb);
dump1.write(h_file,"spline0000");
H5Fclose(h_file);
//double lap;
//TinyVector<double,3> grad;
//for(int k=0; k<nptY-1; k++) {
// //TinyVector<double,3> p(xcut,ycut,gridZ(k)+0.11*gridZ.dr(k));
// TinyVector<double,3> p(xcut,gridY(k)+0.11*gridY.dr(k),ycut);
// aorb.setgrid(p);
// double y=aorb.evaluate(p,grad,lap);
// dfile << std::setw(30) << p[1] << std::setw(30) << infunc.f(p) << std::setw(30) << y << std::setw(30) << infunc.d2f(p) << std::setw(30) << lap << std::endl;
//}
return 0;
}

102
src/QMCTools/spline3_io.cpp
View File

@ -1,102 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <sstream>
#include "OhmmsData/OhmmsElementBase.h"
#include "Numerics/HDFSTLAttrib.h"
#include "Numerics/TriCubicSplineT.h"
#include "Numerics/HDFTriCubicSpline.h"
#include "Utilities/Clock.h"
int main(int argc, char** argv)
{
Pooma::Clock timer;
double ri = 0.0;
double rf = 1.0;
std::vector<int> npts(3);
npts[0]=51;
npts[1]=51;
npts[2]=51;
double xcut=0.23;
double ycut=0.67;
const int nk0=1;
const int nk1=1;
const int nk2=1;
//Create one-dimensional grids for three orthogonal directions
typedef LinearGrid<double> GridType;
GridType gridX, gridY, gridZ;
gridX.set(ri,rf,npts[0]);
gridY.set(ri,rf,npts[1]);
gridZ.set(ri,rf,npts[2]);
//Write to an array
std::vector<double> inData(npts[0]*npts[1]*npts[2]);
timer.start();
hid_t h_file = H5Fopen("spline3d.h5",H5F_ACC_RDWR,H5P_DEFAULT);
HDFAttribIO<std::vector<double> > dummy(inData,npts);
dummy.read(h_file,"orb0000");
H5Fclose(h_file);
timer.stop();
std::cout << "Time to read image data " << timer.cpu_time() << std::endl;
//Create XYZCubicGrid
XYZCubicGrid<double> grid3(&gridX,&gridY,&gridZ);
//Create a TriCubicSpline with PBC: have to think more about fixed-boundary conditions
TriCubicSplineT<double> aorb(&grid3);
//Reset the coefficients
timer.start();
aorb.reset(inData.begin(), inData.end());
timer.stop();
std::cout << "Time to set up spline coefficients " << timer.cpu_time() << std::endl;
double lap,val;
TinyVector<double,3> grad;
//vector<double>::iterator it(inData.begin());
timer.start();
for(int ix=0; ix<npts[0]-1; ix++)
{
double x(gridX(ix));
for(int iy=0; iy<npts[1]-1; iy++)
{
double y(gridY(iy));
int ng = npts[2]*(iy+ix*npts[1]);
for(int iz=0; iz<npts[2]-1; iz++)
{
TinyVector<double,3> p(x,y,gridZ(iz));
//aorb.setgrid(p);
inData[ng++]=aorb.evaluate(p,grad,lap);
//(*it) = aorb.evaluate(p,grad,lap); ++it;
}
}
}
timer.stop();
std::cout << "Time to evaluate the values " << timer.cpu_time() << std::endl;
timer.start();
h_file = H5Fcreate("spline3d_writeback.h5",H5F_ACC_TRUNC,H5P_DEFAULT,H5P_DEFAULT);
HDFAttribIO<std::vector<double> > dump(inData,npts);
dump.write(h_file,"orb0000");
HDFAttribIO<TriCubicSplineT<double> > dump1(aorb);
dump1.write(h_file,"spline0000");
H5Fclose(h_file);
timer.stop();
std::cout << "Time to write to hdf5 " << timer.cpu_time() << std::endl;
return 0;
}

1434
src/QMCTools/trace_density.cpp
View File

File diff suppressed because it is too large Load Diff

60
src/QMCTools/ylmtest.cpp
View File

@ -1,60 +0,0 @@
//////////////////////////////////////////////////////////////////////////////////////
// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
//
// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
//////////////////////////////////////////////////////////////////////////////////////
#include <vector>
#include <string>
#include <iostream>
#include <iomanip>
#include "OhmmsPETE/TinyVector.h"
#include "Numerics/SphericalTensor.h"
int main(int argc, char** argv)
{
int l=5;
TinyVector<double,3> pos(0.1,0.3,-0.45), deltax(0.0001,0.0,0.0), deltay(0.0,0.0001,0.0), deltaz(0.0,0.0,0.0001), tpos, g;
if(argc>1)
l = atoi(argv[1]);
if(argc>4)
{
pos[0] = atof(argv[2]);
pos[1] = atof(argv[3]);
pos[2] = atof(argv[4]);
}
SphericalTensor<double,TinyVector<double,3> > Y1(l,true), Y2(l,true), Yp(l,true), Ym(l,true);
Y1.evaluate(pos);
std::cout.setf(std::ios::scientific, std::ios::floatfield);
for(int lm=0; lm<Y1.size(); lm++)
{
std::cout << lm << std::endl;
std::cout << std::setw(20) << std::setprecision(12) << Y1.Ylm[lm]
<< std::setprecision(12) << Y1.gradYlm[lm] << std::endl;
//std::cout << std::setw(20) << std::setprecision(12) << Y2.Ylm[lm]
// << std::setprecision(12) << Y2.gradYlm[lm] << std::endl;
}
// for(int lm=0; lm<Y1.size(); lm++) {
// tpos = pos+deltax; Yp.evaluate(tpos);
// tpos = pos-deltax; Ym.evaluate(tpos);
// g[0] = (Yp.Ylm[lm]-Ym.Ylm[lm])/0.0002;
// tpos = pos+deltay; Yp.evaluate(tpos);
// tpos = pos-deltay; Ym.evaluate(tpos);
// g[1] = (Yp.Ylm[lm]-Ym.Ylm[lm])/0.0002;
// tpos = pos+deltaz; Yp.evaluate(tpos);
// tpos = pos-deltaz; Ym.evaluate(tpos);
// g[2] = (Yp.Ylm[lm]-Ym.Ylm[lm])/0.0002;
// std::cout << lm << std::endl;
// std::cout << std::setw(20) << std::setprecision(12) << Y1.Ylm[lm]
// << std::setprecision(12) << Y1.gradYlm[lm] - g << std::endl;
// }
}