Mantid
Loading...
Searching...
No Matches
CrystalElectricField.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//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
12
13#include <algorithm>
14#include <array>
15#include <cmath>
16
18
19namespace {
20
21// The missing braces warning is a false positive -
22// https://llvm.org/bugs/show_bug.cgi?id=21629
23GNU_DIAG_OFF("missing-braces")
24
25// number of rare earth ions (trivalent rare earths with unfilled f-shell)
26const int maxNre = 13;
27
28// define some rare earth constants (ggj is the Lande g-factor, ddimj=2J+1)
29const std::array<double, maxNre> ggj = {6.0 / 7., 4.0 / 5., 8.0 / 11., 3.0 / 5., 2.0 / 7., 0.0, 2.0,
30 3.0 / 2., 4.0 / 3., 5.0 / 4., 6.0 / 5., 7.0 / 6., 8.0 / 7.};
31
32const std::array<double, maxNre> ddimj = {6.0, 9.0, 10.0, 9.0, 6.0, 1.0, 8.0, 13.0, 16.0, 17.0, 16.0, 13.0, 8.0};
33//---------------------------
34// set some natural constants
35//---------------------------
36const double pi = 4.0 * atan(1.0);
37const double kb = 1.38062; // x 10**(-23) J/K, Boltzmann constant k_B
38const double hh = 6.626075540; // x 10**(-34) J*sec, Planks constant h
39const double hq = hh / 2 / pi;
40const double ee = 1.6021773349; // x 10**(-19) Coulomb, electric charge
41const double me = 9.109389754; // x 10**(-31) kg, electron mass
42// within the above choose, the factor which connects
43// 1meV and 1 Kelvin is fixed and given by:
44// 10*ee/kb =: fmevkelvin = 11.6047...
45// this means 1 meV is nearly 11.6 K
46const double c_fmevkelvin = 10 * ee / kb;
47// magneton of Bohr in meV per tesla
48const double c_myb = hq / me / 2;
49
50//--------------------------------
51// define the delta function
52//--------------------------------
53double delta(double mj, double nj, double j) {
54 double res = 0.0;
55 if (mj == nj && fabs(mj) <= j && fabs(nj) <= j) {
56 res = 1.0;
57 }
58 return res;
59}
60
61//--------------------
62// jp(n) = <n+1|j+|n>
63//--------------------
64double jp(double nj, double j) {
65 if (fabs(nj) <= j) {
66 return sqrt(j * (j + 1) - nj * (nj + 1));
67 } else {
68 return 0.0;
69 }
70}
71
72//--------------------
73// jm(n) = <n-1|j-|n>
74//--------------------
75double jm(double nj, double j) {
76 if (fabs(nj) <= j) {
77 return sqrt(j * (j + 1) - nj * (nj - 1));
78 } else {
79 return 0.0;
80 }
81}
82
83//-----------------------
84// jp2(n) = <n+2|j+^2|n>
85//-----------------------
86double jp2(double nj, double j) { return jp(nj + 1, j) * jp(nj, j); }
87
88//-----------------------
89// jm2(n) = <n-2|j-^2|n>
90//-----------------------
91double jm2(double nj, double j) { return jm(nj - 1, j) * jm(nj, j); }
92
93//-----------------------
94// jp3(n) = <n+3|j+^3|n>
95//-----------------------
96double jp3(double nj, double j) { return jp(nj + 2, j) * jp(nj + 1, j) * jp(nj, j); }
97
98//-----------------------
99// jm3(n) = <n-3|j-^3|n>
100//-----------------------
101double jm3(double nj, double j) { return jm(nj - 2, j) * jm(nj - 1, j) * jm(nj, j); }
102
103//-----------------------
104// jp4(n) = <n+4|j+^4|n>
105//-----------------------
106double jp4(double nj, double j) { return jp(nj + 3, j) * jp(nj + 2, j) * jp(nj + 1, j) * jp(nj, j); }
107
108//-----------------------
109// jm4(n) = <n-4|j-^4|n>
110//-----------------------
111double jm4(double nj, double j) { return jm(nj - 3, j) * jm(nj - 2, j) * jm(nj - 1, j) * jm(nj, j); }
112
113//-----------------------
114// jp5(n) = <n+5|j+^5|n>
115//-----------------------
116double jp5(double nj, double j) { return jp(nj + 4, j) * jp(nj + 3, j) * jp(nj + 2, j) * jp(nj + 1, j) * jp(nj, j); }
117
118//-----------------------
119// jm5(n) = <n-5|j-^5|n>
120//-----------------------
121double jm5(double nj, double j) { return jm(nj - 4, j) * jm(nj - 3, j) * jm(nj - 2, j) * jm(nj - 1, j) * jm(nj, j); }
122
123//-----------------------
124// jp6(n) = <n+6|j+^6|n>
125//-----------------------
126double jp6(double nj, double j) {
127 return jp(nj + 5, j) * jp(nj + 4, j) * jp(nj + 3, j) * jp(nj + 2, j) * jp(nj + 1, j) * jp(nj, j);
128}
129
130//-----------------------
131// jm6(n) = <n-6|j-^6|n>
132//-----------------------
133double jm6(double nj, double j) {
134 return jm(nj - 5, j) * jm(nj - 4, j) * jm(nj - 3, j) * jm(nj - 2, j) * jm(nj - 1, j) * jm(nj, j);
135}
136
137//-----------------------------
138double f20(double nj, double j) { return 3 * pow(nj, 2) - j * (j + 1); }
139
140//------------------------------
141double f21(double nj, double /*unused*/) { return nj; }
142
143//------------------------------
144double f22(double /*unused*/, double /*unused*/) { return 1.0; }
145
146//------------------------------
147double f40(double nj, double j) {
148 return 35 * pow(nj, 4) - 30 * j * (j + 1) * pow(nj, 2) + 25 * pow(nj, 2) - 6 * j * (j + 1) +
149 3 * pow(j, 2) * pow((j + 1), 2);
150}
151
152//------------------------------
153double f41(double nj, double j) { return 7 * pow(nj, 3) - 3 * j * (j + 1) * nj - nj; }
154
155//------------------------------
156double f42(double nj, double j) { return 7 * pow(nj, 2) - j * (j + 1) - 5; }
157
158//------------------------------
159double f43(double nj, double /*unused*/) { return nj; }
160
161//------------------------------
162double f44(double /*unused*/, double /*unused*/) { return 1.0; }
163
164//------------------------------
165double f60(double nj, double j) {
166 return 231 * pow(nj, 6) - 315 * j * (j + 1) * pow(nj, 4) + 735 * pow(nj, 4) +
167 105 * pow(j, 2) * pow((j + 1), 2) * pow(nj, 2) - 525 * j * (j + 1) * pow(nj, 2) + 294 * pow(nj, 2) -
168 5 * pow(j, 3) * pow((j + 1), 3) + 40 * pow(j, 2) * pow((j + 1), 2) - 60 * j * (j + 1);
169}
170
171//------------------------------
172double f61(double nj, double j) {
173 return 33 * pow(nj, 5) - (30 * j * (j + 1) - 15) * pow(nj, 3) - 10 * j * (j + 1) * nj +
174 5 * pow(j, 2) * pow((j + 1), 2) * nj + 12 * nj;
175}
176
177//------------------------------
178double f62(double nj, double j) {
179 return 33 * pow(nj, 4) - 18 * j * (j + 1) * pow(nj, 2) - 123 * pow(nj, 2) + pow(j, 2) * pow(j + 1, 2) +
180 10 * j * (j + 1) + 102;
181}
182
183//------------------------------
184double f63(double nj, double j) { return 11 * pow(nj, 3) - 3 * j * (j + 1) * nj - 59 * nj; }
185
186//------------------------------
187double f64(double nj, double j) { return 11 * pow(nj, 2) - j * (j + 1) - 38; }
188
189//------------------------------
190double f65(double nj, double /*unused*/) { return nj; }
191
192//------------------------------
193double f66(double /*unused*/, double /*unused*/) { return 1.0; }
194
195//-------------------------------
196// ff(k,q,nj,j) := fkq(nj,j)
197//-------------------------------
198double ff(int k, int q, double nj, double j) {
199 int qq = abs(q);
200 if (k == 2 && qq == 0) {
201 return f20(nj, j);
202 } else if (k == 2 && qq == 1) {
203 return f21(nj, j);
204 } else if (k == 2 && qq == 2) {
205 return f22(nj, j);
206 } else if (k == 4 && qq == 0) {
207 return f40(nj, j);
208 } else if (k == 4 && qq == 1) {
209 return f41(nj, j);
210 } else if (k == 4 && qq == 2) {
211 return f42(nj, j);
212 } else if (k == 4 && qq == 3) {
213 return f43(nj, j);
214 } else if (k == 4 && qq == 4) {
215 return f44(nj, j);
216 } else if (k == 6 && qq == 0) {
217 return f60(nj, j);
218 } else if (k == 6 && qq == 1) {
219 return f61(nj, j);
220 } else if (k == 6 && qq == 2) {
221 return f62(nj, j);
222 } else if (k == 6 && qq == 3) {
223 return f63(nj, j);
224 } else if (k == 6 && qq == 4) {
225 return f64(nj, j);
226 } else if (k == 6 && qq == 5) {
227 return f65(nj, j);
228 } else if (k == 6 && qq == 6) {
229 return f66(nj, j);
230 }
231 throw std::runtime_error("Unknown case in ff function.");
232}
233
234// c----------------------------------
235// c testing if the real number r is 0
236// c----------------------------------
237void ifnull(double &r) {
238 const double macheps = 1.0e-14;
239 if (fabs(r) <= macheps)
240 r = 0.0;
241}
242
243//---------------------------------------------
244// calculates the neutron scattering radius r0
245//---------------------------------------------
246double c_r0() {
247 return -1.91 * pow(ee, 2) / me; // * 10**(-12) cm
248}
249
250// c---------------------------------------
251double epsilon(int k, int q) {
252 const std::array<double, 49> eps = {
253 1.00000000000000, 0.707106781186547, 1.00000000000000, -0.707106781186547, 0.612372435695795,
254 1.22474487139159, 0.500000000000000, -1.22474487139159, 0.612372435695795, 0.559016994374947,
255 1.36930639376292, 0.433012701892219, 0.500000000000000, -0.433012701892219, 1.36930639376292,
256 -0.559016994374947, 0.522912516583797, 1.47901994577490, 0.395284707521047, 0.559016994374947,
257 0.125000000000000, -0.559016994374947, 0.395284707521047, -1.47901994577490, 0.522912516583797,
258 0.496078370824611, 1.56873754975139, 0.369754986443726, 0.603807364424560, 0.114108866146910,
259 0.125000000000000, -0.114108866146910, 0.603807364424560, -0.369754986443726, 1.56873754975139,
260 -0.496078370824611, 0.474958879799083, 1.64530582263602, 0.350780380010057, 0.640434422872475,
261 0.106739070478746, 0.135015431216830, 0.062500000000000, -0.135015431216830, 0.106739070478746,
262 -0.640434422872475, 0.350780380010057, -1.64530582263602, 0.474958879799083};
263 return eps[k * (k + 1) + q];
264}
265
266// c---------------------------------------
267double omega(int k, int q) {
268 const std::array<double, 49> oma = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
269 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
270 1.0, 1.0, 3.0, 3.0, 1.0, 3.0, 3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
271 1.0, 3.0, 3.0, 1.0, 3.0, 3.0, 1.0, 1.0, 1.0, 1.0};
272
273 return oma[k * (k + 1) + q];
274}
275
276//--------------------------------------------------------
277// Function to calculate factorial
278//--------------------------------------------------------
279double fac(double n) {
280 if (n < 0.0)
281 return 0.0;
282 if (n == 0.0)
283 return 1.0;
284 double f = 1.0;
285 auto m = static_cast<int>(std::floor(n));
286 for (int i = 1; i <= m; ++i) {
287 f *= i;
288 }
289 return f;
290}
291
292//--------------------------------------------------------
293// binom (n over k)
294//--------------------------------------------------------
295double binom(int n, int k) { return fac(double(n)) / fac(double(k)) / fac(double(n - k)); }
296
297//--------------------------------------------------------
298// (k)
299// calculates D (a,b,c) [ The Wigner D-matrix ]
300// ms m
301//
302// see Lindner A, 'Drehimpulse in der Quantenmechanik',
303// ISBN 3-519-03061-6 (j)
304// Stuttgart: Teubner, 1984 page 86 (d (beta)=...)
305// for equation (1) m ms
306//
307// see Buckmaster phys. stat. sol. (a) 13 (1972) 9
308// for equation (2)
309//
310//--------------------------------------------------------
311ComplexType ddrot(int j, int m, int ms, double a, double b, double c) {
312 // c equation (1)
313 double d = delta(double(ms), double(m), double(j));
314 ifnull(b);
315 if (b != 0.0) {
316 d = 0.0;
317 for (int n = std::max(0, -(m + ms)); n <= std::min(j - m, j - ms); ++n) { // do n=max(0,-(m+ms)),min(j-m,j-ms)
318 d = d + pow(-1.0, (j - ms - n)) * binom(j - ms, n) * binom(j + ms, j - m - n) *
319 pow(cos(0.5 * b), (2 * n + m + ms)) * pow(sin(0.5 * b), (2 * j - 2 * n - m - ms));
320 }
321 d = d * sqrt(fac(double(j + m)) / fac(double(j + ms)) * fac(double(j - m)) / fac(double(j - ms)));
322 }
323 // c equation (2)
324 return ComplexType(cos(m * a), -sin(m * a)) * d * ComplexType(cos(ms * c), -sin(ms * c));
325}
326//--------------------------------------------------
327// <jm| jp^|q| or jm^|q| |nj>
328//--------------------------------------------------
329double jop(int q, double mj, double nj, double j) {
330 switch (q) {
331 case 0:
332 return delta(mj, nj, j);
333 case 1:
334 return delta(mj, nj + 1, j) * jp(nj, j);
335 case 2:
336 return delta(mj, nj + 2, j) * jp2(nj, j);
337 case 3:
338 return delta(mj, nj + 3, j) * jp3(nj, j);
339 case 4:
340 return delta(mj, nj + 4, j) * jp4(nj, j);
341 case 5:
342 return delta(mj, nj + 5, j) * jp5(nj, j);
343 case 6:
344 return delta(mj, nj + 6, j) * jp6(nj, j);
345 case -1:
346 return delta(mj, nj - 1, j) * jm(nj, j);
347 case -2:
348 return delta(mj, nj - 2, j) * jm2(nj, j);
349 case -3:
350 return delta(mj, nj - 3, j) * jm3(nj, j);
351 case -4:
352 return delta(mj, nj - 4, j) * jm4(nj, j);
353 case -5:
354 return delta(mj, nj - 5, j) * jm5(nj, j);
355 case -6:
356 return delta(mj, nj - 6, j) * jm6(nj, j);
357 default:
358 throw std::runtime_error("Cannot calculate jop with this q value.");
359 }
360}
361
362//--------------------------------------------------
363// calculate the full stevens operators
364//--------------------------------------------------
365double full_okq(int k, int q, double mj, double nj, double j) {
366 return 0.5 * jop(q, mj, nj, j) * (ff(k, q, mj, j) + ff(k, q, nj, j));
367}
368
369//-------------------------------
370// calculates <i|j+|k>
371//-------------------------------
372ComplexType matjp(const ComplexFortranMatrix &ev, int i, int k, int dim) {
373 ComplexFortranVector v(1, dim);
374 auto j = 0.5 * (double(dim) - 1.0);
375 for (int s = dim; s >= 2; --s) { // do 10 s=dim,2,-1
376 auto sj = double(s) - j - 1.0;
377 v(s) = ev(s - 1, k) * jp(sj - 1, j);
378 }
379 v(1) = 0.0;
380 ComplexType res = 0.0;
381 for (int s = 1; s <= dim; ++s) { // do 20 s=1,dim
382 res += std::conj(ev(s, i)) * v(s);
383 }
384 return res;
385}
386//--------------------
387// calculates <i|j-|k>
388//--------------------
389ComplexType matjm(const ComplexFortranMatrix &ev, int i, int k, int dim) {
390 ComplexFortranVector v(1, dim);
391 auto j = 0.5 * (double(dim) - 1.0);
392 for (int s = 1; s <= dim - 1; ++s) { // do 10 s=1,dim-1
393 auto sj = double(s) - j - 1.0;
394 v(s) = ev(s + 1, k) * jm(sj + 1, j);
395 }
396 v(dim) = 0.0;
397 ComplexType res = 0.0;
398 for (int s = 1; s <= dim; ++s) { // do 20 s=1,dim
399 res += std::conj(ev(s, i)) * v(s);
400 }
401 return res;
402}
403//--------------------
404// calculates <i|jx|k>
405//--------------------
406ComplexType matjx(const ComplexFortranMatrix &ev, int i, int k, int dim) {
407 return 0.5 * (matjm(ev, i, k, dim) + matjp(ev, i, k, dim));
408}
409//-------------------------
410// calculates |<i|jx|k>|**2
411//-------------------------
412double matjx2(const ComplexFortranMatrix &ev, int i, int k, int dim) { return std::norm(matjx(ev, i, k, dim)); }
413//--------------------
414// calculates <i|jy|k>
415//--------------------
416ComplexType matjy(const ComplexFortranMatrix &ev, int i, int k, int dim) {
417 auto ci = ComplexType(0.0, 1.0);
418 return 0.5 * ci * (matjm(ev, i, k, dim) - matjp(ev, i, k, dim));
419}
420//-------------------------
421// calculates |<i|jy|k>|**2
422//-------------------------
423double matjy2(const ComplexFortranMatrix &ev, int i, int k, int dim) { return std::norm(matjy(ev, i, k, dim)); }
424//--------------------
425// calculates <i|jz|k>
426//--------------------
427ComplexType matjz(const ComplexFortranMatrix &ev, int i, int k, int dim) {
428 ComplexFortranVector v(1, dim);
429 auto j = 0.5 * (double(dim) - 1.0);
430 for (int s = 1; s <= dim; ++s) { // do 10 s=1,dim
431 auto sj = double(s) - j - 1.0;
432 v(s) = ev(s, k) * sj;
433 }
434 ComplexType res = 0.0;
435 for (int s = 1; s <= dim; ++s) { // do 20 s=1,dim
436 res += std::conj(ev(s, i)) * v(s);
437 }
438 return res;
439}
440//-------------------------
441// calculates |<i|jz|k>|**2
442//-------------------------
443double matjz2(const ComplexFortranMatrix &ev, int i, int k, int dim) { return std::norm(matjz(ev, i, k, dim)); }
444
445//---------------------------------------------------------------
446// calculates all transition matrix elements for a single crystal
447// and a polycrystalline sample (powder)
448//---------------------------------------------------------------
449void matcalc(const ComplexFortranMatrix &ev, const int dim, DoubleFortranMatrix &jx2, DoubleFortranMatrix &jy2,
451 for (int i = 1; i <= dim; ++i) { // do 10 i=1,dim
452 for (int k = 1; k <= dim; ++k) { // do 20 k=1,dim
453 jx2(i, k) = matjx2(ev, i, k, dim);
454 jy2(i, k) = matjy2(ev, i, k, dim);
455 jz2(i, k) = matjz2(ev, i, k, dim);
456 jt2(i, k) = 2.0 / 3 * (jx2(i, k) + jy2(i, k) + jz2(i, k));
457 }
458 }
459}
460
461//------------------------------------------------------------
462// make sure that no underflow and overflow appears within exp
463//------------------------------------------------------------
464double exp_(double z) {
465 const double zmax = 71.0;
466 if (z < -zmax) {
467 return 0.0;
468 }
469 if (z < zmax) {
470 return exp(z);
471 }
472 return exp(zmax);
473}
474
475//------------------------------------------
476// calculates all transition intensities for
477// a polycrystalline sample (powder)
478//------------------------------------------
479void intcalc(const double r0, const double gj, const double z, const DoubleFortranMatrix &jt2,
480 const DoubleFortranVector &e, DoubleFortranMatrix &inten, const int dim, double temp) {
481 // Original code from FOCUS calculated integrated intensity in barn
482 // auto constant = 4.0 * pi * pow(0.5 * r0 * gj, 2);
483 // ISIS normalised data is in milibarn/steradian - need to multiply
484 // by 1000 / 4 / PI
485 auto constant = pow(0.5 * r0 * gj, 2) * 1000.;
486 if (temp == 0.0) {
487 temp = 1.0;
488 }
489 // convert temperature to meV
490 temp /= c_fmevkelvin;
491
492 for (int i = 1; i <= dim; ++i) { // do 10 i=1,dim
493 auto coeff = exp_(-e(i) / temp) / z * constant;
494 for (int k = 1; k <= dim; ++k) { // do 20 k=1,dim
495 inten(i, k) = coeff * jt2(i, k);
496 }
497 }
498}
499//-------------------------------------
500// calculation of the occupation factor
501//-------------------------------------
502double c_occupation_factor(const DoubleFortranVector &energy, double dimj, double temp) {
503 auto dim = static_cast<int>(dimj);
504 double occupation_factor = 0.0;
505 if (temp == 0.0) {
506 temp = 1.0;
507 }
508 // convert temperature to meV
509 temp /= c_fmevkelvin;
510 for (int s = 1; s <= dim; ++s) { // do 10 s=1,dim
511 occupation_factor += exp_(-energy(s) / temp);
512 }
513 return occupation_factor;
514}
515
516//--------------------------------------
517// calculation of the zeeman hamiltonian
518//--------------------------------------
519void zeeman(ComplexFortranMatrix &hamiltonian, const int nre, const DoubleFortranVector &bext,
520 const DoubleFortranVector &bmol) {
521 auto i = ComplexType(0.0, 1.0);
522 auto bmolp = bmol(1) + i * bmol(2);
523 auto bmolm = bmol(1) - i * bmol(2);
524 auto bmolz = bmol(3);
525 auto bextp = bext(1) + i * bext(2);
526 auto bextm = bext(1) - i * bext(2);
527 ComplexType bextz = bext(3);
528 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
529 auto facmol = 2 * (gj - 1) * c_myb;
530 auto facext = gj * c_myb;
531 // Negative nre means arbitrary J, with abs(nre) = 2J. dimj=2J+1
532 auto dimj = (nre > 0) ? ddimj[nre - 1] : (abs(nre) + 1);
533 auto j = 0.5 * (dimj - 1.0);
534 auto dim = static_cast<int>(dimj);
535 hamiltonian.allocate(1, dim, 1, dim);
536 hamiltonian.zero();
537 //-------------------------------------------------------------------
538 // define only the lower triangle of h(m,n)
539 //-------------------------------------------------------------------
540 for (int m = 1; m <= dim; ++m) { // do 10 m=1,dim
541 auto mj = double(m) - j - 1.0;
542 for (int n = 1; n <= m; ++n) { // do 20 n=1,m
543 auto nj = double(n) - j - 1.0;
544 // add the molecular field
545 // f*J*B = f*( 1/2*(J+ * B- + J- * B+) + Jz*Bz )
546 hamiltonian(m, n) =
547 hamiltonian(m, n) + 0.5 * facmol * bmolm * delta(mj, nj + 1, j) * jp(nj, j) +
548 0.5 * facmol * bmolp * delta(mj, nj - 1, j) * jm(nj, j) + facmol * bmolz * delta(mj, nj, j) * nj +
549 // c add an external magnetic field
550 0.5 * facext * bextm * delta(mj, nj + 1, j) * jp(nj, j) +
551 0.5 * facext * bextp * delta(mj, nj - 1, j) * jm(nj, j) + facext * bextz * delta(mj, nj, j) * nj;
552 hamiltonian(n, m) = std::conj(hamiltonian(m, n));
553 }
554 }
555}
556
557//---------------------------------------
558// Calculation of the eigenvalues/vectors
559//---------------------------------------
560void diagonalise(const ComplexFortranMatrix &hamiltonian, DoubleFortranVector &eigenvalues,
561 ComplexFortranMatrix &eigenvectors) {
562 // Diagonalisation of the hamiltonian
563 auto dim = hamiltonian.len1();
564 eigenvalues.allocate(1, dim);
565 eigenvectors.allocate(1, dim, 1, dim);
566 ComplexFortranMatrix h = hamiltonian;
567
568 h.eigenSystemHermitian(eigenvalues, eigenvectors);
569
570 // Get the indicies of the eigenvalues sorted in ascending order
571 auto sortedIndices = eigenvalues.sortIndices();
572
573 // Shift the lowest energy level to 0
574 auto indexMin = static_cast<int>(sortedIndices[0] + 1); // default fortran matrix has base 1.
575 auto eshift = eigenvalues[indexMin];
576 eigenvalues += -eshift;
577
578 eigenvalues.sort(sortedIndices);
579 // Eigenvectors are in columns. Sort the columns
580 // to match the sorted eigenvalues.
581 eigenvectors.sortColumns(sortedIndices);
582}
583
584GNU_DIAG_ON("missing-braces")
585
586} // anonymous namespace
587
601 const ComplexFortranMatrix &hamiltonian, int nre, const DoubleFortranVector &bext) {
602 ComplexFortranMatrix h = hamiltonian;
603 DoubleFortranVector bmol(1, 3);
604 bmol.zero();
605 // Adds the external and molecular fields
607 zeeman(hz, nre, bext, bmol);
608 h -= hz;
609 // Now run the actual diagonalisation
610 diagonalise(h, eigenvalues, eigenvectors);
611}
612
631 ComplexFortranMatrix &hamiltonian, ComplexFortranMatrix &hzeeman, int nre,
632 const DoubleFortranVector &bmol, const DoubleFortranVector &bext,
633 const ComplexFortranMatrix &bkq, double alpha_euler, double beta_euler, double gamma_euler) {
634 if (nre > maxNre) {
635 throw std::out_of_range("nre is out of range");
636 }
637
638 // initialize some rare earth constants
639 auto dimj = (nre > 0) ? ddimj[nre - 1] : (abs(nre) + 1);
640
641 //------------------------------------------------------------
642 // transform the Bkq with
643 // H = sum_k=0 Bk0 Ok0 + sum_k>0_q>0 ReBkq ReOkq + ImBkq ImOkq
644 // to a representation with
645 // *
646 // H = sum_kq dkq Okq
647 //
648 // one finds: dk0 = Bk0 and for q<>0: dkq = Bkq/2
649 //-------------------------------------------------------
650 ComplexFortranMatrix dkq_star(1, 6, -6, 6);
651 dkq_star.zero();
652 auto i = ComplexType(0.0, 1.0);
653 for (int k = 2; k <= 6; k += 2) { // do k=2,6,2
654 for (int q = 0; q <= k; ++q) { // do q=0,k
655 ComplexType b_kq = bkq(k, q);
656 dkq_star(k, q) = b_kq;
657 if (q != 0) {
658 dkq_star(k, q) = dkq_star(k, q) / 2.0;
659 }
660 dkq_star(k, -q) = std::conj(dkq_star(k, q));
661 }
662 }
663 //-------------------------------------------------------------------
664 // the parameters are conjg(D_kq)
665 //-------------------------------------------------------------------
666 for (int k = 2; k <= 6; k += 2) { // do k=2,6,2
667 for (int q = -k; q <= k; ++q) { // do q=-k,k
668 dkq_star(k, q) = std::conj(dkq_star(k, q));
669 }
670 }
671 //-------------------------------------------------------------------
672 // Rotate the crystal field quantisation axis by the specified
673 // Euler angles. In some cases the number of CF parameters can be
674 // reduced by chosing the quantisation axis along a high symmetry
675 // rotation axis, rather than along a crystallographic axis.
676 // The eigenvalues should remain the same, but the eigenvectors will
677 // change.
678 // The rotation is done using a Wigner D-matrix. As noted by
679 // Buckmaster, the Stevens operators cannot be rotated as is, because
680 // the have an inconsistent normalisation between differen k, q terms
681 // Thus they have to be converted to and from the "Wybourne"
682 // normalisation using the epsilon / omega values.
683 // There was a bug in the original FOCUS code. Multiplying by
684 // omega*epsilon converts Wybourne parameters to Stevens parameters.
685 // Thus to convert the original dkq_star(k,qs) from Steven to Wybourn
686 // we should divide by epsilon(k,qs)*omega(k,qs) and then multiply by
687 // the new dkq_star(k,q) by epsilon(k,q)*omega(k,q).
688 //-------------------------------------------------------------------
689 ComplexFortranMatrix rdkq_star(1, 6, -6, 6);
690 for (int k = 2; k <= 6; k += 2) { // do k=2,6,2
691 for (int q = -k; q <= k; ++q) { // do q=-k,k
692 rdkq_star(k, q) = ComplexType(0.0, 0.0);
693 for (int qs = -k; qs <= k; ++qs) { // do qs=-k,k
694 rdkq_star(k, q) = rdkq_star(k, q) + dkq_star(k, qs) * epsilon(k, q) / epsilon(k, qs) * omega(k, q) /
695 omega(k, qs) * ddrot(k, q, qs, alpha_euler, beta_euler, gamma_euler);
696 }
697 }
698 }
699 // -------------------------------------------------------------------
700 // rotate the external magnetic field
701 // -------------------------------------------------------------------
702 // aj = sqrt(3/j/(j+1)/(2*j+1))
703 //
704 // the following hamiltonian will be rotated
705 //
706 // cry ----- cry cry
707 // h = f * > Bex V with <j|| V1 ||j> = sqrt(3)
708 // mag ----- 1q 1q
709 // q
710 // -------------------------------------------------------------------
711
712 ComplexFortranMatrix bex(1, 1, -1, 1);
713 bex(1, 1) = -(bext(1) - i * bext(2)) * M_SQRT1_2;
714 bex(1, -1) = (bext(1) + i * bext(2)) * M_SQRT1_2;
715 bex(1, 0) = bext(3);
716
717 // calculates Bex(1,q) for a canted moment
718 //
719 // moment cry cry cry +
720 // h = R (abc) h R (abc)
721 // mag mag
722 //
723 ComplexFortranMatrix rbex(1, 1, -1, 1);
724 for (int q = -1; q <= 1; ++q) { // do q=-1,1
725 rbex(1, q) = ComplexType(0.0, 0.0);
726 for (int qs = -1; qs <= 1; ++qs) { // do qs=-1,1
727 rbex(1, q) = rbex(1, q) + bex(1, qs) * ddrot(1, q, qs, alpha_euler, beta_euler, gamma_euler);
728 }
729 }
730 DoubleFortranVector rbext(1, 3);
731 rbext(1) = real(static_cast<ComplexType>(rbex(1, -1))) * M_SQRT2;
732 rbext(2) = imag(static_cast<ComplexType>(rbex(1, -1))) * M_SQRT2;
733 rbext(3) = real(static_cast<ComplexType>(rbex(1, 0)));
734
735 auto dim = static_cast<int>(dimj);
736 hamiltonian.allocate(1, dim, 1, dim);
737 hamiltonian.zero();
738 //-------------------------------------------------------------------
739 // calculate the crystal field hamiltonian
740 // define only the lower triangle of h(m,n)
741 //-------------------------------------------------------------------
742 // total momentum J of R3+
743 auto j = 0.5 * (dimj - 1.0);
744 for (int m = 1; m <= dim; ++m) { // do m=1,dim
745 auto mj = double(m) - j - 1.0;
746 for (int n = 1; n <= m; ++n) { // do n=1,m
747 auto nj = double(n) - j - 1.0;
748 for (int k = 2; k <= 6; k += 2) { // do k=2,6,2
749 for (int q = -k; q <= k; ++q) { // do q=-k,k
750 hamiltonian(m, n) = hamiltonian(m, n) + rdkq_star(k, q) * full_okq(k, q, mj, nj, j);
751 }
752 }
753 hamiltonian(n, m) = std::conj(hamiltonian(m, n));
754 }
755 }
756
757 // Adds the external and molecular fields
758 zeeman(hzeeman, nre, rbext, bmol);
759 hamiltonian -= hzeeman;
760
761 // Now run the actual diagonalisation
762 diagonalise(hamiltonian, eigenvalues, eigenvectors);
763}
764
765//-------------------------
766// transforms the indices
767//-------------------------
768int no(int i, const IntFortranVector &d, int n) {
769 // implicit none
770 // integer i,d(17*17),n,up,lo
771 int up = 0;
772 for (int nno = 1; nno <= n; ++nno) { // do no=1,n
773 int lo = up + 1;
774 up = lo + d(nno) - 1;
775 if (i >= lo && i <= up) {
776 return nno;
777 }
778 }
779 return n;
780}
781
793 DoubleFortranVector &e_energies, DoubleFortranMatrix &i_energies, double de) {
794 // real*8 energy(17) ! already defined in CF_FABI.INC
795 // real*8 mat(17,17)
796 // integer degeneration(17*17) ! stores the degeneration of a level
797 // integer n_energies ! no. of degenerated energy levels
798 // real*8 e_energies(17*17) ! energy values of the degenerated energy
799 // levels
800 // real*8 i_energies(17,17) ! intensities of the degenerated energy
801 // levels
802 // real*8 de,di
803
804 // energy levels which are closer than de are assumed to be degenerated
805 // only those excitations are taken into account whose intensities are
806 // greater equal than di
807
808 const auto dim = static_cast<int>(energy.size());
809 IntFortranVector tempDegeneration(dim);
810 DoubleFortranVector tempEnergies(dim);
811
812 // find out how many degenerated energy levels exists
813 tempEnergies(1) = 0.0;
814 tempDegeneration(1) = 1;
815 int n_energies = 1;
816 for (int i = 2; i <= dim; ++i) {
817 if (std::fabs(tempEnergies(n_energies) - energy(i)) >= de) {
818 n_energies += 1;
819 tempEnergies(n_energies) = energy(i);
820 tempDegeneration(n_energies) = 1;
821 } else {
822 tempDegeneration(n_energies) = tempDegeneration(n_energies) + 1;
823 }
824 }
825
826 // Resize the output arrays
827 degeneration.allocate(n_energies);
828 e_energies.allocate(n_energies);
829 i_energies.allocate(n_energies, n_energies);
830
831 // store the intensities of the degenarated levels
832 i_energies.zero();
833 for (int i = 1; i <= dim; ++i) {
834 int ii = no(i, tempDegeneration, n_energies);
835 degeneration(ii) = tempDegeneration(ii);
836 e_energies(ii) = tempEnergies(ii);
837 for (int k = 1; k <= dim; ++k) {
838 int kk = no(k, tempDegeneration, n_energies);
839 i_energies(ii, kk) = i_energies(ii, kk) + mat(i, k);
840 }
841 }
842}
843
854void calculateIntensities(int nre, const DoubleFortranVector &energies, const ComplexFortranMatrix &wavefunctions,
855 double temperature, double de, IntFortranVector &degeneration,
856 DoubleFortranVector &e_energies, DoubleFortranMatrix &i_energies) {
857 auto dim = static_cast<int>(energies.size());
858 auto dimj = (nre > 0) ? ddimj[nre - 1] : (abs(nre) + 1);
859 if (static_cast<double>(dim) != dimj) {
860 throw std::runtime_error("calculateIntensities was called for a wrong ion");
861 }
862
863 // calculates the transition matrixelements for a single crystal and
864 // a powdered sample
865 DoubleFortranMatrix jx2mat(1, dim, 1, dim);
866 DoubleFortranMatrix jy2mat(1, dim, 1, dim);
867 DoubleFortranMatrix jz2mat(1, dim, 1, dim);
868 DoubleFortranMatrix jt2mat(1, dim, 1, dim);
869 matcalc(wavefunctions, dim, jx2mat, jy2mat, jz2mat, jt2mat);
870
871 // calculates the sum over all occupation_factor
872 auto occupation_factor = c_occupation_factor(energies, dimj, temperature);
873
874 // calculates the transition intensities for a powdered sample
875 auto r0 = c_r0();
876 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
877 DoubleFortranMatrix mat(1, dim, 1, dim);
878 intcalc(r0, gj, occupation_factor, jt2mat, energies, mat, dim, temperature);
879
880 deg_on(energies, mat, degeneration, e_energies, i_energies, de);
881}
882
893void calculateExcitations(const DoubleFortranVector &e_energies, const DoubleFortranMatrix &i_energies, double de,
894 double di, DoubleFortranVector &e_excitations, DoubleFortranVector &i_excitations) {
895 auto n_energies = static_cast<int>(e_energies.size());
896 auto dimj = n_energies;
897 // Calculate transition energies (excitations) and corresponding
898 // intensities.
899 DoubleFortranVector eex(n_energies * n_energies);
900 DoubleFortranVector iex(n_energies * n_energies);
901 int nex = 0;
902 for (int i = 1; i <= n_energies; ++i) { // do i=1,n_energies
903 for (int k = 1; k <= n_energies; ++k) { // do k=1,n_energies
904 nex = nex + 1;
905 eex(nex) = e_energies(k) - e_energies(i);
906 iex(nex) = i_energies(i, k); // !I(i->k)
907 ifnull(eex(nex));
908 }
909 }
910 // nex at this point is the largest possible number
911 // of transitions
912
913 // Sort excitations in ascending order.
914 std::vector<size_t> ind = eex.sortIndices();
915 eex.sort(ind);
916
917 // Define lambda that transforms C++ indices to fortran indices
918 auto index = [&ind](int i) { return static_cast<int>(ind[i - 1] + 1); };
919
920 DoubleFortranVector tempEex(nex);
921 IntFortranVector degeneration(nex);
922 tempEex(1) = eex(1);
923 degeneration(1) = 1;
924
925 int k = 1;
926 // Check if there are any overlapping (degenerate) excitations
927 for (int i = 2; i <= nex; ++i) { // do i=2,nex
928 if (std::fabs(tempEex(k) - eex(i)) >= de) {
929 k = k + 1;
930 tempEex(k) = eex(i);
931 degeneration(k) = 1;
932 } else {
933 degeneration(k) = degeneration(k) + 1;
934 }
935 }
936 // Number of non-overlapping transitions
937 int n_excitation = k;
938
939 DoubleFortranVector tempIex(n_excitation);
940 for (int i = 1; i <= nex; ++i) { // do i=1,nex
941 auto ii = no(i, degeneration, n_excitation);
942 tempIex(ii) = tempIex(ii) + iex(index(i));
943 }
944
945 ind = tempIex.sortIndices(false);
946 tempIex.sort(ind);
947
948 // Keep only transitions that are strong enough
949 // i >= di
950 e_excitations.allocate(n_excitation);
951 i_excitations.allocate(n_excitation);
952 k = 0;
953 for (int i = 1; i <= n_excitation; ++i) { // do i=1,n_excitation
954 if (tempIex(i) >= di || dimj == 1) {
955 k = k + 1;
956 i_excitations(k) = tempIex(i);
957 e_excitations(k) = tempEex(index(i));
958 }
959 }
960 // nex now is the actual number of excitations that will
961 // be outputted.
962 nex = k;
963 if (nex != n_excitation) {
964 e_excitations.allocate(nex);
965 i_excitations.allocate(nex);
966 }
967}
968
975void calculateMagneticMoment(const ComplexFortranMatrix &ev, const DoubleFortranVector &Hdir, const int nre,
976 DoubleFortranVector &moment) {
977 int dim = (nre > 0) ? (int)ddimj[nre - 1] : (abs(nre) + 1);
978 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
979 moment.allocate(dim);
980 for (auto i = 1; i <= dim; ++i) {
981 moment(i) = real(matjx(ev, i, i, dim)) * Hdir(1) + // <ev|jx|ev>
982 real(matjy(ev, i, i, dim)) * Hdir(2) + // <ev|jy|ev>
983 real(matjz(ev, i, i, dim)) * Hdir(3); // <ev|jz|ev>
984 }
985 moment *= gj;
986}
987
993void calculateMagneticMomentMatrix(const ComplexFortranMatrix &ev, const std::vector<double> &Hdir, const int nre,
994 ComplexFortranMatrix &mumat) {
995 int dim = (nre > 0) ? (int)ddimj[nre - 1] : (abs(nre) + 1);
996 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
997 mumat.allocate(1, dim, 1, dim);
998 for (auto i = 1; i <= dim; ++i) {
999 for (auto j = 1; j <= dim; ++j) {
1000 mumat(i, j) = matjx(ev, i, j, dim) * Hdir[0] + // <ev|jx|ev'>
1001 matjy(ev, i, j, dim) * Hdir[1] + // <ev|jy|ev'>
1002 matjz(ev, i, j, dim) * Hdir[2]; // <ev|jz|ev'>
1003 }
1004 }
1005 mumat *= gj;
1006}
1007
1008} // namespace Mantid::CurveFitting::Functions
double energy
Definition: GetAllEi.cpp:157
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
#define GNU_DIAG_ON(x)
#define GNU_DIAG_OFF(x)
This is a collection of macros for turning compiler warnings off in a controlled manner.
void zero()
Set all elements to zero.
void eigenSystemHermitian(EigenVector &eigenValues, ComplexMatrix &eigenVectors)
Calculate the eigensystem of a Hermitian matrix.
void zero()
Set all elements to zero.
void sort(const std::vector< size_t > &indices)
Sort this vector in order defined by an index array.
std::vector< size_t > sortIndices(bool ascending=true) const
Create an index array that would sort this vector.
size_t size() const
Size of the vector.
void allocate(const int iFrom, const int iTo, const int jFrom, const int jTo)
Resize the matrix.
void allocate(int firstIndex, int lastIndex)
Resize the vector.
void MANTID_CURVEFITTING_DLL calculateMagneticMoment(const ComplexFortranMatrix &ev, const DoubleFortranVector &Hmag, const int nre, DoubleFortranVector &moment)
Calculate the diagonal matrix elements of the magnetic moment operator in a particular eigenvector ba...
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.
void deg_on(const DoubleFortranVector &energy, const DoubleFortranMatrix &mat, IntFortranVector &degeneration, DoubleFortranVector &e_energies, DoubleFortranMatrix &i_energies, double de)
Find out how many degenerated energy levels exists.
void MANTID_CURVEFITTING_DLL calculateExcitations(const DoubleFortranVector &e_energies, const DoubleFortranMatrix &i_energies, double de, double di, DoubleFortranVector &e_excitations, DoubleFortranVector &i_excitations)
Calculate the excitations (transition energies) and their intensities.
void MANTID_CURVEFITTING_DLL calculateZeemanEigensystem(DoubleFortranVector &eigenvalues, ComplexFortranMatrix &eigenvectors, const ComplexFortranMatrix &hamiltonian, int nre, const DoubleFortranVector &bext)
Calculates the eigenvalues/vectors of a crystal field Hamiltonian in a specified external magnetic fi...
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.
int no(int i, const IntFortranVector &d, int n)
void MANTID_CURVEFITTING_DLL calculateIntensities(int nre, const DoubleFortranVector &energies, const ComplexFortranMatrix &wavefunctions, double temperature, double de, IntFortranVector &degeneration, DoubleFortranVector &e_energies, DoubleFortranMatrix &i_energies)
Calculate the intensities of transitions.
FortranMatrix< EigenMatrix > DoubleFortranMatrix
FortranVector< ComplexVector > ComplexFortranVector
std::complex< double > ComplexType
FortranMatrix< ComplexMatrix > ComplexFortranMatrix
FortranVector< EigenVector > DoubleFortranVector
STL namespace.