Skip to content

Commit 2d94530

Browse files
mduclemducle
committed
Added cached ham, evec/val. Tests evec/val. Removed invstevfact.
1 parent 3fa3014 commit 2d94530

File tree

5 files changed

+84
-24
lines changed

5 files changed

+84
-24
lines changed

‎include/cf1ion.hpp

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,21 @@ namespace libMcPhase {
1616

1717
class cf1ion: public cfpars {
1818

19+
protected:
20+
bool m_ham_calc = false; // Flag to indicate if Hamiltonian calculated
21+
bool m_ev_calc = false; // Flag to indicate if eigenvectors/values calculated
22+
RowMatrixXcd m_hamiltonian; // Cached Hamiltonian
23+
RowMatrixXcd m_eigenvectors; // Cached eigenvectors
24+
VectorXd m_eigenvalues; // Cached eigenvalues
25+
1926
public:
27+
// Setters
28+
virtual void set_unit(const Units newunit);
29+
virtual void set_type(const Type newtype);
30+
virtual void set_name(const std::string &ionname);
31+
virtual void set_J(const double J);
32+
virtual void set(const Blm blm, double val);
33+
virtual void set(int l, int m, double val);
2034
// Constructors
2135
cf1ion() : cfpars() {};
2236
cf1ion(const int J2) : cfpars(J2) {};

‎include/cfpars.hpp

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,6 @@ class cfpars {
3838
std::array<double, 27> m_Bo{}; // Output array of parameters.
3939
std::array<double, 3> m_rk{}; // Radial integrals <r^k> (if constructed from ionname)
4040
std::array<double, 3> m_stevfact = {{1., 1., 1.}}; // Stevens factors \theta_k
41-
std::array<double, 3> m_invstevfact = {{1., 1., 1.}}; // Inverse Stevens factors
4241
std::array<double, 27> m_convfact; // Conversion factor from internal to external parameters
4342
double m_econv = 1.; // Conversion factor from internal to external energy
4443
Units m_unit = Units::meV; // Energy units of parameters
@@ -62,12 +61,12 @@ class cfpars {
6261
const double gamma() { return m_stevfact[2]; }
6362
const double get_J() { return (double)(m_J2 / 2.); }
6463
// Setters
65-
void set_unit(const Units newunit);
66-
void set_type(const Type newtype);
67-
void set_name(const std::string &ionname);
68-
void set_J(const double J);
69-
void set(const Blm blm, double val);
70-
void set(int l, int m, double val);
64+
virtual void set_unit(const Units newunit);
65+
virtual void set_type(const Type newtype);
66+
virtual void set_name(const std::string &ionname);
67+
virtual void set_J(const double J);
68+
virtual void set(const Blm blm, double val);
69+
virtual void set(int l, int m, double val);
7170
// Constructors
7271
cfpars();
7372
cfpars(const int J2);

‎singleion/cf1ion.cpp

Lines changed: 58 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -7,27 +7,69 @@
77
*/
88

99
#include "cf1ion.hpp"
10-
#include <iostream>
1110

1211
namespace libMcPhase {
1312

1413
static const std::array<std::array<int, 4>, 12> idq = { {{2,2,0,4}, {2,1,1,3}, {4,4,5,13}, {4,3,6,12}, {4,2,7,11}, {4,1,8,10},
1514
{6,6,14,26}, {6,5,15,25}, {6,4,16,24}, {6,3,17,23}, {6,2,18,22}, {6,1,19,21}} };
1615
static const std::array<std::array<int, 2>, 3> idq0 = { {{2,2}, {4,9}, {6,20}} };
1716

17+
// --------------------------------------------------------------------------------------------------------------- //
18+
// Setter/getter methods for cfpars class
19+
// --------------------------------------------------------------------------------------------------------------- //
20+
void cf1ion::set(const Blm blm, double val) {
21+
cfpars::set(blm, val);
22+
m_ham_calc = false;
23+
m_ev_calc = false;
24+
}
25+
26+
void cf1ion::set(int l, int m, double val) {
27+
cfpars::set(l, m, val);
28+
m_ham_calc = false;
29+
m_ev_calc = false;
30+
}
31+
32+
void cf1ion::set_unit(cfpars::Units const newunit) {
33+
cfpars::set_unit(newunit);
34+
m_ham_calc = false;
35+
m_ev_calc = false;
36+
}
37+
38+
void cf1ion::set_type(const cfpars::Type newtype) {
39+
cfpars::set_type(newtype);
40+
m_ham_calc = false;
41+
m_ev_calc = false;
42+
}
43+
44+
void cf1ion::set_name(const std::string &ionname) {
45+
cfpars::set_name(ionname);
46+
m_ham_calc = false;
47+
m_ev_calc = false;
48+
}
49+
50+
void cf1ion::set_J(const double J) {
51+
cfpars::set_J(J);
52+
m_ham_calc = false;
53+
m_ev_calc = false;
54+
}
55+
1856
// --------------------------------------------------------------------------------------------------------------- //
1957
// General methods for the cf1ion class
2058
// --------------------------------------------------------------------------------------------------------------- //
2159
RowMatrixXcd cf1ion::hamiltonian(bool upper) {
60+
if (m_ham_calc) {
61+
return m_hamiltonian;
62+
}
2263
if (m_J2 <= 0) {
2364
throw std::runtime_error("Invalid value of J - must be strictly greater than zero");
2465
}
2566

2667
int dimj = m_J2 + 1;
27-
RowMatrixXcd ham = RowMatrixXcd::Zero(dimj, dimj);
68+
m_hamiltonian = RowMatrixXcd::Zero(dimj, dimj);
2869
// For J=1/2, all Stevens operators are zero
2970
if (dimj == 2) {
30-
return ham;
71+
m_ham_calc = true;
72+
return m_hamiltonian;
3173
}
3274

3375
// The crystal field potential operator is given by:
@@ -84,11 +126,10 @@ RowMatrixXcd cf1ion::hamiltonian(bool upper) {
84126
// First the diagonal elements (q=0 terms)
85127
for (auto iq: idq0) {
86128
int k = iq[0], m = iq[1];
87-
std::cout << "k=" << k << ", m=" << m << ", m_Bi[m]=" << m_Bi[m] << ", rme[k]=" << rme[k] << ", m_econv=" << m_econv << "\n";
88129
if (std::fabs(m_Bi[m]) > 1e-12) {
89130
for (int i=0; i<dimj; i++) {
90131
int mj = 2*i - m_J2;
91-
ham(i,i) += pow(-1., (m_J2-mj)/2.) * m_racah.threej(m_J2, 2*k, m_J2, -mj, 0, mj) * rme[k] * m_Bi[m] * m_econv;
132+
m_hamiltonian(i,i) += pow(-1., (m_J2-mj)/2.) * m_racah.threej(m_J2, 2*k, m_J2, -mj, 0, mj) * rme[k] * m_Bi[m] * m_econv;
92133
}
93134
}
94135
}
@@ -100,14 +141,14 @@ RowMatrixXcd cf1ion::hamiltonian(bool upper) {
100141
for (int i=0; i<(dimj-q); i++) {
101142
int mj = 2*i - m_J2, mjp = 2*(i+q) - m_J2;
102143
double tjp = m_racah.threej(m_J2, 2*k, m_J2, -mj, 2*q, mjp) - m_racah.threej(m_J2, 2*k, m_J2, -mj, -2*q, mjp);
103-
ham(i+q,i) += std::complex<double>(0., pow(-1., (m_J2-mj)/2.) * tjp * rme[k] * m_Bi[m] * m_econv);
144+
m_hamiltonian(i+q,i) += std::complex<double>(0., pow(-1., (m_J2-mj)/2.) * tjp * rme[k] * m_Bi[m] * m_econv);
104145
}
105146
}
106147
if (std::fabs(m_Bi[p]) > 1e-12) {
107148
for (int i=0; i<(dimj-q); i++) {
108149
int mj = 2*i - m_J2, mjp = 2*(i+q) - m_J2;
109150
double tjp = m_racah.threej(m_J2, 2*k, m_J2, -mj, 2*q, mjp) + m_racah.threej(m_J2, 2*k, m_J2, -mj, -2*q, mjp);
110-
ham(i+q,i) += pow(-1., (m_J2-mj)/2.) * tjp * rme[k] * m_Bi[p] * m_econv;
151+
m_hamiltonian(i+q,i) += pow(-1., (m_J2-mj)/2.) * tjp * rme[k] * m_Bi[p] * m_econv;
111152
}
112153
}
113154
}
@@ -116,19 +157,23 @@ RowMatrixXcd cf1ion::hamiltonian(bool upper) {
116157
if (upper) {
117158
for (int i=0; i<dimj; i++) {
118159
for (int j=i+1; j<dimj; j++) {
119-
ham(i,j) = std::conj(ham(j,i));
160+
m_hamiltonian(i,j) = std::conj(m_hamiltonian(j,i));
120161
}
121162
}
122163
}
123164

124-
return ham;
165+
m_ham_calc = true;
166+
return m_hamiltonian;
125167
}
126168

127169
std::tuple<RowMatrixXcd, VectorXd> cf1ion::eigensystem() {
128-
SelfAdjointEigenSolver<RowMatrixXcd> es(hamiltonian(false));
129-
RowMatrixXcd eigenvectors = es.eigenvectors();
130-
VectorXd eigenvalues = es.eigenvalues();
131-
return std::tuple<RowMatrixXcd, VectorXd>(eigenvectors, eigenvalues);
170+
if (!m_ev_calc) {
171+
SelfAdjointEigenSolver<RowMatrixXcd> es(hamiltonian(false));
172+
m_eigenvectors = es.eigenvectors();
173+
m_eigenvalues = es.eigenvalues();
174+
m_ev_calc = true;
175+
}
176+
return std::tuple<RowMatrixXcd, VectorXd>(m_eigenvectors, m_eigenvalues);
132177
}
133178

134179
} // namespace libMcPhase

‎singleion/cfpars.cpp

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -238,7 +238,6 @@ void cfpars::set_name(const std::string &ionname) {
238238
m_Bi[id] = m_Bo[id] / m_econv;
239239
}
240240
m_stevfact = {alpha, beta, gamma};
241-
m_invstevfact = {1./alpha, 1./beta, 1./gamma};
242241
const Map3 &rktable = RKTABLE();
243242
auto rk_table = rktable.find(ion);
244243
m_rk = rk_table->second;
@@ -256,7 +255,6 @@ void cfpars::set_J(const double J) {
256255
}
257256
m_ionname = "";
258257
m_stevfact = {1., 1., 1.};
259-
m_invstevfact = {1., 1., 1.};
260258
m_rk = {0., 0., 0.};
261259
m_convertible = false;
262260
m_convfact = lambda;
@@ -304,7 +302,6 @@ cfpars::cfpars(const std::string &ionname) {
304302
m_J2 = J2[n];
305303
m_convfact = lambda;
306304
m_stevfact = {alpha, beta, gamma};
307-
m_invstevfact = {1./alpha, 1./beta, 1./gamma};
308305
const Map3 &rktable = RKTABLE();
309306
auto rk_table = rktable.find(ion);
310307
m_rk = rk_table->second;

‎tests/test_cf1ion.py

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,9 @@ class cf1ionTests(unittest.TestCase):
1515
[ 0, 1.0034e+00, -3.8676e-01, -9.5181e-01, 7.3185e-01, 3.0767e-01, -8.2511e-01, 1.3833e+00, 1.8754e+00],
1616
[ 0, 0, 7.5851e-01, -7.7351e-01, 9.7844e-01, 1.0278e+00, -6.8087e-02, 1.8754e+00, -6.9491e-01]])
1717

18-
def test_creation(self):
18+
eval_ref = np.array([-4.0496e+00, -3.4501e+00, -1.8251e+00, -1.7653e+00, 7.7521e-01, 1.3800e+00, 1.6579e+00, 3.5307e+00, 3.7464e+00])
19+
20+
def test_creation_and_diagonalisation(self):
1921
with self.assertRaises(RuntimeError): # J must be integer or half-integer
2022
cfp = libMcPhase.cf1ion(2.3)
2123
with self.assertRaises(RuntimeError): # Unknown ion name (must be a rare earth 3+)
@@ -25,3 +27,6 @@ def test_creation(self):
2527
B40=-7.2952e-04, B41=6.4828e-03, B42=-3.3059e-03, B43=1.5821e-02, B44=-4.5648e-03,
2628
B60=-2.1369e-05, B61=4.3678e-04, B62=8.8995e-05, B63=1.0380e-04, B64=4.7701e-04, B65=-1.1485e-03, B66=3.9818e-04)
2729
nptest.assert_allclose(cfp.hamiltonian(), self.ham_ref, rtol=1e-4)
30+
V, E = cfp.eigensystem()
31+
nptest.assert_allclose(E, self.eval_ref, rtol=1e-4)
32+
nptest.assert_allclose(V @ np.diag(E) @ V.conj().T, cfp.hamiltonian(), atol=1e-4)

0 commit comments

Comments
 (0)