diff --git a/src/action.hpp b/src/action.hpp
index 1066e57835a2050cdc004897e4f90a527f406021..d4b5fcb46453e2e9a36ac7e53ffb3c1b396b8f12 100644
--- a/src/action.hpp
+++ b/src/action.hpp
@@ -83,7 +83,14 @@ inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
const int N_i[2] = {1, 0};
const int N_j[2] = {0, 1};
- /* Directions cycle (left->right, down->up, right->left, up->down) */
+ /*
+ * Directions cycle (left->right, down->up, right->left, up->down)
+ *
+ * U(i,j,k) = U_{(i2,j2),(i,j)} where i2 = i + N_i[k], j2 = j + N_i[k].
+ *
+ * Hence, spatial index of the link matrices U indicate the starting point,
+ * and the direction index the direction away from the starting point.
+ */
const size_t invdir[4] = {2, 3, 0, 1};
Scalar S{0.0};
@@ -108,8 +115,8 @@ inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
continue;
auto Q2 = Qmat.part(i2,j2).matrix();
- auto U12 = U.part(i,j,k).matrix();
- auto U21 = U.part(i2,j2,invdir[k]).matrix();
+ auto U21 = U.part(i,j,k).matrix();
+ auto U12 = U.part(i2,j2,invdir[k]).matrix();
auto M = nambu_diag_mul(U12, Q2, U21) - Q1;
diff --git a/src/core.cpp b/src/core.cpp
index 068c42124dc584a202650c6cc59f8ea75ea7f16d..2ce489a90a9d3a24aa4342a51dee058fba2deea9 100644
--- a/src/core.cpp
+++ b/src/core.cpp
@@ -69,8 +69,8 @@ private:
std::unordered_map<std::pair<size_t,size_t>, bool, pair_hash> hes_pattern_map_;
- CppAD::ADFun<Complex> J_f_; /**< Derivatives for currents */
- bool J_computed_;
+ CppAD::ADFun<Complex> A_f_; /**< Derivatives for currents */
+ bool A_computed_;
CppAD::ADFun<Complex> S_f_; /**< Derivatives for densities */
bool S_computed_;
@@ -91,7 +91,7 @@ private:
public:
Action(double Lx, double Ly, Complex alpha, Complex D, MaskType&& mask, UType&& U, OmegaType&& Omega)
: U_(std::move(U)), Omega_(std::move(Omega)), mask_(std::move(mask)), nx_(U_.dim(0)), ny_(U_.dim(1)), Lx_(Lx), Ly_(Ly), alpha_(alpha), D_(D), omega_(0),
- hes_pattern_map_(), J_computed_(false), S_computed_(false)
+ hes_pattern_map_(), A_computed_(false), S_computed_(false)
{
if (Omega.dim(0) != nx_ || Omega.dim(1) != ny_)
throw std::domain_error("U and Omega have incompatible shape");
@@ -310,9 +310,9 @@ public:
* Computing currents and densities.
*
* Current flowing to i is obtained as derivatives vs. xi with
- * U_{ij} -> exp(xi J) * U_{ij} =~ U_{ij} + (xi J) U_ij
- * where J is the generator matrix. We implement this as derivatives
- * vs. the elements of the generator J.
+ * U_{ij} -> exp(xi A(i)) * U_{ij} exp(xi A(j)) =~ U_{ij} + xi A(i) U_ij + xi U_ij A(j)
+ * where A is the generator matrix. We implement this as derivatives
+ * vs. the elements of the generator A.
*
* Densities are derivatives vs. the 0-component of the gauge
* potential, i.e., derivatives vs Omega. We implement this
@@ -348,27 +348,28 @@ public:
return dynamic;
}
- CppAD::ADFun<Complex> eval_dS_dJ(const QType& Q)
+ CppAD::ADFun<Complex> eval_dS_dA(const QType& Q)
{
if (Q.dim(0) != nx_ || Q.dim(1) != ny_)
throw std::out_of_range("Q array wrong size");
Array<ADComplex, Shape<Dynamic,Dynamic,4,4,4> >
- J({nx_, ny_, 4, 4, 4});
+ A({nx_, ny_, 4, 4, 4});
- J.fill(Complex(0));
+ A.fill(Complex(0));
std::vector<ADComplex> dynamic = dynamic_Omega_Q<ADComplex>(Q);
ArrayView<ADComplex, Shape<Dynamic,Dynamic,4,4> > Omegadyn(dynamic.data()+1, dynamic.size()-1, Omega_.shape());
ArrayView<ADComplex, Shape<Dynamic,Dynamic,2,2,2> > Qdyn(dynamic.data()+1+Omega_.size(), dynamic.size()-1-Omega_.size(), Q.shape());
- CppAD::Independent(J.storage(), dynamic);
+ CppAD::Independent(A.storage(), dynamic);
Array<ADComplex, Shape<Dynamic,Dynamic,4,4,4> >
U({nx_, ny_, 4, 4, 4});
const size_t invdir[4] = {2, 3, 0, 1};
const int neigh[4][2] = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};
+ const Complex sgn_A[4] = {1, 1, -1, -1};
for (size_t i = 0; i < nx_; ++i) {
for (size_t j = 0; j < ny_; ++j) {
@@ -381,8 +382,8 @@ public:
continue;
U.part(i,j,k).matrix() +=
- J.part(i,j,k).matrix() * U_.part(i,j,k).matrix()
- - U_.part(i,j,k).matrix() * J.part(i2,j2,invdir[k]).matrix();
+ sgn_A[k] * U_.part(i,j,k).matrix() * A.part(i,j,k).matrix()
+ + sgn_A[k] * A.part(i2,j2,invdir[k]).matrix() * U_.part(i,j,k).matrix();
}
}
}
@@ -390,7 +391,7 @@ public:
Array<ADComplex, Shape<1> > res;
res(0u) = S(mask_, Qdyn, dynamic[0], U, Omegadyn, Lx_, Ly_, alpha_, D_);
- return {J.storage(), res.storage()};
+ return {A.storage(), res.storage()};
}
CppAD::ADFun<Complex> eval_dS_dOmega(const QType& Q)
@@ -417,11 +418,11 @@ public:
return {Omega.storage(), res.storage()};
}
- void compute_dS_dJ(const QType& Q)
+ void compute_dS_dA(const QType& Q)
{
- J_f_ = eval_dS_dJ(Q);
- J_f_.optimize();
- J_computed_ = true;
+ A_f_ = eval_dS_dA(Q);
+ A_f_.optimize();
+ A_computed_ = true;
}
void compute_dS_dOmega(const QType& Q)
@@ -431,16 +432,16 @@ public:
S_computed_ = true;
}
- auto grad_J(const QType& Q)
+ auto grad_A(const QType& Q)
{
Vector<Complex> x(U_.size());
- if (!J_computed_)
- compute_dS_dJ(Q);
+ if (!A_computed_)
+ compute_dS_dA(Q);
else
- J_f_.new_dynamic(dynamic_Omega_Q<Complex>(Q));
+ A_f_.new_dynamic(dynamic_Omega_Q<Complex>(Q));
- auto jac = J_f_.Jacobian(x);
+ auto jac = A_f_.Jacobian(x);
return Array<Complex, Shape<Dynamic,Dynamic,4,4,4>, Vector<Complex> >(std::move(jac), {nx_, ny_, 4, 4, 4});
}
@@ -448,7 +449,7 @@ public:
{
Vector<Complex> x = ravel(Omega_);
- if (!J_computed_)
+ if (!S_computed_)
compute_dS_dOmega(Q);
else
S_f_.new_dynamic(dynamic_Q<Complex>(Q));
@@ -479,7 +480,7 @@ PYBIND11_MODULE(_core, m)
.def("compute", &Action::compute, "Precompute for AD", py::arg("Q"))
.def_property("omega", &Action::get_omega, &Action::set_omega)
.def_property_readonly("hess_computed", &Action::get_hess_computed)
- .def("grad_J", &Action::grad_J, "Evaluate gradient of action vs. current generators", py::arg("Q"))
+ .def("grad_A", &Action::grad_A, "Evaluate gradient of action vs. current generators", py::arg("Q"))
.def("grad_Omega", &Action::grad_Omega, "Evaluate gradient of action vs. Omega", py::arg("Q"))
;
diff --git a/tests/test_solver.py b/tests/test_solver.py
index ecbb70b1ef6a0852c9c97d5c338e18fe090d3909..212def5a09493ec168318b69725dbe8de4133d20 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -312,7 +312,8 @@ def test_soc_analytic(n=10):
Jtot = 0
for ww, aa in zip(w, a):
res0 = sol0.solve(omega=ww, finite_difference=False)
- J = sol._core.grad_J(res0.Phi.copy()).transpose(0, 1, 2, 4, 3) * (-1 / 8)
+ sol.omega = ww
+ J = sol._core.grad_A(res0.Phi.copy()).transpose(0, 1, 2, 4, 3) * (-1 / 8)
Jtot += aa * J
diff --git a/usadelndsoc/core.py b/usadelndsoc/core.py
index 7ab47fb2980e2a7972228fb9f96b8cb0d4981621..ff89b460333a21355c004d8c85a52bc910c035a4 100644
--- a/usadelndsoc/core.py
+++ b/usadelndsoc/core.py
@@ -105,9 +105,9 @@ class Core:
self._ensure_action()
return self._action.hess_mul(Q, dQ)
- def grad_J(self, Q):
+ def grad_A(self, Q):
self._ensure_action()
- return self._action.grad_J(Q)
+ return self._action.grad_A(Q)
def grad_Omega(self, Q):
self._ensure_action()
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index 780db9edb981cbffb50a65c48431ce52766352fd..22e933713df93e9e5727dab867ea3b86c8659d33 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -597,7 +597,7 @@ class Result:
@property
def J(self):
- return self._core.grad_J(self.Phi).transpose(0, 1, 2, 4, 3) * (-1 / 16)
+ return self._core.grad_A(self.Phi).transpose(0, 1, 2, 4, 3) * (-1 / 16)
@property
def S(self):