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