Merge pull request #4800 from PDoakORNL/reference_points

Class for handling reference points in EnergyDensityEstimator

src/Estimators/CMakeLists.txt

@@ -33,7 +33,9 @@ set(QMCEST_SRC
     MagnetizationDensity.cpp
     MagnetizationDensityInput.cpp
     PerParticleHamiltonianLoggerInput.cpp
-    PerParticleHamiltonianLogger.cpp)
+    PerParticleHamiltonianLogger.cpp
+    ReferencePointsInput.cpp
+    NEReferencePoints.cpp)
 
 ####################################
 # create libqmcestimators

src/Estimators/NEReferencePoints.cpp

@@ -0,0 +1,116 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2023 QMCPACK developers.
+//
+// File developed by: Jaron T. Krogel, [email protected], Oak Ridge National Laboratory
+//                    Mark A. Berrill, [email protected], Oak Ridge National Laboratory
+//                    Peter W. Doak. [email protected], Oak Ridge National Laboratory
+//
+// File refactored from: QMCHamiltonian/ReferencePoints.cpp
+//////////////////////////////////////////////////////////////////////////////////////
+
+
+#include "NEReferencePoints.h"
+#include "Utilities/string_utils.h"
+#include "OhmmsData/AttributeSet.h"
+#include "QMCHamiltonians/ObservableHelper.h"
+
+namespace qmcplusplus
+{
+
+NEReferencePoints::NEReferencePoints(const ReferencePointsInput& rp_input,
+                                     const ParticleSet& pset,
+                                     RefVector<ParticleSet>& ref_psets)
+    : input_(rp_input)
+{
+  processParticleSets(pset, ref_psets);
+  for (int i = 0; i < OHMMS_DIM; i++)
+    for (int d = 0; d < OHMMS_DIM; d++)
+      axes(d, i) = pset.getLattice().a(i)[d];
+  Axes crd;
+  // no need to handle error here rp_input will have a valid value for coord_form
+  switch (input_.get_coord_form())
+  {
+  case Coord::CELL:
+    crd = axes;
+    break;
+  case Coord::CARTESIAN:
+    for (int i = 0; i < OHMMS_DIM; i++)
+      for (int d = 0; d < OHMMS_DIM; d++)
+        if (d == i)
+          crd(i, i) = 1.0;
+        else
+          crd(d, i) = 0.0;
+    break;
+  }
+
+  for (const auto& [key, value] : input_.get_points())
+    points_[key] = dot(crd, value);
+}
+
+void NEReferencePoints::processParticleSets(const ParticleSet& P, RefVector<ParticleSet>& Psets)
+{
+  //get axes and origin information from the ParticleSet
+  points_["zero"] = 0 * P.getLattice().a(0);
+  points_["a1"]   = P.getLattice().a(0);
+  points_["a2"]   = P.getLattice().a(1);
+  points_["a3"]   = P.getLattice().a(2);
+  //points_["center"]= .5*(P.getLattice().a(0)+P.getLattice().a(1)+P.Lattice.a(2))
+  //set points_ on face centers
+  points_["f1p"] = points_["zero"] + .5 * points_["a1"];
+  points_["f1m"] = points_["zero"] - .5 * points_["a1"];
+  points_["f2p"] = points_["zero"] + .5 * points_["a2"];
+  points_["f2m"] = points_["zero"] - .5 * points_["a2"];
+  points_["f3p"] = points_["zero"] + .5 * points_["a3"];
+  points_["f3m"] = points_["zero"] - .5 * points_["a3"];
+  //set points_ on cell corners
+  points_["cmmm"] = points_["zero"] + .5 * (-1 * points_["a1"] - points_["a2"] - points_["a3"]);
+  points_["cpmm"] = points_["zero"] + .5 * (points_["a1"] - points_["a2"] - points_["a3"]);
+  points_["cmpm"] = points_["zero"] + .5 * (-1 * points_["a1"] + points_["a2"] - points_["a3"]);
+  points_["cmmp"] = points_["zero"] + .5 * (-1 * points_["a1"] - points_["a2"] + points_["a3"]);
+  points_["cmpp"] = points_["zero"] + .5 * (-1 * points_["a1"] + points_["a2"] + points_["a3"]);
+  points_["cpmp"] = points_["zero"] + .5 * (points_["a1"] - points_["a2"] + points_["a3"]);
+  points_["cppm"] = points_["zero"] + .5 * (points_["a1"] + points_["a2"] - points_["a3"]);
+  points_["cppp"] = points_["zero"] + .5 * (points_["a1"] + points_["a2"] + points_["a3"]);
+  //get points from requested particle sets
+  int cshift = 1;
+  for (ParticleSet& pset : Psets)
+  {
+    for (int p = 0; p < pset.getTotalNum(); p++)
+    {
+      std::stringstream ss;
+      ss << p + cshift;
+      points_[pset.getName() + ss.str()] = pset.R[p];
+    }
+  }
+}
+
+void NEReferencePoints::write_description(std::ostream& os, const std::string& indent) const
+{
+  os << indent + "reference_points" << std::endl;
+  std::map<std::string, Point>::const_iterator it, end = points_.end();
+  for (it = points_.begin(); it != end; ++it)
+  {
+    os << indent + "  " << it->first << ": " << it->second << std::endl;
+  }
+  os << indent + "end reference_points" << std::endl;
+  return;
+}
+
+void NEReferencePoints::write(hdf_archive& file) const
+{
+  file.push(std::string_view("reference_points"));
+  for (auto it = points_.cbegin(); it != points_.cend(); ++it)
+    file.write(const_cast<Point&>(it->second), it->first);
+  file.pop();
+}
+
+std::ostream& operator<<(std::ostream& out, const NEReferencePoints& rhs)
+{
+  rhs.write_description(out, "");
+  return out;
+}
+
+} // namespace qmcplusplus

src/Estimators/NEReferencePoints.h

@@ -0,0 +1,82 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2023 QMCPACK developers.
+//
+// File developed by: Jaron T. Krogel, [email protected], Oak Ridge National Laboratory
+//                    Mark A. Berrill, [email protected], Oak Ridge National Laboratory
+//                    Peter W. Doak. [email protected], Oak Ridge National Laboratory
+//
+// File refactored from: QMCHamiltonian/ReferencePoints.h
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifndef QMCPLUSPLUS_NEREFERENCE_POINTS_H
+#define QMCPLUSPLUS_NEREFERENCE_POINTS_H
+
+/** @file
+ *  \todo When QMCHamiltonian/ReferencePoints is removed rename this class to ReferencePoints
+ */
+
+#include <Configuration.h>
+#include "OhmmsData/OhmmsElementBase.h"
+#include "Particle/ParticleSet.h"
+#include "QMCHamiltonians/ObservableHelper.h"
+#include "OhmmsPETE/Tensor.h"
+#include "ReferencePointsInput.h"
+
+namespace qmcplusplus
+{
+
+/** This class creates, contains, and writes both user and machine readable referencepoints.
+ *  they are derived from the lattice of the pset passed and the particle positions in the ref_psets
+ *  an arbitrary number of additional points can be defined in the input that ReferencePointsInput
+ *  presents as native input.
+ *  It is a dependency of Estimators/NESpaceGrid and Estimatorss/EnergyDensityEstimator
+ */
+class NEReferencePoints
+{
+public:
+  using Real   = QMCTraits::RealType;
+  using Axes   = Tensor<Real, OHMMS_DIM>;
+  using Point  = TinyVector<Real, OHMMS_DIM>;
+  using Points = std::map<std::string, Point>;
+  using Coord  = typename ReferencePointsInput::Coord;
+  /** Usual constructor
+   *  \param[in] rp_input   Input object for reference points which can contain and arbitrary set of points beyond
+   *                        those take from the pset, and ref_psets
+   *  \param[in] pset       pset that supplies the lattice information for reference points
+   *  \param[in] ref_psets  pset reference vector the particle points in this/these psets are reference points
+   */
+  NEReferencePoints(const ReferencePointsInput& rp_input, const ParticleSet& pset, RefVector<ParticleSet>& ref_psets);
+  NEReferencePoints(const NEReferencePoints& nerp) = default;
+
+  /** writes a human readable representation of the reference points.
+   *  \param[inout] os       ostream to write description to
+   *  \param[in]    indent   spaces or other text to preface each line of output with.  needed to preserve
+   *                         legacy output format.
+   */
+  void write_description(std::ostream& os, const std::string& indent) const;
+
+  /** machine readable output
+   *  \param[inout]  file  hdf5 file to write to.  Respects current state of file.
+   */
+  void write(hdf_archive& file) const;
+
+  /** return const ref to map of reference points.
+   *  labeling scheme unchanged from legacy
+   */
+  const Points& get_points() const { return points_; }
+
+protected:
+  Points points_;
+private:
+  void processParticleSets(const ParticleSet& P, RefVector<ParticleSet>& Pref);
+  Axes axes;
+  ReferencePointsInput input_;
+};
+
+std::ostream& operator<<(std::ostream& out, const NEReferencePoints& rhs);
+
+}
+#endif

src/Estimators/ReferencePointsInput.cpp

@@ -0,0 +1,88 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2022 QMCPACK developers.
+//
+// File developed by: Jaron T. Krogel, [email protected], Oak Ridge National Laboratory
+//                    Mark A. Berrill, [email protected], Oak Ridge National Laboratory
+//                    Peter W. Doak, [email protected], Oak Ridge National Laboratory
+//
+// Some code refactored from:  QMCHamiltonians/ReferencePoints.cpp
+//////////////////////////////////////////////////////////////////////////////////////
+
+#include "ReferencePointsInput.h"
+
+#include <sstream>
+
+#include "EstimatorInput.h"
+#include "ModernStringUtils.hpp"
+
+namespace qmcplusplus
+{
+
+ReferencePointsInput::ReferencePointsInput(xmlNodePtr cur)
+{
+  input_section_.readXML(cur);
+  auto setIfInInput = LAMBDA_setIfInInput;
+  setIfInInput(coord_form_, "coord");
+  readRefPointsXML(cur);
+}
+
+// ReferencePointsInput::ReferencePointsInput(const Points& points, const CoordForm coord_form)
+//     : points_(points), coord_form_(coord_form)
+// {}
+
+void ReferencePointsInput::readRefPointsXML(xmlNodePtr cur)
+{
+  using modernstrutil::split;
+  using modernstrutil::strip;
+
+  // read refpoints values they are some sequence of value nodes
+  std::string node_str{XMLNodeString{cur}};
+  std::vector<std::string_view> lines = split(strip(node_str), "\n");
+  for (int i = 0; i < lines.size(); i++)
+  {
+    auto stripped_line                   = strip(lines[i]);
+    std::vector<std::string_view> tokens = split(stripped_line, " ");
+    if (tokens.size() != OHMMS_DIM + 1)
+    {
+      std::ostringstream error;
+      error << error_tag << "reference point has 4 entries, given " << tokens.size() << ": " << lines[i];
+      throw UniformCommunicateError(error.str());
+    }
+    else
+    {
+      Point rp;
+      for (int d = 0; d < OHMMS_DIM; d++)
+      {
+        try
+        {
+          rp[d] = std::stod(std::string(tokens[d + 1].begin(), tokens[d + 1].size()));
+        }
+        catch (const std::invalid_argument& ia)
+        {
+          throw UniformCommunicateError(ia.what());
+        }
+      }
+
+      // This must be done in constructor of ReferencePoints
+      // rp                = dot(crd, rp);
+      points_[std::string(tokens[0].begin(), tokens[0].size())] = rp;
+    }
+  }
+}
+
+std::any ReferencePointsInput::ReferencePointsInputSection::assignAnyEnum(const std::string& name) const
+{
+  return lookupAnyEnum(name, get<std::string>(name), lookup_input_enum_value);
+}
+
+std::any makeReferencePointsInput(xmlNodePtr cur, std::string& value_label)
+{
+  ReferencePointsInput rpi{cur};
+  value_label = "referencepoints";
+  return rpi;
+}
+
+} // namespace qmcplusplus