Merge pull request #4032 from ye-luo/input-llocal

Allow overriding NLPP local channel via input.

docs/hamiltonianobservable.rst

@@ -433,6 +433,8 @@ attributes:
   +-----------------------------------+--------------+-----------------+-------------+---------------------------+
   | ``nrule``:math:`^o`               | integer      |                 |             | Integration grid order    |
   +-----------------------------------+--------------+-----------------+-------------+---------------------------+
+  | ``l-local``:math:`^o`             | integer      |                 |             | Override local channel    |
+  +-----------------------------------+--------------+-----------------+-------------+---------------------------+
 
 .. code-block::
   :caption: QMCPXML element for pseudopotential of single ionic species.

src/QMCHamiltonians/ECPComponentBuilder.2.cpp

@@ -54,33 +54,34 @@ void ECPComponentBuilder::buildSemiLocalAndLocal(std::vector<xmlNodePtr>& semiPt
   int ndown = 1;
   int nup   = 0;
   int nso   = 0;
-  Llocal    = -1;
   OhmmsAttributeSet aAttrib;
-  int quadRule = -1;
+  int quad_rule = -1;
+  int local_channel = -1;
   aAttrib.add(eunits, "units");
   aAttrib.add(format, "format");
   aAttrib.add(ndown, "npots-down");
   aAttrib.add(nup, "npots-up");
-  aAttrib.add(Llocal, "l-local");
-  aAttrib.add(quadRule, "nrule");
+  aAttrib.add(local_channel, "l-local");
+  aAttrib.add(quad_rule, "nrule");
   aAttrib.add(Srule, "srule");
   aAttrib.add(nso, "npots-so");
 
   xmlNodePtr cur_semilocal = semiPtr[0];
   aAttrib.put(cur_semilocal);
 
-  if (quadRule > -1 && Nrule > -1)
+  // settle Nrule. Priority: current value (from input file) > PP XML file > lmax derived
+  if (quad_rule > -1 && Nrule > -1)
   {
     app_warning() << " Nrule setting found in both qmcpack input (Nrule = " << Nrule
-                  << ") and pseudopotential file (Nrule = " << quadRule << ")."
+                  << ") and pseudopotential file (Nrule = " << quad_rule << ")."
                   << " Using nrule setting in qmcpack input file." << std::endl;
   }
-  else if (quadRule > -1 && Nrule == -1)
+  else if (quad_rule > -1 && Nrule == -1)
   {
     app_log() << " Nrule setting found in pseudopotential file and used." << std::endl;
-    Nrule = quadRule;
+    Nrule = quad_rule;
   }
-  else if (quadRule == -1 && Nrule > -1)
+  else if (quad_rule == -1 && Nrule > -1)
     app_log() << " Nrule setting found in qmcpack input file and used." << std::endl;
   else
   {
@@ -197,13 +198,35 @@ void ECPComponentBuilder::buildSemiLocalAndLocal(std::vector<xmlNodePtr>& semiPt
     }
     cur_vps = cur_vps->next;
   }
+
   if (rmax < 0)
     rmax = 1.8;
-  if (angList.size() == 1)
+
+  // settle Llocal. Priority: current value (from input file) > PP XML file
+  if (local_channel > -1 && Llocal > -1)
+  {
+    app_warning() << " l-local setting found in both qmcpack input (l-local = " << Llocal
+                  << ") and pseudopotential file (l-local = " << local_channel << ")."
+                  << " Using l-local setting in qmcpack input file." << std::endl;
+  }
+  else if (local_channel > -1 && Llocal == -1)
+  {
+    app_log() << " l-local setting found in pseudopotential file and used." << std::endl;
+    Llocal = local_channel;
+  }
+  else if (local_channel == -1 && Llocal > -1)
+    app_log() << " l-local setting found in qmcpack input file and used." << std::endl;
+  else if (angList.size() == 1)
   {
     Llocal = Lmax;
     app_log() << "    Only one vps is found. Set the local component=" << Lmax << std::endl;
   }
+  else
+  {
+    app_error() << "The local channel is specified in neither the pseudopotential file nor the input file.\n"
+                << "Please add \'l-local=\"n\"\' attribute to either file.\n";
+    myComm->barrier_and_abort("ECPComponentBuilder::doBreakUp");
+  }
 
   if (angListSO.size() != nso)
   {
@@ -280,7 +303,7 @@ void ECPComponentBuilder::buildSemiLocalAndLocal(std::vector<xmlNodePtr>& semiPt
         vnnso[i][j] *= grid_global->r(j);
   }
   app_log() << "   Number of angular momentum channels " << angList.size() << std::endl;
-  app_log() << "   Maximum angular momentum channel " << Lmax << std::endl;
+  app_log() << "   Maximum angular momentum channel (Lmax) " << Lmax << std::endl;
   doBreakUp(angList, vnn, rmax, Vprefactor);
 
   //If any spinorbit terms are found...
@@ -487,13 +510,8 @@ void ECPComponentBuilder::doBreakUp(const std::vector<int>& angList,
   // If d is not reset, we generate an error in the interpolated PP!
   d        = agrid->Delta;
   int ngIn = vnn.cols() - 2;
-  if (Llocal == -1 && Lmax > 0)
-  {
-    app_error() << "The local channel is not specified in the pseudopotential file.\n"
-                << "Please add \'l-local=\"n\"\' attribute the semilocal section of the fsatom XML file.\n";
-    myComm->barrier_and_abort("ECPComponentBuilder::doBreakUp");
-    // Llocal = Lmax;
-  }
+
+  assert(angList.size() > 0 && "ECPComponentBuilder::doBreakUp angList cannot be empty!");
   //find the index of local
   int iLlocal = -1;
   for (int l = 0; l < angList.size(); l++)

src/QMCHamiltonians/ECPComponentBuilder.cpp

@@ -27,10 +27,11 @@
 
 namespace qmcplusplus
 {
-ECPComponentBuilder::ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule)
+ECPComponentBuilder::ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule, int llocal)
     : MPIObjectBase(c),
       NumNonLocal(0),
       Lmax(0),
+      Llocal(llocal),
       Nrule(nrule),
       Srule(8),
       AtomicNumber(0),

src/QMCHamiltonians/ECPComponentBuilder.h

@@ -48,7 +48,7 @@ struct ECPComponentBuilder : public MPIObjectBase, public QMCTraits
   std::unique_ptr<L2RadialPotential> pp_L2;
   std::map<std::string, int> angMon;
 
-  ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule = -1);
+  ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule = -1, int llocal = -1);
 
   bool parse(const std::string& fname, xmlNodePtr cur);
   bool put(xmlNodePtr cur);

src/QMCHamiltonians/ECPotentialBuilder.cpp

@@ -212,13 +212,15 @@ void ECPotentialBuilder::useXmlFormat(xmlNodePtr cur)
       std::string ionName("none");
       std::string format("xml");
       int nrule = -1;
+      int llocal = -1;
       //RealType rc(2.0);//use 2 Bohr
       OhmmsAttributeSet hAttrib;
       hAttrib.add(href, "href");
       hAttrib.add(ionName, "elementType");
       hAttrib.add(ionName, "symbol");
       hAttrib.add(format, "format");
       hAttrib.add(nrule, "nrule");
+      hAttrib.add(llocal, "l-local");
       //hAttrib.add(rc,"cutoff");
       hAttrib.put(cur);
       SpeciesSet& ion_species(IonConfig.getSpeciesSet());
@@ -232,7 +234,7 @@ void ECPotentialBuilder::useXmlFormat(xmlNodePtr cur)
       {
         app_log() << std::endl << "  Adding pseudopotential for " << ionName << std::endl;
 
-        ECPComponentBuilder ecp(ionName, myComm, nrule);
+        ECPComponentBuilder ecp(ionName, myComm, nrule, llocal);
         if (format == "xml")
         {
           if (href == "none")

src/QMCHamiltonians/tests/test_ecp.cpp

@@ -83,7 +83,7 @@ TEST_CASE("ReadFileBuffer_ecp", "[hamiltonian]")
 {
   Communicate* c = OHMMS::Controller;
 
-  ECPComponentBuilder ecp("test_read_ecp", c);
+  ECPComponentBuilder ecp("test_read_ecp", c, 4, 1);
 
   bool okay = ecp.read_pp_file("C.BFD.xml");
   REQUIRE(okay);

src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp

@@ -779,7 +779,7 @@ void MultiDiracDeterminant::mw_evaluateDetsAndGradsForPtclMove(
     for (size_t iw = 0; iw < nw; iw++)
     {
       GradType ratioGradRef_local(0);
-      //PRAGMA_OFFLOAD("omp parallel for reduction(+ : ratioGradRef_local)")
+      PRAGMA_OFFLOAD("omp parallel for reduction(+ : ratioGradRef_local)")
       for (size_t i = 0; i < NumPtcls; i++)
       {
         const size_t J            = confgListOccup_ptr[i];