Refactor physical properties calculations

src/CMakeLists.txt

@@ -11,6 +11,7 @@ set(SINGLEION_SOURCES
     singleion/ic_socf.cpp
     singleion/ic_tensorops.cpp
     singleion/ic1ion.cpp
+    singleion/physprop.cpp
 )
 
 # Adds the source modules as object libraries

src/include/cf1ion.hpp

@@ -11,34 +11,43 @@
 
 #include "eigen.hpp"
 #include "cfpars.hpp"
+#include "physprop.hpp"
 
 namespace libMcPhase {
 
-class cf1ion: public cfpars {
+class cf1ion: public cfpars, public physprop {
 
     protected:
         bool m_ham_calc = false;                              // Flag to indicate if Hamiltonian calculated
+        bool m_upper = true;                                  // Flag to indicate if upper triangle of Ham is calc
         bool m_ev_calc = false;                               // Flag to indicate if eigenvectors/values calculated
+        bool m_magops_calc = false;                           // Flag to indicate if magnetic operators calculated
         RowMatrixXcd m_hamiltonian;                           // Cached Hamiltonian
         RowMatrixXcd m_eigenvectors;                          // Cached eigenvectors
         VectorXd m_eigenvalues;                               // Cached eigenvalues
+        std::vector<RowMatrixXcd> m_magops;                   // Cached magnetic operators
+        RowMatrixXcd _hamiltonian(bool upper=true);
+        void calc_mag_ops();
+        void fill_upper();
 
     public:
         // Setters
-        virtual void set_unit(const Units newunit);
-        virtual void set_type(const Type newtype);
-        virtual void set_name(const std::string &ionname);
-        virtual void set_J(const double J);
-        virtual void set(const Blm blm, double val);
-        virtual void set(int l, int m, double val);
+        void set_unit(const Units newunit);
+        void set_type(const Type newtype);
+        void set_name(const std::string &ionname);
+        void set_J(const double J);
+        void set(const Blm blm, double val);
+        void set(int l, int m, double val);
         // Constructors
         cf1ion() : cfpars() {};
         cf1ion(const int J2) : cfpars(J2) {};
         cf1ion(const double J) : cfpars(J) {};
         cf1ion(const std::string &ionname) : cfpars(ionname) {};
         // Methods
-        RowMatrixXcd hamiltonian(bool upper=true);
+        RowMatrixXcd hamiltonian();
         std::tuple<RowMatrixXcd, VectorXd> eigensystem();
+        RowMatrixXcd zeeman_hamiltonian(double H, std::vector<double> Hdir);
+        std::vector<RowMatrixXcd> calculate_moments_matrix(RowMatrixXcd ev);
 
 }; // class cf1ion
 

src/include/cfpars.hpp

@@ -24,10 +24,8 @@
 
 namespace libMcPhase {
 
-// Boltzmann constant and bohr magneton defined in cfpars.cpp and ic1ion.cpp because they depend on energy
-// and the default energy unit in cfpars is meV whilst in ic1ion it is cm
-static const double NAMUB = 5.5849397;         // N_A*mu_B - J/T/mol - product of Bohr magneton and Avogadro's number
-static const double GS = 2.0023193043622;      // The electron gyromagnetic ratio
+// Conversion factors for different energy units[from][to], order: [meV, cm, K].
+static const std::array<double, 3> ENERGYCONV = { {1., 8.065544005, 11.6045221} };
 
 class cfpars {
 
@@ -39,7 +37,6 @@ class cfpars {
                     B44S = 5, B43S = 6, B42S = 7, B41S = 8, B40 = 9, B41 = 10, B42 = 11, B43 = 12, B44 = 13,
                     B66S = 14, B65S = 15, B64S = 16, B63S = 17, B62S = 18, B61S = 19, 
                     B60 = 20, B61 = 21, B62 = 22, B63 = 23, B64 = 24, B65 = 25, B66 = 26};
-    enum class MagUnits {bohr = 0, cgs = 1, SI = 2};
 
     protected:
         std::array<double, 27> m_Bi{};                        // Internal array of values (in Wybourne/theta_k in meV)
@@ -53,6 +50,7 @@ class cfpars {
         Type m_type = Type::Blm;                              // Type of crystal field parameters
         std::string m_ionname;                                // Name of ion
         int m_J2 = 0;                                         // 2*J == twice the total angular momentum
+        double m_GJ = -1.;                                    // Lande g-factor
         bool m_convertible = false;                           // True if can convert between types and normalisations
         racah m_racah{};                                      // Class to calc n-j symbols and cache factorials
         int m_n = 1;                                          // Number of open shell electrons in this configuration
@@ -61,21 +59,23 @@ class cfpars {
         // Methods
         virtual void getfromionname(const std::string &ionname);
         // Getters
-        const Units get_unit() const { return m_unit; }
-        const Normalisation get_normalisation() const { return m_norm; }
-        const Type get_type() const { return m_type; }
-        const std::string get_name() const { return m_ionname; }
-        const double get(const Blm blm) const { return m_Bo[(int)blm]; }
-        const double get(int l, int m) const;
-        const double alpha() { return m_stevfact[0]; }
-        const double beta() { return m_stevfact[1]; }
-        const double gamma() { return m_stevfact[2]; }
-        const double get_J() { return (double)(m_J2 / 2.); }
+        Units get_unit() const { return m_unit; }
+        Normalisation get_normalisation() const { return m_norm; }
+        Type get_type() const { return m_type; }
+        std::string get_name() const { return m_ionname; }
+        double get(const Blm blm) const { return m_Bo[(int)blm]; }
+        double get(int l, int m) const;
+        double get_GJ() const { return m_GJ; }
+        double alpha() const { return m_stevfact[0]; }
+        double beta() const { return m_stevfact[1]; }
+        double gamma() const { return m_stevfact[2]; }
+        double get_J() const { return (double)(m_J2 / 2.); }
         // Setters
         virtual void set_unit(const Units newunit);
         virtual void set_type(const Type newtype);
         virtual void set_name(const std::string &ionname);
         virtual void set_J(const double J);
+        virtual void set_GJ(const double GJ) { m_GJ = GJ; }
         virtual void set(const Blm blm, double val);
         virtual void set(int l, int m, double val);
         // Constructors

src/include/ic1ion.hpp

@@ -13,17 +13,15 @@
 #define IC1ION_H
 
 #include "cfpars.hpp"
+#include "physprop.hpp"
 #include "eigen.hpp"
 #include "ic_states.hpp"
 #include <algorithm>
 #include <vector>
-#include <unordered_map>
 
 namespace libMcPhase {
 
-// Declare K_B and MU_B in cm here (meV version declared in cfpars.cpp)
-static const double K_B = 0.6950348004;    // cm/K - Boltzmann constant
-static const double MU_B = 0.46686447783;  // cm/T - Bohr magneton
+static const double GS = 2.0023193043622;  // The electron gyromagnetic ratio
 
 struct pair_hash
 {
@@ -33,7 +31,7 @@ struct pair_hash
     }
 };
 
-class ic1ion : public cfpars {
+class ic1ion : public cfpars, public physprop {
 
     public:
         enum class CoulombType { Slater = 0, CondonShortley = 1, Racah = 2 };
@@ -80,18 +78,17 @@ class ic1ion : public cfpars {
         void calculate_eigensystem();                          // Diagonalises the Hamiltonian
         // Declarations for functions in ic_tensoropts.cpp
         void calc_tensorops(int num);                          // Populates m_tensorops vector
-        std::vector< RowMatrixXcd > calculate_moments_matrix(RowMatrixXcd ev);
 
     public:
         // Setters
-        virtual void set_unit(const Units newunit);
-        virtual void set_type(const Type newtype);
-        virtual void set_name(const std::string &ionname);
-        virtual void set(const Blm blm, double val);
-        virtual void set(int l, int m, double val);
-        virtual void set_coulomb(std::vector<double> val, CoulombType type = CoulombType::Slater);
-        virtual void set_ci(std::vector<double> val);
-        virtual void set_spinorbit(double val, SpinOrbType type = SpinOrbType::Zeta);
+        void set_unit(const Units newunit);
+        void set_type(const Type newtype);
+        void set_name(const std::string &ionname);
+        void set(const Blm blm, double val);
+        void set(int l, int m, double val);
+        void set_coulomb(std::vector<double> val, CoulombType type = CoulombType::Slater);
+        void set_ci(std::vector<double> val);
+        void set_spinorbit(double val, SpinOrbType type = SpinOrbType::Zeta);
         // Getters
         std::array<double, 4> get_coulomb() const { return m_F; };
         double get_spinorbit() const { return m_xi; };
@@ -104,13 +101,9 @@ class ic1ion : public cfpars {
         // Methods
         RowMatrixXcd hamiltonian();
         std::tuple<RowMatrixXcd, VectorXd> eigensystem();
-        std::vector<double> magnetisation(std::vector<double> H, std::vector<double> Hdir, double T, MagUnits type);
-        std::vector<double> susceptibility(std::vector<double> T, std::vector<double> Hdir, MagUnits type);
         RowMatrixXcd zeeman_hamiltonian(double H,              // Calculates the Zeeman Hamiltonian for applied
             std::vector<double> Hdir);                         //   field H in direction Hdir
-        std::vector< std::vector<double> > calculate_moments(RowMatrixXcd ev);
-        std::vector<double> calculate_boltzmann(               // Calculates the Boltzmann factor exp(-E/kT)
-            VectorXd en, double T);
+        std::vector< RowMatrixXcd > calculate_moments_matrix(RowMatrixXcd ev);
         std::vector<fstates_t> get_states();
 
 }; // class ic1ion

src/include/ic_states.hpp

@@ -72,7 +72,7 @@ class qR7
 
       // Other member functions //
       int set(int w1_, int w2_, int w3_);					// Set value
-      bool isequal (const char * Wstr);						// (isequal to string)
+      bool isequal (int w);
 
       // Overloaded operators //
       bool operator == (const qR7 & wp) const;					// (isequal)
@@ -102,7 +102,7 @@ class qG2
 
       // Other member functions //
       int set(int u1_, int u2_);						// Set value
-      bool isequal (const char * Ustr);						// (isequal to string)
+      bool isequal (int u);
 
       // Overloaded operators //
       bool operator == (const qG2 & up) const;					// (isequal)

src/include/physprop.hpp

@@ -0,0 +1,69 @@
+/* physprop.hpp
+ *
+ * Header file for the physical properties base class
+ *
+ * (c) 2024 Duc Le - [email protected]
+ * This program is licensed under the GNU General Purpose License, version 3. Please see the COPYING file
+ */
+
+#pragma once
+
+#include "eigen.hpp"
+#include <vector>
+#include <unordered_map>
+
+namespace libMcPhase {
+
+// Basic physical constants (needs cm constants as internal energy units in ic1ion is cm)
+static const double K_B = 0.08617343183;       // meV/K - Boltzmann constant
+static const double MU_B = 0.0578838263;       // meV/T - Bohr magneton
+static const double K_Bc = 0.6950348004;       // cm/K - Boltzmann constant
+static const double MU_Bc = 0.46686447783;     // cm/T - Bohr magneton
+
+// Conversion factors for different magnetic units for magnetisation. Order: [bohr, cgs, SI].
+// NAMUB is N_A * MU_B in J/T/mol == Am^2/mol is the SI unit. 
+// The cgs unit is N_A * MU_B in erg/G/mol == emu/mol is different only by a factor of 1000 larger
+static const double NAMUB = 5.5849397;  // N_A*mu_B - J/T/mol - product of Bohr magneton and Avogadro's number
+static const std::array<double, 3> MAGCONV = {1., NAMUB*1000, NAMUB};
+
+// Note these constants are strange because the default energy unit in this module is cm-1
+static const double NAMUBSQc_ERG = 0.26074098;     // N_A * muB[erg/G] * muB[cm/T] * [Tesla_to_Gauss=1e-4]
+static const double NAMUBSQc_JOULE = 3.276568e-06; // N_A * muB[J/T] * muB[cm/T] * mu0
+static const double NAMUBSQ_ERG = 0.03232776;      // N_A * muB[erg/G] * muB[meV/T] * [Tesla_to_Gauss=1e-4]
+static const double NAMUBSQ_JOULE = 4.062426e-07;  // N_A * muB[J/T] * muB[meV/T] * mu0
+// Factor of mu0 is needed in SI due to different definition of the magnetisation, B and H fields
+
+// Conversion factors for different magnetic units for magnetic susceptibility. Order: [bohr, cgs, SI].
+// The susceptibility prefactor is (in principle) N_A * MU_B^2 but we need to account for various units...
+// The atomic (bohr) susceptibility is in uB/T/ion; cgs is in erg/G^2/mol==cm^3/mol; SI in J/T^2/mol==m^3/mol
+// Note that chi_SI = (4pi*10^-6)chi_cgs [*not* 4pi*10-7!]
+static const std::array<double, 3> SUSCCONV = {MU_B, NAMUBSQ_ERG, NAMUBSQ_JOULE};
+
+// Conversion factor for heat capacity calculations
+static const double NAMEV = 96.48533212;          // J/mol = N_A * meV
+
+// EPSILON to determine if energy levels are degenerate or not
+static const double DELTA_EPS = 1e-6;
+
+// Base class for physical properties calculations. Must derive and implement zeeman_hamiltonian
+class physprop {
+    protected:
+        double m_meVconv = 1.0;   // Conversion factor from user energy units to meV
+
+    public:
+        enum class MagUnits {bohr = 0, cgs = 1, SI = 2};
+        virtual RowMatrixXcd hamiltonian() = 0;
+        virtual RowMatrixXcd zeeman_hamiltonian(double H, std::vector<double> Hdir) = 0;
+        virtual std::tuple<RowMatrixXcd, VectorXd> eigensystem() = 0;
+        virtual std::vector<RowMatrixXcd> calculate_moments_matrix(RowMatrixXcd ev) = 0;
+        VectorXd calculate_boltzmann(VectorXd en, double T);
+        VectorXd heatcapacity(std::vector<double> Tvec);
+        RowMatrixXd magnetisation(std::vector<double> H, std::vector<double> Hdir, std::vector<double> Tvec, MagUnits type);
+        VectorXd susceptibility(std::vector<double> T, std::vector<double> Hdir, MagUnits type);
+};
+
+// Mapping for Python binding to map string to enum
+static const std::unordered_map<std::string, physprop::MagUnits> mag_unit_names = {
+    {"bohr", physprop::MagUnits::bohr}, {"cgs", physprop::MagUnits::cgs}, {"SI", physprop::MagUnits::SI} };
+
+} // namespace libMcPhase

src/libmcphase/pycf1ion.cpp

@@ -9,14 +9,13 @@
 #include "cf1ion.hpp"
 #include <pybind11/pybind11.h>
 #include <pybind11/eigen.h>
+#include <pybind11/stl.h>
 #include "pycfpars.hpp"
 
 namespace py = pybind11;
 using namespace libMcPhase;
 using namespace pybind11::literals;
 
-void cf_parse(cfpars *cls, py::args args, py::kwargs kwargs);
-
 cf1ion *cf1ion_init(py::args args, py::kwargs kwargs) {
     cf1ion *cls = new cf1ion;
     cf_parse(static_cast<cfpars*>(cls), args, kwargs);
@@ -31,8 +30,16 @@ void wrap_cf1ion(py::module &m) {
         .def(py::init<const std::string &>(), py::arg("ionname"))
         .def(py::init<const double &>(), py::arg("J"))
         .def(py::init(&cf1ion_init), cfpars_init_str)
-        .def("hamiltonian", &cf1ion::hamiltonian, "the crystal field Hamiltonian", "upper"_a=true)
-        .def("eigensystem", &cf1ion::eigensystem, "the eigenvectors and eigenvalues of the crystal field Hamiltonian");
+        .def_property("GJ", [](cf1ion const &self) { return self.get_GJ(); }, [](cf1ion &self, double v) { self.set_GJ(v); })
+        .def("hamiltonian", &cf1ion::hamiltonian, "the crystal field Hamiltonian")
+        .def("eigensystem", &cf1ion::eigensystem, "the eigenvectors and eigenvalues of the crystal field Hamiltonian")
+        .def("zeeman_hamiltonian", &cf1ion::zeeman_hamiltonian, "the Zeeman Hamiltonian")
+        .def("calculate_boltzmann", &cf1ion::calculate_boltzmann, "")
+        .def("heatcapacity", &cf1ion::heatcapacity, "the heat capacity of the crystal field Hamiltonian in J/mol/K")
+        .def("magnetisation", [](cf1ion &self, std::vector<double> H, std::vector<double> Hdir, std::vector<double> T, std::string unit) { return self.magnetisation(H, Hdir, T,
+             set_enum(unit, mag_unit_names, "Invalid magnetic unit, must be one of: 'bohr', 'cgs', or 'SI'")); })
+        .def("susceptibility", [](cf1ion &self, std::vector<double> T, std::vector<double> Hdir, std::string unit) { return self.susceptibility(T, Hdir,
+             set_enum(unit, mag_unit_names, "Invalid magnetic unit, must be one of: 'bohr', 'cgs', or 'SI'")); });
 }
 
 

src/libmcphase/pycfpars.cpp

@@ -6,8 +6,6 @@
  * This program is licensed under the GNU General Purpose License, version 3. Please see the LICENSE file
  */
 
-#include "cfpars.hpp"
-#include <pybind11/pybind11.h>
 #include "pycfpars.hpp"
 
 namespace py = pybind11;
@@ -29,25 +27,28 @@ static const std::unordered_map<std::string, cfpars::Units> unit_names = {
     {"meV", cfpars::Units::meV}, {"cm", cfpars::Units::cm}, {"K", cfpars::Units::K} };
 static const std::string unit_err = "Invalid unit, must be one of 'meV', 'cm', or 'K'";
 
-void cf_parse(cfpars *cls, py::args args, py::kwargs kwargs) {
+void cf_parse(cfpars *cls, py::args args, py::kwargs kwargs, bool is_ic1ion) {
     if (!args && !kwargs) {
         return;
     }
     if (args && args.size() > 0) {
         try {
-            double J = args[0].cast<double>();
-            cls->set_J(J);
+            std::string ionname = args[0].cast<std::string>();
+            cls->set_name(ionname);
         } catch (py::cast_error) {
+            if (is_ic1ion) {
+                throw std::runtime_error("ic1ion module can only be initialised by an ion name");
+            }
             try {
-                std::string ionname = args[0].cast<std::string>();
-                cls->set_name(ionname);
+                double J = args[0].cast<double>();
+                cls->set_J(J);
             } catch (py::cast_error) {
                 throw std::runtime_error("Invalid first argument: must be the ion name as a string "
                                          "or the total angular momentum quantum number J");
             }
         }
     }
-    else if (kwargs.contains("J")) {
+    else if (!is_ic1ion && kwargs.contains("J")) {
         cls->set_J(kwargs["J"].cast<double>());
     }
     else if (kwargs.contains("ionname")) {

src/libmcphase/pycfpars.hpp

@@ -8,6 +8,11 @@
 
 #pragma once
 
+#include "cfpars.hpp"
+#include <pybind11/pybind11.h>
+namespace py = pybind11;
+using namespace libMcPhase;
+
 template <typename T> T set_enum(std::string key, std::unordered_map<std::string, T> enum_map, std::string errmsg) {
     auto it = enum_map.find(key);
     if (it == enum_map.end()) {
@@ -27,3 +32,4 @@ static const char* cfpars_init_str = "Construct a cfpars object\n"
                                      "               see online documentation or McPhase webpage for definitions\n"
                                      "        Blm - value of parameters, e.g. cfp = cfpars('Pr3+', B20=0.1, B22=-0.01, B40=0.001)\n";
 
+void cf_parse(cfpars *cls, py::args args, py::kwargs kwargs, bool is_ic1ion=false);