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[1])
else
return E2(p.pos[1])
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[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
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 