How to deal with sparse Hamiltonians ?

[devmessias]

Hi, I saw the tutorial about the Keppler problem.
https://github.com/SciML/SciMLTutorials.jl/blob/d5ece1515327f5931618e48c3a06eec42461e5ce/notebook/models/05-kepler_problem.ipynb

I'm tyring to study a Hamiltonian where $N~2500$ rotors interact in a periodic 2d grid,

I've implemented the Verlet method using the following work as reference
A Study of XY Model by Spin-dynamic simulation
Unfortunately, my code is really slow.

I’m trying to understand if it’s possible to attack this kind of problem using your package. It is?

I could'nt find any point in the documentation which could help me.

[ChrisRackauckas]

You'll want to use reverse mode AD (like Zygote.gradient or ReverseDiff.gradient) to define the dynamics, or use ModelingToolkit. For the latter, specify the Hamiltonian in ModelingToolkit (https://mtk.sciml.ai/dev/) and then calculate the sparse derivatives analytically and build a fast ODE function. Your issue is that forward mode AD for the gradients here won't scale.

[devmessias]

There is a reason to use AD in this case? Because we can get the equations of motions analytically.
For example, this is my implementation for the Verlet-algorithm applied in this kind of Hamiltonian

import cupy as cp

def cpgetTorque(thetas, A, J =1):

    z = thetas
    diffZsin = cp.sin(z.T-A*z)
    diffZsin = cp.multiply(diffZsin, A)
    torque = -J*diffZsin.sum(axis=1)
    return torque

def evolveThetas(thetas, momentum, torque, I, dt):

    nextThetas = thetas+momentum*dt + torque*dt**2/(2*I)
    return nextThetas

def evolveMomentum(momentum, torque, nextTorque, I, dt):
    
    nextMomentum  = momentum + (torque +nextTorque)*dt/(2*I)
    return nextMomentum

[ChrisRackauckas]

Oh, then just do a straight definition of the DynamicalODEProblem. That would be faster than using AD.

[devmessias]

Ah ok!
The equations are,

then I just to need to specify the port and qdot function in

prob = DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)

It's my first contact with Julia and your package. Thanks for your help.

[ChrisRackauckas]

Yup that would be the way to do it. Then use VelocityVerlet at first to make sure you got it right, then use the other methods to get higher order and such.