# 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
pos::Array{Float64,1}
k::Array{Float64,1}
c1::Number
c2::Number
surface::Int
m::Int
end
``````

and then I append methods like

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

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

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

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 += -force*dt/2
# Complete step position
p.pos += p.k*dt/p.m
# Half step momentum
force = der_surface(p)
p.k += -force*dt/2
end
``````

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 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.

4 Likes

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.

2 Likes

Welcome aboard!

On my reading of the code above, `k` and `pos` 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).