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< mGridType > createGrid (xmlNodePtr cur, bool useLinear=false)
 
RadialPotentialType * createVrWithBasisGroup (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...
 
Communicate * getCommunicator () const
 return myComm More...
 
Communicate & getCommRef () 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< mGridType > grid_global
 
std::map< std::string, std::unique_ptr< mGridType > > grid_inp
 
std::unique_ptr< RadialPotentialType > pp_loc
 
std::unique_ptr< NonLocalECPComponent > pp_nonloc
 
std::unique_ptr< SOECPComponent > pp_so
 
std::unique_ptr< L2RadialPotential > pp_L2
 
std::map< std::string, int > angMon
 

Additional Inherited Members

- Protected Attributes inherited from MPIObjectBase
Communicate * myComm
 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: