diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 635e2c9..253b964 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/include/cf1ion.hpp b/src/include/cf1ion.hpp index e095212..e834ce9 100644 --- a/src/include/cf1ion.hpp +++ b/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 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 eigensystem(); + RowMatrixXcd zeeman_hamiltonian(double H, std::vector Hdir); + std::vector calculate_moments_matrix(RowMatrixXcd ev); }; // class cf1ion diff --git a/src/include/cfpars.hpp b/src/include/cfpars.hpp index 7acddeb..744579e 100644 --- a/src/include/cfpars.hpp +++ b/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 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 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 diff --git a/src/include/ic1ion.hpp b/src/include/ic1ion.hpp index d92d887..8efe9cc 100644 --- a/src/include/ic1ion.hpp +++ b/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 #include -#include 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 val, CoulombType type = CoulombType::Slater); - virtual void set_ci(std::vector 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 val, CoulombType type = CoulombType::Slater); + void set_ci(std::vector val); + void set_spinorbit(double val, SpinOrbType type = SpinOrbType::Zeta); // Getters std::array 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 eigensystem(); - std::vector magnetisation(std::vector H, std::vector Hdir, double T, MagUnits type); - std::vector susceptibility(std::vector T, std::vector Hdir, MagUnits type); RowMatrixXcd zeeman_hamiltonian(double H, // Calculates the Zeeman Hamiltonian for applied std::vector Hdir); // field H in direction Hdir - std::vector< std::vector > calculate_moments(RowMatrixXcd ev); - std::vector calculate_boltzmann( // Calculates the Boltzmann factor exp(-E/kT) - VectorXd en, double T); + std::vector< RowMatrixXcd > calculate_moments_matrix(RowMatrixXcd ev); std::vector get_states(); }; // class ic1ion diff --git a/src/include/ic_states.hpp b/src/include/ic_states.hpp index cb2c647..d8b15f2 100644 --- a/src/include/ic_states.hpp +++ b/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) diff --git a/src/include/physprop.hpp b/src/include/physprop.hpp new file mode 100644 index 0000000..37e1246 --- /dev/null +++ b/src/include/physprop.hpp @@ -0,0 +1,69 @@ +/* physprop.hpp + * + * Header file for the physical properties base class + * + * (c) 2024 Duc Le - duc.le@stfc.ac.uk + * This program is licensed under the GNU General Purpose License, version 3. Please see the COPYING file + */ + +#pragma once + +#include "eigen.hpp" +#include +#include + +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 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 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 Hdir) = 0; + virtual std::tuple eigensystem() = 0; + virtual std::vector calculate_moments_matrix(RowMatrixXcd ev) = 0; + VectorXd calculate_boltzmann(VectorXd en, double T); + VectorXd heatcapacity(std::vector Tvec); + RowMatrixXd magnetisation(std::vector H, std::vector Hdir, std::vector Tvec, MagUnits type); + VectorXd susceptibility(std::vector T, std::vector Hdir, MagUnits type); +}; + +// Mapping for Python binding to map string to enum +static const std::unordered_map mag_unit_names = { + {"bohr", physprop::MagUnits::bohr}, {"cgs", physprop::MagUnits::cgs}, {"SI", physprop::MagUnits::SI} }; + +} // namespace libMcPhase diff --git a/src/libmcphase/pycf1ion.cpp b/src/libmcphase/pycf1ion.cpp index a0aeb1e..1113d7d 100644 --- a/src/libmcphase/pycf1ion.cpp +++ b/src/libmcphase/pycf1ion.cpp @@ -9,14 +9,13 @@ #include "cf1ion.hpp" #include #include +#include #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(cls), args, kwargs); @@ -31,8 +30,16 @@ void wrap_cf1ion(py::module &m) { .def(py::init(), py::arg("ionname")) .def(py::init(), 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 H, std::vector Hdir, std::vector 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 T, std::vector 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'")); }); } diff --git a/src/libmcphase/pycfpars.cpp b/src/libmcphase/pycfpars.cpp index c40bb9a..1589757 100644 --- a/src/libmcphase/pycfpars.cpp +++ b/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 #include "pycfpars.hpp" namespace py = pybind11; @@ -29,25 +27,28 @@ static const std::unordered_map 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(); - cls->set_J(J); + std::string ionname = args[0].cast(); + 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(); - cls->set_name(ionname); + double J = args[0].cast(); + 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()); } else if (kwargs.contains("ionname")) { diff --git a/src/libmcphase/pycfpars.hpp b/src/libmcphase/pycfpars.hpp index 5cc225c..76eb46c 100644 --- a/src/libmcphase/pycfpars.hpp +++ b/src/libmcphase/pycfpars.hpp @@ -8,6 +8,11 @@ #pragma once +#include "cfpars.hpp" +#include +namespace py = pybind11; +using namespace libMcPhase; + template T set_enum(std::string key, std::unordered_map 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); diff --git a/src/libmcphase/pyic1ion.cpp b/src/libmcphase/pyic1ion.cpp index e3dafed..42fdfd9 100644 --- a/src/libmcphase/pyic1ion.cpp +++ b/src/libmcphase/pyic1ion.cpp @@ -6,8 +6,8 @@ * This program is licensed under the GNU General Purpose License, version 3. Please see the LICENSE file */ -#include "cfpars.hpp" #include "ic1ion.hpp" +#include #include #include #include @@ -16,20 +16,16 @@ namespace py = pybind11; using namespace libMcPhase; -static const std::unordered_map mag_unit_names = { - {"bohr", cfpars::MagUnits::bohr}, {"cgs", cfpars::MagUnits::cgs}, {"SI", cfpars::MagUnits::SI} }; - static const std::unordered_map coulomb_names = { {"Slater", ic1ion::CoulombType::Slater}, {"CondonShortley", ic1ion::CoulombType::CondonShortley}, {"Racah", ic1ion::CoulombType::Racah} }; static const std::unordered_map spinorb_names = { {"Zeta", ic1ion::SpinOrbType::Zeta}, {"Lambda", ic1ion::SpinOrbType::Lambda} }; -void cf_parse(cfpars *cls, py::args args, py::kwargs kwargs); - ic1ion *ic1ion_init(py::args args, py::kwargs kwargs) { ic1ion *cls = new ic1ion; - cf_parse(static_cast(cls), args, kwargs); + bool is_ic1ion = true; + cf_parse(static_cast(cls), args, kwargs, is_ic1ion); if (kwargs.contains("zeta")) { cls->set_spinorbit(kwargs["zeta"].cast(), ic1ion::SpinOrbType::Zeta); } @@ -56,12 +52,13 @@ void wrap_ic1ion(py::module &m) { .def("get_ci", &ic1ion::get_ci) .def("hamiltonian", &ic1ion::hamiltonian, "the crystal field Hamiltonian") .def("eigensystem", &ic1ion::eigensystem, "the eigenvectors and eigenvalues of the crystal field Hamiltonian") + .def_property("GJ", [](ic1ion const &self) { return self.get_GJ(); }, [](ic1ion &self, double v) { self.set_GJ(v); }) .def_property("zeta", [](ic1ion const &self) { return self.get_spinorbit(); }, [](ic1ion &self, double v) { self.set_spinorbit(v, ic1ion::SpinOrbType::Zeta); }) .def_property("slater", [](ic1ion const &self) { return self.get_coulomb(); }, [](ic1ion &self, std::vector v) { self.set_coulomb(v, ic1ion::CoulombType::Slater); }) .def("zeeman_hamiltonian", &ic1ion::zeeman_hamiltonian, "the Zeeman Hamiltonian") .def("calculate_boltzmann", &ic1ion::calculate_boltzmann, "") - .def("calculate_moments", &ic1ion::calculate_moments, "") - .def("magnetisation", [](ic1ion &self, std::vector H, std::vector Hdir, double T, std::string unit) { return self.magnetisation(H, Hdir, T, + .def("heatcapacity", &ic1ion::heatcapacity, "the heat capacity of the crystal field Hamiltonian in J/mol/K") + .def("magnetisation", [](ic1ion &self, std::vector H, std::vector Hdir, std::vector 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", [](ic1ion &self, std::vector T, std::vector 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'")); }) diff --git a/src/singleion/cf1ion.cpp b/src/singleion/cf1ion.cpp index 831e83a..8a8d4b6 100644 --- a/src/singleion/cf1ion.cpp +++ b/src/singleion/cf1ion.cpp @@ -31,6 +31,7 @@ void cf1ion::set(int l, int m, double val) { void cf1ion::set_unit(cfpars::Units const newunit) { cfpars::set_unit(newunit); + physprop::m_meVconv = m_econv; m_ham_calc = false; m_ev_calc = false; } @@ -45,19 +46,33 @@ void cf1ion::set_name(const std::string &ionname) { cfpars::set_name(ionname); m_ham_calc = false; m_ev_calc = false; + m_magops_calc = false; } void cf1ion::set_J(const double J) { cfpars::set_J(J); m_ham_calc = false; m_ev_calc = false; + m_magops_calc = false; } // --------------------------------------------------------------------------------------------------------------- // // General methods for the cf1ion class // --------------------------------------------------------------------------------------------------------------- // -RowMatrixXcd cf1ion::hamiltonian(bool upper) { +void cf1ion::fill_upper() { + int dimj = m_J2 + 1; + for (int i=0; i cf1ion::eigensystem() { if (!m_ev_calc) { - SelfAdjointEigenSolver es(hamiltonian(false)); + SelfAdjointEigenSolver es(_hamiltonian(false)); m_eigenvectors = es.eigenvectors(); m_eigenvalues = es.eigenvalues(); m_ev_calc = true; @@ -176,4 +192,66 @@ std::tuple cf1ion::eigensystem() { return std::tuple(m_eigenvectors, m_eigenvalues); } +void cf1ion::calc_mag_ops() { + // Calculates the magnetic operators Jx, Jy, Jz + if (m_magops_calc) + return; + int dimj = m_J2 + 1; + m_magops = std::vector(3, RowMatrixXcd::Zero(dimj, dimj)); + // RM1 = (1/2) * sqrt( factorial(2*J+1+1) / factorial(2*J-1) ); + double rm1 = sqrt( m_racah.f(m_J2 + 2) / m_racah.f(m_J2 - 1) ) / 2.; + double rm1_sq2 = rm1 / sqrt(2.); + for (size_t i=0; i(0., Om + Op); // Jy + } + } + // Fill the upper triangle + for (int i=0; i Hdir) { + if (m_GJ < 0) { + throw std::runtime_error("Lande g-factor not defined. Either set the ion name or set GJ."); + } + if (Hdir.size() != 3) { + throw std::runtime_error("cf1ion::zeeman_hamiltonian(): Hdir must be a 3-vector"); + } + // Normalise the field direction vector + double Hnorm = sqrt(Hdir[0] * Hdir[0] + Hdir[1] * Hdir[1] + Hdir[2] * Hdir[2]); + if (fabs(Hnorm) < 1.e-6) { + throw std::runtime_error("cf1ion::zeeman_hamiltonian(): Direction vector cannot be zero"); + } + std::vector nHdir; + std::transform(Hdir.begin(), Hdir.end(), std::back_inserter(nHdir), [Hnorm](double Hd){ return Hd / Hnorm; }); + // Calculates the Jx, Jy, Jz operators if not done already. + calc_mag_ops(); + RowMatrixXcd zeeman = RowMatrixXcd::Zero(m_J2+1, m_J2+1); + zeeman = (m_magops[0] * nHdir[0]) + (m_magops[1] * nHdir[1]) + (m_magops[2] * nHdir[2]); + double H_in_energy = m_GJ * H * MU_B * m_econv; // want H in user requested external energy units + return (zeeman * H_in_energy); +} + +std::vector cf1ion::calculate_moments_matrix(RowMatrixXcd ev) { + calc_mag_ops(); + std::vector moments; + for (size_t ii=0; ii<3; ii++) { + moments.push_back((ev.adjoint()) * (m_magops[ii] * ev) * m_GJ); + } + return moments; +} + } // namespace libMcPhase diff --git a/src/singleion/cfpars.cpp b/src/singleion/cfpars.cpp index 7d8ea7a..d345943 100644 --- a/src/singleion/cfpars.cpp +++ b/src/singleion/cfpars.cpp @@ -10,9 +10,6 @@ namespace libMcPhase { -static const double K_B = 0.08617343183; // meV/K - Boltzmann constant -static const double MU_B = 0.0578838263; // meV/T - Bohr magneton - // --------------------------------------------------------------------------------------------------------------- // // Reference tables (values taken from program cfield, by Peter Fabi, FZ Juelich, file theta.c) // --------------------------------------------------------------------------------------------------------------- // @@ -136,9 +133,6 @@ const Map3 &RKTABLE() { return rk_table; } -// Conversion factors for different energy units[from][to], order: [meV, cm, K]. -static const std::array ENERGYCONV = { {1., 8.065544005, 11.6045221} }; - // --------------------------------------------------------------------------------------------------------------- // // General methods for cfpars class // --------------------------------------------------------------------------------------------------------------- // @@ -179,7 +173,7 @@ void cfpars::set(int l, int m, double val) { m_Bi[id] = val / m_convfact[id] / m_econv; } -const double cfpars::get(int l, int m) const { +double cfpars::get(int l, int m) const { switch(l) { case 2: return m_Bo[2 + m]; case 4: return m_Bo[9 + m]; @@ -250,6 +244,7 @@ void cfpars::set_name(const std::string &ionname) { } m_stevfact = {alpha, beta, gamma}; m_convertible = true; + m_GJ = GJ[m_n]; // Now reset the conversion table (from internal to external parameters) this->set_type(m_type); } diff --git a/src/singleion/ic1ion.cpp b/src/singleion/ic1ion.cpp index 4e8cd71..8622845 100644 --- a/src/singleion/ic1ion.cpp +++ b/src/singleion/ic1ion.cpp @@ -25,25 +25,6 @@ static const double F625 = 0.016104; // Ratios of F_6/F_2 slater integrals // Conversion factors for different energy units[from][to], order: [meV, cm, K]. static const std::array ICENERGYCONV = {0.1239841973, 1., 1.4387773587}; -// 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 std::array MAGCONV = {1., NAMUB*1000, NAMUB}; - -// Note these constants are strange because the default energy unit in this module is cm-1 -static const double NAMUBSQ_ERG = 0.26074098; // N_A * muB[erg/G] * muB[cm/T] * [Tesla_to_Gauss=1e-4] -static const double NAMUBSQ_JOULE = 3.276568e-06; // N_A * muB[J/T] * muB[cm/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 SUSCCONV = {MU_B, NAMUBSQ_ERG, NAMUBSQ_JOULE}; - -// EPSILON to determine if energy levels are degenerate or not -static const double DELTA_EPS = 1e-6; - // Helper vectors for indexing into CF parameters array static const std::array, 12> idq = { {{2,2,0,4}, {2,1,1,3}, {4,4,5,13}, {4,3,6,12}, {4,2,7,11}, {4,1,8,10}, {6,6,14,26}, {6,5,15,25}, {6,4,16,24}, {6,3,17,23}, {6,2,18,22}, {6,1,19,21}} }; @@ -90,6 +71,7 @@ void ic1ion::set_unit(cfpars::Units const newunit) { if (ii < 3) m_alpha[ii] = m_alpha_i[ii] * m_econv; } + physprop::m_meVconv = ENERGYCONV[(int)m_unit]; } void ic1ion::set_type(const cfpars::Type newtype) { @@ -637,120 +619,6 @@ RowMatrixXcd ic1ion::hamiltonian() { return m_hamiltonian; } -// --------------------------------------------------------------------------------------------------------------- // -// Calculates bulk properties (magnetisation, susceptibility) -// --------------------------------------------------------------------------------------------------------------- // -std::vector ic1ion::calculate_boltzmann(VectorXd en, double T) -{ - std::vector boltzmann, en_meV; - // Need kBT in external energy units. K_B is in meV/K - double kBT = K_B * T * m_econv; - double Emin = std::numeric_limits::max(); - for (size_t i=0; i < (size_t)en.size(); i++) { - Emin = (en(i) < Emin) ? en(i) : Emin; - } - for (size_t i=0; i < (size_t)en.size(); i++) { - const double expi = exp(-(en(i) - Emin) / kBT); - boltzmann.push_back((fabs(expi) > DELTA_EPS) ? expi : 0.); - } - return boltzmann; -} - -std::vector ic1ion::magnetisation(std::vector Hvec, std::vector Hdir, double T, MagUnits unit_type) -{ - // Normalise the field direction vector - double Hnorm = sqrt(Hdir[0] * Hdir[0] + Hdir[1] * Hdir[1] + Hdir[2] * Hdir[2]); - if (fabs(Hnorm) < 1.e-6) { - throw std::runtime_error("ic1ion::magnetisation(): Direction vector cannot be zero"); - } - std::vector nHdir; - std::transform(Hdir.begin(), Hdir.end(), std::back_inserter(nHdir), [Hnorm](double Hd){ return Hd / Hnorm; }); - // Calculates Magnetisation M(H) at specified T - if (!m_ham_calc) - calculate_hamiltonian(); - std::vector M; - M.reserve(Hvec.size()); - // Loops through all the input field magnitudes and calculates the magnetisation - for (auto H: Hvec) { - if (unit_type == MagUnits::cgs) { - H /= 1e4; // For cgs, input field is in Gauss, need to convert to Tesla for Zeeman calculation - } - RowMatrixXcd ham = m_hamiltonian - zeeman_hamiltonian(H, Hdir); - SelfAdjointEigenSolver es(ham); - // calculate_moments returns a vector of 3 moments *squared* vectors, in the x, y, z directions - std::vector< std::vector > moments_vec = calculate_moments(es.eigenvectors()); - std::vector boltzmann = calculate_boltzmann(es.eigenvalues(), T); - std::vector Mdir; - for (auto moments: moments_vec) { - double Mexp = 0., Z = 0.; - //std::inner_product(moments.begin(), moments.end(), boltzmann.begin(), Mexp); - //std::accumulate(boltzmann.begin(), boltzmann.end(), Z); - for (int ii=0; ii ic1ion::susceptibility(std::vector Tvec, std::vector Hdir, MagUnits unit_type) -{ - // Normalise the field direction vector - double Hnorm = sqrt(Hdir[0] * Hdir[0] + Hdir[1] * Hdir[1] + Hdir[2] * Hdir[2]); - if (fabs(Hnorm) < 1.e-6) { - throw std::runtime_error("ic1ion::magnetisation(): Direction vector cannot be zero"); - } - std::vector nHdir; - std::transform(Hdir.begin(), Hdir.end(), std::back_inserter(nHdir), [Hnorm](double Hd){ return Hd / Hnorm; }); - // Calculates the susceptibility chi(T) - if (!m_ev_calc) - calculate_eigensystem(); - std::vector chi; - chi.reserve(Tvec.size()); - // Calculates the moments matrices in the x, y, z directions, and get the resultant - std::vector moments_mat_vec = calculate_moments_matrix(m_eigenvectors); - RowMatrixXcd moments_mat = moments_mat_vec[0] * nHdir[0] - + moments_mat_vec[1] * nHdir[1] - + moments_mat_vec[2] * nHdir[2]; - // Now calculate the first and second order terms in the Van Vleck equation - size_t nlev = m_eigenvectors.cols(); - std::vector mu(nlev, 0.); - std::vector mu2(nlev, 0.); - for (size_t ii=0; ii --- ] - // chi(T) = --- > [ ------------ - 2 > ------------ ] exp(-E/k_BT) - // Z --- [ k_B T --- En - Em ] - // n m!=n - - for (auto T: Tvec) { - std::vector boltzmann = calculate_boltzmann(m_eigenvalues, T); - const double beta = 1. / (K_B * T); - double U = 0., Z = 0.; - for (size_t ii=0; ii0) { id = (L-1) + ( L>3 ? 1 : 0 ) + ( L>5 ? 1 : 0 ) + ( L>6 ? 1 : 0 ) + ( L>7 ? 1 : 0 ); } else { id = (L==-3 ? 3 : 0) + (L!=-3 ? -2*L-4 : 0 ); } @@ -151,16 +151,16 @@ double racah_chi(orbital Lp, orbital L, qG2 U, qG2 Up) 0, 0, 0, 0, 0, 1672., 0, // N 10 13 9 4 0, 0, 0, 0, 0, 220., 0};// O 11 14 10 4 - if(U.isequal("10")) { retval = t[id*7]; } - else if(U.isequal("11")) { retval = t[id*7+1]; } - else if(U.isequal("20")) { retval = t[id*7+2]; } - else if(U.isequal("21")) { retval = t[id*7+3]; } - else if(U.isequal("30")) { retval = t[id*7+4]; } - else if(U.isequal("31")) { + if(U.isequal(10)) { retval = t[id*7]; } + else if(U.isequal(11)) { retval = t[id*7+1]; } + else if(U.isequal(20)) { retval = t[id*7+2]; } + else if(U.isequal(21)) { retval = t[id*7+3]; } + else if(U.isequal(30)) { retval = t[id*7+4]; } + else if(U.isequal(31)) { if(L==Lp) { retval = t[id*7+5]; } else { retval = t[id*7+6]; } } } - else if(Up.isequal("40")) // Use table VIc + else if(Up.isequal(40)) // Use table VIc { #define CD(c,d) case c: id = d; break switch (L) { CD(0,0); CD(2,1); CD(3,2); CD(4,3); CD(-4,4); CD(5,5); CD(6,6); @@ -182,7 +182,7 @@ double racah_chi(orbital Lp, orbital L, qG2 U, qG2 Up) 0, 0, 0, 0, 528., 0, // N 10 12 0, 0, 0, 0, 22., 0};// Q 12 13 - if(U.isequal("40")) { if(L==Lp) { retval = t[id*6+4]; } else { retval = t[id*6+5]; } } + if(U.isequal(40)) { if(L==Lp) { retval = t[id*6+4]; } else { retval = t[id*6+5]; } } else { retval = t[id*6+U.u1]; } } else // Use table VIa @@ -204,21 +204,21 @@ double racah_chi(orbital Lp, orbital L, qG2 U, qG2 Up) 260., 0, -25., 0, 94., 104., -181., 0, -36., 0, 40.}; // 22 22 if(L>=0) { - if(Up.isequal("20") && U.isequal("20")) { retval = t[L]; } - else if(Up.isequal("21")) { - if(U.isequal("11")) { retval = t[L+11]; } - else if(U.isequal("20")) { retval = t[L+22]; } - else if(U.isequal("21")) { retval = 0; } } - else if(Up.isequal("30")) { - if(U.isequal("10")) { retval = t[L+55]; } - else if(U.isequal("11")) { retval = t[L+66]; } - else if(U.isequal("20")) { retval = t[L+77]; } - else if(U.isequal("21")) { retval = t[L+88]; } - else if(U.isequal("30")) { retval = t[L+99]; } } - else if(Up.isequal("22")) { - if(U.isequal("20")) { retval = t[L+110]; } - else if(U.isequal("21")) { retval = t[L+121]; } - else if(U.isequal("22")) { retval = t[L+132]; } + if(Up.isequal(20) && U.isequal(20)) { retval = t[L]; } + else if(Up.isequal(21)) { + if(U.isequal(11)) { retval = t[L+11]; } + else if(U.isequal(20)) { retval = t[L+22]; } + else if(U.isequal(21)) { retval = 0; } } + else if(Up.isequal(30)) { + if(U.isequal(10)) { retval = t[L+55]; } + else if(U.isequal(11)) { retval = t[L+66]; } + else if(U.isequal(20)) { retval = t[L+77]; } + else if(U.isequal(21)) { retval = t[L+88]; } + else if(U.isequal(30)) { retval = t[L+99]; } } + else if(Up.isequal(22)) { + if(U.isequal(20)) { retval = t[L+110]; } + else if(U.isequal(21)) { retval = t[L+121]; } + else if(U.isequal(22)) { retval = t[L+132]; } }} } return retval; @@ -251,12 +251,12 @@ double racah_e2prod(qR7 W, qG2 U, qG2 Up, orbital L, orbital Lp) { double x1=0,x2=0,c1,c2; double retval = 0; - if(U.isequal("21") && Up.isequal("21")) + if(U.isequal(21) && Up.isequal(21)) { - if(W.isequal("210")) { x1 = (3./7); x2 = 0.; } - else if(W.isequal("211")) { x1 = (4./7); x2 = 3.; } - else if(W.isequal("220")) { x1 = (-6./7); x2 = -3.; } - else if(W.isequal("221")) { x1 = (-1./7); x2 = (12./11); } + if(W.isequal(210)) { x1 = (3./7); x2 = 0.; } + else if(W.isequal(211)) { x1 = (4./7); x2 = 3.; } + else if(W.isequal(220)) { x1 = (-6./7); x2 = -3.; } + else if(W.isequal(221)) { x1 = (-1./7); x2 = (12./11); } switch(L) { case 0: c1= 0; c2= 0; break; @@ -319,7 +319,7 @@ double racah_yfn(int n, int v, int S2, qG2 U, int vp, qG2 Up) switch(n) { - case 2: if(S2==0 && vp==2 && v==2 && U.isequal("20") && Up.isequal("20")) { y = 2; } break; + case 2: if(S2==0 && vp==2 && v==2 && U.isequal(20) && Up.isequal(20)) { y = 2; } break; case 3: if(S2==1 && vp==3) { // Table XV from Racah IV. // Up = (11) (20) (21) U @@ -348,8 +348,8 @@ double racah_yfn(int n, int v, int S2, qG2 U, int vp, qG2 Up) 0, 0, 0, 0, 0, -sqrt(7./80), 0, 1./4}; // 4 0 (22) if(S2==2) { if(v==2) { id = U.u2*8 + 2*(Up.u1-1)+Up.u2; y = t[id]; } else if(v==4) { id = (2*(U.u1-1)+U.u2+2)*8 + 2*(Up.u1-1)+Up.u2; y = t[id]; } } - else if(S2==0) { if(v==0 && U.isequal("00")) { y = t[Up.u2+61]; } - else if(v==2 && U.isequal("20")) { y = t[Up.u2+69]; } + else if(S2==0) { if(v==0 && U.isequal(00)) { y = t[Up.u2+61]; } + else if(v==2 && U.isequal(20)) { y = t[Up.u2+69]; } else if(v==4) { id = (U.u2+9)*8 + Up.u2+5; y = t[id]; } } } break; case 5: if(vp==5) @@ -448,7 +448,7 @@ double racah_phi(qG2 U, qG2 Up, orbital Lp, orbital L) double retval=0; int id; - if(Up.isequal("31")) + if(Up.isequal(31)) { // Table XIVb from Racah IV. U' = (31) // U = (10) (11) (21) (30) (31) (31)' L Lnum id double t[]={0,sqrt(330.), 0, 17*sqrt(143.), 209., 0, // P 1 0 @@ -468,13 +468,13 @@ double racah_phi(qG2 U, qG2 Up, orbital Lp, orbital L) 0, 0, 0, 0, 352., 0}; // O 11 14 switch (L) { CD(1,0); CD(2,1); CD(3,2); CD(-3,3); CD(4,4); CD(5,5); CD(-5,6); CD(6,7); CD(-6,8); CD(7,9); CD(-7,10); CD(8,11); CD(9,12); CD(10,13); CD(11,14); default: return retval; } - if(U.isequal("31")) { if(L==Lp) { retval = t[id*6+4]; } else { retval = t[id*6+5]; } } - else if(U.isequal("30")) { retval = t[id*6+3]; } - else if(U.isequal("21")) { retval = t[id*6+2]; } - else if(U.isequal("11")) { retval = t[id*6+1]; } - else if(U.isequal("10")) { retval = t[id*6]; } + if(U.isequal(31)) { if(L==Lp) { retval = t[id*6+4]; } else { retval = t[id*6+5]; } } + else if(U.isequal(30)) { retval = t[id*6+3]; } + else if(U.isequal(21)) { retval = t[id*6+2]; } + else if(U.isequal(11)) { retval = t[id*6+1]; } + else if(U.isequal(10)) { retval = t[id*6]; } } - else if(Up.isequal("40")) + else if(Up.isequal(40)) { // Table XIVc from Racah IV. Up=(40) // U = (20) (21) (22) L Lnum id double t[] = { 0, 0, 2*sqrt(2145.), // S 0 0 @@ -491,9 +491,9 @@ double racah_phi(qG2 U, qG2 Up, orbital Lp, orbital L) 0, -84*sqrt(19./31), sqrt(2530.)}; // N 10 11 switch (L) { CD(0,0); CD(2,1); CD(3,2); CD(4,3); CD(-4,4); CD(5,5); CD(6,6); CD(-6,7); CD(7,8); CD(8,9); CD(-8,10); CD(10,11); default: return retval; } - if(U.isequal("20")) { retval = t[id*3]; } - else if(U.isequal("21")) { retval = t[id*3+1]; } - else if(U.isequal("22")) { retval = t[id*3+2]; } + if(U.isequal(20)) { retval = t[id*3]; } + else if(U.isequal(21)) { retval = t[id*3+1]; } + else if(U.isequal(22)) { retval = t[id*3+2]; } } else if(L>=0) { // Table XIVa from Racah IV. @@ -513,31 +513,31 @@ double racah_phi(qG2 U, qG2 Up, orbital Lp, orbital L) 1., 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // (00) (22) 0, 0,3*sqrt(429.), 0, 4*sqrt(65.), 0, 3*sqrt(85.), 0, 0, 0, 0, // (20) (22) 144., 0, 69., 0, -148., 72., 39., 0, -96., 0, 56.}; // (22) (22) - if(U.isequal("11") && Up.isequal("11")) { retval = t[L]; } - else if(Up.isequal("20")) + if(U.isequal(11) && Up.isequal(11)) { retval = t[L]; } + else if(Up.isequal(20)) { - if (U.isequal("20")) { retval = t[L+11]; } - else if(U.isequal("21")) { retval = t[L+22]; } + if (U.isequal(20)) { retval = t[L+11]; } + else if(U.isequal(21)) { retval = t[L+22]; } } - else if(Up.isequal("10") && U.isequal("21")) { retval = t[L+33]; } - else if(Up.isequal("21")) + else if(Up.isequal(10) && U.isequal(21)) { retval = t[L+33]; } + else if(Up.isequal(21)) { - if (U.isequal("10")) { retval = t[L+33]; } - else if(U.isequal("20")) { retval = t[L+44]; } - else if(U.isequal("21")) { retval = t[L+55]; } + if (U.isequal(10)) { retval = t[L+33]; } + else if(U.isequal(20)) { retval = t[L+44]; } + else if(U.isequal(21)) { retval = t[L+55]; } } - else if(Up.isequal("30")) + else if(Up.isequal(30)) { - if (U.isequal("11")) { retval = t[L+66]; } - else if(U.isequal("20")) { retval = t[L+77]; } - else if(U.isequal("21")) { retval = t[L+88]; } - else if(U.isequal("30")) { retval = t[L+99]; } + if (U.isequal(11)) { retval = t[L+66]; } + else if(U.isequal(20)) { retval = t[L+77]; } + else if(U.isequal(21)) { retval = t[L+88]; } + else if(U.isequal(30)) { retval = t[L+99]; } } - else if(Up.isequal("22")) + else if(Up.isequal(22)) { - if (U.isequal("00")) { retval = t[L+110]; } - else if(U.isequal("20")) { retval = t[L+121]; } - else if(U.isequal("22")) { retval = t[L+132]; } + if (U.isequal(00)) { retval = t[L+110]; } + else if(U.isequal(20)) { retval = t[L+121]; } + else if(U.isequal(22)) { retval = t[L+132]; } } } return retval; @@ -844,7 +844,7 @@ RowMatrixXd racah_ci(int n, double alpha) // For //Up.set(1,1); std::vector xwu2 = racah_xwu(W,U,Up); //xwu = xwu2; std::cout << xwu[0] << ", " << xwu[1] << "\n"; - //qG2 U(0,0); std::cout << U.isequal("00") << "\n"; + //qG2 U(0,0); std::cout << U.isequal(00) << "\n"; //qG2 Up(2,2); //std::cout << racah_phi(U,Up,S,S) << "\n"; //std::vector< std::vector > e2 = racah_e2(2); diff --git a/src/singleion/ic_states.cpp b/src/singleion/ic_states.cpp index e69c9d8..78622fc 100644 --- a/src/singleion/ic_states.cpp +++ b/src/singleion/ic_states.cpp @@ -116,27 +116,27 @@ int qR7::set(int w1_, int w2_, int w3_) bool qR7::operator == (const qR7 & wp) const { - if (w1==wp.w1 && w2==wp.w2 && w3==wp.w3) { return true; } - else { return false; } + return (w1==wp.w1 && w2==wp.w2 && w3==wp.w3); } bool qR7::operator != (const qR7 & wp) const { - if (w1!=wp.w1 || w2!=wp.w2 || w3!=wp.w3) { return true; } - else { return false; } + return (w1!=wp.w1 || w2!=wp.w2 || w3!=wp.w3); } -bool qR7::isequal (const char * Wstr) +bool qR7::isequal (int w) { - char ws=0; - int w1n=0,w2n=1,w3n=2; - - ws = Wstr[0]; w1n = atoi(&ws); - ws = Wstr[1]; w2n = atoi(&ws); - ws = Wstr[2]; w3n = atoi(&ws); - - if (w1==w1n && w2==w2n && w3==w3n) { return true; } - else { return false; } + // In C++ a literal beginning with 0 is octal + if (w > 99) { // Decimal literal + std::div_t w1d = std::div(w, 100); + std::div_t w2d = std::div(w1d.rem, 10); + return (w1==w1d.quot && w2==w2d.quot && w3==w2d.rem); + } else if (w > 0) { // Octal literal, w1=0 + std::div_t w2d = std::div(w, 8); + return (w1==0 && w2==w2d.quot && w3==w2d.rem); + } else { + return (w1==0 && w2==0 && w3==0); + } } // --------------------------------------------------------------------------------------------------------------- // @@ -149,26 +149,24 @@ int qG2::set(int u1_, int u2_) bool qG2::operator == (const qG2 & up) const { - if (u1==up.u1 && u2==up.u2) { return true; } - else { return false; } + return (u1==up.u1 && u2==up.u2); } bool qG2::operator != (const qG2 & up) const { - if (u1!=up.u1 || u2!=up.u2) { return true; } - else { return false; } + return (u1!=up.u1 || u2!=up.u2); } -bool qG2::isequal (const char * Ustr) +bool qG2::isequal (int u) { - char us=0; - int u1n=0,u2n=1; - - us = Ustr[0]; u1n = atoi(&us); - us = Ustr[1]; u2n = atoi(&us); - - if (u1==u1n && u2==u2n) { return true; } - else { return false; } + if (u > 9) { // Decimal literal + std::div_t u1d = std::div(u, 10); + return (u1==u1d.quot && u2==u1d.rem); + } else if (u > 0) { // Octal literal + return (u1==0 && u2==u); + } else { + return (u1==0 && u2==0); + } } // --------------------------------------------------------------------------------------------------------------- // @@ -861,7 +859,7 @@ void fconf::set(int n, bool mJflag, orbital l) << ", U=[" << conf.states[0].U.u1 << " " << conf.states[0].U.u2 << "], W=[" << conf.states[0].W.w1 << " " << conf.states[0].W.w2 << " " << conf.states[0].W.w3 << "]\n"; - WisW = conf.states[0].W.isequal("100"); + WisW = conf.states[0].W.isequal(100); //fprintf(stdout,"W=[%i %i %i] == [1 0 0] is %s\n",conf.states[0].W.w1,conf.states[0].W.w2,conf.states[0].W.w3, // WisW ? "true" : "false"); std::cout << "W=[" << conf.states[0].W.w1 << " " << conf.states[0].W.w2 << " " << conf.states[0].W.w3 diff --git a/src/singleion/ic_tensorops.cpp b/src/singleion/ic_tensorops.cpp index 258c12d..a87d286 100644 --- a/src/singleion/ic_tensorops.cpp +++ b/src/singleion/ic_tensorops.cpp @@ -22,7 +22,7 @@ static const std::array q = {0,0,0,0,0,0,-2,-1,0,1,2,-3,-2,-1,0,1,2,3,- // Calculates the tensor operators up to // --------------------------------------------------------------------------------------------------------------- // void ic1ion::calc_tensorops(int num_op) { - if (num_op < m_tensorops.size() && num_op <= 0) + if (num_op < m_tensorops.size() || num_op <= 0) return; // Calculates the L and S operator matrix for the each directions // Note operators up to 6 are the magnetic moment operators in x, y, z basis @@ -33,6 +33,7 @@ void ic1ion::calc_tensorops(int num_op) { racah_mumat(-1, Lm1, Sm1); // Note also that m_tensorops is a vector of _real_ matrices. // BUT negative order (-q) operators are purely imaginary, as are the Sy and Ly matrices. + m_tensorops.clear(); m_tensorops.push_back((Sm1 - Sp1) / sqrt(2.)); // Sx m_tensorops.push_back((Lm1 - Lp1) / sqrt(2.)); // Lx m_tensorops.push_back((Sm1 + Sp1) / sqrt(2.)); // Sy @@ -56,7 +57,7 @@ RowMatrixXcd ic1ion::zeeman_hamiltonian(double H, std::vector Hdir) { // Normalise the field direction vector double Hnorm = sqrt(Hdir[0] * Hdir[0] + Hdir[1] * Hdir[1] + Hdir[2] * Hdir[2]); if (fabs(Hnorm) < 1.e-6) { - throw std::runtime_error("ic1ion::magnetisation(): Direction vector cannot be zero"); + throw std::runtime_error("ic1ion::zeeman_hamiltonian(): Direction vector cannot be zero"); } std::vector nHdir; std::transform(Hdir.begin(), Hdir.end(), std::back_inserter(nHdir), [Hnorm](double Hd){ return Hd / Hnorm; }); @@ -66,7 +67,7 @@ RowMatrixXcd ic1ion::zeeman_hamiltonian(double H, std::vector Hdir) { zeeman.real() = (GS * m_tensorops[0] + m_tensorops[1]) * nHdir[0] // gSx + Lx +(GS * m_tensorops[4] + m_tensorops[5]) * nHdir[2]; // gSz + Lz zeeman.imag() = (GS * m_tensorops[2] + m_tensorops[3]) * nHdir[1]; // gSy + Ly - double H_in_energy = H * MU_B * m_econv; // MU_B in internal units (cm/T); want H in external energy units + double H_in_energy = H * MU_Bc * m_econv; // MU_B in internal units (cm/T); want H in external energy units return (zeeman * H_in_energy); } @@ -85,20 +86,6 @@ std::vector ic1ion::calculate_moments_matrix(RowMatrixXcd ev) { return moments; } -std::vector< std::vector > ic1ion::calculate_moments(RowMatrixXcd ev) { - std::vector mommat = calculate_moments_matrix(ev); - std::vector< std::vector > moments_diag; - for (auto mom: mommat) { - std::vector moment_vec; - for (int ii=0; ii::max(); + for (size_t i=0; i < (size_t)en.size(); i++) { + Emin = (en(i) < Emin) ? en(i) : Emin; + } + for (size_t i=0; i < (size_t)en.size(); i++) { + const double expi = exp(-(en(i) - Emin) * beta); + boltzmann(i) = (fabs(expi) > DELTA_EPS) ? expi : 0.; + } + return boltzmann; +} + +VectorXd physprop::heatcapacity(std::vector Tvec) { + auto es = eigensystem(); + VectorXd out = VectorXd::Zero(Tvec.size()); + size_t sz = std::get<1>(es).size(); + VectorXd en = VectorXd::Zero(sz); + double Emin = std::numeric_limits::max(); + for (size_t i=0; i < sz; i++) { + en(i) = std::get<1>(es)(i) / m_meVconv; + Emin = (en[i] < Emin) ? en[i] : Emin; + } + for (size_t i=0; i < sz; i++) { + en[i] -= Emin; + } + for (size_t tt=0; tt(es), T); + for (size_t i=0; i < sz; i++) { + Z += expfact[i]; + U += en[i] * expfact[i]; + U2 += en[i] * en[i] * expfact[i]; + } + U /= Z; + U2 /= Z; + out(tt) = ((U2 - U * U) / (K_B * T * T)) * NAMEV; + } + return out; +} + +RowMatrixXd physprop::magnetisation(std::vector Hvec, std::vector Hdir, std::vector Tvec, MagUnits unit_type) +{ + // Normalise the field direction vector + double Hnorm = sqrt(Hdir[0] * Hdir[0] + Hdir[1] * Hdir[1] + Hdir[2] * Hdir[2]); + if (fabs(Hnorm) < 1.e-6) { + throw std::runtime_error("physprop::magnetisation(): Direction vector cannot be zero"); + } + std::vector nHdir; + std::transform(Hdir.begin(), Hdir.end(), std::back_inserter(nHdir), [Hnorm](double Hd){ return Hd / Hnorm; }); + // Calculates Magnetisation M(H) at specified T + RowMatrixXcd ham0 = hamiltonian(); + RowMatrixXd M = RowMatrixXd::Zero(Hvec.size(), Tvec.size()); + // Loops through all the input field magnitudes and calculates the magnetisation + std::vector mag_ops = calculate_moments_matrix(RowMatrixXcd::Identity(ham0.rows(), ham0.cols())); + RowMatrixXcd Jmat = nHdir[0] * mag_ops[0] + nHdir[1] * mag_ops[1] + nHdir[2] * mag_ops[2]; + for (size_t hh=0; hh es(ham); + for (size_t tt=0; tt Tvec, std::vector Hdir, MagUnits unit_type) +{ + // Normalise the field direction vector + double Hnorm = sqrt(Hdir[0] * Hdir[0] + Hdir[1] * Hdir[1] + Hdir[2] * Hdir[2]); + if (fabs(Hnorm) < 1.e-6) { + throw std::runtime_error("physprop::magnetisation(): Direction vector cannot be zero"); + } + std::vector nHdir; + std::transform(Hdir.begin(), Hdir.end(), std::back_inserter(nHdir), [Hnorm](double Hd){ return Hd / Hnorm; }); + // Calculates the susceptibility chi(T) + std::tuple es = eigensystem(); + VectorXd chi = VectorXd::Zero(Tvec.size()); + // Calculates the moments matrices in the x, y, z directions, and get the resultant + std::vector moments_mat_vec = calculate_moments_matrix(std::get<0>(es)); + RowMatrixXcd moments_mat = moments_mat_vec[0] * nHdir[0] + + moments_mat_vec[1] * nHdir[1] + + moments_mat_vec[2] * nHdir[2]; + // Now calculate the first and second order terms in the Van Vleck equation + size_t nlev = std::get<0>(es).cols(); + std::vector mu(nlev, 0.); + std::vector mu2(nlev, 0.); + std::vector eigenvalues(nlev, 0.); + for (size_t ii=0; ii(es)[ii] / m_meVconv; } // Convert energies to meV + for (size_t ii=0; ii --- ] + // chi(T) = --- > [ ------------ - 2 > ------------ ] exp(-E/k_BT) + // Z --- [ k_B T --- En - Em ] + // n m!=n + + for (size_t tt=0; tt(es), Tvec[tt]); + const double beta = 1. / (K_B * Tvec[tt]); + double U = 0., Z = 0.; + for (size_t ii=0; ii