Mantid
Loading...
Searching...
No Matches
CrystalFieldPeaksBase.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 +
7#include "MantidCurveFitting/Functions/CrystalFieldPeaksBase.h"
8#include "MantidCurveFitting/Functions/CrystalElectricField.h"
9#include "MantidKernel/Exception.h"
10
11#include <cctype>
12#include <functional>
13#include <map>
14
15namespace Mantid::CurveFitting::Functions {
16
17namespace {
18
19// Maps ion name to its int code.
20const std::map<std::string, int> ION_2_NRE{{"Ce", 1}, {"Pr", 2}, {"Nd", 3}, {"Pm", 4}, {"Sm", 5},
21 {"Eu", 6}, {"Gd", 7}, {"Tb", 8}, {"Dy", 9}, {"Ho", 10},
22 {"Er", 11}, {"Tm", 12}, {"Yb", 13}};
23
24const bool REAL_PARAM_PART = true;
25const bool IMAG_PARAM_PART = false;
26
27void fixx(API::IFunction &fun, const std::string &par) {
28 fun.setParameter(par, 0.0);
29 fun.fixParameter(par);
30 const std::string ipar = "I" + par;
31 fun.setParameter(ipar, 0.0);
32 fun.fixParameter(ipar);
33}
34
35void free(API::IFunction &fun, const std::string &par, bool realOnly) {
36 fun.unfixParameter(par);
37 const std::string ipar = "I" + par;
38 if (realOnly) {
39 fun.setParameter(ipar, 0.0);
40 fun.fixParameter(ipar);
41 } else {
42 fun.unfixParameter(ipar);
43 }
44}
45
46// Set symmetry C1 or Ci
47void setSymmetryC1(API::IFunction &fun) {
48 fun.clearTies();
49 auto i = fun.parameterIndex("B20");
50 for (; i < fun.nParams(); ++i) {
51 fun.unfix(i);
52 }
53}
54
55// Set symmetry C2, Cs or C2h
56void setSymmetryC2(API::IFunction &fun) {
57 fun.clearTies();
58 fun.unfixParameter("B20");
59 fun.unfixParameter("B40");
60 fun.unfixParameter("B60");
61 fixx(fun, "B21");
62 free(fun, "B22", REAL_PARAM_PART);
63 fixx(fun, "B41");
64 free(fun, "B42", IMAG_PARAM_PART);
65 fixx(fun, "B43");
66 free(fun, "B44", IMAG_PARAM_PART);
67 fixx(fun, "B61");
68 free(fun, "B62", IMAG_PARAM_PART);
69 fixx(fun, "B63");
70 free(fun, "B64", IMAG_PARAM_PART);
71 fixx(fun, "B65");
72 free(fun, "B66", IMAG_PARAM_PART);
73}
74
75// Set symmetry C2v, D2 or D2h
76void setSymmetryC2v(API::IFunction &fun) {
77 fun.clearTies();
78 fun.unfixParameter("B20");
79 fun.unfixParameter("B40");
80 fun.unfixParameter("B60");
81 fixx(fun, "B21");
82 free(fun, "B22", REAL_PARAM_PART);
83 fixx(fun, "B41");
84 free(fun, "B42", REAL_PARAM_PART);
85 fixx(fun, "B43");
86 free(fun, "B44", REAL_PARAM_PART);
87 fixx(fun, "B61");
88 free(fun, "B62", REAL_PARAM_PART);
89 fixx(fun, "B63");
90 free(fun, "B64", REAL_PARAM_PART);
91 fixx(fun, "B65");
92 free(fun, "B66", REAL_PARAM_PART);
93}
94
95// Set symmetry C4, S4 or C4h
96void setSymmetryC4(API::IFunction &fun) {
97 fun.clearTies();
98 fun.unfixParameter("B20");
99 fun.unfixParameter("B40");
100 fun.unfixParameter("B60");
101 fixx(fun, "B21");
102 fixx(fun, "B22");
103 fixx(fun, "B41");
104 fixx(fun, "B42");
105 fixx(fun, "B43");
106 free(fun, "B44", REAL_PARAM_PART);
107 fixx(fun, "B61");
108 fixx(fun, "B62");
109 fixx(fun, "B63");
110 free(fun, "B64", IMAG_PARAM_PART);
111 fixx(fun, "B65");
112 fixx(fun, "B66");
113}
114
115// Set symmetry D4, C4v, D2d or D4h
116void setSymmetryD4(API::IFunction &fun) {
117 fun.clearTies();
118 fun.unfixParameter("B20");
119 fun.unfixParameter("B40");
120 fun.unfixParameter("B60");
121 fixx(fun, "B21");
122 fixx(fun, "B22");
123 fixx(fun, "B41");
124 fixx(fun, "B42");
125 fixx(fun, "B43");
126 free(fun, "B44", REAL_PARAM_PART);
127 fixx(fun, "B61");
128 fixx(fun, "B62");
129 fixx(fun, "B63");
130 free(fun, "B64", REAL_PARAM_PART);
131 fixx(fun, "B65");
132 fixx(fun, "B66");
133}
134
135// Set symmetry C3 or S6
136void setSymmetryC3(API::IFunction &fun) {
137 fun.clearTies();
138 fun.unfixParameter("B20");
139 fun.unfixParameter("B40");
140 fun.unfixParameter("B60");
141 fixx(fun, "B21");
142 fixx(fun, "B22");
143 fixx(fun, "B41");
144 fixx(fun, "B42");
145 free(fun, "B43", REAL_PARAM_PART);
146 fixx(fun, "B44");
147 fixx(fun, "B61");
148 fixx(fun, "B62");
149 free(fun, "B63", IMAG_PARAM_PART);
150 fixx(fun, "B64");
151 fixx(fun, "B65");
152 free(fun, "B66", IMAG_PARAM_PART);
153}
154
155// Set symmetry D3, C3v or D3d
156void setSymmetryD3(API::IFunction &fun) {
157 fun.clearTies();
158 fun.unfixParameter("B20");
159 fun.unfixParameter("B40");
160 fun.unfixParameter("B60");
161 fixx(fun, "B21");
162 fixx(fun, "B22");
163 fixx(fun, "B41");
164 fixx(fun, "B42");
165 free(fun, "B43", REAL_PARAM_PART);
166 fixx(fun, "B44");
167 fixx(fun, "B61");
168 fixx(fun, "B62");
169 free(fun, "B63", REAL_PARAM_PART);
170 fixx(fun, "B64");
171 fixx(fun, "B65");
172 free(fun, "B66", REAL_PARAM_PART);
173}
174
175// Set symmetry C6, C3h, C6h, D6, C6v, D3h, or D6h
176void setSymmetryC6(API::IFunction &fun) {
177 fun.clearTies();
178 fun.unfixParameter("B20");
179 fun.unfixParameter("B40");
180 fun.unfixParameter("B60");
181 fixx(fun, "B21");
182 fixx(fun, "B22");
183 fixx(fun, "B41");
184 fixx(fun, "B42");
185 fixx(fun, "B43");
186 fixx(fun, "B44");
187 fixx(fun, "B61");
188 fixx(fun, "B62");
189 fixx(fun, "B63");
190 fixx(fun, "B64");
191 fixx(fun, "B65");
192 free(fun, "B66", REAL_PARAM_PART);
193}
194
195// Set symmetry T, Td, Th, O, or Oh
196void setSymmetryT(API::IFunction &fun) {
197 fun.clearTies();
198 fun.setParameter("B20", 0.0);
199 fun.fixParameter("B20");
200 fun.unfixParameter("B40");
201 fun.unfixParameter("B60");
202 fixx(fun, "B21");
203 fixx(fun, "B22");
204 fixx(fun, "B41");
205 fixx(fun, "B42");
206 fixx(fun, "B43");
207 free(fun, "B44", REAL_PARAM_PART);
208 fixx(fun, "B61");
209 fixx(fun, "B62");
210 fixx(fun, "B63");
211 free(fun, "B64", REAL_PARAM_PART);
212 fixx(fun, "B65");
213 fixx(fun, "B66");
214 fun.tie("B44", "5*B40");
215 fun.tie("B64", "-21*B60");
216}
217
219const std::map<std::string, std::function<void(API::IFunction &)>> SYMMETRY_MAP{
220 // Set symmetry C1 or Ci
221 {"C1", setSymmetryC1},
222 {"Ci", setSymmetryC1},
223 // Set symmetry C2, Cs or C2h
224 {"C2", setSymmetryC2},
225 {"Cs", setSymmetryC2},
226 {"C2h", setSymmetryC2},
227 // Set symmetry C2v, D2 or D2h
228 {"C2v", setSymmetryC2v},
229 {"D2", setSymmetryC2v},
230 {"D2h", setSymmetryC2v},
231 // Set symmetry C4, S4 or C4h
232 {"C4", setSymmetryC4},
233 {"S4", setSymmetryC4},
234 {"C4h", setSymmetryC4},
235 // Set symmetry D4, C4v, D2d or D4h
236 {"D4", setSymmetryD4},
237 {"C4v", setSymmetryD4},
238 {"D2d", setSymmetryD4},
239 {"D4h", setSymmetryD4},
240 // Set symmetry C3 or S6
241 {"C3", setSymmetryC3},
242 {"S6", setSymmetryC3},
243 // Set symmetry D3, C3v or D3d
244 {"D3", setSymmetryD3},
245 {"C3v", setSymmetryD3},
246 {"D3d", setSymmetryD3},
247 // Set symmetry C6, C3h, C6h, D6, C6v, D3h, or D6h
248 {"C6", setSymmetryC6},
249 {"C3h", setSymmetryC6},
250 {"C6h", setSymmetryC6},
251 {"D6", setSymmetryC6},
252 {"C6v", setSymmetryC6},
253 {"D3h", setSymmetryC6},
254 {"D6h", setSymmetryC6},
255 // Set symmetry T, Td, Th, O, or Oh
256 {"T", setSymmetryT},
257 {"Td", setSymmetryT},
258 {"Th", setSymmetryT},
259 {"O", setSymmetryT},
260 {"Oh", setSymmetryT}};
261
262} // anonymous namespace
263
265CrystalFieldPeaksBase::CrystalFieldPeaksBase() : API::ParamFunction(), m_defaultDomainSize(0) {
266
267 declareAttribute("Ion", Attribute("Ce"));
268 declareAttribute("Symmetry", Attribute("Ci"));
269 declareAttribute("ToleranceEnergy", Attribute(1.0e-10));
270 declareAttribute("ToleranceIntensity", Attribute(1.0e-1));
271 declareAttribute("MaxPeakCount", Attribute(0));
272
273 declareParameter("BmolX", 0.0, "The x-component of the molecular field.");
274 declareParameter("BmolY", 0.0, "The y-component of the molecular field.");
275 declareParameter("BmolZ", 0.0, "The z-component of the molecular field.");
276
277 declareParameter("BextX", 0.0, "The x-component of the external field.");
278 declareParameter("BextY", 0.0, "The y-component of the external field.");
279 declareParameter("BextZ", 0.0, "The z-component of the external field.");
280
281 declareParameter("B20", 0.0, "Real part of the B20 field parameter.");
282 declareParameter("B21", 0.0, "Real part of the B21 field parameter.");
283 declareParameter("B22", 0.0, "Real part of the B22 field parameter.");
284 declareParameter("B40", 0.0, "Real part of the B40 field parameter.");
285 declareParameter("B41", 0.0, "Real part of the B41 field parameter.");
286 declareParameter("B42", 0.0, "Real part of the B42 field parameter.");
287 declareParameter("B43", 0.0, "Real part of the B43 field parameter.");
288 declareParameter("B44", 0.0, "Real part of the B44 field parameter.");
289 declareParameter("B60", 0.0, "Real part of the B60 field parameter.");
290 declareParameter("B61", 0.0, "Real part of the B61 field parameter.");
291 declareParameter("B62", 0.0, "Real part of the B62 field parameter.");
292 declareParameter("B63", 0.0, "Real part of the B63 field parameter.");
293 declareParameter("B64", 0.0, "Real part of the B64 field parameter.");
294 declareParameter("B65", 0.0, "Real part of the B65 field parameter.");
295 declareParameter("B66", 0.0, "Real part of the B66 field parameter.");
296
297 declareParameter("IB21", 0.0, "Imaginary part of the B21 field parameter.");
298 declareParameter("IB22", 0.0, "Imaginary part of the B22 field parameter.");
299 declareParameter("IB41", 0.0, "Imaginary part of the B41 field parameter.");
300 declareParameter("IB42", 0.0, "Imaginary part of the B42 field parameter.");
301 declareParameter("IB43", 0.0, "Imaginary part of the B43 field parameter.");
302 declareParameter("IB44", 0.0, "Imaginary part of the B44 field parameter.");
303 declareParameter("IB61", 0.0, "Imaginary part of the B61 field parameter.");
304 declareParameter("IB62", 0.0, "Imaginary part of the B62 field parameter.");
305 declareParameter("IB63", 0.0, "Imaginary part of the B63 field parameter.");
306 declareParameter("IB64", 0.0, "Imaginary part of the B64 field parameter.");
307 declareParameter("IB65", 0.0, "Imaginary part of the B65 field parameter.");
308 declareParameter("IB66", 0.0, "Imaginary part of the B66 field parameter.");
309
310 setSymmetryC1(*this);
311}
312
319void CrystalFieldPeaksBase::calculateEigenSystem(DoubleFortranVector &en, ComplexFortranMatrix &wf,
320 ComplexFortranMatrix &ham, ComplexFortranMatrix &hz, int &nre) const {
321
322 auto ion = getAttribute("Ion").asString();
323 if (ion.empty()) {
324 throw std::runtime_error("Ion name must be specified.");
325 }
326
327 auto ionIter = ION_2_NRE.find(ion);
328 if (ionIter == ION_2_NRE.end()) {
329 // If 'Ion=S2', or 'Ion=J2.5' etc, interpret as arbitrary J values with gJ=2
330 // Allow lower case, but values must be half-integral. E.g. 'Ion=S2.4' fails
331 switch (ion[0]) {
332 case 'S':
333 case 's':
334 case 'J':
335 case 'j': {
336 if (ion.size() > 1 && std::isdigit(static_cast<unsigned char>(ion[1]))) {
337 // Need to store as 2J to allow half-integer values
338 try {
339 auto J2 = std::stof(ion.substr(1)) * 2.;
340 if (J2 > 99.) {
341 throw std::out_of_range("");
342 }
343 if (fabs(J2 - (int)J2) < 0.001) {
344 nre = -(int)J2;
345 break;
346 }
347 // Catch exceptions thrown by stof so we get a more meaningful error
348 } catch (const std::invalid_argument &) {
349 throw std::runtime_error("Invalid value '" + ion.substr(1) + "' of J passed to CrystalFieldPeaks.");
350 } catch (const std::out_of_range &) {
351 throw std::runtime_error("Value of J: '" + ion.substr(1) + "' passed to CrystalFieldPeaks is too big.");
352 }
353 }
354 // fall through
355 }
356 default:
357 throw std::runtime_error("Unknown ion name '" + ion + "' passed to CrystalFieldPeaks.");
358 }
359 } else {
360 nre = ionIter->second;
361 }
362
363 DoubleFortranVector bmol(1, 3);
364 bmol(1) = getParameter("BmolX");
365 bmol(2) = getParameter("BmolY");
366 bmol(3) = getParameter("BmolZ");
367
368 // For CrystalFieldSusceptibility and CrystalFieldMagnetisation we need
369 // to be able to override the external field set here, since in these
370 // measurements, a different external field is applied.
371 DoubleFortranVector bext(1, 3);
372 bext(1) = getParameter("BextX");
373 bext(2) = getParameter("BextY");
374 bext(3) = getParameter("BextZ");
375
376 double B20 = getParameter("B20");
377 double B21 = getParameter("B21");
378 double B22 = getParameter("B22");
379 double B40 = getParameter("B40");
380 double B41 = getParameter("B41");
381 double B42 = getParameter("B42");
382 double B43 = getParameter("B43");
383 double B44 = getParameter("B44");
384 double B60 = getParameter("B60");
385 double B61 = getParameter("B61");
386 double B62 = getParameter("B62");
387 double B63 = getParameter("B63");
388 double B64 = getParameter("B64");
389 double B65 = getParameter("B65");
390 double B66 = getParameter("B66");
391
392 double IB21 = getParameter("IB21");
393 double IB22 = getParameter("IB22");
394 double IB41 = getParameter("IB41");
395 double IB42 = getParameter("IB42");
396 double IB43 = getParameter("IB43");
397 double IB44 = getParameter("IB44");
398 double IB61 = getParameter("IB61");
399 double IB62 = getParameter("IB62");
400 double IB63 = getParameter("IB63");
401 double IB64 = getParameter("IB64");
402 double IB65 = getParameter("IB65");
403 double IB66 = getParameter("IB66");
404
405 ComplexFortranMatrix bkq(0, 6, 0, 6);
406 bkq(2, 0) = ComplexType(B20, 0.0);
407 bkq(2, 1) = ComplexType(B21, IB21);
408 bkq(2, 2) = ComplexType(B22, IB22);
409 bkq(4, 0) = ComplexType(B40, 0.0);
410 bkq(4, 1) = ComplexType(B41, IB41);
411 bkq(4, 2) = ComplexType(B42, IB42);
412 bkq(4, 3) = ComplexType(B43, IB43);
413 bkq(4, 4) = ComplexType(B44, IB44);
414 bkq(6, 0) = ComplexType(B60, 0.0);
415 bkq(6, 1) = ComplexType(B61, IB61);
416 bkq(6, 2) = ComplexType(B62, IB62);
417 bkq(6, 3) = ComplexType(B63, IB63);
418 bkq(6, 4) = ComplexType(B64, IB64);
419 bkq(6, 5) = ComplexType(B65, IB65);
420 bkq(6, 6) = ComplexType(B66, IB66);
421
422 calculateEigensystem(en, wf, ham, hz, nre, bmol, bext, bkq);
423 // MaxPeakCount is a read-only "mutable" attribute.
424 const_cast<CrystalFieldPeaksBase *>(this)->setAttributeValue("MaxPeakCount", static_cast<int>(en.size()));
425}
426
428void CrystalFieldPeaksBase::setAttribute(const std::string &name, const IFunction::Attribute &attr) {
429 if (name == "Symmetry") {
430 auto symmIter = SYMMETRY_MAP.find(attr.asString());
431 if (symmIter == SYMMETRY_MAP.end()) {
432 throw std::runtime_error("Unknown symmetry passed to CrystalFieldPeaks: " + attr.asString());
433 }
434 symmIter->second(*this);
435 }
436 IFunction::setAttribute(name, attr);
437}
438
439std::string CrystalFieldPeaksBaseImpl::name() const { return "CrystalFieldPeaksBaseImpl"; }
440
441void CrystalFieldPeaksBaseImpl::function(const API::FunctionDomain & /*domain*/,
442 API::FunctionValues & /*values*/) const {
443 throw Kernel::Exception::NotImplementedError("Method is intentionally not implemented.");
444}
445
446} // namespace Mantid::CurveFitting::Functions
name
std::string name
Definition Run.cpp:60
fabs
#define fabs(x)
Definition Matrix.cpp:22
Mantid::API::FunctionDomain
Base class that represents the domain of a function.
Definition FunctionDomain.h:31
Mantid::API::FunctionValues
A class to store values calculated by a function.
Definition FunctionValues.h:25
Mantid::API::IFunction::Attribute
Attribute is a non-fitting parameter.
Definition IFunction.h:285
Mantid::API::IFunction::Attribute::asString
std::string asString() const
Returns string value if attribute is a string, throws exception otherwise.
Definition IFunction.cpp:735
Mantid::API::IFunction::getAttribute
virtual Attribute getAttribute(const std::string &name) const
Return a value of attribute attName.
Definition IFunction.cpp:1469
Mantid::API::IFunction::declareAttribute
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
Definition IFunction.cpp:1493
Mantid::API::IFunction::setAttribute
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
Definition IFunction.cpp:1484
Mantid::API::IFunction::name
virtual std::string name() const =0
Returns the function's name.
Mantid::API::IFunction::setAttributeValue
void setAttributeValue(const std::string &attName, const T &value)
Set an attribute value.
Definition IFunction.h:601
Mantid::API::ParamFunction
Implements the part of IFunction interface dealing with parameters.
Definition ParamFunction.h:33
Mantid::API::ParamFunction::declareParameter
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
Definition ParamFunction.cpp:239
Mantid::API::ParamFunction::getParameter
double getParameter(size_t i) const override
Get i-th parameter.
Definition ParamFunction.cpp:62
Mantid::CurveFitting::EigenVector::size
size_t size() const
Size of the vector.
Definition EigenVector.cpp:108
Mantid::CurveFitting::Functions::CrystalFieldPeaksBaseImpl::function
void function(const API::FunctionDomain &, API::FunctionValues &) const override
Evaluates the function for all arguments in the domain.
Definition CrystalFieldPeaksBase.cpp:441
Mantid::CurveFitting::Functions::CrystalFieldPeaksBaseImpl::name
std::string name() const override
Returns the function's name.
Definition CrystalFieldPeaksBase.cpp:439
Mantid::CurveFitting::Functions::CrystalFieldPeaksBase
CrystalFieldPeaks is a function that calculates crystal field peak positions and intensities.
Definition CrystalFieldPeaksBase.h:22
Mantid::CurveFitting::Functions::CrystalFieldPeaksBase::calculateEigenSystem
void calculateEigenSystem(DoubleFortranVector &en, ComplexFortranMatrix &wf, ComplexFortranMatrix &ham, ComplexFortranMatrix &hz, int &nre) const
Calculate the crystal field eigensystem.
Definition CrystalFieldPeaksBase.cpp:319
Mantid::CurveFitting::Functions::CrystalFieldPeaksBase::setAttribute
void setAttribute(const std::string &name, const Attribute &) override
Perform a castom action when an attribute is set.
Definition CrystalFieldPeaksBase.cpp:428
Mantid::Kernel::Exception::NotImplementedError
Marks code as not implemented yet.
Definition Exception.h:138
Mantid::CurveFitting::Functions::calculateEigensystem
void MANTID_CURVEFITTING_DLL calculateEigensystem(DoubleFortranVector &eigenvalues, ComplexFortranMatrix &eigenvectors, ComplexFortranMatrix &hamiltonian, ComplexFortranMatrix &hzeeman, int nre, const DoubleFortranVector &bmol, const DoubleFortranVector &bext, const ComplexFortranMatrix &bkq, double alpha_euler=0.0, double beta_euler=0.0, double gamma_euler=0.0)
Calculate eigenvalues and eigenvectors of the crystal field hamiltonian.
Definition CrystalElectricField.cpp:633
Mantid::CurveFitting::ComplexType
std::complex< double > ComplexType
Definition EigenComplexVector.h:19
Generated by doxygen 1.9.8