# Best implementation for semiclassical dynamics in Julia

**URL:** <https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423>\
**Category:** Performance\
**Tags:** question\
**Created:** [March 17, 2021, 10:51pm UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423 "2021-03-17T22:51:05Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![ErikPillon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erikpillon/32/4068_2.png) [@ErikPillon](https://discourse.julialang.org/u/ErikPillon)\
**Post date:** [March 17, 2021, 10:51pm UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423/1 "2021-03-17T22:51:06Z")

</div>

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:

```julia
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

```julia
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 😉

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [March 17, 2021, 11:13pm UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423/2 "2021-03-17T23:13:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [March 17, 2021, 11:15pm UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423/3 "2021-03-17T23:15:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 17, 2021, 11:16pm UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423/4 "2021-03-17T23:16:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![thisrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thisrod/32/10642_2.png) [@thisrod](https://discourse.julialang.org/u/thisrod)\
**Post date:** [March 18, 2021, 12:22am UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423/5 "2021-03-18T00:22:28Z")

</div>

Welcome aboard!

> [@ErikPillon](#):
>
> nonadiabatic semi-classical dynamics

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?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 18, 2021, 10:20am UTC](https://discourse.julialang.org/t/best-implementation-for-semiclassical-dynamics-in-julia/57423/6 "2021-03-18T10:20:53Z")

</div>

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