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