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.};
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};
36const double pi = 4.0 * atan(1.0);
37const double kb = 1.38062;
38const double hh = 6.626075540;
39const double hq = hh / 2 / pi;
40const double ee = 1.6021773349;
41const double me = 9.109389754;
46const double c_fmevkelvin = 10 * ee / kb;
48const double c_myb = hq / me / 2;
53double delta(
double mj,
double nj,
double j) {
55 if (mj == nj &&
fabs(mj) <= j &&
fabs(nj) <= j) {
64double jp(
double nj,
double j) {
66 return sqrt(j * (j + 1) - nj * (nj + 1));
75double jm(
double nj,
double j) {
77 return sqrt(j * (j + 1) - nj * (nj - 1));
86double jp2(
double nj,
double j) {
return jp(nj + 1, j) * jp(nj, j); }
91double jm2(
double nj,
double j) {
return jm(nj - 1, j) * jm(nj, j); }
96double jp3(
double nj,
double j) {
return jp(nj + 2, j) * jp(nj + 1, j) * jp(nj, j); }
101double jm3(
double nj,
double j) {
return jm(nj - 2, j) * jm(nj - 1, j) * jm(nj, j); }
106double jp4(
double nj,
double j) {
return jp(nj + 3, j) * jp(nj + 2, j) * jp(nj + 1, j) * jp(nj, j); }
111double jm4(
double nj,
double j) {
return jm(nj - 3, j) * jm(nj - 2, j) * jm(nj - 1, j) * jm(nj, j); }
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); }
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); }
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);
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);
138double f20(
double nj,
double j) {
return 3 * pow(nj, 2) - j * (j + 1); }
141double f21(
double nj,
double ) {
return nj; }
144double f22(
double ,
double ) {
return 1.0; }
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);
153double f41(
double nj,
double j) {
return 7 * pow(nj, 3) - 3 * j * (j + 1) * nj - nj; }
156double f42(
double nj,
double j) {
return 7 * pow(nj, 2) - j * (j + 1) - 5; }
159double f43(
double nj,
double ) {
return nj; }
162double f44(
double ,
double ) {
return 1.0; }
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);
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;
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;
184double f63(
double nj,
double j) {
return 11 * pow(nj, 3) - 3 * j * (j + 1) * nj - 59 * nj; }
187double f64(
double nj,
double j) {
return 11 * pow(nj, 2) - j * (j + 1) - 38; }
190double f65(
double nj,
double ) {
return nj; }
193double f66(
double ,
double ) {
return 1.0; }
198double ff(
int k,
int q,
double nj,
double j) {
200 if (k == 2 && qq == 0) {
202 }
else if (k == 2 && qq == 1) {
204 }
else if (k == 2 && qq == 2) {
206 }
else if (k == 4 && qq == 0) {
208 }
else if (k == 4 && qq == 1) {
210 }
else if (k == 4 && qq == 2) {
212 }
else if (k == 4 && qq == 3) {
214 }
else if (k == 4 && qq == 4) {
216 }
else if (k == 6 && qq == 0) {
218 }
else if (k == 6 && qq == 1) {
220 }
else if (k == 6 && qq == 2) {
222 }
else if (k == 6 && qq == 3) {
224 }
else if (k == 6 && qq == 4) {
226 }
else if (k == 6 && qq == 5) {
228 }
else if (k == 6 && qq == 6) {
231 throw std::runtime_error(
"Unknown case in ff function.");
237void ifnull(
double &r) {
238 const double macheps = 1.0e-14;
239 if (
fabs(r) <= macheps)
247 return -1.91 * pow(ee, 2) / me;
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];
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};
273 return oma[k * (k + 1) + q];
279double fac(
double n) {
285 auto m =
static_cast<int>(std::floor(
n));
286 for (
int i = 1; i <=
m; ++i) {
295double binom(
int n,
int k) {
return fac(
double(
n)) / fac(
double(k)) / fac(
double(
n - k)); }
311ComplexType ddrot(
int j,
int m,
int ms,
double a,
double b,
double c) {
313 double d =
delta(
double(ms),
double(
m),
double(j));
317 for (
int n = std::max(0, -(
m + ms));
n <= std::min(j -
m, j - ms); ++
n) {
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));
321 d =
d * sqrt(fac(
double(j +
m)) / fac(
double(j + ms)) * fac(
double(j -
m)) / fac(
double(j - ms)));
329double jop(
int q,
double mj,
double nj,
double j) {
332 return delta(mj, nj, j);
334 return delta(mj, nj + 1, j) * jp(nj, j);
336 return delta(mj, nj + 2, j) * jp2(nj, j);
338 return delta(mj, nj + 3, j) * jp3(nj, j);
340 return delta(mj, nj + 4, j) * jp4(nj, j);
342 return delta(mj, nj + 5, j) * jp5(nj, j);
344 return delta(mj, nj + 6, j) * jp6(nj, j);
346 return delta(mj, nj - 1, j) * jm(nj, j);
348 return delta(mj, nj - 2, j) * jm2(nj, j);
350 return delta(mj, nj - 3, j) * jm3(nj, j);
352 return delta(mj, nj - 4, j) * jm4(nj, j);
354 return delta(mj, nj - 5, j) * jm5(nj, j);
356 return delta(mj, nj - 6, j) * jm6(nj, j);
358 throw std::runtime_error(
"Cannot calculate jop with this q value.");
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));
374 auto j = 0.5 * (double(dim) - 1.0);
375 for (
int s = dim; s >= 2; --s) {
376 auto sj = double(s) - j - 1.0;
377 v(s) = ev(s - 1, k) * jp(sj - 1, j);
381 for (
int s = 1; s <= dim; ++s) {
382 res += std::conj(ev(s, i)) * v(s);
391 auto j = 0.5 * (double(dim) - 1.0);
392 for (
int s = 1; s <= dim - 1; ++s) {
393 auto sj = double(s) - j - 1.0;
394 v(s) = ev(s + 1, k) * jm(sj + 1, j);
398 for (
int s = 1; s <= dim; ++s) {
399 res += std::conj(ev(s, i)) * v(s);
407 return 0.5 * (matjm(ev, i, k, dim) + matjp(ev, i, k, dim));
412double matjx2(
const ComplexFortranMatrix &ev,
int i,
int k,
int dim) {
return std::norm(matjx(ev, i, k, dim)); }
418 return 0.5 * ci * (matjm(ev, i, k, dim) - matjp(ev, i, k, dim));
423double matjy2(
const ComplexFortranMatrix &ev,
int i,
int k,
int dim) {
return std::norm(matjy(ev, i, k, dim)); }
429 auto j = 0.5 * (double(dim) - 1.0);
430 for (
int s = 1; s <= dim; ++s) {
431 auto sj = double(s) - j - 1.0;
432 v(s) = ev(s, k) * sj;
435 for (
int s = 1; s <= dim; ++s) {
436 res += std::conj(ev(s, i)) * v(s);
443double matjz2(
const ComplexFortranMatrix &ev,
int i,
int k,
int dim) {
return std::norm(matjz(ev, i, k, dim)); }
451 for (
int i = 1; i <= dim; ++i) {
452 for (
int k = 1; k <= dim; ++k) {
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));
464double exp_(
double z) {
465 const double zmax = 71.0;
485 auto constant = pow(0.5 * r0 * gj, 2) * 1000.;
490 temp /= c_fmevkelvin;
492 for (
int i = 1; i <= dim; ++i) {
493 auto coeff = exp_(-e(i) / temp) /
z * constant;
494 for (
int k = 1; k <= dim; ++k) {
495 inten(i, k) = coeff * jt2(i, k);
503 auto dim =
static_cast<int>(dimj);
504 double occupation_factor = 0.0;
509 temp /= c_fmevkelvin;
510 for (
int s = 1; s <= dim; ++s) {
511 occupation_factor += exp_(-
energy(s) / temp);
513 return occupation_factor;
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);
528 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
529 auto facmol = 2 * (gj - 1) * c_myb;
530 auto facext = gj * c_myb;
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);
540 for (
int m = 1;
m <= dim; ++
m) {
541 auto mj = double(
m) - j - 1.0;
542 for (
int n = 1;
n <=
m; ++
n) {
543 auto nj = double(
n) - j - 1.0;
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 +
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));
563 auto dim = hamiltonian.len1();
564 eigenvalues.allocate(1, dim);
565 eigenvectors.allocate(1, dim, 1, dim);
571 auto sortedIndices = eigenvalues.sortIndices();
574 auto indexMin =
static_cast<int>(sortedIndices[0] + 1);
575 auto eshift = eigenvalues[indexMin];
576 eigenvalues += -eshift;
578 eigenvalues.sort(sortedIndices);
581 eigenvectors.sortColumns(sortedIndices);
607 zeeman(hz, nre, bext, bmol);
610 diagonalise(h, eigenvalues, eigenvectors);
635 throw std::out_of_range(
"nre is out of range");
639 auto dimj = (nre > 0) ? ddimj[nre - 1] : (abs(nre) + 1);
653 for (
int k = 2; k <= 6; k += 2) {
654 for (
int q = 0; q <= k; ++q) {
656 dkq_star(k, q) = b_kq;
658 dkq_star(k, q) = dkq_star(k, q) / 2.0;
660 dkq_star(k, -q) = std::conj(dkq_star(k, q));
666 for (
int k = 2; k <= 6; k += 2) {
667 for (
int q = -k; q <= k; ++q) {
668 dkq_star(k, q) = std::conj(dkq_star(k, q));
690 for (
int k = 2; k <= 6; k += 2) {
691 for (
int q = -k; q <= k; ++q) {
693 for (
int qs = -k; qs <= k; ++qs) {
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);
713 bex(1, 1) = -(bext(1) - i * bext(2)) * M_SQRT1_2;
714 bex(1, -1) = (bext(1) + i * bext(2)) * M_SQRT1_2;
724 for (
int q = -1; q <= 1; ++q) {
726 for (
int qs = -1; qs <= 1; ++qs) {
727 rbex(1, q) = rbex(1, q) + bex(1, qs) * ddrot(1, q, qs, alpha_euler, beta_euler, gamma_euler);
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)));
735 auto dim =
static_cast<int>(dimj);
736 hamiltonian.
allocate(1, dim, 1, dim);
743 auto j = 0.5 * (dimj - 1.0);
744 for (
int m = 1;
m <= dim; ++
m) {
745 auto mj = double(
m) - j - 1.0;
746 for (
int n = 1;
n <=
m; ++
n) {
747 auto nj = double(
n) - j - 1.0;
748 for (
int k = 2; k <= 6; k += 2) {
749 for (
int q = -k; q <= k; ++q) {
750 hamiltonian(
m,
n) = hamiltonian(
m,
n) + rdkq_star(k, q) * full_okq(k, q, mj, nj, j);
753 hamiltonian(
n,
m) = std::conj(hamiltonian(
m,
n));
758 zeeman(hzeeman, nre, rbext, bmol);
759 hamiltonian -= hzeeman;
762 diagonalise(hamiltonian, eigenvalues, eigenvectors);
772 for (
int nno = 1; nno <=
n; ++nno) {
774 up = lo +
d(nno) - 1;
775 if (i >= lo && i <= up) {
808 const auto dim =
static_cast<int>(
energy.size());
813 tempEnergies(1) = 0.0;
814 tempDegeneration(1) = 1;
816 for (
int i = 2; i <= dim; ++i) {
817 if (std::fabs(tempEnergies(n_energies) -
energy(i)) >= de) {
819 tempEnergies(n_energies) =
energy(i);
820 tempDegeneration(n_energies) = 1;
822 tempDegeneration(n_energies) = tempDegeneration(n_energies) + 1;
829 i_energies.
allocate(n_energies, n_energies);
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);
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");
869 matcalc(wavefunctions, dim, jx2mat, jy2mat, jz2mat, jt2mat);
872 auto occupation_factor = c_occupation_factor(energies, dimj, temperature);
876 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
878 intcalc(r0, gj, occupation_factor, jt2mat, energies, mat, dim, temperature);
880 deg_on(energies, mat, degeneration, e_energies, i_energies, de);
895 auto n_energies =
static_cast<int>(e_energies.
size());
896 auto dimj = n_energies;
902 for (
int i = 1; i <= n_energies; ++i) {
903 for (
int k = 1; k <= n_energies; ++k) {
905 eex(nex) = e_energies(k) - e_energies(i);
906 iex(nex) = i_energies(i, k);
918 auto index = [&ind](
int i) {
return static_cast<int>(ind[i - 1] + 1); };
927 for (
int i = 2; i <= nex; ++i) {
928 if (std::fabs(tempEex(k) - eex(i)) >= de) {
933 degeneration(k) = degeneration(k) + 1;
937 int n_excitation = k;
940 for (
int i = 1; i <= nex; ++i) {
941 auto ii =
no(i, degeneration, n_excitation);
942 tempIex(ii) = tempIex(ii) + iex(
index(i));
950 e_excitations.
allocate(n_excitation);
951 i_excitations.
allocate(n_excitation);
953 for (
int i = 1; i <= n_excitation; ++i) {
954 if (tempIex(i) >= di || dimj == 1) {
956 i_excitations(k) = tempIex(i);
957 e_excitations(k) = tempEex(
index(i));
963 if (nex != n_excitation) {
977 int dim = (nre > 0) ? (
int)ddimj[nre - 1] : (abs(nre) + 1);
978 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
980 for (
auto i = 1; i <= dim; ++i) {
981 moment(i) = real(matjx(ev, i, i, dim)) * Hdir(1) +
982 real(matjy(ev, i, i, dim)) * Hdir(2) +
983 real(matjz(ev, i, i, dim)) * Hdir(3);
995 int dim = (nre > 0) ? (
int)ddimj[nre - 1] : (abs(nre) + 1);
996 auto gj = (nre > 0) ? ggj[nre - 1] : 2.;
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] +
1001 matjy(ev, i, j, dim) * Hdir[1] +
1002 matjz(ev, i, j, dim) * Hdir[2];
std::map< DeltaEMode::Type, std::string > index
#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 °eneration, 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 °eneration, 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