Refactored qvs

fedluc · Sep 8, 2024 · 55f0530 · 55f0530
1 parent ab7dfdf
commit 55f0530
Show file tree

Hide file tree

Showing 9 changed files with 447 additions and 425 deletions.
diff --git a/include/qvs.hpp b/include/qvs.hpp
@@ -9,190 +9,180 @@
 #include <cmath>
 #include <memory>
 
-// // -----------------------------------------------------------------
-// // Class to handle simultaneous state point calculations
-// // -----------------------------------------------------------------
-
-// class QStlsCSR : public CSR<Vector2D, Qstls, QVSStlsInput> {
-
-// public:
-
-//   // Constructor
-//   explicit QStlsCSR(const QVSStlsInput &in_)
-//       : CSR(in_, Qstls(in_.toQstlsInput(), false, false)) {}
-//   // Compute auxiliary density response
-//   void computeAdrStls();
-//   void computeAdr();
-//   // Update the static structure factor
-//   void updateSsf() { ssf = ssfOld; };
-//   // Initialize the scheme
-//   void init();
-//   // Compute Q
-//   double getQAdder() const;
-
-// private:
-
-//   // Helper methods to compute the derivatives
-//   double getDerivative(const std::shared_ptr<Vector2D> &f,
-//                        const int &l,
-//                        const size_t &idx,
-//                        const Derivative &type);
-// };
-
-// // -----------------------------------------------------------------
-// // Class to handle the Q-adder in the free parameter expression
-// // -----------------------------------------------------------------
-
-// class QAdder {
-
-// public:
-
-//   // Constructor
-//   QAdder(const double &Theta_,
-//          const double &mu_,
-//          const double &limitMin,
-//          const double &limitMax,
-//          const std::vector<double> &itgGrid_,
-//          Integrator1D &itg1_,
-//          Integrator2D &itg2_,
-//          const Interpolator1D &interp_)
-//       : Theta(Theta_),
-//         mu(mu_),
-//         limits(limitMin, limitMax),
-//         itgGrid(itgGrid_),
-//         itg1(itg1_),
-//         itg2(itg2_),
-//         interp(interp_) {}
-//   // Get Q-adder
-//   double get() const;
-
-// private:
-
-//   const double lambda = pow(4.0 / (9.0 * M_PI), 1.0 / 3.0);
-//   // Degeneracy parameter
-//   const double Theta;
-//   // Chemical potential
-//   const double mu;
-//   // Integration limits
-//   const std::pair<double, double> limits;
-//   // Grid for 2D integration
-//   const std::vector<double> &itgGrid;
-//   // Integrator objects
-//   Integrator1D &itg1;
-//   Integrator2D &itg2;
-//   // Interpolator 1D class instance
-//   const Interpolator1D &interp;
-
-//   // SSF interpolation
-//   double ssf(const double &y) const;
-//   // Integrands
-//   double integrandDenominator(const double q) const;
-//   double integrandNumerator1(const double q) const;
-//   double integrandNumerator2(const double w) const;
-//   // Get Integral denominator
-//   void getIntDenominator(double &res) const;
-// };
-
-// // -----------------------------------------------------------------
-// // Class to handle the state point derivatives
-// // -----------------------------------------------------------------
-
-// class QStructProp : public StructPropBase<QStlsCSR, QVSStlsInput> {
-
-// public:
-
-//   // Constructor
-//   explicit QStructProp(const QVSStlsInput &in)
-//       : StructPropBase(in) {}
-//   // Get Q term
-//   std::vector<double> getQ() const;
-
-// private:
-
-//   using StructPropBase = StructPropBase<QStlsCSR, QVSStlsInput>;
-//   // Setup input for the CSR objects
-//   std::vector<QVSStlsInput> setupCSRInput(const QVSStlsInput &in);
-//   // Setup dependencies in the CSR objects
-//   void setupCSRDependencies();
-//   // Perform iterations to compute structural properties
-//   void doIterations();
-// };
-
-// // -----------------------------------------------------------------
-// // Class to handle the thermodynamic properties and derivatives
-// // -----------------------------------------------------------------
-
-// class QThermoProp : public ThermoPropBase<QStructProp, QVSStlsInput> {
-
-// public:
-
-//   // Typdef
-//   using ThermoPropBase = ThermoPropBase<QStructProp, QVSStlsInput>;
-//   // Constructors
-//   explicit QThermoProp(const QVSStlsInput &in)
-//       : ThermoPropBase(in) {}
-//   // Get internal energy and internal energy derivatives
-//   std::vector<double> getQData() const;
-//   // Get structural properties
-//   Vector2D getAdr();
-// };
-
-// // -----------------------------------------------------------------
-// // Solver for the qVS-STLS scheme
-// // -----------------------------------------------------------------
-
-// class QVSStls : public VSBase<QThermoProp, Rpa, QVSStlsInput> {
-
-// public:
-
-//   // Typedef
-//   using VSBase = VSBase<QThermoProp, Rpa, QVSStlsInput>;
-//   // Constructor from initial data
-//   explicit QVSStls(const QVSStlsInput &in_)
-//       : VSBase(in_) {}
-//   // Constructor for recursive calculations
-//   QVSStls(const QVSStlsInput &in_, const QThermoProp &thermoProp_)
-//       : VSBase(in_, thermoProp_) {}
-//   // Getters
-//   Vector2D getAdr() const { return adr; }
-
-// private:
-
-//   // Compute free parameter
-//   double computeAlpha();
-//   // Iterations to solve the qvsstls-scheme
-//   void updateSolution();
-//   // Setup free energy integrand
-//   void initFreeEnergyIntegrand();
-//   // Auxiliary density response
-//   Vector2D adr;
-// };
+class QThermoProp;
+class QStructProp;
+class QstlsCSR;
+class QAdder;
 
 // -----------------------------------------------------------------
 // Solver for the qVS-STLS scheme
 // -----------------------------------------------------------------
 
-class QVSStls : public VSBase {
+class QVSStls : public VSBase, public Qstls {
 
 public:
 
-  // Constructor
-  QVSStls(const QVSStlsInput &in_)
-      : VSBase(in_) {}
-  // Getters
-  Vector2D getAdr() const { return adr; }
+  // Constructor from initial data
+  explicit QVSStls(const QVSStlsInput &in_);
+  // Constructor for recursive calculations
+  QVSStls(const QVSStlsInput &in_, const QThermoProp &thermoProp_);
+  // Solve the scheme
+  using VSBase::compute;
 
 private:
 
+  // Input
+  QVSStlsInput in;
+  // Verbosity
+  using VSBase::verbose;
+  // Thermodyanmic properties
+  std::shared_ptr<QThermoProp> thermoProp;
+  // Initialize
+  void initScheme();
+  void initFreeEnergyIntegrand();
   // Compute free parameter
-  double computeAlpha() { return 0.0; }
+  double computeAlpha();
   // Iterations to solve the qvsstls-scheme
-  void updateSolution() {};
-  // Setup free energy integrand
-  void initScheme() {};
-  void initFreeEnergyIntegrand() {};
-  // Auxiliary density response
-  Vector2D adr;
+  void updateSolution();
+};
+
+// -----------------------------------------------------------------
+// Class to handle the thermodynamic properties
+// -----------------------------------------------------------------
+
+class QThermoProp : public ThermoPropBase {
+
+public:
+
+  // Constructors
+  explicit QThermoProp(const QVSStlsInput &in_);
+  // Get internal energy and internal energy derivatives
+  std::vector<double> getQData() const;
+  // Get structural properties
+  Vector2D getAdr();
+
+private:
+
+  std::shared_ptr<QStructProp> structProp;
+};
+
+// -----------------------------------------------------------------
+// Class to handle the structural properties
+// -----------------------------------------------------------------
+
+class QStructProp : public StructPropBase {
+
+public:
+
+  // Constructor
+  explicit QStructProp(const QVSStlsInput &in_);
+  // Get structural properties for output
+  const QstlsCSR &getCsr(const Idx &idx) const;
+  // Get Q term
+  std::vector<double> getQ() const;
+
+private:
+
+  // Vector containing NPOINTS state points to be solved simultaneously
+  std::vector<std::shared_ptr<QstlsCSR>> csr;
+  // Setup dependencies in the CSR objects
+  std::vector<QVSStlsInput> setupCSRInput(const QVSStlsInput &in);
+  void setupCSR(const QVSStlsInput &in_);
+  // Perform iterations to compute structural properties
+  void doIterations();
+};
+
+// -----------------------------------------------------------------
+// Class to solve one state point
+// -----------------------------------------------------------------
+
+class QstlsCSR : public CSR, public Qstls {
+
+public:
+
+  // Constructor
+  explicit QstlsCSR(const QVSStlsInput &in_)
+      : CSR(in_),
+        Qstls(in_.toQstlsInput(), false, false),
+        in(in_) {}
+  // Compute auxiliary density response
+  void computeAdrQStls();
+  void computeAdr();
+  // Publicly esposed private stls methods
+  void init();
+  void initialGuess() { Qstls::initialGuess(); }
+  void computeSsf() { Qstls::computeSsf(); }
+  double computeError() { return Qstls::computeError(); }
+  void updateSolution() { Qstls::updateSolution(); }
+  // Update the static structure factor
+  void updateSsf() { ssf = ssfOld; };
+  // Compute Q
+  double getQAdder() const;
+  // Getters
+  std::vector<double> getSsf() const { return Qstls::getSsf(); }
+  std::vector<double> getSlfc() const { return Qstls::getSlfc(); }
+  std::vector<double> getWvg() const { return Qstls::getWvg(); }
+
+private:
+
+  // Input
+  QVSStlsInput in;
+  // Compute derivative contribution to the auxiliary density response
+  Vector2D getDerivativeContribution() const;
+};
+
+// -----------------------------------------------------------------
+// Class to handle the Q-adder in the free parameter expression
+// -----------------------------------------------------------------
+
+class QAdder {
+
+public:
+
+  // Constructor
+  QAdder(const double &Theta_,
+         const double &mu_,
+         const double &limitMin,
+         const double &limitMax,
+         const std::vector<double> &itgGrid_,
+         Integrator1D &itg1_,
+         Integrator2D &itg2_,
+         const Interpolator1D &interp_)
+      : Theta(Theta_),
+        mu(mu_),
+        limits(limitMin, limitMax),
+        itgGrid(itgGrid_),
+        itg1(itg1_),
+        itg2(itg2_),
+        interp(interp_) {}
+  // Get Q-adder
+  double get() const;
+
+private:
+
+  const double lambda = pow(4.0 / (9.0 * M_PI), 1.0 / 3.0);
+  // Degeneracy parameter
+  const double Theta;
+  // Chemical potential
+  const double mu;
+  // Integration limits
+  const std::pair<double, double> limits;
+  // Grid for 2D integration
+  const std::vector<double> &itgGrid;
+  // Integrator objects
+  Integrator1D &itg1;
+  Integrator2D &itg2;
+  // Interpolator 1D class instance
+  const Interpolator1D &interp;
+
+  // SSF interpolation
+  double ssf(const double &y) const;
+  // Integrands
+  double integrandDenominator(const double q) const;
+  double integrandNumerator1(const double q) const;
+  double integrandNumerator2(const double w) const;
+  // Get Integral denominator
+  void getIntDenominator(double &res) const;
 };
 
 #endif
diff --git a/include/vector2D.hpp b/include/vector2D.hpp
@@ -80,6 +80,9 @@ class Vector2D {
   // Element-wise division between this vector and v_
   void div(const Vector2D &v_);
 
+  // Element-wise linear combination v += num*v_
+  void linearCombination(const Vector2D &v_, const double &num_);
+
 private:
 
   // Underlying vector structure with s1*s2 entries

diff --git a/include/vector3D.hpp b/include/vector3D.hpp
@@ -81,6 +81,9 @@ class Vector3D {
   // Element-wise division between this vector and v_
   void div(const Vector3D &v_);
 
+  // Element-wise linear combination v += num*v_
+  void linearCombination(const Vector3D &v_, const double &num_);
+
 private:
 
   // Underlying vector structure with s1*s2*s3 entries

diff --git a/include/vector_util.hpp b/include/vector_util.hpp
@@ -28,6 +28,12 @@ namespace vecUtil {
   // Element-wise multiplication of a vector and a scalar
   std::vector<double> mult(const std::vector<double> &v, const double a);
 
+  // Linear combination of two vectors
+  std::vector<double> linearCombination(const std::vector<double> &v1,
+                                        const double a,
+                                        const std::vector<double> &v2,
+                                        const double b);
+
   // Root mean square difference between two vectors
   double rms(const std::vector<double> &v1,
              const std::vector<double> &v2,