Mantid
Loading...
Searching...
No Matches
CrystalFieldSusceptibility.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
13#include "MantidAPI/Jacobian.h"
19#include <boost/algorithm/string/predicate.hpp>
20#include <cmath>
21
23
24namespace {
25
26// Does the actual calculation of the susceptibility
27void calculate(double *out, const double *xValues, const size_t nData, const DoubleFortranVector &en,
28 const ComplexFortranMatrix &wf, const int nre, const std::vector<double> &H, const double convfact) {
29 // Some constants
30 const double k_B = PhysicalConstants::BoltzmannConstant; // in meV/K
31 const double eps = 1.e-6; // for degeneracy calculations
32 // First calculates the matrix elements of the magnetic moment operator
34 calculateMagneticMomentMatrix(wf, H, nre, mumat);
35 int nlevels = en.len();
36 DoubleFortranVector mu(nlevels);
37 DoubleFortranVector mu2(nlevels);
38 mu.zero();
39 mu2.zero();
40 for (auto i = 1; i <= nlevels; i++) {
41 for (auto j = 1; j <= nlevels; j++) {
42 const double den = en(i) - en(j);
43 if (fabs(den) < eps) { // First order term
44 mu(i) += real(mumat(i, j) * conj(mumat(i, j)));
45 } else { // Second order term
46 mu2(i) += real(mumat(i, j) * conj(mumat(i, j))) / den;
47 }
48 }
49 }
50 // Now calculates the temperature dependence.
51 mu *= convfact;
52 mu2 *= convfact;
53 for (size_t iT = 0; iT < nData; iT++) {
54 double Z = 0.;
55 double U = 0.;
56 const double beta = 1 / (k_B * xValues[iT]);
57 for (auto i = 1; i <= nlevels; i++) {
58 double expfact = exp(-beta * en(i));
59 Z += expfact;
60 U += ((mu(i) * beta) - (2 * mu2(i))) * expfact;
61 }
62 out[iT] = U / Z;
63 }
64}
65
66// Calculate powder average - Mpowder = (Mx + My + Mz)/3
67void calculate_powder(double *out, const double *xValues, const size_t nData, const DoubleFortranVector &en,
68 const ComplexFortranMatrix &wf, const int nre, const double convfact) {
69 for (size_t j = 0; j < nData; j++) {
70 out[j] = 0.;
71 }
72 // Loop over the x, y, z directions
73 std::vector<double> H;
74 std::vector<double> tmp(nData, 0.);
75 for (int i = 0; i < 3; i++) {
76 H.assign(3, 0.);
77 H[i] = 1.;
78 calculate(&tmp[0], xValues, nData, en, wf, nre, H, convfact);
79 for (size_t j = 0; j < nData; j++) {
80 out[j] += tmp[j];
81 }
82 }
83 for (size_t j = 0; j < nData; j++) {
84 out[j] /= 3.;
85 }
86}
87} // namespace
88
90 // Declare attributes inherited by subclasses
91 declareAttribute("Hdir", Attribute(std::vector<double>{0., 0., 1.}));
92 declareAttribute("Unit", Attribute("cgs"));
93 declareAttribute("inverse", Attribute(false));
94 declareAttribute("powder", Attribute(false));
95 declareAttribute("ScaleFactor", Attribute(1.0)); // Only for multi-site use
96 // Cannot declare fitting parameters as this is not a ParamFunction
97}
98
99void CrystalFieldSusceptibilityBase::function1D(double *out, const double *xValues, const size_t nData) const {
100 auto H = getAttribute("Hdir").asVector();
101 if (H.size() != 3) {
102 throw std::invalid_argument("Hdir must be a three-element vector.");
103 }
104 auto powder = getAttribute("powder").asBool();
105 double Hnorm = sqrt(H[0] * H[0] + H[1] * H[1] + H[2] * H[2]);
106 if (fabs(Hnorm) > 1.e-6) {
107 for (auto i = 0; i < 3; i++) {
108 H[i] /= Hnorm;
109 }
110 }
111 // Get the unit conversion factor.
112 auto unit = getAttribute("Unit").asString();
113 const double NAMUB2cgs = 0.03232776; // N_A * muB(erg/G) * muB(meV/G)
114 // The constant above is in strange units because we need the output
115 // to be in erg/G^2/mol==cm^3/mol, but in the calculation all energies
116 // are in meV and the expression for chi has a energy denominator.
117 const double NAMUB2si = 4.062426e-7; // N_A * muB(J/T) * muB(meV/T) * mu0
118 // Again, for SI units, we need to have one of the muB in meV not J.
119 // The additional factor of mu0 is due to the different definitions of
120 // the magnetisation, B- and H-fields in the SI and cgs systems.
121 double convfact = boost::iequals(unit, "bohr") ? 0.057883818 : (boost::iequals(unit, "SI") ? NAMUB2si : NAMUB2cgs);
122 // Note the constant for "bohr" is the bohr magneton in meV/T, this will
123 // give the susceptibility in "atomic" units of uB/T/ion.
124 // Note that chi_SI = (4pi*10^-6)chi_cgs
125 // Default unit is cgs (cm^3/mol).
126 // Use stored values
127 if (powder) {
128 calculate_powder(out, xValues, nData, m_en, m_wf, m_nre, convfact);
129 } else {
130 calculate(out, xValues, nData, m_en, m_wf, m_nre, H, convfact);
131 }
132 const double lambda = getParameter("Lambda");
133 const double EPS = 1.e-6;
134 if (fabs(lambda) > EPS) {
135 for (size_t i = 0; i < nData; i++) {
136 out[i] /= (1. - lambda * out[i]); // chi = chi_cf/(1 - lambda.chi_cf)
137 }
138 }
139 const double chi0 = getParameter("Chi0");
140 if (fabs(lambda) > 1.e-6) {
141 for (size_t i = 0; i < nData; i++) {
142 out[i] += chi0;
143 }
144 }
145 if (getAttribute("inverse").asBool()) {
146 for (size_t i = 0; i < nData; i++) {
147 out[i] = 1. / out[i];
148 }
149 }
150 auto fact = getAttribute("ScaleFactor").asDouble();
151 if (fact != 1.0) {
152 for (size_t i = 0; i < nData; i++) {
153 out[i] *= fact;
154 }
155 }
156}
157
159
161 : CrystalFieldPeaksBase(), CrystalFieldSusceptibilityBase(), m_setDirect(false) {
162 declareParameter("Lambda", 0.0, "Effective exchange interaction");
163 declareParameter("Chi0", 0.0, "Background or remnant susceptibility");
164}
165
166// Sets the eigenvectors / values directly
168 const int nre) {
169 m_en = en;
170 m_wf = wf;
171 m_nre = nre;
172 m_setDirect = true;
173}
174
175void CrystalFieldSusceptibility::function1D(double *out, const double *xValues, const size_t nData) const {
176 if (!m_setDirect) {
178 }
180}
181
183 declareParameter("Lambda", 0.0, "Effective exchange interaction");
184 declareParameter("Chi0", 0.0, "Background or remnant susceptibility");
185}
186
187// Sets the eigenvectors / values directly
189 const ComplexFortranMatrix &wf, const int nre) {
190 m_en = en;
191 m_wf = wf;
192 m_nre = nre;
193}
194
195} // namespace Mantid::CurveFitting::Functions
const std::vector< double > * lambda
gsl_vector * tmp
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
#define fabs(x)
Definition: Matrix.cpp:22
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
std::vector< double > asVector() const
Returns vector<double> if attribute is vector<double>, throws exception otherwise.
Definition: IFunction.cpp:765
std::string asString() const
Returns string value if attribute is a string, throws exception otherwise.
Definition: IFunction.cpp:660
double asDouble() const
Returns double value if attribute is a double, throws exception otherwise.
Definition: IFunction.cpp:739
bool asBool() const
Returns bool value if attribute is a bool, throws exception otherwise.
Definition: IFunction.cpp:752
virtual Attribute getAttribute(const std::string &name) const
Return a value of attribute attName.
Definition: IFunction.cpp:1394
virtual double getParameter(size_t i) const =0
Get i-th parameter.
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
Definition: IFunction.cpp:1418
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
CrystalFieldPeaks is a function that calculates crystal field peak positions and intensities.
void calculateEigenSystem(DoubleFortranVector &en, ComplexFortranMatrix &wf, ComplexFortranMatrix &ham, ComplexFortranMatrix &hz, int &nre) const
Calculate the crystal field eigensystem.
CrystalFieldSusceptibility is a function that calculates the molar magnetic susceptibility (in cm^3/m...
void function1D(double *out, const double *xValues, const size_t nData) const override
Function you want to fit to.
void setEigensystem(const DoubleFortranVector &en, const ComplexFortranMatrix &wf, const int nre)
void function1D(double *out, const double *xValues, const size_t nData) const override
Function you want to fit to.
void setEigensystem(const DoubleFortranVector &en, const ComplexFortranMatrix &wf, const int nre)
void MANTID_CURVEFITTING_DLL calculateMagneticMomentMatrix(const ComplexFortranMatrix &ev, const std::vector< double > &Hdir, const int nre, ComplexFortranMatrix &mumat)
Calculate the full magnetic moment matrix in a particular eigenvector basis.
FortranMatrix< ComplexMatrix > ComplexFortranMatrix
FortranVector< EigenVector > DoubleFortranVector
MANTID_GEOMETRY_DLL Raster calculate(const Kernel::V3D &beamDirection, const IObject &shape, const double cubeSizeInMetre)
Definition: Rasterize.cpp:181
static constexpr double BoltzmannConstant
Boltzmann Constant in meV/K Taken from http://physics.nist.gov/cuu/Constants on 10/07/2012.