Symplectic midpoint solver for spin systems on a product of spheres.
Switch branches/tags
Nothing to show
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Type Name Latest commit message Commit time
Failed to load latest commit information.

spinsys - Spherical midpoint for spin systems

Build Status Coverage Status ![Python version](, 3.4-blue.svg?style=flat-square)


This is a symplectic midpoint solver for spin systems. These systems are symplectic differential equations derived from the Lie–Poisson structure of a product of spheres. The papers "Symplectic integrators for spin systems", "A minimal-coordinate symplectic integrator on spheres" and "Geometry of discrete-time spin systems" gives the details of the setting and the method, as well as a theoretical background.

The method and the implementation are a joint work of R. I. McLachlan, K. Modin and O. Verdier.

How to use it

What you need

  • m0: An initial condition: an Nx3 matrix where each row is a vector of length one.
  • gradient: A gradient which returns the gradient of a Hamiltonian at each such points.

You may also optionally need

  • strengths: A N vector of strength which determines the symplectic form on the product of spheres.

If such a strenghts vector is not provided, the strengths are assumed to be one for every sphere.

Run the simulation

The code is then as follows:

import midpoint
import sphere

dt = .01 # time step
nb_steps = 10 # number of steps

vector = sphere.radially_constant(gradient, strengths)
generator =, m0=m0.ravel(), dt=dt, nb_steps=nb_steps)

# run the solver
ms = [m.reshape(-1,3) for m in generator]

Now the solution ms is a list of Nx3 matrices which correspond to the solutions at each time step.

Available systems

There are two common spin systems which are already available in this package.

Spin chain

The Heisenberg spin chain system for N points is initialized as follows:

from spinsys import spinchain
gradient = spinchain.get_gradient(1/N**2)

The strength vector may be omitted in the call of radially_constant:

vector = sphere.radially_constant(gradient)

Point vortices

The gradient is obtained from the strengths vector:

from spinsys import vortex
import numpy as np

strengths = np.array([1., 1., -1.])
gradient = vortex.get_gradient(strengths)