QMCPACK
AOBasisBuilder.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2020 QMCPACK developers.
6 //
7 // File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
9 // Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
10 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
11 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
12 //
13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
14 //////////////////////////////////////////////////////////////////////////////////////
15 
16 
17 #include "AOBasisBuilder.h"
18 #include "Utilities/ProgressReportEngine.h"
19 #include "OhmmsData/AttributeSet.h"
20 #include "RadialOrbitalSetBuilder.h"
21 #include "SoaAtomicBasisSet.h"
22 #include "MultiQuinticSpline1D.h"
23 #include "MultiFunctorAdapter.h"
24 #include "Numerics/SoaCartesianTensor.h"
25 #include "Numerics/SoaSphericalTensor.h"
26 
27 namespace qmcplusplus
28 {
29 template<typename COT>
30 AOBasisBuilder<COT>::AOBasisBuilder(const std::string& eName, Communicate* comm)
31  : MPIObjectBase(comm),
32  addsignforM(false),
33  expandlm(GAUSSIAN_EXPAND),
34  Morder("gaussian"),
35  sph("default"),
36  basisType("Numerical"),
37  elementType(eName),
38  Normalized("yes")
39 {
40  // mmorales: for "Cartesian Gaussian", m is an integer that maps
41  // the component to Gamess notation, see Numerics/CartesianTensor.h
42  nlms_id["n"] = q_n;
43  nlms_id["l"] = q_l;
44  nlms_id["m"] = q_m;
45  nlms_id["s"] = q_s;
46 }
47 
48 template<class COT>
49 bool AOBasisBuilder<COT>::put(xmlNodePtr cur)
50 {
51  ReportEngine PRE("AtomicBasisBuilder", "put(xmlNodePtr)");
52  //Register valid attributes attributes
53  OhmmsAttributeSet aAttrib;
54  aAttrib.add(basisType, "type");
55  aAttrib.add(sph, "angular");
56  aAttrib.add(addsignforM, "expM");
57  aAttrib.add(Morder, "expandYlm");
58  aAttrib.add(Normalized, "normalized");
59  aAttrib.put(cur);
60  PRE.echo(cur);
61  if (sph == "spherical")
62  addsignforM = 1; //include (-1)^m
63 
64  if (Morder == "gaussian")
65  expandlm = GAUSSIAN_EXPAND;
66  else if (Morder == "natural")
67  expandlm = NATURAL_EXPAND;
68  else if (Morder == "no")
69  expandlm = DONOT_EXPAND;
70  else if (Morder == "pyscf")
71  {
72  expandlm = MOD_NATURAL_EXPAND;
73  addsignforM = 1;
74  if (sph != "spherical")
75  {
76  myComm->barrier_and_abort(" Error: expandYlm='pyscf' only compatible with angular='spherical'. Aborting.\n");
77  }
78  }
79 
80  if (sph == "cartesian" || Morder == "Gamess")
81  {
82  expandlm = CARTESIAN_EXPAND;
83  addsignforM = 0;
84  }
85 
86  if (Morder == "Dirac")
87  {
88  expandlm = DIRAC_CARTESIAN_EXPAND;
89  addsignforM = 0;
90  if (sph != "cartesian")
91  myComm->barrier_and_abort(" Error: expandYlm='Dirac' only compatible with angular='cartesian'. Aborting\n");
92  }
93 
94  // Numerical basis is a special case
95  if (basisType == "Numerical")
96  myComm->barrier_and_abort("Purely numerical atomic orbitals are not supported any longer.");
97 
98  return true;
99 }
100 
101 template<class COT>
102 bool AOBasisBuilder<COT>::putH5(hdf_archive& hin)
103 {
104  ReportEngine PRE("AtomicBasisBuilder", "putH5(hin)");
105  std::string CenterID, basisName;
106 
107  if (myComm->rank() == 0)
108  {
109  hin.read(sph, "angular");
110  hin.read(CenterID, "elementType");
111  hin.read(Normalized, "normalized");
112  hin.read(Morder, "expandYlm");
113  hin.read(basisName, "name");
114  }
115 
116  myComm->bcast(sph);
117  myComm->bcast(Morder);
118  myComm->bcast(CenterID);
119  myComm->bcast(Normalized);
120  myComm->bcast(basisName);
121  myComm->bcast(basisType);
122  myComm->bcast(addsignforM);
123 
124  if (sph == "spherical")
125  addsignforM = 1; //include (-1)^m
126 
127  if (Morder == "gaussian")
128  expandlm = GAUSSIAN_EXPAND;
129  else if (Morder == "natural")
130  expandlm = NATURAL_EXPAND;
131  else if (Morder == "no")
132  expandlm = DONOT_EXPAND;
133  else if (Morder == "pyscf")
134  {
135  expandlm = MOD_NATURAL_EXPAND;
136  addsignforM = 1;
137  if (sph != "spherical")
138  {
139  myComm->barrier_and_abort(" Error: expandYlm='pyscf' only compatible with angular='spherical'. Aborting.\n");
140  }
141  }
142 
143  if (sph == "cartesian" || Morder == "Gamess")
144  {
145  expandlm = CARTESIAN_EXPAND;
146  addsignforM = 0;
147  }
148 
149  if (Morder == "Dirac")
150  {
151  expandlm = DIRAC_CARTESIAN_EXPAND;
152  addsignforM = 0;
153  if (sph != "cartesian")
154  myComm->barrier_and_abort(" Error: expandYlm='Dirac' only compatible with angular='cartesian'. Aborting\n");
155  }
156  app_log() << R"(<input node="atomicBasisSet" name=")" << basisName << "\" expandYlm=\"" << Morder << "\" angular=\""
157  << sph << "\" elementType=\"" << CenterID << "\" normalized=\"" << Normalized << "\" type=\"" << basisType
158  << "\" expM=\"" << addsignforM << "\" />" << std::endl;
159 
160  return true;
161 }
162 
163 
164 template<typename COT>
165 std::unique_ptr<COT> AOBasisBuilder<COT>::createAOSet(xmlNodePtr cur)
166 {
167  ReportEngine PRE("AtomicBasisBuilder", "createAOSet(xmlNodePtr)");
168  app_log() << " AO BasisSet for " << elementType << "\n";
169 
170  if (expandlm != CARTESIAN_EXPAND)
171  {
172  if (addsignforM)
173  app_log() << " Spherical Harmonics contain (-1)^m factor" << std::endl;
174  else
175  app_log() << " Spherical Harmonics DO NOT contain (-1)^m factor" << std::endl;
176  }
177 
178  switch (expandlm)
179  {
180  case (GAUSSIAN_EXPAND):
181  app_log() << " Angular momentum m expanded according to Gaussian" << std::endl;
182  break;
183  case (NATURAL_EXPAND):
184  app_log() << " Angular momentum m expanded as -l, ... ,l" << std::endl;
185  break;
186  case (MOD_NATURAL_EXPAND):
187  app_log() << " Angular momentum m expanded as -l, ... ,l, with the exception of L=1 (1,-1,0)" << std::endl;
188  break;
189  case (CARTESIAN_EXPAND):
190  app_log() << " Angular momentum expanded in cartesian functions x^lx y^ly z^lz according to Gamess" << std::endl;
191  break;
192  case (DIRAC_CARTESIAN_EXPAND):
193  app_log() << " Angular momentum expanded in cartesian functions in DIRAC ordering" << std::endl;
194  break;
195  default:
196  app_log() << " Angular momentum m is explicitly given." << std::endl;
197  }
198 
199  QuantumNumberType nlms;
200  std::string rnl;
201  int Lmax(0); //maxmimum angular momentum of this center
202  int num(0); //the number of localized basis functions of this center
203  //process the basic property: maximun angular momentum, the number of basis functions to be added
204  std::vector<xmlNodePtr> radGroup;
205  xmlNodePtr cur1 = cur->xmlChildrenNode;
206  xmlNodePtr gptr = 0;
207  while (cur1 != NULL)
208  {
209  std::string cname1((const char*)(cur1->name));
210  if (cname1 == "basisGroup")
211  {
212  radGroup.push_back(cur1);
213  const int l = std::stoi(getXMLAttributeValue(cur1, "l"));
214  Lmax = std::max(Lmax, l);
215  //expect that only Rnl is given
216  if (expandlm == CARTESIAN_EXPAND || expandlm == DIRAC_CARTESIAN_EXPAND)
217  num += (l + 1) * (l + 2) / 2;
218  else if (expandlm)
219  num += 2 * l + 1;
220  else
221  num++;
222  }
223  else if (cname1 == "grid")
224  {
225  gptr = cur1;
226  }
227  cur1 = cur1->next;
228  }
229 
230  //create a new set of atomic orbitals sharing a center with (Lmax, num)
231  //if(addsignforM) the basis function has (-1)^m sqrt(2)Re(Ylm)
232  auto aos = std::make_unique<COT>(Lmax, addsignforM);
233  aos->LM.resize(num);
234  aos->NL.resize(num);
235 
236  //Now, add distinct Radial Orbitals and (l,m) channels
237  RadialOrbitalSetBuilder<COT> radFuncBuilder(myComm, *aos);
238  radFuncBuilder.Normalized = (Normalized == "yes");
239  radFuncBuilder.addGrid(gptr, basisType); //assign a radial grid for the new center
240  std::vector<xmlNodePtr>::iterator it(radGroup.begin());
241  std::vector<xmlNodePtr>::iterator it_end(radGroup.end());
242  std::vector<int> all_nl;
243  while (it != it_end)
244  {
245  cur1 = (*it);
246  xmlAttrPtr att = cur1->properties;
247  while (att != NULL)
248  {
249  std::string aname((const char*)(att->name));
250  if (aname == "rid" || aname == "id")
251  //accept id/rid
252  {
253  rnl = (const char*)(att->children->content);
254  }
255  else
256  {
257  std::map<std::string, int>::iterator iit = nlms_id.find(aname);
258  if (iit != nlms_id.end())
259  //valid for n,l,m,s
260  {
261  nlms[(*iit).second] = atoi((const char*)(att->children->content));
262  }
263  }
264  att = att->next;
265  }
266  //add Ylm channels
267  app_log() << " R(n,l,m,s) " << nlms[0] << " " << nlms[1] << " " << nlms[2] << " " << nlms[3] << std::endl;
268  std::map<std::string, int>::iterator rnl_it = RnlID.find(rnl);
269  if (rnl_it == RnlID.end())
270  {
271  int nl = aos->RnlID.size();
272  if (radFuncBuilder.addRadialOrbital(cur1, basisType, nlms))
273  RnlID[rnl] = nl;
274  all_nl.push_back(nl);
275  }
276  else
277  {
278  all_nl.push_back((*rnl_it).second);
279  }
280  ++it;
281  }
282 
283  if (expandYlm(aos.get(), all_nl, expandlm) != num)
284  myComm->barrier_and_abort("expandYlm doesn't match the number of basis.");
285  radFuncBuilder.finalize();
286  aos->finalize();
287  app_log() << " Maximum Angular Momentum = " << aos->Ylm.lmax() << std::endl
288  << " Number of Radial functors = " << aos->RnlID.size() << std::endl
289  << " Basis size = " << aos->getBasisSetSize() << "\n\n";
290  return aos;
291 }
292 
293 
294 template<typename COT>
295 std::unique_ptr<COT> AOBasisBuilder<COT>::createAOSetH5(hdf_archive& hin)
296 {
297  ReportEngine PRE("AOBasisBuilder:", "createAOSetH5(std::string)");
298  app_log() << " AO BasisSet for " << elementType << "\n";
299 
300  if (expandlm != CARTESIAN_EXPAND)
301  {
302  if (addsignforM)
303  app_log() << " Spherical Harmonics contain (-1)^m factor" << std::endl;
304  else
305  app_log() << " Spherical Harmonics DO NOT contain (-1)^m factor" << std::endl;
306  }
307 
308  switch (expandlm)
309  {
310  case (GAUSSIAN_EXPAND):
311  app_log() << " Angular momentum m expanded according to Gaussian" << std::endl;
312  break;
313  case (NATURAL_EXPAND):
314  app_log() << " Angular momentum m expanded as -l, ... ,l" << std::endl;
315  break;
316  case (MOD_NATURAL_EXPAND):
317  app_log() << " Angular momentum m expanded as -l, ... ,l, with the exception of L=1 (1,-1,0)" << std::endl;
318  break;
319  case (CARTESIAN_EXPAND):
320  app_log() << " Angular momentum expanded in cartesian functions x^lx y^ly z^lz according to Gamess" << std::endl;
321  break;
322  case (DIRAC_CARTESIAN_EXPAND):
323  app_log() << " Angular momentum expanded in cartesian functions in DIRAC ordering" << std::endl;
324  break;
325  default:
326  app_log() << " Angular momentum m is explicitly given." << std::endl;
327  }
328 
329  QuantumNumberType nlms;
330  std::string rnl;
331  int Lmax(0); //maxmimum angular momentum of this center
332  int num(0); //the number of localized basis functions of this center
333 
334  int numbasisgroups(0);
335  if (myComm->rank() == 0)
336  {
337  if (!hin.readEntry(numbasisgroups, "NbBasisGroups"))
338  PRE.error("Could not read NbBasisGroups in H5; Probably Corrupt H5 file", true);
339  }
340  myComm->bcast(numbasisgroups);
341 
342  for (int i = 0; i < numbasisgroups; i++)
343  {
344  std::string basisGroupID = "basisGroup" + std::to_string(i);
345  int l(0);
346  if (myComm->rank() == 0)
347  {
348  hin.push(basisGroupID);
349  hin.read(l, "l");
350  hin.pop();
351  }
352  myComm->bcast(l);
353 
354  Lmax = std::max(Lmax, l);
355  //expect that only Rnl is given
356  if (expandlm == CARTESIAN_EXPAND || expandlm == DIRAC_CARTESIAN_EXPAND)
357  num += (l + 1) * (l + 2) / 2;
358  else if (expandlm)
359  num += 2 * l + 1;
360  else
361  num++;
362  }
363 
364  //create a new set of atomic orbitals sharing a center with (Lmax, num)
365  //if(addsignforM) the basis function has (-1)^m sqrt(2)Re(Ylm)
366  auto aos = std::make_unique<COT>(Lmax, addsignforM);
367  aos->LM.resize(num);
368  aos->NL.resize(num);
369 
370  //Now, add distinct Radial Orbitals and (l,m) channels
371  RadialOrbitalSetBuilder<COT> radFuncBuilder(myComm, *aos);
372  radFuncBuilder.Normalized = (Normalized == "yes");
373  radFuncBuilder.addGridH5(hin); //assign a radial grid for the new center
374  std::vector<int> all_nl;
375  for (int i = 0; i < numbasisgroups; i++)
376  {
377  std::string basisGroupID = "basisGroup" + std::to_string(i);
378  if (myComm->rank() == 0)
379  {
380  hin.push(basisGroupID);
381  hin.read(rnl, "rid");
382  hin.read(nlms[0], "n");
383  hin.read(nlms[1], "l");
384  }
385  myComm->bcast(rnl);
386  myComm->bcast(nlms[0]);
387  myComm->bcast(nlms[1]);
388 
389  //add Ylm channels
390  app_log() << " R(n,l,m,s) " << nlms[0] << " " << nlms[1] << " " << nlms[2] << " " << nlms[3] << std::endl;
391  std::map<std::string, int>::iterator rnl_it = RnlID.find(rnl);
392  if (rnl_it == RnlID.end())
393  {
394  int nl = aos->RnlID.size();
395  if (radFuncBuilder.addRadialOrbitalH5(hin, basisType, nlms))
396  RnlID[rnl] = nl;
397  all_nl.push_back(nl);
398  }
399  else
400  {
401  all_nl.push_back((*rnl_it).second);
402  }
403 
404  if (myComm->rank() == 0)
405  hin.pop();
406  }
407 
408  if (expandYlm(aos.get(), all_nl, expandlm) != num)
409  myComm->barrier_and_abort("expandYlm doesn't match the number of basis.");
410  radFuncBuilder.finalize();
411  aos->finalize();
412  app_log() << " Maximum Angular Momentum = " << aos->Ylm.lmax() << std::endl
413  << " Number of Radial functors = " << aos->RnlID.size() << std::endl
414  << " Basis size = " << aos->getBasisSetSize() << "\n\n";
415  return aos;
416 }
417 
418 
419 template<typename COT>
420 int AOBasisBuilder<COT>::expandYlm(COT* aos, std::vector<int>& all_nl, int expandlm)
421 {
422  int num = 0;
423  if (expandlm == GAUSSIAN_EXPAND)
424  {
425  app_log() << "Expanding Ylm according to Gaussian98" << std::endl;
426  for (int nl = 0; nl < aos->RnlID.size(); nl++)
427  {
428  int l = aos->RnlID[nl][q_l];
429  app_log() << "Adding " << 2 * l + 1 << " spherical orbitals for l= " << l << std::endl;
430  switch (l)
431  {
432  case (0):
433  aos->LM[num] = aos->Ylm.index(0, 0);
434  aos->NL[num] = nl;
435  num++;
436  break;
437  case (1): //px(1),py(-1),pz(0)
438  aos->LM[num] = aos->Ylm.index(1, 1);
439  aos->NL[num] = nl;
440  num++;
441  aos->LM[num] = aos->Ylm.index(1, -1);
442  aos->NL[num] = nl;
443  num++;
444  aos->LM[num] = aos->Ylm.index(1, 0);
445  aos->NL[num] = nl;
446  num++;
447  break;
448  default: //0,1,-1,2,-2,...,l,-l
449  aos->LM[num] = aos->Ylm.index(l, 0);
450  aos->NL[num] = nl;
451  num++;
452  for (int tm = 1; tm <= l; tm++)
453  {
454  aos->LM[num] = aos->Ylm.index(l, tm);
455  aos->NL[num] = nl;
456  num++;
457  aos->LM[num] = aos->Ylm.index(l, -tm);
458  aos->NL[num] = nl;
459  num++;
460  }
461  break;
462  }
463  }
464  }
465  else if (expandlm == MOD_NATURAL_EXPAND)
466  {
467  app_log() << "Expanding Ylm as L=1 as (1,-1,0) and L>1 as -l,-l+1,...,l-1,l" << std::endl;
468  for (int nl = 0; nl < aos->RnlID.size(); nl++)
469  {
470  int l = aos->RnlID[nl][q_l];
471  app_log() << " Adding " << 2 * l + 1 << " spherical orbitals" << std::endl;
472  if (l == 1)
473  {
474  //px(1),py(-1),pz(0)
475  aos->LM[num] = aos->Ylm.index(1, 1);
476  aos->NL[num] = nl;
477  num++;
478  aos->LM[num] = aos->Ylm.index(1, -1);
479  aos->NL[num] = nl;
480  num++;
481  aos->LM[num] = aos->Ylm.index(1, 0);
482  aos->NL[num] = nl;
483  num++;
484  }
485  else
486  {
487  for (int tm = -l; tm <= l; tm++, num++)
488  {
489  aos->LM[num] = aos->Ylm.index(l, tm);
490  aos->NL[num] = nl;
491  }
492  }
493  }
494  }
495  else if (expandlm == NATURAL_EXPAND)
496  {
497  app_log() << "Expanding Ylm as -l,-l+1,...,l-1,l" << std::endl;
498  for (int nl = 0; nl < aos->RnlID.size(); nl++)
499  {
500  int l = aos->RnlID[nl][q_l];
501  app_log() << " Adding " << 2 * l + 1 << " spherical orbitals" << std::endl;
502  for (int tm = -l; tm <= l; tm++, num++)
503  {
504  aos->LM[num] = aos->Ylm.index(l, tm);
505  aos->NL[num] = nl;
506  }
507  }
508  }
509  else if (expandlm == CARTESIAN_EXPAND)
510  {
511  app_log() << "Expanding Ylm (angular function) according to Gamess using cartesian gaussians" << std::endl;
512  for (int nl = 0; nl < aos->RnlID.size(); nl++)
513  {
514  int l = aos->RnlID[nl][q_l];
515  app_log() << "Adding " << (l + 1) * (l + 2) / 2 << " cartesian gaussian orbitals for l= " << l << std::endl;
516  int nbefore = 0;
517  for (int i = 0; i < l; i++)
518  nbefore += (i + 1) * (i + 2) / 2;
519  for (int i = 0; i < (l + 1) * (l + 2) / 2; i++)
520  {
521  aos->LM[num] = nbefore + i;
522  aos->NL[num] = nl;
523  num++;
524  }
525  }
526  }
527  else if (expandlm == DIRAC_CARTESIAN_EXPAND)
528  {
529  app_log() << "Expanding Ylm (angular function) according to DIRAC using cartesian gaussians" << std::endl;
530  for (int nl = 0; nl < aos->RnlID.size(); nl++)
531  {
532  int l = aos->RnlID[nl][q_l];
533  app_log() << "Adding " << (l + 1) * (l + 2) / 2 << " cartesian gaussian orbitals for l= " << l << std::endl;
534  int nbefore = 0;
535  for (int i = 0; i < l; i++)
536  nbefore += (i + 1) * (i + 2) / 2;
537  switch (l)
538  {
539  case (0):
540  aos->LM[num] = nbefore + 0;
541  aos->NL[num] = nl;
542  num++;
543  break;
544  case (1):
545  aos->LM[num] = nbefore + 0;
546  aos->NL[num] = nl;
547  num++;
548  aos->LM[num] = nbefore + 1;
549  aos->NL[num] = nl;
550  num++;
551  aos->LM[num] = nbefore + 2;
552  aos->NL[num] = nl;
553  num++;
554  break;
555  case (2):
556  aos->LM[num] = nbefore + 0; //xx
557  aos->NL[num] = nl;
558  num++;
559  aos->LM[num] = nbefore + 3; //xy
560  aos->NL[num] = nl;
561  num++;
562  aos->LM[num] = nbefore + 4; //xz
563  aos->NL[num] = nl;
564  num++;
565  aos->LM[num] = nbefore + 1; //yy
566  aos->NL[num] = nl;
567  num++;
568  aos->LM[num] = nbefore + 5; //yz
569  aos->NL[num] = nl;
570  num++;
571  aos->LM[num] = nbefore + 2; //zz
572  aos->NL[num] = nl;
573  num++;
574  break;
575  case (3):
576  aos->LM[num] = nbefore + 0; //xxx
577  aos->NL[num] = nl;
578  num++;
579  aos->LM[num] = nbefore + 3; //xxy
580  aos->NL[num] = nl;
581  num++;
582  aos->LM[num] = nbefore + 4; //xxz
583  aos->NL[num] = nl;
584  num++;
585  aos->LM[num] = nbefore + 5; //xyy
586  aos->NL[num] = nl;
587  num++;
588  aos->LM[num] = nbefore + 9; //xyz
589  aos->NL[num] = nl;
590  num++;
591  aos->LM[num] = nbefore + 7; //xzz
592  aos->NL[num] = nl;
593  num++;
594  aos->LM[num] = nbefore + 1; //yyy
595  aos->NL[num] = nl;
596  num++;
597  aos->LM[num] = nbefore + 6; //yyz
598  aos->NL[num] = nl;
599  num++;
600  aos->LM[num] = nbefore + 8; //yzz
601  aos->NL[num] = nl;
602  num++;
603  aos->LM[num] = nbefore + 2; //zzz
604  aos->NL[num] = nl;
605  num++;
606  break;
607  case (4):
608  aos->LM[num] = nbefore + 0; //400
609  aos->NL[num] = nl;
610  num++;
611  aos->LM[num] = nbefore + 3; //310
612  aos->NL[num] = nl;
613  num++;
614  aos->LM[num] = nbefore + 4; //301
615  aos->NL[num] = nl;
616  num++;
617  aos->LM[num] = nbefore + 9; //220
618  aos->NL[num] = nl;
619  num++;
620  aos->LM[num] = nbefore + 12; //211
621  aos->NL[num] = nl;
622  num++;
623  aos->LM[num] = nbefore + 10; //202
624  aos->NL[num] = nl;
625  num++;
626  aos->LM[num] = nbefore + 5; //130
627  aos->NL[num] = nl;
628  num++;
629  aos->LM[num] = nbefore + 13; //121
630  aos->NL[num] = nl;
631  num++;
632  aos->LM[num] = nbefore + 14; //112
633  aos->NL[num] = nl;
634  num++;
635  aos->LM[num] = nbefore + 7; //103
636  aos->NL[num] = nl;
637  num++;
638  aos->LM[num] = nbefore + 1; //040
639  aos->NL[num] = nl;
640  num++;
641  aos->LM[num] = nbefore + 6; //031
642  aos->NL[num] = nl;
643  num++;
644  aos->LM[num] = nbefore + 11; //022
645  aos->NL[num] = nl;
646  num++;
647  aos->LM[num] = nbefore + 8; //013
648  aos->NL[num] = nl;
649  num++;
650  aos->LM[num] = nbefore + 2; //004
651  aos->NL[num] = nl;
652  num++;
653  break;
654  case (5):
655  aos->LM[num] = nbefore + 0; //500
656  aos->NL[num] = nl;
657  num++;
658  aos->LM[num] = nbefore + 3; //410
659  aos->NL[num] = nl;
660  num++;
661  aos->LM[num] = nbefore + 4; //401
662  aos->NL[num] = nl;
663  num++;
664  aos->LM[num] = nbefore + 9; //320
665  aos->NL[num] = nl;
666  num++;
667  aos->LM[num] = nbefore + 15; //311
668  aos->NL[num] = nl;
669  num++;
670  aos->LM[num] = nbefore + 10; //302
671  aos->NL[num] = nl;
672  num++;
673  aos->LM[num] = nbefore + 11; //230
674  aos->NL[num] = nl;
675  num++;
676  aos->LM[num] = nbefore + 18; //221
677  aos->NL[num] = nl;
678  num++;
679  aos->LM[num] = nbefore + 19; //212
680  aos->NL[num] = nl;
681  num++;
682  aos->LM[num] = nbefore + 13; //203
683  aos->NL[num] = nl;
684  num++;
685  aos->LM[num] = nbefore + 5; //140
686  aos->NL[num] = nl;
687  num++;
688  aos->LM[num] = nbefore + 16; //131
689  aos->NL[num] = nl;
690  num++;
691  aos->LM[num] = nbefore + 20; //122
692  aos->NL[num] = nl;
693  num++;
694  aos->LM[num] = nbefore + 17; //113
695  aos->NL[num] = nl;
696  num++;
697  aos->LM[num] = nbefore + 7; //104
698  aos->NL[num] = nl;
699  num++;
700  aos->LM[num] = nbefore + 1; //050
701  aos->NL[num] = nl;
702  num++;
703  aos->LM[num] = nbefore + 6; //041
704  aos->NL[num] = nl;
705  num++;
706  aos->LM[num] = nbefore + 12; //032
707  aos->NL[num] = nl;
708  num++;
709  aos->LM[num] = nbefore + 14; //023
710  aos->NL[num] = nl;
711  num++;
712  aos->LM[num] = nbefore + 8; //014
713  aos->NL[num] = nl;
714  num++;
715  aos->LM[num] = nbefore + 2; //005
716  aos->NL[num] = nl;
717  num++;
718  break;
719  case (6):
720  aos->LM[num] = nbefore + 0; //600
721  aos->NL[num] = nl;
722  num++;
723  aos->LM[num] = nbefore + 3; //510
724  aos->NL[num] = nl;
725  num++;
726  aos->LM[num] = nbefore + 4; //501
727  aos->NL[num] = nl;
728  num++;
729  aos->LM[num] = nbefore + 9; //420
730  aos->NL[num] = nl;
731  num++;
732  aos->LM[num] = nbefore + 15; //411
733  aos->NL[num] = nl;
734  num++;
735  aos->LM[num] = nbefore + 10; //402
736  aos->NL[num] = nl;
737  num++;
738  aos->LM[num] = nbefore + 18; //330
739  aos->NL[num] = nl;
740  num++;
741  aos->LM[num] = nbefore + 21; //321
742  aos->NL[num] = nl;
743  num++;
744  aos->LM[num] = nbefore + 22; //312
745  aos->NL[num] = nl;
746  num++;
747  aos->LM[num] = nbefore + 19; //303
748  aos->NL[num] = nl;
749  num++;
750  aos->LM[num] = nbefore + 11; //240
751  aos->NL[num] = nl;
752  num++;
753  aos->LM[num] = nbefore + 23; //231
754  aos->NL[num] = nl;
755  num++;
756  aos->LM[num] = nbefore + 27; //222
757  aos->NL[num] = nl;
758  num++;
759  aos->LM[num] = nbefore + 25; //213
760  aos->NL[num] = nl;
761  num++;
762  aos->LM[num] = nbefore + 13; //204
763  aos->NL[num] = nl;
764  num++;
765  aos->LM[num] = nbefore + 5; //150
766  aos->NL[num] = nl;
767  num++;
768  aos->LM[num] = nbefore + 16; //141
769  aos->NL[num] = nl;
770  num++;
771  aos->LM[num] = nbefore + 24; //132
772  aos->NL[num] = nl;
773  num++;
774  aos->LM[num] = nbefore + 26; //123
775  aos->NL[num] = nl;
776  num++;
777  aos->LM[num] = nbefore + 17; //114
778  aos->NL[num] = nl;
779  num++;
780  aos->LM[num] = nbefore + 7; //105
781  aos->NL[num] = nl;
782  num++;
783  aos->LM[num] = nbefore + 1; //060
784  aos->NL[num] = nl;
785  num++;
786  aos->LM[num] = nbefore + 6; //051
787  aos->NL[num] = nl;
788  num++;
789  aos->LM[num] = nbefore + 12; //042
790  aos->NL[num] = nl;
791  num++;
792  aos->LM[num] = nbefore + 20; //033
793  aos->NL[num] = nl;
794  num++;
795  aos->LM[num] = nbefore + 14; //024
796  aos->NL[num] = nl;
797  num++;
798  aos->LM[num] = nbefore + 8; //015
799  aos->NL[num] = nl;
800  num++;
801  aos->LM[num] = nbefore + 2; //006
802  aos->NL[num] = nl;
803  num++;
804  break;
805  default:
806  myComm->barrier_and_abort("Cartesian Tensor only defined up to Lmax=6. Aborting\n");
807  break;
808  }
809  }
810  }
811  else
812  {
813  for (int ind = 0; ind < all_nl.size(); ind++)
814  {
815  int nl = all_nl[ind];
816  int l = aos->RnlID[nl][q_l];
817  int m = aos->RnlID[nl][q_m];
818  //assign the index for real Spherical Harmonic with (l,m)
819  aos->LM[num] = aos->Ylm.index(l, m);
820  //assign the index for radial orbital with (n,l)
821  aos->NL[num] = nl;
822  //increment number of basis functions
823  num++;
824  }
825  }
826  return num;
827 }
828 
829 template class AOBasisBuilder<
830  SoaAtomicBasisSet<MultiQuinticSpline1D<QMCTraits::RealType>, SoaCartesianTensor<QMCTraits::RealType>>>;
831 template class AOBasisBuilder<
832  SoaAtomicBasisSet<MultiQuinticSpline1D<QMCTraits::RealType>, SoaSphericalTensor<QMCTraits::RealType>>>;
833 template class AOBasisBuilder<SoaAtomicBasisSet<MultiFunctorAdapter<GaussianCombo<QMCTraits::RealType>>,
834  SoaCartesianTensor<QMCTraits::RealType>>>;
835 template class AOBasisBuilder<SoaAtomicBasisSet<MultiFunctorAdapter<GaussianCombo<QMCTraits::RealType>>,
836  SoaSphericalTensor<QMCTraits::RealType>>>;
837 template class AOBasisBuilder<
838  SoaAtomicBasisSet<MultiFunctorAdapter<SlaterCombo<QMCTraits::RealType>>, SoaCartesianTensor<QMCTraits::RealType>>>;
839 template class AOBasisBuilder<
840  SoaAtomicBasisSet<MultiFunctorAdapter<SlaterCombo<QMCTraits::RealType>>, SoaSphericalTensor<QMCTraits::RealType>>>;
841 
842 } // namespace qmcplusplus
qmcplusplus::ReportEngine::echo
void echo(xmlNodePtr cur, bool recursive=false)
Definition: ProgressReportEngine.cpp:20
qmcplusplus::MPIObjectBase
Base class for any object which needs to know about a MPI communicator.
Definition: MPIObjectBase.h:26
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
qmcplusplus::app_log
std::ostream & app_log()
Definition: OutputManager.h:65
OhmmsAttributeSet::put
bool put(xmlNodePtr cur)
assign attributes to the set
Definition: AttributeSet.h:55
qmcplusplus::RadialOrbitalSetBuilder
Build a set of radial orbitals at the origin.
Definition: MultiFunctorAdapter.h:187
qmcplusplus::hdf_archive
class to handle hdf file
Definition: hdf_archive.h:51
qmcplusplus::SoaSphericalTensor
SoaSphericalTensor that evaluates the Real Spherical Harmonics.
Definition: SoaSphericalTensor.h:40
qmcplusplus::RadialOrbitalSetBuilder::addGrid
bool addGrid(xmlNodePtr cur, const std::string &rad_type)
add a grid
Definition: RadialOrbitalSetBuilder.h:160
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
qmcplusplus::AOBasisBuilder
atomic basisset builder
Definition: AOBasisBuilder.h:33
qmcplusplus::ReportEngine
Final class and should not be derived.
Definition: ProgressReportEngine.h:33
qmcplusplus::RadialOrbitalSetBuilder::finalize
void finalize()
This is when the radial orbitals are actually created.
Definition: RadialOrbitalSetBuilder.h:328
OhmmsAttributeSet
class to handle a set of attributes of an xmlNode
Definition: AttributeSet.h:24
ProgressReportEngine.h
declaration of ProgressReportEngine
qmcplusplus::RadialOrbitalSetBuilder::addRadialOrbital
bool addRadialOrbital(xmlNodePtr cur, const std::string &m_infunctype, const QuantumNumberType &nlms)
add a radial functor
Definition: RadialOrbitalSetBuilder.h:237
getXMLAttributeValue
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...
Definition: XMLParsingString.cpp:19
qmcplusplus::AOBasisBuilder::AOBasisBuilder
AOBasisBuilder(const std::string &eName, Communicate *comm)
Definition: AOBasisBuilder.cpp:30
qmcplusplus::hdf_archive::push
void push(const std::string &gname, bool createit=true)
push a group to the group stack
Definition: hdf_archive.cpp:198
qmcplusplus::AOBasisBuilder::createAOSet
std::unique_ptr< COT > createAOSet(xmlNodePtr cur)
Definition: AOBasisBuilder.cpp:165
qmcplusplus::ReportEngine::error
void error(const std::string &msg, bool fatal=false)
Definition: ProgressReportEngine.h:62
qmcplusplus::AOBasisBuilder::nlms_id
std::map< std::string, int > nlms_id
map for (n,l,m,s) to its quantum number index
Definition: AOBasisBuilder.h:59
qmcplusplus::hdf_archive::read
void read(T &data, const std::string &aname)
read the data from the group aname and check status runtime error is issued on I/O error ...
Definition: hdf_archive.h:306
qmcplusplus::AOBasisBuilder::createAOSetH5
std::unique_ptr< COT > createAOSetH5(hdf_archive &hin)
Definition: AOBasisBuilder.cpp:295
qmcplusplus::hdf_archive::readEntry
bool readEntry(T &data, const std::string &aname)
read the data from the group aname and return status use read() for inbuilt error checking ...
Definition: hdf_archive.h:293
qmcplusplus::AOBasisBuilder::expandYlm
int expandYlm(COT *aos, std::vector< int > &all_nl, int expandlm=DONOT_EXPAND)
Definition: AOBasisBuilder.cpp:420
qmcplusplus::SoaCartesianTensor
CartesianTensor according to Gamess order.
Definition: SoaCartesianTensor.h:40
OhmmsAttributeSet::add
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
qmcplusplus::RadialOrbitalSetBuilder::Normalized
bool Normalized
true, if the RadialOrbitalType is normalized
Definition: RadialOrbitalSetBuilder.h:99
qmcplusplus::AOBasisBuilder::put
bool put(xmlNodePtr cur)
Definition: AOBasisBuilder.cpp:49
qmcplusplus::RadialOrbitalSetBuilder::addRadialOrbitalH5
bool addRadialOrbitalH5(hdf_archive &hin, const std::string &radtype, const QuantumNumberType &nlms)
Definition: RadialOrbitalSetBuilder.h:259
qmcplusplus::AOBasisBuilder::putH5
bool putH5(hdf_archive &hin)
Definition: AOBasisBuilder.cpp:102
MultiQuinticSpline1D.h
Assume that coeffs.D1 and the LogLightGrid r_values.size() are equal Therefore r must be < r_max...