Hi, I'm trying to add a biquadratic term of the form K*(S_i . S_j)^2 to the Hamiltonian of a classical spin system consisting of Heisenberg spins. Also, I'm working on a frustrated lattice, so I'm only using local updates.

My thinking at the moment is that I need to primarly make changes to the functions "bond_energy" and "bond_energy_change" in the on.h file. However, that leaves me with two questions.

1) What would be the best (i.e. least obtrusive) way to pass the parameter K into those functions? At the moment they only seem to be able to handle two-spin interactions. Could I perhaps make the code in the on.h file aware of the simulation parameters directly (i.e. such that the code double K = parms["K"] works)? Or, should I make K a field/property of the ONMoment class? Or some other approach (function overloading, then passing K in through a call in localupdate.h)?

2) Could someone please explain what in fact the function bond_energy_change is doing? Here is the function definition: template <class M> inline double bond_energy_change(const ONMoment<M::dim>& m1, const ONMoment<M::dim>& m2, M J, typename ONMoment<M::dim>::update_type dir) { TinyVector<double,M::dim> tmp = dir*dot(dir,m1.state_); return 2.*J.vec_mat_vec(tmp, m2.state_); } I am confused about the above definition in two ways. One is that ONMoment is initialzed as ONMoment<3> for the Heisenberg model, yet here the template calls for a class "M" and refers to "M::dim", though surely an int (i.e. 3) doesn't have a property "dim". Also, my understanding is that the argument "dir" in the above code is a random unit vector and independent of m1.state_. If so, how exactly does the above function give the energy change? Wouldn't the energy change be equal to J.vec_mat_vec((dir - m1.state_),m2.state_) ? That is, whence the factor of 2 and the definition of tmp?

Thanks!