QMCPACK
ECPComponentBuilder Struct Reference
+ Inheritance diagram for ECPComponentBuilder:
+ Collaboration diagram for ECPComponentBuilder:

Public Types

using GridType = LocalECPotential::GridType
 
using mRealType = ParticleSet::Scalar_t
 
using mGridType = OneDimGridBase< mRealType >
 
using RadialPotentialType = LocalECPotential::RadialPotentialType
 
- Public Types inherited from MPIObjectBase
using mpi_comm_type = Communicate::mpi_comm_type
 
- Public Types inherited from QMCTraits
enum  { DIM = OHMMS_DIM, DIM_VGL = OHMMS_DIM + 2 }
 
using QTBase = QMCTypes< OHMMS_PRECISION, DIM >
 
using QTFull = QMCTypes< OHMMS_PRECISION_FULL, DIM >
 
using RealType = QTBase::RealType
 
using ComplexType = QTBase::ComplexType
 
using ValueType = QTBase::ValueType
 
using PosType = QTBase::PosType
 
using GradType = QTBase::GradType
 
using TensorType = QTBase::TensorType
 
using IndexType = OHMMS_INDEXTYPE
 define other types More...
 
using FullPrecRealType = QTFull::RealType
 
using FullPrecValueType = QTFull::ValueType
 
using PropertySetType = RecordNamedProperty< FullPrecRealType >
 define PropertyList_t More...
 
using PtclGrpIndexes = std::vector< std::pair< int, int > >
 

Public Member Functions

 ECPComponentBuilder (const std::string &aname, Communicate *c, int nrule=-1, int llocal=-1, int srule=8)
 constructor spin grid used for numerical integration. More...
 
bool parse (const std::string &fname, xmlNodePtr cur)
 
bool put (xmlNodePtr cur)
 
void addSemiLocal (xmlNodePtr cur)
 
void buildLocal (xmlNodePtr cur)
 build a Local Pseudopotential More...
 
void buildSemiLocalAndLocal (std::vector< xmlNodePtr > &semiPtr)
 
void buildL2 (xmlNodePtr cur)
 
bool parseCasino (const std::string &fname, xmlNodePtr cur)
 
void SetQuadratureRule (int rule)
 
std::unique_ptr< mGridTypecreateGrid (xmlNodePtr cur, bool useLinear=false)
 
RadialPotentialTypecreateVrWithBasisGroup (xmlNodePtr cur, mGridType *agrid)
 
void doBreakUp (const std::vector< int > &angList, const Matrix< mRealType > &vnn, RealType rmax, mRealType Vprefactor=1.0)
 Separate local from non-local potentials. More...
 
void buildSO (const std::vector< int > &angList, const Matrix< mRealType > &vnnso, RealType rmax, mRealType Vprefactor=1.0)
 brief buildSO - takes the previously parsed angular momenta and spin-orbit tabulated potentials and uses them to construct SOECPComponent* pp_so. More...
 
void printECPTable ()
 
bool read_pp_file (const std::string &fname)
 
- Public Member Functions inherited from MPIObjectBase
 MPIObjectBase (Communicate *c)
 constructor with communicator More...
 
int rank () const
 return the rank of the communicator More...
 
int getGroupID () const
 return the group id of the communicator More...
 
CommunicategetCommunicator () const
 return myComm More...
 
CommunicategetCommRef () const
 return a TEMPORARY reference to Communicate More...
 
mpi_comm_type getMPI () const
 return MPI communicator if one wants to use MPI directly More...
 
bool is_manager () const
 return true if the rank == 0 More...
 
const std::string & getName () const
 return the name More...
 
void setName (const std::string &aname)
 

Public Attributes

int NumNonLocal
 
int Lmax
 
int Llocal
 
int Nrule
 
int Srule
 
int NumSO
 
int LmaxSO
 
int AtomicNumber
 
RealType Zeff
 
RealType RcutMax
 
std::string Species
 
std::unique_ptr< mGridTypegrid_global
 
std::map< std::string, std::unique_ptr< mGridType > > grid_inp
 
std::unique_ptr< RadialPotentialTypepp_loc
 
std::unique_ptr< NonLocalECPComponentpp_nonloc
 
std::unique_ptr< SOECPComponentpp_so
 
std::unique_ptr< L2RadialPotentialpp_L2
 
std::map< std::string, int > angMon
 

Additional Inherited Members

- Protected Attributes inherited from MPIObjectBase
CommunicatemyComm
 pointer to Communicate More...
 
std::string ClassName
 class Name More...
 
std::string myName
 name of the object More...
 

Detailed Description

Definition at line 28 of file ECPComponentBuilder.h.

Member Typedef Documentation

◆ GridType

Definition at line 30 of file ECPComponentBuilder.h.

◆ mGridType

Definition at line 32 of file ECPComponentBuilder.h.

◆ mRealType

Definition at line 31 of file ECPComponentBuilder.h.

◆ RadialPotentialType

Constructor & Destructor Documentation

◆ ECPComponentBuilder()

ECPComponentBuilder ( const std::string &  aname,
Communicate c,
int  nrule = -1,
int  llocal = -1,
int  srule = 8 
)

constructor spin grid used for numerical integration.

use 0 for exact integration.

Definition at line 31 of file ECPComponentBuilder.cpp.

References ECPComponentBuilder::angMon.

32  : MPIObjectBase(c),
33  NumNonLocal(0),
34  Lmax(0),
35  Llocal(llocal),
36  Nrule(nrule),
37  Srule(srule),
38  AtomicNumber(0),
39  Zeff(0),
40  RcutMax(-1),
41  Species(aname)
42 {
43  angMon["s"] = 0;
44  angMon["p"] = 1;
45  angMon["d"] = 2;
46  angMon["f"] = 3;
47  angMon["g"] = 4;
48  angMon["h"] = 5;
49  angMon["i"] = 6;
50  angMon["j"] = 7;
51  angMon["k"] = 8;
52  angMon["0"] = 0;
53  angMon["1"] = 1;
54  angMon["2"] = 2;
55  angMon["3"] = 3;
56  angMon["4"] = 4;
57  angMon["5"] = 5;
58  angMon["6"] = 6;
59  angMon["7"] = 7;
60  angMon["8"] = 8;
61 }
std::map< std::string, int > angMon
MPIObjectBase(Communicate *c)
constructor with communicator

Member Function Documentation

◆ addSemiLocal()

void addSemiLocal ( xmlNodePtr  cur)

Definition at line 22 of file ECPComponentBuilder.1.cpp.

References ECPComponentBuilder::angMon, ECPComponentBuilder::createGrid(), ECPComponentBuilder::createVrWithBasisGroup(), getXMLAttributeValue(), ECPComponentBuilder::Lmax, ECPComponentBuilder::NumNonLocal, and ECPComponentBuilder::pp_nonloc.

Referenced by ECPComponentBuilder::put().

23 {
24  std::unique_ptr<mGridType> grid_semilocal;
25  RealType rmax = pp_nonloc->Rmax;
26  cur = cur->children;
27  while (cur != NULL)
28  {
29  std::string cname((const char*)cur->name);
30  if (cname == "grid")
31  {
32  grid_semilocal = createGrid(cur);
33  rmax = grid_semilocal->rmax();
34  }
35  else if (cname == "vps")
36  {
37  //should be able to overwrite rmax
38  int l = angMon[getXMLAttributeValue(cur, "l")];
39  Lmax = std::max(l, Lmax);
40  xmlNodePtr cur1 = cur->children;
41  while (cur1 != NULL)
42  {
43  std::string cname1((const char*)cur1->name);
44  if (cname1 == "basisGroup")
45  {
46  pp_nonloc->add(l, createVrWithBasisGroup(cur1, grid_semilocal.get()));
47  }
48  cur1 = cur1->next;
49  }
50  NumNonLocal++;
51  }
52  cur = cur->next;
53  }
54  pp_nonloc->lmax = Lmax;
55  pp_nonloc->Rmax = rmax;
56 }
std::map< std::string, int > angMon
std::unique_ptr< NonLocalECPComponent > pp_nonloc
std::string getXMLAttributeValue(const xmlNodePtr cur, const std::string_view name)
get the value string for attribute name if name is unfound in cur you get an empty string back this i...
QMCTraits::RealType RealType
RadialPotentialType * createVrWithBasisGroup(xmlNodePtr cur, mGridType *agrid)
std::unique_ptr< mGridType > createGrid(xmlNodePtr cur, bool useLinear=false)

◆ buildL2()

void buildL2 ( xmlNodePtr  cur)

Definition at line 16 of file ECPComponentBuilder_L2.cpp.

References OhmmsAttributeSet::add(), qmcplusplus::app_error(), qmcplusplus::app_log(), Communicate::barrier_and_abort(), qmcplusplus::ceil(), qmcplusplus::Units::charge::e, ECPComponentBuilder::grid_global, LINEAR_1DGRID, omptarget::min(), MPIObjectBase::myComm, ECPComponentBuilder::pp_L2, OhmmsAttributeSet::put(), putContent(), and OneDimCubicSpline< Td, Tg, CTd, CTg >::spline().

Referenced by ECPComponentBuilder::put().

17 {
18  // report to log and check for global grid (pseudo/grid)
19  app_log() << " ECPComponentBuilder::buildL2 " << std::endl;
20  if (grid_global == 0)
21  {
22  app_error() << " Global grid needs to be defined." << std::endl;
23  myComm->barrier_and_abort("ECPComponentBuilder::buildL2");
24  }
25 
26  // read <L2> attributes
27  std::string eunits = "hartree";
28  std::string format = "r*V";
29  RealType rcut = -1.0;
30  OhmmsAttributeSet attrib;
31  attrib.add(eunits, "units");
32  attrib.add(format, "format");
33  attrib.add(rcut, "cutoff");
34  attrib.put(cur);
35 
36  // check validity of <L2> attributes
37  RealType Vprefactor = 1.0;
38  if (eunits.find("ydberg") < eunits.size())
39  {
40  app_log() << " Input energy unit = Rydberg " << std::endl;
41  Vprefactor = 0.5;
42  }
43  else
44  {
45  app_log() << " Assuming Hartree unit" << std::endl;
46  }
47  bool is_r_times_V(true);
48  if (format == "r*V")
49  is_r_times_V = true;
50  else if (format == "V")
51  is_r_times_V = false;
52  else
53  {
54  app_error() << " Unrecognized format \"" << format << "\" in PP file." << std::endl;
55  myComm->barrier_and_abort("ECPComponentBuilder::buildL2");
56  }
57  if (rcut < 0.0)
58  {
59  app_error() << " L2 attribute cutoff is missing or negative. Cutoff must be a positive real number." << std::endl;
60  myComm->barrier_and_abort("ECPComponentBuilder::buildL2");
61  }
62  int npts = grid_global->size();
63 
64  // read in the data
65  std::vector<RealType> vL2in(npts);
66  xmlNodePtr c = cur->children;
67  while (c != NULL)
68  {
69  if (xmlStrEqual(c->name, (const xmlChar*)"radfunc"))
70  {
71  xmlNodePtr c1 = c->children;
72  while (c1 != NULL)
73  {
74  if (xmlStrEqual(c1->name, (const xmlChar*)"data"))
75  putContent(vL2in, c1);
76  c1 = c1->next;
77  }
78  }
79  c = c->next;
80  }
81 
82  // convert to r*V if needed
83  if (!is_r_times_V)
84  {
85  app_log() << " Input pseudopotential is converted into r*V" << std::endl;
86  for (int i = 0; i < npts; i++)
87  vL2in[i] *= grid_global->r(i);
88  }
89 
90  // create a new grid with maximum extent set to the cutoff radius
91  RealType rmax = rcut;
92  const int max_points = 100000;
93  app_log() << " Creating a Linear Grid Rmax=" << rmax << std::endl;
94  auto grid = std::make_unique<LinearGrid<RealType>>();
95  RealType d = 1e-4;
96  int ng;
97  if (grid_global->getGridTag() == LINEAR_1DGRID)
98  {
99  ng = (int)std::ceil(rmax * grid_global->DeltaInv) + 1;
100  if (ng <= max_points)
101  {
102  app_log() << " Using global grid with delta = " << grid_global->Delta << std::endl;
103  rmax = grid_global->Delta * (ng - 1);
104  grid->set(0.0, rmax, ng);
105  }
106  else
107  grid->set(0.0, rmax, max_points);
108  }
109  else
110  {
111  ng = std::min(max_points, static_cast<int>(rmax / d) + 1);
112  grid->set(0, rmax, ng);
113  }
114  d = grid->Delta;
115 
116  // follow doBreakUp and reinterpolate the data
117  // likely necessary for data inputted on other than a linear grid
118  // but seems extraneous for data on linear grid (most common case)
119  int ngIn = npts - 2;
120  std::vector<RealType> new_vL2(ng);
121  std::vector<RealType> new_vL2in(ngIn);
122  for (int i = 0; i < ngIn; i++)
123  new_vL2in[i] = Vprefactor * vL2in[i];
124  OneDimCubicSpline<mRealType> infunc(grid_global->makeClone(), new_vL2in);
125  infunc.spline(0, 0.0, ngIn - 1, 0.0);
126  for (int i = 1; i < ng - 1; i++)
127  {
128  RealType r = d * i;
129  new_vL2[i] = infunc.splint(r) / r;
130  }
131  new_vL2[0] = new_vL2[1];
132  new_vL2[ng - 1] = 0.0;
133  auto vL2 = std::make_unique<RadialPotentialType>(std::move(grid), new_vL2);
134  vL2->spline();
135 
136  // save the splined L2 potential
137  pp_L2 = std::make_unique<L2RadialPotential>();
138  pp_L2->vL2 = std::move(vL2);
139  pp_L2->rcut = rcut;
140 
141  //app_log()<<std::endl;
142  //for(int i=0;i<ng;i++)
143  // app_log()<<d*i<<" "<<pp_L2->vL2->splint(d*i)*d*i<<std::endl;
144  //app_log()<<std::endl;
145  //app_log()<<pp_L2->rcut;
146  //APP_ABORT("here");
147 }
std::ostream & app_log()
Definition: OutputManager.h:65
bool put(xmlNodePtr cur)
assign attributes to the set
Definition: AttributeSet.h:55
std::ostream & app_error()
Definition: OutputManager.h:67
T min(T a, T b)
std::unique_ptr< mGridType > grid_global
MakeReturn< UnaryNode< FnCeil, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t ceil(const Vector< T1, C1 > &l)
Communicate * myComm
pointer to Communicate
Definition: MPIObjectBase.h:62
class to handle a set of attributes of an xmlNode
Definition: AttributeSet.h:24
bool putContent(T &a, xmlNodePtr cur)
replaces a&#39;s value with the first "element" in the "string" returned by XMLNodeString{cur}.
Definition: libxmldefs.h:88
QMCTraits::RealType RealType
std::unique_ptr< L2RadialPotential > pp_L2
void barrier_and_abort(const std::string &msg) const
void add(PDT &aparam, const std::string &aname, std::vector< PDT > candidate_values={}, TagStatus status=TagStatus::OPTIONAL)
add a new attribute
Definition: AttributeSet.h:42

◆ buildLocal()

void buildLocal ( xmlNodePtr  cur)

build a Local Pseudopotential

For numerical stability, the radial function of the local pseudopotential is stored as $rV(r)/Z_{eff}/Z_{e}$ for the distance r and $Z_e = -1$. The coulomb factor $Z_{e}Z_{eff}/r$ is applied by LocalECPotential (see LocalECPotential::evaluate).

Definition at line 104 of file ECPComponentBuilder.1.cpp.

References qmcplusplus::abs(), qmcplusplus::app_log(), Communicate::barrier_and_abort(), ECPComponentBuilder::createGrid(), CUSTOM_1DGRID, qmcplusplus::Units::charge::e, getXMLAttributeValue(), ECPComponentBuilder::grid_global, MPIObjectBase::myComm, ECPComponentBuilder::pp_loc, and ECPComponentBuilder::Zeff.

Referenced by ECPComponentBuilder::put().

105 {
106  if (pp_loc)
107  return; //something is wrong
108 
109  std::string vFormat("V");
110  const std::string v_str(getXMLAttributeValue(cur, "format"));
111  if (!v_str.empty())
112  vFormat = v_str;
113 
114  int vPowerCorrection = 1;
115  if (vFormat == "r*V")
116  {
117  app_log() << " Local pseudopotential format = r*V" << std::endl;
118  vPowerCorrection = 0;
119  }
120  else
121  {
122  app_log() << " Local pseudopotential format = V" << std::endl;
123  }
124  using InFuncType = GaussianTimesRN<RealType>;
125  std::unique_ptr<GridType> grid_local;
126  std::unique_ptr<mGridType> grid_local_inp;
127  InFuncType vr;
128  bool bareCoulomb = true;
129  cur = cur->children;
130  while (cur != NULL)
131  {
132  std::string cname((const char*)cur->name);
133  if (cname == "grid")
134  {
135  grid_local_inp = createGrid(cur, true);
136  }
137  else if (cname == "basisGroup")
138  {
139  vr.putBasisGroup(cur, vPowerCorrection);
140  bareCoulomb = false;
141  }
142  cur = cur->next;
143  }
144  if (grid_local_inp == nullptr)
145  {
146  if (grid_global == nullptr)
147  myComm->barrier_and_abort("ECPComponentBuilder::buildLocal Missing grid information. ");
148  grid_local = std::make_unique<LinearGrid<RealType>>();
149  grid_local->set(grid_global->rmin(), grid_global->rmax(), grid_global->size());
150  }
151  else
152  {
153  grid_local = std::make_unique<LinearGrid<RealType>>();
154  grid_local->set(grid_local_inp->rmin(), grid_local_inp->rmax(), grid_local_inp->size());
155  }
156  if (grid_local->GridTag == CUSTOM_1DGRID)
157  myComm->barrier_and_abort("ECPComponentBuilder::buildLocal Custom grid is used. Need to recast to the linear grid");
158  else
159  {
160  std::vector<RealType> v;
161  if (bareCoulomb)
162  {
163  app_log() << " Bare Coulomb potential is used." << std::endl;
164  grid_local->set(0.0, 1., 3);
165  v.resize(3);
166  for (int ig = 0; ig < 3; ig++)
167  v[ig] = 1.0;
168  pp_loc = std::make_unique<RadialPotentialType>(std::move(grid_local), v);
169  pp_loc->spline(0, 0.0, 2, 0.0);
170  }
171  else
172  {
173  app_log() << " Guassian basisGroup is used: base power " << vr.basePower << std::endl;
174  const RealType eps = 1e-12; //numeric_limits<RealType>::epsilon();//1e-12;
175  RealType zinv = 1.0 / Zeff;
176  RealType r = 10.;
177  bool ignore = true;
178  int last = grid_local->size() - 1;
179  while (ignore && last)
180  {
181  r = (*grid_local)[last];
182  ignore = (std::abs(zinv * vr.f(r)) < eps);
183  --last;
184  }
185  if (last == 0)
186  myComm->barrier_and_abort("ECPComponentBuilder::buildLocal. Illegal Local Pseudopotential");
187  //Add the reset values here
188  int ng = static_cast<int>(r / 1e-3) + 1;
189  app_log() << " Use a Linear Grid: [0," << r << "] Number of points = " << ng << std::endl;
190  grid_local->set(0.0, r, ng);
191  v.resize(ng);
192  for (int ig = 1; ig < ng - 1; ig++)
193  {
194  double r = (*grid_local)[ig];
195  v[ig] = 1.0 - zinv * vr.f(r);
196  }
197  v[0] = 2.0 * v[1] - v[2];
198  v[ng - 1] = 1.0;
199  pp_loc = std::make_unique<RadialPotentialType>(std::move(grid_local), v);
200  pp_loc->spline(); //use the fixed conditions
201  }
202  }
203 }
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
std::ostream & app_log()
Definition: OutputManager.h:65
std::unique_ptr< mGridType > grid_global
Communicate * myComm
pointer to Communicate
Definition: MPIObjectBase.h:62
std::unique_ptr< RadialPotentialType > pp_loc
std::string getXMLAttributeValue(const xmlNodePtr cur, const std::string_view name)
get the value string for attribute name if name is unfound in cur you get an empty string back this i...
QMCTraits::RealType RealType
void barrier_and_abort(const std::string &msg) const
std::unique_ptr< mGridType > createGrid(xmlNodePtr cur, bool useLinear=false)

◆ buildSemiLocalAndLocal()

void buildSemiLocalAndLocal ( std::vector< xmlNodePtr > &  semiPtr)

Definition at line 27 of file ECPComponentBuilder.2.cpp.

References OhmmsAttributeSet::add(), ECPComponentBuilder::angMon, qmcplusplus::app_error(), qmcplusplus::app_log(), qmcplusplus::app_warning(), Communicate::barrier_and_abort(), ECPComponentBuilder::buildSO(), copy(), ECPComponentBuilder::doBreakUp(), ECPComponentBuilder::grid_global, ECPComponentBuilder::Llocal, ECPComponentBuilder::Lmax, ECPComponentBuilder::LmaxSO, MPIObjectBase::myComm, ECPComponentBuilder::Nrule, ECPComponentBuilder::pp_nonloc, ECPComponentBuilder::pp_so, OhmmsAttributeSet::put(), putContent(), and ECPComponentBuilder::Srule.

Referenced by ECPComponentBuilder::put().

28 {
29  app_log() << " ECPComponentBuilder::buildSemiLocalAndLocal " << std::endl;
30  if (grid_global == 0)
31  myComm->barrier_and_abort("ECPComponentBuilder::buildSemiLocalAndLocal. Global grid needs to be defined.\n");
32  // There should only be one semilocal tag
33  if (semiPtr.size() > 1)
34  {
35  std::stringstream err_msg;
36  err_msg << "ECPComponentBuilder::buildSemiLocalAndLocal. "
37  << "We have more than one semilocal sections in the PP xml file";
38  myComm->barrier_and_abort(err_msg.str());
39  }
40  RealType rmax = -1;
41  //attributes: initailize by defaults
42  std::string eunits("hartree");
43  std::string format("r*V");
44  std::string lloc;
45  int ndown = 1;
46  int nup = 0;
47  int nso = 0;
48  OhmmsAttributeSet aAttrib;
49  int quad_rule = -1;
50  int local_channel = -1;
51  aAttrib.add(eunits, "units");
52  aAttrib.add(format, "format");
53  aAttrib.add(ndown, "npots-down");
54  aAttrib.add(nup, "npots-up");
55  aAttrib.add(local_channel, "l-local");
56  aAttrib.add(quad_rule, "nrule");
57  aAttrib.add(Srule, "srule");
58  aAttrib.add(nso, "npots-so");
59 
60  xmlNodePtr cur_semilocal = semiPtr[0];
61  aAttrib.put(cur_semilocal);
62 
63  // settle Nrule. Priority: current value (from input file) > PP XML file > lmax derived
64  if (quad_rule > -1 && Nrule > -1)
65  {
66  app_warning() << " Nrule setting found in both qmcpack input (Nrule = " << Nrule
67  << ") and pseudopotential file (Nrule = " << quad_rule << ")."
68  << " Using nrule setting in qmcpack input file." << std::endl;
69  }
70  else if (quad_rule > -1 && Nrule == -1)
71  {
72  app_log() << " Nrule setting found in pseudopotential file and used." << std::endl;
73  Nrule = quad_rule;
74  }
75  else if (quad_rule == -1 && Nrule > -1)
76  app_log() << " Nrule setting found in qmcpack input file and used." << std::endl;
77  else
78  {
79  //Sperical quadrature rules set by exact integration up to lmax of
80  //nonlocal channels.
81  //From J. Chem. Phys. 95 (3467) (1991)
82  //Keeping Nrule = 4 as default for lmax <= 5.
83  switch (pp_nonloc->lmax)
84  {
85  case 0:
86  case 1:
87  case 2:
88  case 3:
89  case 4:
90  case 5:
91  Nrule = 4;
92  break;
93  case 6:
94  case 7:
95  Nrule = 6;
96  break;
97  case 8:
98  case 9:
99  case 10:
100  case 11:
101  Nrule = 7;
102  break;
103  default:
104  myComm->barrier_and_abort("Default value for pseudopotential nrule not determined.");
105  break;
106  }
107  app_warning() << "Nrule was not determined from qmcpack input or pseudopotential file. Setting sensible default."
108  << std::endl;
109  }
110 
111 
112  RealType Vprefactor = 1.0;
113  if (eunits.find("ydberg") < eunits.size())
114  {
115  app_log() << " Input energy unit = Rydberg " << std::endl;
116  Vprefactor = 0.5;
117  }
118  else
119  {
120  app_log() << " Assuming Hartree unit" << std::endl;
121  }
122  bool is_r_times_V(true);
123  if (format == "r*V")
124  is_r_times_V = true;
125  else if (format == "V")
126  is_r_times_V = false;
127  else
128  {
129  std::stringstream err_msg;
130  err_msg << "ECPComponentBuilder::buildSemiLocalAndLocal."
131  << "Unrecognized format \"" << format << "\" in PP file.\n";
132  myComm->barrier_and_abort(err_msg.str());
133  }
134  // We cannot construct the potentials as we construct them since
135  // we may not know which one is local yet.
136 
137  std::vector<int> angList;
138  std::vector<int> angListSO; //For spin-orbit, if it exists
139  std::vector<xmlNodePtr> vpsPtr;
140  std::vector<xmlNodePtr> vpsoPtr; //For spin-orbit, if it exists.
141  Lmax = -1;
142  LmaxSO = -1;
143  // Now read vps sections
144  xmlNodePtr cur_vps = cur_semilocal->children;
145  while (cur_vps != NULL)
146  {
147  std::string vname((const char*)cur_vps->name);
148  if (vname == "vps")
149  {
150  OhmmsAttributeSet aAttrib;
151  std::string lstr("s");
152  RealType rc = -1.0;
153  aAttrib.add(lstr, "l");
154  aAttrib.add(rc, "cutoff");
155  aAttrib.put(cur_vps);
156  rmax = std::max(rmax, rc);
157  if (angMon.find(lstr) == angMon.end())
158  {
159  std::stringstream err_msg;
160  err_msg << "ECPComponentBuilder::buildSemiLocalAndLocal. "
161  << "Requested angular momentum " << lstr << " not available.\n";
162  myComm->barrier_and_abort(err_msg.str());
163  }
164  int l = angMon[lstr];
165  angList.push_back(l);
166  vpsPtr.push_back(cur_vps);
167  Lmax = std::max(Lmax, l); //count the maximum L
168  }
169  else if (vname == "vps_so") //This accumulates the spin-orbit corrections, if defined.
170  {
171  OhmmsAttributeSet aAttrib;
172  std::string lstr("s");
173  RealType rc = -1.0;
174  aAttrib.add(lstr, "l");
175  aAttrib.add(rc, "cutoff");
176  aAttrib.put(cur_vps);
177  rmax = std::max(rmax, rc);
178  if (angMon.find(lstr) == angMon.end())
179  {
180  std::stringstream err_msg;
181  err_msg << "ECPComponentBuilder::buildSemiLocalAndLocal. "
182  << "Requested angular momentum " << lstr << " not available for SO.\n";
183  myComm->barrier_and_abort(err_msg.str());
184  }
185  int l = angMon[lstr];
186  angListSO.push_back(l);
187  vpsoPtr.push_back(cur_vps);
188  LmaxSO = std::max(LmaxSO, l); //count the maximum L
189  }
190  cur_vps = cur_vps->next;
191  }
192 
193  if (rmax < 0)
194  rmax = 1.8;
195 
196  // settle Llocal. Priority: current value (from input file) > PP XML file
197  if (local_channel > -1 && Llocal > -1)
198  {
199  app_warning() << " l-local setting found in both qmcpack input (l-local = " << Llocal
200  << ") and pseudopotential file (l-local = " << local_channel << ")."
201  << " Using l-local setting in qmcpack input file." << std::endl;
202  }
203  else if (local_channel > -1 && Llocal == -1)
204  {
205  app_log() << " l-local setting found in pseudopotential file and used." << std::endl;
206  Llocal = local_channel;
207  }
208  else if (local_channel == -1 && Llocal > -1)
209  app_log() << " l-local setting found in qmcpack input file and used." << std::endl;
210  else if (angList.size() == 1)
211  {
212  Llocal = Lmax;
213  app_log() << " Only one vps is found. Set the local component=" << Lmax << std::endl;
214  }
215  else
216  {
217  app_error() << "The local channel is specified in neither the pseudopotential file nor the input file.\n"
218  << "Please add \'l-local=\"n\"\' attribute to either file.\n";
219  myComm->barrier_and_abort("ECPComponentBuilder::doBreakUp");
220  }
221 
222  if (angListSO.size() != nso)
223  {
224  std::stringstream ssout;
225  ssout << "Error. npots-so=" << angListSO.size() << " while declared number of SO channels is " << nso << std::endl;
226  std::string outstring("");
227  outstring = ssout.str();
228 
229  myComm->barrier_and_abort(outstring.c_str());
230  }
231  int npts = grid_global->size();
232  Matrix<mRealType> vnn(angList.size(), npts);
233  for (int l = 0; l < angList.size(); l++)
234  {
235  std::vector<mRealType> vt(npts);
236  xmlNodePtr c = vpsPtr[l]->children;
237  while (c != NULL)
238  {
239  if (xmlStrEqual(c->name, (const xmlChar*)"radfunc"))
240  {
241  xmlNodePtr c1 = c->children;
242  while (c1 != NULL)
243  {
244  if (xmlStrEqual(c1->name, (const xmlChar*)"data"))
245  putContent(vt, c1);
246  c1 = c1->next;
247  }
248  }
249  c = c->next;
250  }
251  //copy the numerical data with the correct map
252  copy(vt.begin(), vt.end(), vnn[angList[l]]);
253  }
254 
255  //Grabbing the spin-orbit functions from XML.
256  Matrix<mRealType> vnnso(angListSO.size(), npts);
257  for (int l = 0; l < angListSO.size(); l++)
258  {
259  std::vector<mRealType> vtso(npts);
260  xmlNodePtr c = vpsoPtr[l]->children;
261  while (c != NULL)
262  {
263  if (xmlStrEqual(c->name, (const xmlChar*)"radfunc"))
264  {
265  xmlNodePtr c1 = c->children;
266  while (c1 != NULL)
267  {
268  if (xmlStrEqual(c1->name, (const xmlChar*)"data"))
269  putContent(vtso, c1);
270  c1 = c1->next;
271  }
272  }
273  c = c->next;
274  }
275  //copy the numerical data with the correct map
276  //So this is weird, but I feel like l should be the proper index for vnnso,
277  //with angListSO[l] being the actual angular momentum channel referred to by l.
278  //This differs from the parsing of the nonlocal pseudopotential piece, but whatever.
279  copy(vtso.begin(), vtso.end(), vnnso[l]);
280  }
281 
282 
283  ////rather stupid to do this but necessary
284  //vector<RealType> temp(npts);
285  //for(int i=0; i<npts; i++) temp[i]=grid_global->r(i);
286  if (!is_r_times_V)
287  {
288  app_log() << " Input pseudopotential is converted into r*V" << std::endl;
289  for (int i = 0; i < vnn.rows(); i++)
290  for (int j = 0; j < npts; j++)
291  vnn[i][j] *= grid_global->r(j);
292  for (int i = 0; i < vnnso.rows(); i++)
293  for (int j = 0; j < npts; j++)
294  vnnso[i][j] *= grid_global->r(j);
295  }
296  app_log() << " Number of angular momentum channels " << angList.size() << std::endl;
297  app_log() << " Maximum angular momentum channel (Lmax) " << Lmax << std::endl;
298  doBreakUp(angList, vnn, rmax, Vprefactor);
299 
300  //If any spinorbit terms are found...
301  if (nso > 0)
302  buildSO(angListSO, vnnso, rmax, 1.0);
303  else
304  {
305  //No SO channels found. Delete pp_so
306  pp_so.reset();
307  }
308 }
std::ostream & app_warning()
Definition: OutputManager.h:69
std::map< std::string, int > angMon
std::ostream & app_log()
Definition: OutputManager.h:65
bool put(xmlNodePtr cur)
assign attributes to the set
Definition: AttributeSet.h:55
std::ostream & app_error()
Definition: OutputManager.h:67
std::unique_ptr< NonLocalECPComponent > pp_nonloc
void copy(const Array< T1, 3 > &src, Array< T2, 3 > &dest)
Definition: Blitz.h:639
std::unique_ptr< mGridType > grid_global
std::unique_ptr< SOECPComponent > pp_so
Communicate * myComm
pointer to Communicate
Definition: MPIObjectBase.h:62
class to handle a set of attributes of an xmlNode
Definition: AttributeSet.h:24
bool putContent(T &a, xmlNodePtr cur)
replaces a&#39;s value with the first "element" in the "string" returned by XMLNodeString{cur}.
Definition: libxmldefs.h:88
QMCTraits::RealType RealType
void doBreakUp(const std::vector< int > &angList, const Matrix< mRealType > &vnn, RealType rmax, mRealType Vprefactor=1.0)
Separate local from non-local potentials.
void barrier_and_abort(const std::string &msg) const
void add(PDT &aparam, const std::string &aname, std::vector< PDT > candidate_values={}, TagStatus status=TagStatus::OPTIONAL)
add a new attribute
Definition: AttributeSet.h:42
void buildSO(const std::vector< int > &angList, const Matrix< mRealType > &vnnso, RealType rmax, mRealType Vprefactor=1.0)
brief buildSO - takes the previously parsed angular momenta and spin-orbit tabulated potentials and u...

◆ buildSO()

void buildSO ( const std::vector< int > &  angList,
const Matrix< mRealType > &  vnnso,
RealType  rmax,
mRealType  Vprefactor = 1.0 
)

brief buildSO - takes the previously parsed angular momenta and spin-orbit tabulated potentials and uses them to construct SOECPComponent* pp_so.

This is called in "doBreakUp".

param std::vector<int>& angList The angular momentum for each SO potential. param Matrix<mRealType>& vnnso (npot x ngrid) matrix storing all tabulated SO potentials. param RealType rmax max r on the specified grid. param mRealType Vprefactor optional scale factor.

return void

Definition at line 313 of file ECPComponentBuilder.2.cpp.

References qmcplusplus::app_log(), qmcplusplus::ceil(), Matrix< T, Alloc >::cols(), qmcplusplus::Units::charge::e, ECPComponentBuilder::grid_global, LINEAR_1DGRID, omptarget::min(), ECPComponentBuilder::NumSO, ECPComponentBuilder::pp_so, and OneDimCubicSpline< Td, Tg, CTd, CTg >::spline().

Referenced by ECPComponentBuilder::buildSemiLocalAndLocal().

317 {
318  const int max_points = 100000;
319  app_log() << " Creating a Linear Grid Rmax=" << rmax << std::endl;
320  //this is a new grid
321  mRealType d = 1e-4;
322  auto agrid = std::make_unique<LinearGrid<RealType>>();
323  // If the global grid is already linear, do not interpolate the data
324  int ng;
325  if (grid_global->getGridTag() == LINEAR_1DGRID)
326  {
327  ng = (int)std::ceil(rmax * grid_global->DeltaInv) + 1;
328  if (ng <= max_points)
329  {
330  app_log() << " Using global grid with delta = " << grid_global->Delta << std::endl;
331  rmax = grid_global->Delta * (ng - 1);
332  agrid->set(0.0, rmax, ng);
333  }
334  else
335  agrid->set(0.0, rmax, max_points);
336  }
337  else
338  {
339  ng = std::min(max_points, static_cast<int>(rmax / d) + 1);
340  agrid->set(0, rmax, ng);
341  }
342  // This is critical!!!
343  // If d is not reset, we generate an error in the interpolated PP!
344  d = agrid->Delta;
345  int ngIn = vnnso.cols() - 2;
346  std::vector<RealType> newP(ng);
347  std::vector<mRealType> newPin(ngIn);
348  for (int l = 0; l < angList.size(); l++)
349  {
350  const mRealType* restrict vp = vnnso[l];
351  for (int i = 0; i < ngIn; i++)
352  newPin[i] = Vprefactor * vp[i];
353 
354  OneDimCubicSpline<mRealType> infunc(grid_global->makeClone(), newPin);
355  infunc.spline(0, 0.0, ngIn - 1, 0.0);
356  for (int i = 1; i < ng - 1; i++)
357  {
358  mRealType r = d * i;
359  newP[i] = infunc.splint(r) / r;
360  }
361  newP[0] = newP[1];
362  newP[ng - 1] = 0.0;
363  RadialPotentialType* app = new RadialPotentialType(agrid->makeClone(), newP);
364  app->spline();
365  pp_so->add(angList[l], app);
366  }
367  NumSO = angList.size();
368  pp_so->setRmax(rmax);
369 }
std::ostream & app_log()
Definition: OutputManager.h:65
void spline(int imin, value_type yp1, int imax, value_type ypn) override
Evaluate the 2nd derivative on the grid points.
T min(T a, T b)
size_type cols() const
Definition: OhmmsMatrix.h:78
std::unique_ptr< mGridType > grid_global
std::unique_ptr< SOECPComponent > pp_so
MakeReturn< UnaryNode< FnCeil, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t ceil(const Vector< T1, C1 > &l)
LocalECPotential::RadialPotentialType RadialPotentialType

◆ createGrid()

std::unique_ptr< ECPComponentBuilder::mGridType > createGrid ( xmlNodePtr  cur,
bool  useLinear = false 
)

Definition at line 205 of file ECPComponentBuilder.1.cpp.

References OhmmsAttributeSet::add(), qmcplusplus::app_log(), qmcplusplus::Units::charge::e, ECPComponentBuilder::grid_inp, OhmmsAttributeSet::put(), and putContent().

Referenced by ECPComponentBuilder::addSemiLocal(), ECPComponentBuilder::buildLocal(), and ECPComponentBuilder::put().

206 {
207  mRealType ri = 1e-6;
208  mRealType rf = 100.0;
209  mRealType ascale = -1.0e0;
210  mRealType astep = -1.0;
211  //mRealType astep = 1.25e-2;
212  int npts = 1001;
213  std::string gridType("log");
214  std::string gridID("global");
215  OhmmsAttributeSet radAttrib;
216  radAttrib.add(gridType, "type");
217  radAttrib.add(gridID, "grid_id");
218  radAttrib.add(gridID, "grid_def");
219  radAttrib.add(gridID, "name");
220  radAttrib.add(gridID, "id");
221  radAttrib.add(npts, "npts");
222  radAttrib.add(ri, "ri");
223  radAttrib.add(rf, "rf");
224  radAttrib.add(ascale, "ascale");
225  radAttrib.add(astep, "astep");
226  radAttrib.add(ascale, "scale");
227  radAttrib.add(astep, "step");
228  radAttrib.put(cur);
229  auto git = grid_inp.find(gridID);
230  if (git != grid_inp.end())
231  {
232  return (*git).second->makeClone(); //use the same grid
233  }
234  //overwrite the grid type to linear starting at 0.0
235  if (useLinear)
236  {
237  gridType = "linear";
238  ri = 0.0;
239  }
240  std::unique_ptr<mGridType> agrid;
241  if (gridType == "log")
242  {
243  if (ascale > 0.0)
244  {
245  app_log() << " Log grid scale=" << ascale << " step=" << astep << " npts=" << npts << std::endl;
246  agrid = std::make_unique<LogGridZero<mRealType>>();
247  agrid->set(astep, ascale, npts);
248  }
249  else
250  {
251  if (ri < std::numeric_limits<mRealType>::epsilon())
252  {
253  ri = std::numeric_limits<mRealType>::epsilon();
254  }
255  agrid = std::make_unique<LogGrid<mRealType>>();
256  agrid->set(ri, rf, npts);
257  }
258  }
259  else if (gridType == "linear")
260  {
261  agrid = std::make_unique<LinearGrid<mRealType>>();
262  if (astep > 0.0)
263  {
264  npts = static_cast<int>((rf - ri) / astep) + 1;
265  }
266  agrid->set(ri, rf, npts);
267  app_log() << " Linear grid ri=" << ri << " rf=" << rf << " npts = " << npts << std::endl;
268  }
269  else
270  //accept numerical data
271  {
272  xmlNodePtr cur1 = cur->children;
273  while (cur1 != NULL)
274  {
275  std::string cname((const char*)cur1->name);
276  if (cname == "data")
277  {
278  std::vector<double> gIn(npts);
279  putContent(gIn, cur1);
280  agrid = std::make_unique<NumericalGrid<mRealType>>(gIn);
281  app_log() << " Numerical grid npts = " << gIn.size() << std::endl;
282  }
283  cur1 = cur1->next;
284  }
285  }
286  return agrid;
287 }
std::ostream & app_log()
Definition: OutputManager.h:65
bool put(xmlNodePtr cur)
assign attributes to the set
Definition: AttributeSet.h:55
class to handle a set of attributes of an xmlNode
Definition: AttributeSet.h:24
std::map< std::string, std::unique_ptr< mGridType > > grid_inp
bool putContent(T &a, xmlNodePtr cur)
replaces a&#39;s value with the first "element" in the "string" returned by XMLNodeString{cur}.
Definition: libxmldefs.h:88
void add(PDT &aparam, const std::string &aname, std::vector< PDT > candidate_values={}, TagStatus status=TagStatus::OPTIONAL)
add a new attribute
Definition: AttributeSet.h:42

◆ createVrWithBasisGroup()

ECPComponentBuilder::RadialPotentialType * createVrWithBasisGroup ( xmlNodePtr  cur,
mGridType agrid 
)

Definition at line 58 of file ECPComponentBuilder.1.cpp.

References qmcplusplus::abs(), qmcplusplus::app_log(), qmcplusplus::Units::charge::e, GaussianTimesRN< T >::putBasisGroup(), OneDimGridBase< T, CT >::rmax(), OneDimGridBase< T, CT >::size(), and OneDimCubicSpline< Td, Tg, CTd, CTg >::spline().

Referenced by ECPComponentBuilder::addSemiLocal().

59 {
60  //todo rcut should be reset if necessary
61  using InFuncType = GaussianTimesRN<RealType>;
62  InFuncType a;
63  a.putBasisGroup(cur);
64  bool ignore = true;
65  const RealType eps = 1e-4;
66  mRealType rout = agrid->rmax() * 2;
67  while (ignore && rout > agrid->rmax())
68  {
69  ignore = (std::abs(a.f(rout)) < eps);
70  rout -= 0.01;
71  }
72  rout += 0.01;
73  app_log() << " cutoff for non-local pseudopotential = " << agrid->rmax() << std::endl;
74  app_log() << " calculated cutoff for " << eps << " = " << rout << std::endl;
75  int ng = agrid->size();
76  //all ways reset grid. Ye Luo
77  std::unique_ptr<GridType> agrid_new;
78  //if(agrid->GridTag != LINEAR_1DGRID)// || rout > agrid->rmax())
79  {
80  RealType ri = 0.0;
81  RealType rf = std::max(rout, agrid->rmax());
82  agrid_new = std::make_unique<LinearGrid<RealType>>();
83  agrid_new->set(ri, rf, ng);
84  app_log() << " Reset the grid for SemiLocal component to a LinearGrid. " << std::endl;
85  }
86  std::vector<RealType> v(ng);
87  for (int ig = 0; ig < ng; ig++)
88  {
89  v[ig] = a.f((*agrid_new)[ig]);
90  }
91  v[ng - 1] = 0.0;
92  RadialPotentialType* app = new RadialPotentialType(std::move(agrid_new), v);
93  app->spline();
94  return app;
95 }
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
std::ostream & app_log()
Definition: OutputManager.h:65
LocalECPotential::RadialPotentialType RadialPotentialType
QMCTraits::RealType RealType

◆ doBreakUp()

void doBreakUp ( const std::vector< int > &  angList,
const Matrix< mRealType > &  vnn,
RealType  rmax,
mRealType  Vprefactor = 1.0 
)

Separate local from non-local potentials.

Parameters
angListangular momentum list
vnnsemilocal tables of size (angList.size(),global_grid->size())
rmaxcutoff radius
Vprefactorconversion factor to Hartree

Note that local pseudopotential is r*V !!!

Definition at line 463 of file ECPComponentBuilder.2.cpp.

References qmcplusplus::app_log(), qmcplusplus::ceil(), Matrix< T, Alloc >::cols(), qmcplusplus::Units::charge::e, qmcplusplus::floor(), ECPComponentBuilder::grid_global, LINEAR_1DGRID, ECPComponentBuilder::Llocal, ECPComponentBuilder::Lmax, qmcplusplus::Units::distance::m, omptarget::min(), ECPComponentBuilder::NumNonLocal, ECPComponentBuilder::pp_loc, ECPComponentBuilder::pp_nonloc, OneDimCubicSpline< Td, Tg, CTd, CTg >::spline(), and ECPComponentBuilder::Zeff.

Referenced by ECPComponentBuilder::buildSemiLocalAndLocal(), and ECPComponentBuilder::parseCasino().

467 {
468  //ALERT magic number
469  const int max_points = 100000;
470  app_log() << " Creating a Linear Grid Rmax=" << rmax << std::endl;
471  //this is a new grid
472  mRealType d = 1e-4;
473  auto agrid = std::make_unique<LinearGrid<RealType>>();
474  // If the global grid is already linear, do not interpolate the data
475  int ng;
476  if (grid_global->getGridTag() == LINEAR_1DGRID)
477  {
478  ng = (int)std::ceil(rmax * grid_global->DeltaInv) + 1;
479  if (ng <= max_points)
480  {
481  app_log() << " Using global grid with delta = " << grid_global->Delta << std::endl;
482  rmax = grid_global->Delta * (ng - 1);
483  agrid->set(0.0, rmax, ng);
484  }
485  else
486  agrid->set(0.0, rmax, max_points);
487  }
488  else
489  {
490  ng = std::min(max_points, static_cast<int>(rmax / d) + 1);
491  agrid->set(0, rmax, ng);
492  }
493  // This is critical!!!
494  // If d is not reset, we generate an error in the interpolated PP!
495  d = agrid->Delta;
496  int ngIn = vnn.cols() - 2;
497 
498  assert(angList.size() > 0 && "ECPComponentBuilder::doBreakUp angList cannot be empty!");
499  //find the index of local
500  int iLlocal = -1;
501  for (int l = 0; l < angList.size(); l++)
502  if (angList[l] == Llocal)
503  iLlocal = l;
504  std::vector<RealType> newP(ng);
505  std::vector<mRealType> newPin(ngIn);
506  for (int l = 0; l < angList.size(); l++)
507  {
508  if (angList[l] == Llocal)
509  continue;
510  const mRealType* restrict vp = vnn[angList[l]];
511  const mRealType* restrict vpLoc = vnn[iLlocal];
512  int ll = angList[l];
513  for (int i = 0; i < ngIn; i++)
514  newPin[i] = Vprefactor * (vp[i] - vpLoc[i]);
515  OneDimCubicSpline<mRealType> infunc(grid_global->makeClone(), newPin);
516  infunc.spline(0, 0.0, ngIn - 1, 0.0);
517  for (int i = 1; i < ng - 1; i++)
518  {
519  mRealType r = d * i;
520  newP[i] = infunc.splint(r) / r;
521  }
522  newP[0] = newP[1];
523  newP[ng - 1] = 0.0;
524  auto app = new RadialPotentialType(agrid->makeClone(), newP);
525  app->spline();
526  pp_nonloc->add(angList[l], app);
527  }
528  NumNonLocal = Lmax;
529  if (Llocal == Lmax)
530  Lmax--;
531  if (NumNonLocal)
532  {
533  pp_nonloc->lmax = Lmax;
534  pp_nonloc->Rmax = rmax;
535  }
536  else
537  {
538  //only one component, remove non-local potentials
539  pp_nonloc.reset();
540  }
541  {
542  // Spline local potential on original grid
543  newPin[0] = 0.0;
544  mRealType vfac = -Vprefactor / Zeff;
545  const mRealType* restrict vpLoc = vnn[iLlocal];
546  for (int i = 0; i < ngIn; i++)
547  newPin[i] = vfac * vpLoc[i];
548  double dy0 = (newPin[1] - newPin[0]) / ((*grid_global)[1] - (*grid_global)[0]);
549  OneDimCubicSpline<mRealType> infunc(grid_global->makeClone(), newPin);
550  infunc.spline(0, dy0, ngIn - 1, 0.0);
551  int m = grid_global->size();
552  double loc_max = grid_global->r(m - 1);
553  int nloc = (int)std::floor(loc_max / d);
554  loc_max = (nloc - 1) * d;
555  auto grid_loc = std::make_unique<LinearGrid<RealType>>();
556  grid_loc->set(0.0, loc_max, nloc);
557  app_log() << " Making L=" << Llocal << " a local potential with a radial cutoff of " << loc_max << std::endl;
558  std::vector<RealType> newPloc(nloc);
559  for (int i = 1; i < nloc - 1; i++)
560  {
561  mRealType r = d * i;
562  newPloc[i] = infunc.splint(r);
563  }
564  newPloc[0] = 0.0;
565  newPloc[nloc - 1] = 1.0;
566  pp_loc = std::make_unique<RadialPotentialType>(std::move(grid_loc), newPloc);
567  pp_loc->spline(0, dy0, nloc - 1, 0.0);
568  // for (double r=0.0; r<3.50001; r+=0.001)
569  // fprintf (stderr, "%10.5f %10.5f\n", r, pp_loc->splint(r));
570  }
571 }
std::ostream & app_log()
Definition: OutputManager.h:65
void spline(int imin, value_type yp1, int imax, value_type ypn) override
Evaluate the 2nd derivative on the grid points.
T min(T a, T b)
std::unique_ptr< NonLocalECPComponent > pp_nonloc
size_type cols() const
Definition: OhmmsMatrix.h:78
std::unique_ptr< mGridType > grid_global
MakeReturn< UnaryNode< FnCeil, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t ceil(const Vector< T1, C1 > &l)
std::unique_ptr< RadialPotentialType > pp_loc
LocalECPotential::RadialPotentialType RadialPotentialType
MakeReturn< UnaryNode< FnFloor, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t floor(const Vector< T1, C1 > &l)

◆ parse()

bool parse ( const std::string &  fname,
xmlNodePtr  cur 
)

Definition at line 63 of file ECPComponentBuilder.cpp.

References getXMLAttributeValue(), ECPComponentBuilder::RcutMax, and ECPComponentBuilder::read_pp_file().

Referenced by ECPotentialBuilder::useXmlFormat().

64 {
65  const std::string cutoff_str(getXMLAttributeValue(cur, "cutoff"));
66  if (!cutoff_str.empty())
67  RcutMax = std::stod(cutoff_str);
68 
69  return read_pp_file(fname);
70 }
std::string getXMLAttributeValue(const xmlNodePtr cur, const std::string_view name)
get the value string for attribute name if name is unfound in cur you get an empty string back this i...
bool read_pp_file(const std::string &fname)

◆ parseCasino()

bool parseCasino ( const std::string &  fname,
xmlNodePtr  cur 
)

Definition at line 371 of file ECPComponentBuilder.2.cpp.

References qmcplusplus::abs(), OhmmsAttributeSet::add(), qmcplusplus::app_error(), qmcplusplus::app_log(), ECPComponentBuilder::AtomicNumber, Communicate::barrier_and_abort(), ECPComponentBuilder::doBreakUp(), ECPComponentBuilder::grid_global, ECPComponentBuilder::Llocal, ECPComponentBuilder::Lmax, MPIObjectBase::myComm, ECPComponentBuilder::Nrule, ECPComponentBuilder::pp_nonloc, OhmmsAttributeSet::put(), ECPComponentBuilder::SetQuadratureRule(), and ECPComponentBuilder::Zeff.

Referenced by ECPotentialBuilder::useXmlFormat().

372 {
373  app_log() << " Start ECPComponentBuilder::parseCasino" << std::endl;
374  RealType rmax = 2.0;
375  Llocal = -1;
376  Lmax = -1;
377  OhmmsAttributeSet aAttrib;
378  aAttrib.add(rmax, "cutoff");
379  aAttrib.add(Llocal, "l-local");
380  aAttrib.add(Lmax, "lmax");
381  aAttrib.add(Nrule, "nrule");
382  aAttrib.put(cur);
383 
384  std::ifstream fin(fname.c_str(), std::ios_base::in);
385  if (!fin)
386  {
387  app_error() << "Could not open file " << fname << std::endl;
388  myComm->barrier_and_abort("ECPComponentBuilder::parseCasino");
389  }
390  if (!pp_nonloc)
391  pp_nonloc = std::make_unique<NonLocalECPComponent>();
392  OhmmsAsciiParser aParser;
393  int npts = 0, idummy;
394  std::string eunits("rydberg");
395  app_log() << " ECPComponentBuilder::parseCasino" << std::endl;
396  aParser.skiplines(fin, 1); //Header
397  aParser.skiplines(fin, 1); //Atomic number and pseudo-charge
398  aParser.getValue(fin, AtomicNumber, Zeff);
399  app_log() << " Atomic number = " << AtomicNumber << " Zeff = " << Zeff << std::endl;
400  aParser.skiplines(fin, 1); //Energy units (rydberg/hartree/ev):
401  aParser.getValue(fin, eunits);
402  app_log() << " Unit of the potentials = " << eunits << std::endl;
403  mRealType Vprefactor = (eunits == "rydberg") ? 0.5 : 1.0;
404  aParser.skiplines(fin, 1); //Angular momentum of local component (0=s,1=p,2=d..)
405  aParser.getValue(fin, idummy);
406  if (Lmax < 0)
407  Lmax = idummy;
408  aParser.skiplines(fin, 1); //NLRULE override (1) VMC/DMC (2) config gen (0 ==> input/default value)
409  aParser.skiplines(fin, 1); //0 0, not sure what to do yet
410  aParser.skiplines(fin, 1); //Number of grid points
411  aParser.getValue(fin, npts);
412  app_log() << " Input Grid size = " << npts << std::endl;
413  std::vector<mRealType> temp(npts);
414  aParser.skiplines(fin, 1); //R(i) in atomic units
415  aParser.getValues(fin, temp.begin(), temp.end());
416  //create a global grid of numerical type
417  grid_global = std::make_unique<NumericalGrid<mRealType>>(temp);
418  Matrix<mRealType> vnn(Lmax + 1, npts);
419  for (int l = 0; l <= Lmax; l++)
420  {
421  aParser.skiplines(fin, 1);
422  aParser.getValues(fin, vnn[l], vnn[l] + npts);
423  }
424  std::vector<int> angList(Lmax + 1);
425  for (int l = 0; l <= Lmax; l++)
426  angList[l] = l;
427  // Now, check to see what maximum cutoff should be
428  if (vnn.size() > 1)
429  {
430  const double tolerance = 1.0e-5;
431  double rc_check = grid_global->r(npts - 1);
432  for (int j = npts - 1; j > 0; j--)
433  {
434  bool closeEnough = true;
435  for (int i = 0; i < vnn.rows(); i++)
436  for (int k = i + 1; k < vnn.rows(); k++)
437  if (std::abs(vnn[i][j] - vnn[k][j]) > tolerance)
438  closeEnough = false;
439  if (!closeEnough)
440  {
441  rc_check = grid_global->r(j);
442  break;
443  }
444  }
445  app_log() << " Maxium cutoff for non-local pseudopotentials " << rc_check << std::endl;
446  }
447  doBreakUp(angList, vnn, rmax, Vprefactor);
449  app_log() << " Non-local pseudopotential parameters" << std::endl;
450  pp_nonloc->print(app_log());
451  return true;
452 }
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
std::ostream & app_log()
Definition: OutputManager.h:65
bool put(xmlNodePtr cur)
assign attributes to the set
Definition: AttributeSet.h:55
std::ostream & app_error()
Definition: OutputManager.h:67
std::unique_ptr< NonLocalECPComponent > pp_nonloc
std::unique_ptr< mGridType > grid_global
Communicate * myComm
pointer to Communicate
Definition: MPIObjectBase.h:62
class to handle a set of attributes of an xmlNode
Definition: AttributeSet.h:24
QMCTraits::RealType RealType
void doBreakUp(const std::vector< int > &angList, const Matrix< mRealType > &vnn, RealType rmax, mRealType Vprefactor=1.0)
Separate local from non-local potentials.
void barrier_and_abort(const std::string &msg) const
void add(PDT &aparam, const std::string &aname, std::vector< PDT > candidate_values={}, TagStatus status=TagStatus::OPTIONAL)
add a new attribute
Definition: AttributeSet.h:42

◆ printECPTable()

void printECPTable ( )

Definition at line 223 of file ECPComponentBuilder.cpp.

References DEBUG, QMCState::io_node, OutputManagerClass::isActive(), ECPComponentBuilder::Lmax, QMCState::mpi_groups, outputManager, ECPComponentBuilder::pp_loc, ECPComponentBuilder::pp_nonloc, qmcplusplus::qmc_common, and ECPComponentBuilder::Species.

Referenced by ECPotentialBuilder::useXmlFormat().

224 {
226  return;
228  return;
229  std::array<char, 12> fname;
230  std::snprintf(fname.data(), fname.size(), "%s.pp.dat", Species.c_str());
231  std::ofstream fout(fname.data());
232  fout.setf(std::ios::scientific, std::ios::floatfield);
233  fout.precision(12);
234  int nl = pp_nonloc ? pp_nonloc->nlpp_m.size() : 0;
235  RealType d = 1.7e-2;
236  RealType rt = 0.13 * d;
237  if (nl)
238  {
239  fout << "# Lmax = " << Lmax + 1 << " nonlocal L channels" << nl << std::endl;
240  fout << "# Units = bohr hartree " << std::endl;
241  fout << "# r -r*V/zeff Vnl ... " << std::endl;
242  while (rt < 5)
243  {
244  fout << rt << std::setw(25) << pp_loc->splint(rt);
245  for (int l = 0; l < nl; l++)
246  fout << std::setw(25) << pp_nonloc->nlpp_m[l]->splint(rt);
247  fout << std::endl;
248  rt += d;
249  }
250  }
251  else
252  {
253  fout << "# Units = bohr hartree " << std::endl;
254  fout << "# r -r*V/zeff " << std::endl;
255  while (rt < 5)
256  {
257  fout << rt << std::setw(25) << pp_loc->splint(rt) << std::endl;
258  ;
259  rt += d;
260  }
261  }
262 }
std::unique_ptr< NonLocalECPComponent > pp_nonloc
OutputManagerClass outputManager(Verbosity::HIGH)
bool io_node
true, print out file
Definition: qmc_common.h:39
std::unique_ptr< RadialPotentialType > pp_loc
int mpi_groups
number of mpi groups
Definition: qmc_common.h:33
QMCTraits::RealType RealType
QMCState qmc_common
a unique QMCState during a run
Definition: qmc_common.cpp:111
bool isActive(Verbosity level) const

◆ put()

bool put ( xmlNodePtr  cur)

Definition at line 166 of file ECPComponentBuilder.cpp.

References ECPComponentBuilder::addSemiLocal(), qmcplusplus::app_log(), ECPComponentBuilder::AtomicNumber, ECPComponentBuilder::buildL2(), ECPComponentBuilder::buildLocal(), ECPComponentBuilder::buildSemiLocalAndLocal(), ECPComponentBuilder::createGrid(), getXMLAttributeValue(), ECPComponentBuilder::grid_global, ECPComponentBuilder::Nrule, ECPComponentBuilder::pp_loc, ECPComponentBuilder::pp_nonloc, ECPComponentBuilder::pp_so, ECPComponentBuilder::SetQuadratureRule(), and ECPComponentBuilder::Zeff.

Referenced by ECPComponentBuilder::read_pp_file(), and ECPotentialBuilder::useXmlFormat().

167 {
168  // int nk=0;
169  //vector<RealType> kpts;
170  std::vector<xmlNodePtr> semiPtr;
171  cur = cur->children;
172  while (cur != NULL)
173  {
174  std::string cname((const char*)cur->name);
175  if (cname == "header")
176  {
177  Zeff = std::stoi(getXMLAttributeValue(cur, "zval"));
178  AtomicNumber = std::stoi(getXMLAttributeValue(cur, "atomic-number"));
179  }
180  else if (cname == "grid")
181  {
182  //capture the global grid
183  grid_global = createGrid(cur);
184  }
185  else if (cname == "L2")
186  {
187  buildL2(cur);
188  }
189  else if (cname == "semilocal")
190  {
191  semiPtr.push_back(cur); //save the pointer
192  }
193  else if (cname == "local")
194  {
195  buildLocal(cur);
196  }
197  cur = cur->next;
198  }
199  if (semiPtr.size())
200  {
201  if (!pp_nonloc)
202  pp_nonloc = std::make_unique<NonLocalECPComponent>();
203  if (pp_so == 0)
204  pp_so = std::make_unique<SOECPComponent>();
205  if (pp_loc)
206  {
207  for (int i = 0; i < semiPtr.size(); i++)
208  addSemiLocal(semiPtr[i]);
209  }
210  else
211  buildSemiLocalAndLocal(semiPtr);
212  }
213  if (pp_nonloc)
214  {
216  app_log() << " Non-local pseudopotential parameters" << std::endl;
217  pp_nonloc->print(app_log());
218  app_log() << " Maximum cutoff radius " << pp_nonloc->Rmax << std::endl;
219  }
220  return true;
221 }
void buildLocal(xmlNodePtr cur)
build a Local Pseudopotential
std::ostream & app_log()
Definition: OutputManager.h:65
std::unique_ptr< NonLocalECPComponent > pp_nonloc
std::unique_ptr< mGridType > grid_global
std::unique_ptr< SOECPComponent > pp_so
std::unique_ptr< RadialPotentialType > pp_loc
std::string getXMLAttributeValue(const xmlNodePtr cur, const std::string_view name)
get the value string for attribute name if name is unfound in cur you get an empty string back this i...
void buildSemiLocalAndLocal(std::vector< xmlNodePtr > &semiPtr)
std::unique_ptr< mGridType > createGrid(xmlNodePtr cur, bool useLinear=false)

◆ read_pp_file()

bool read_pp_file ( const std::string &  fname)

Definition at line 136 of file ECPComponentBuilder.cpp.

References Communicate::barrier_and_abort(), ReadFileBuffer::contents(), ReadFileBuffer::length, MPIObjectBase::myComm, qmcplusplus::okay, ReadFileBuffer::open_file(), ECPComponentBuilder::put(), and ReadFileBuffer::read_contents().

Referenced by qmcplusplus::doSOECPotentialTest(), ECPComponentBuilder::parse(), and qmcplusplus::TEST_CASE().

137 {
138  ReadFileBuffer buf(myComm);
139  bool okay = buf.open_file(fname);
140  if (!okay)
141  myComm->barrier_and_abort("ECPComponentBuilder::read_pp_file Missing PP file " + fname + "\n");
142 
143  okay = buf.read_contents();
144  if (!okay)
145  myComm->barrier_and_abort("ECPComponentBuilder::read_pp_file Unable to read PP file " + fname + "\n");
146 
147  xmlDocPtr m_doc = xmlReadMemory(buf.contents(), buf.length, NULL, NULL, 0);
148 
149  if (m_doc == NULL)
150  {
151  xmlFreeDoc(m_doc);
152  myComm->barrier_and_abort("ECPComponentBuilder::read_pp_file xml file " + fname + " is invalid");
153  }
154  // Check the document is of the right kind
155  xmlNodePtr cur = xmlDocGetRootElement(m_doc);
156  if (cur == NULL)
157  {
158  xmlFreeDoc(m_doc);
159  myComm->barrier_and_abort("Empty document");
160  }
161  bool success = put(cur);
162  xmlFreeDoc(m_doc);
163  return success;
164 }
Communicate * myComm
pointer to Communicate
Definition: MPIObjectBase.h:62
void barrier_and_abort(const std::string &msg) const

◆ SetQuadratureRule()

void SetQuadratureRule ( int  rule)

Definition at line 264 of file ECPComponentBuilder.cpp.

References qmcplusplus::app_log(), ECPComponentBuilder::Lmax, Quadrature3D< T >::nk, ECPComponentBuilder::NumNonLocal, ECPComponentBuilder::NumSO, ECPComponentBuilder::pp_nonloc, ECPComponentBuilder::pp_so, ECPComponentBuilder::Srule, Quadrature3D< T >::weight_m, and Quadrature3D< T >::xyz_m.

Referenced by ECPComponentBuilder::parseCasino(), and ECPComponentBuilder::put().

265 {
266  app_log() << " Quadrature Nrule: " << rule << std::endl;
267  Quadrature3D<RealType> myRule(rule);
268  pp_nonloc->sgridxyz_m = myRule.xyz_m;
269  pp_nonloc->sgridweight_m = myRule.weight_m;
270  // Allocate storage for wave function ratios
271  pp_nonloc->resize_warrays(myRule.nk, NumNonLocal, Lmax);
272  if (pp_so)
273  { //added here bc must have nonlocal terms to have SO contributions
274  pp_so->sgridxyz_m_ = myRule.xyz_m;
275  pp_so->sgridweight_m_ = myRule.weight_m;
276  pp_so->resize_warrays(myRule.nk, NumSO, Srule);
277  }
278 }
std::ostream & app_log()
Definition: OutputManager.h:65
std::unique_ptr< NonLocalECPComponent > pp_nonloc
std::unique_ptr< SOECPComponent > pp_so

Member Data Documentation

◆ angMon

◆ AtomicNumber

◆ grid_global

◆ grid_inp

std::map<std::string, std::unique_ptr<mGridType> > grid_inp

Definition at line 44 of file ECPComponentBuilder.h.

Referenced by ECPComponentBuilder::createGrid().

◆ Llocal

◆ Lmax

◆ LmaxSO

int LmaxSO

Definition at line 38 of file ECPComponentBuilder.h.

Referenced by ECPComponentBuilder::buildSemiLocalAndLocal().

◆ Nrule

◆ NumNonLocal

◆ NumSO

int NumSO

◆ pp_L2

std::unique_ptr<L2RadialPotential> pp_L2

◆ pp_loc

◆ pp_nonloc

◆ pp_so

◆ RcutMax

RealType RcutMax

Definition at line 41 of file ECPComponentBuilder.h.

Referenced by ECPComponentBuilder::parse().

◆ Species

std::string Species

Definition at line 42 of file ECPComponentBuilder.h.

Referenced by ECPComponentBuilder::printECPTable().

◆ Srule

◆ Zeff


The documentation for this struct was generated from the following files: