Best implementation for semiclassical dynamics in Julia

Dear all,

I’m a Julia user since 2017 and after getting familiar with multiple dispatch and mutable structures, I finally found the courage to interact for the first time with the Julia community. Here’s my question.

During my PhD I need to implement a code for nonadiabatic semi-classical dynamics. At the moment I have already have two functioning codes, one is much faster than the other one but is much less easily maintainable.

The first code is essentially made by a series of functions that takes a series of inputs (mass, velocity, force, time, position, complex coefficients of the wavefunction, Hamiltonian matrix, ecc…) and essentially evaluates the time evolution of my pseudo-particle. The code is much faster than the python equivalent that I implemented several months ago, but is generally difficult to extend and maintain.

The other implementation instead defines a mutable structure pseudo particle and the time evolution is performed via all the methods of the structure. Below you can find a minimalistic example:

mutable struct Particle

and then I append methods like

function kinetic_energy(p::Particle)
    return sum(p.k.*p.k)/2/p.m

function potential_energy(p::Particle)
    if p.surface == 1
        return E1(p.pos[1])
        return E2(p.pos[1])

function total_energy(p::Particle)
    return potential_energy(p)+kinetic_energy(p)

function time_evolution(p::Particle, dt::Real)
    # This method defines the time evolution of a particle and
    # update all the variables accordingly.
    # Position is approximated at the first order,
    # momentum is approximated at the second order.
    # ---------------------------------------------------------
    # Half step momentum
    force = der_surface(p)
    p.k[1] += -force*dt/2
    # Complete step position
    p.pos[1] += p.k[1]*dt/p.m
    # Half step momentum
    force = der_surface(p)
    p.k[1] += -force*dt/2

This second implementation is much more readable (IMO also more elegant) but also much slower.

Now my question is: considering that Julia is not a object oriented language, is it “in the Julia spirit” to continue trying to optimize the second implementation (probably there’s some methods that are not really numerically efficient) or go back to the first implementation?

Every help is appreciated :wink:

1 Like

Go ahead: post some MWE (or the entire code). I am sure you will get plenty of tips on how to make it fly.

1 Like

c1::Number for instance is an abstract type. Make it concrete, it will help the compiler to optimize.


The only clear “bad” thing I can see on your code is using Number for some fields, which is an abstract type.

From the rest it is not obvious what your “fast” option could do different. If you can provide a more explicit example that would help.

Ps. Don’t be afraid of asking questions! I have done and still do all the stupid questions possible. Asking and answering is the best way to learn. And this forum is really, really great.


Welcome aboard!

On my reading of the code above, k[1] and pos[1] are the only bits of p that get modifed. If I was optimising this, I would delete k[2:end] and pos[2:end], which the time evolution code never touches. That reduces the problem to solving F=ma for a single classical particle in 1 dimension, which DifferentialEquations can do in the blink of an eye.

What are you really doing, and how does quantum mechanics enter into it?

Broadly, given that pos and k are likely to be small vectors, I would use StaticArrays.SVector and code it in a functional style (not mutating existing object, but making new ones).