Conversation
* Had to modify the antiferro_chain.py model * Fix missing 2pi in spin-spin phase factor * Fix use of lower triangular factor in decomposition We use LDL whereas Matlab code uses Cholesky and the inv() of upper or lower triangles are v. different; * TODO: investigate why maths only work for u-tri
* Change to using adjoint of T matrix instead of just transpose. Now numbers agree exactly with Matlab (previously off by 5%) * Remove sorting in Python code - this swapped the order of the energies (eigenvalues) and intensities (eigenvectors) in subsequent calculations and is incorrect - it was only used because the Python was a copy of the rust code which used a -1 multiply but we can just take the abs for the sqrt(E) calcs.
lucas-wilkins
left a comment
lucas-wilkins
left a comment
There was a problem hiding this comment.
There's some more work to be done on how the q points scale with the supercell, and I'm not 100% sure things will work right treating as you have here. I think we need a layer of converting q in and out of RLU, or things will get very confusing. Not saying that should be done here, in fact, I have previously made an issue for this #102 .
| unique_id_to_index[input_coupling.site_1._unique_id], | ||
| np.array(input_coupling.coupling_matrix.T, **rust_kw), | ||
| -input_coupling.vector(expanded.structure.unit_cell) | ||
| -input_coupling.cell_offset.vector.astype('double') |
There was a problem hiding this comment.
This has changed from cartesian to LU coordinates, were you aware of this? is this the bugfix?
There was a problem hiding this comment.
is this the bugfix?
Yes. I think that we need this in lattice units rather than Cartesian... I'm not 100% sure but in the original code many of the examples fail with nonsensical energies...
|
I'm still very skeptical about this. Fundamentally, I think the calculation should be invariant to the choice of unit cell, and I don't think this solution achieves that, even after correcting for the supercell scaling. |
* Include double counting factor of 1/2 in calculation of reverse couplings * Fix cartesian transformation bug in lattice_distances (no transpose) * Revert limits calculations (-l,l) vs (0,l) in lattice_distances * Allow LaticeSite "moments" to be basisstate (can be complex) This requires taking the real part in supercell.expand()
* Update antiferromagnetic_chain, af_with_field, kagome_ferromagnet * Add kagome_antiferromagnet (tutorial 7) * Add triangular_antiferromagnet (tutorial 12) * Add support for LU (lattice units) units in LatticeSite
* Change rust to make sqrt_hamiltonian mutable and remove kk * Rename variables in hamiltonian.py
* Remove StrEnum as not in Python 3.10 * Fix hamiltonian.py code as called with unnormed transformation * Propagate LatticeSite unit in Structure * Remove old incorrect lattice_distance comment * Fix af_with_field example coupling distance as changed lattice pars
* Previous rust code was not adding delta to diagonal * Handles changes Py3.11 <-> 3.12 in numpy where the LinAlgError was not triggered giving a matrix with NaNs (in 3.11) * Change limits of ferromagnetic chain example to exclude (100) where the intensity diverges and rust and Python disagree
Refactor to remove lu unit from LatticeSite class Add spin length/magnitude to LatticeSite class this means that mi,mj,mk can just indicate direction consistent with Matlab code Update examples code to use genmagstr
Fix matlab docstring code in examples Remove unneeded imports Remove debug_plot Plot only positive energies in some examples to make it clearer
* Remove magnitude from LatticeSite / corresponding Structures * Refactor genmagstr to define magnitude in its own args * Add complex numpy array serialisation functionality
|
@lucas-wilkins I think this is ready for review... There is still a scale factor discrepancy in the intensity calculation but there's so much in this already I want to merge it before doing another PR for that fix (which I would need to investigate in more detail but don't have much time this week). |
alexhroom
left a comment
alexhroom
left a comment
There was a problem hiding this comment.
just reviewing the Rust as you asked - mostly quality code, just a few little bits here and there for performance
alexhroom
left a comment
alexhroom
left a comment
There was a problem hiding this comment.
happy with the Rust now!
3 participants
This PR ensures that output from PySpinW agrees with Matlab SpinW numerically. There were a couple of errors in the computations:
inter_site_vectoris actually the vector between unit cells (G) rather than the full vector between the sites.adjoint(L) \ U * sqrt(E)rather thanL \ U * sqrt(E)- in Matlab they use a Cholesky decomposition which yields an upper-triangularKwhereas we use the LDL decomposition which yields a lower triangularLbut the computation only works for an upper-triangular matrix.hamiltonian.pythe calculation of reverse couplings was missing the double counting factor of1/2e1,e2,e3hade2parallel to[0,1,0]ife3is parallel to[1,0,0]but the Matlab code had it parallel to[0,0,1]- and this caused thekagome_supercellcalculation to fail (not sure why).In addition, I changed the driver routines to correct the above mistakes; and also added more examples from the Matlab test suite.
However, there are two major discrepancies between the theory presented in Toth & Lake and the Matlab implementation:
2in the Matlab Hamiltonian calculation of the single-ion anisotropy terms (calledA20/D20in the code, computed in lines 618-9 compared to the1/(2*prod(nExt))(factor of 1/2 is in line 1066) which is the number of unit cells in the supercell - e.g. the Matlab code computes the intensity per unit cell whereas the theory gives it per magnetic atom.In this PR I've added code to make the Python calculations consistent with the Matlab, but we might want to make it such that it agrees with the paper instead.