# Lammps interface to julia vs. python

**URL:** <https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164>\
**Category:** General Usage\
**Created:** [January 4, 2018, 7:49pm UTC](https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164 "2018-01-04T19:49:01Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)\
**Post date:** [January 4, 2018, 7:49pm UTC](https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164/1 "2018-01-04T19:49:01Z")

</div>

So LAMMPS ([http://lammps.sandia.gov/](http://lammps.sandia.gov/)) has a python interface which I have been trying to use within Julia, with some success. However, I had two points I wanted to query the community about.

First, when I wish to access the positions of the atoms within the environment, I run the following commands:

```julia
using PyCall
@pyimport lammps
@pyimport numpy.ctypeslib as ctlib
lmp = lammps.lammps()
lmp[:file]("setup.lammps") # raw lammps scripting set up file

# n = number of atoms, something I would set beforehand
Xobj = @pycall ctlib.as_array(lmp[:extract_atom]("x",3)[:contents],shape=(n,3))::Any
X = unsafe_wrap(Array{Float64,2}, Ptr{Float64}(Xobj[: __array_interface__]["data"][1]), reverse(Xobj[: __array_interface__]["shape"]))
Fobj = @pycall ctlib.as_array(lmp[:extract_atom]("f",3)[:contents],shape=(n,3))::Any
F = unsafe_wrap(Array{Float64,2}, Ptr{Float64}(Fobj[: __array_interface__]["data"][1]), reverse(Fobj[: __array_interface__]["shape"]))

```

At this point I have (3,n) Julia arrays which gives me access to the underlying positions and forces of the atoms, `X` and `F`. Clearly, it takes two commands for each such array. Is there a way to condense this down to single command? Is there a way to do this more efficiently?

Second, suppose I want to do a very elementary simulation, running the overdamped langevin problem, as

```julia
# set nsteps number of time steps, Dt the step size, and kBT the temperature
for j =1:nsteps
 
   X.+= Dt.*F .+ sqrt(2. * kBT * Dt) .* randn(3,n)
    lmp[:command]("run 1 pre yes post no") # forces the lammps environment to update forces as atoms move
# other commands go here
end

```

While this works, I find that the analogous code within python is about twice as fast. Are there any obvious mistakes and/or ways to get back the Python speed?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 4, 2018, 8:08pm UTC](https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164/2 "2018-01-04T20:08:03Z")

</div>

> [@gideonsimpson](#):
>
> ```julia
> # set nsteps number of time steps, Dt the step size, and kBT the temperature
> for j =1:nsteps
> 
> X.+= Dt.*F .+ sqrt(2. * kBT * Dt) .* randn(3,n)
> lmp[:command]("run 1 pre yes post no") # forces the lammps environment to update forces as atoms move
> # other commands go here
> end
> 
> ```
> 
> While this works, I find that the analogous code within python is about twice as fast. Are there any obvious mistakes and/or ways to get back the Python speed?

Is the limiting factor the `X .+=` line or the `lmp` call? If it’s the latter, I don’t have much to say since I don’t know anything about LAMMPS. The `X .+=` line can probably be made faster, though, in two ways:

- Don’t run performance-critical code in global scope (or using global variables): put your performance-critical code into a function. (It’s not clear whether you are doing this already?) e.g. write a function `timestep(X, F, kBT, Dt, nsteps, lmp)` that performs your loop.

- Your `randn` call is allocating a temporary array on each call. You can instead do `X.+= Dt.*F .+ sqrt(2. * kBT * Dt) .* randn.()`, or just set `c = sqrt(2kBT * Dt)` and do `@. X += Dt*F + c*randn()` to omit all of the dots.

Rather than `X = unsafe_wrap(Array{Float64,2}, Ptr{Float64}(..)`, by the way, have you tried just using:

```julia
X = @pycall ctlib.as_array(lmp[:extract_atom]("x",3)[:contents],shape=(n,3))::PyArray

```

since `PyArray` is a copy-free abstract-array wrapper around a NumPy array? (Possibly PyArray could still use some performance tuning, though, since there have been a lot of array updates since it was first implemented.)

---

<div class="post-metadata">

**Author:** ![gideonsimpson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gideonsimpson/32/1928_2.png) [@gideonsimpson](https://discourse.julialang.org/u/gideonsimpson)\
**Post date:** [January 4, 2018, 9:17pm UTC](https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164/3 "2018-01-04T21:17:37Z")

</div>

The suggestions did seem to give a bit of speedup, so I’ll keep playing with that.

Regarding your suggestion about using

```julia
X = @pycall ctlib.as_array(lmp[:extract_atom]("x",3)[:contents],shape=(n,3))::PyArray

```

While this certainly works, it introduces a new question. The output is now an (n,3) array. Suppose I have a problem that I wish to constrain to 2D by ignoring the 3rd coordinate. I do not seem to have the ability to call X[:,1:2], like I would with a regular julia array. Is this kind of indexing not possible with PyArray’s?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 4, 2018, 10:36pm UTC](https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164/4 "2018-01-04T22:36:24Z")

</div>

> [@gideonsimpson](#):
>
> I do not seem to have the ability to call X[:,1:2], like I would with a regular julia array. Is this kind of indexing not possible with PyArray’s?

It is possible, but slice-indexing hasn’t been implemented yet for `PyArray`. However, realize that expressions like ` X[:,1:2]` make a copy, even with `Array`, unlike NumPy. If you want slicing to create a view, use e.g. `@view(X[:,1:2])`, and this _will_ work with PyArray since it just wraps a `SubArray` around it.

---

<div class="post-metadata">

**Author:** ![ScottPJones](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/scottpjones/32/146_2.png) [@ScottPJones](https://discourse.julialang.org/u/ScottPJones)\
**Post date:** [January 5, 2018, 1:30am UTC](https://discourse.julialang.org/t/lammps-interface-to-julia-vs-python/8164/5 "2018-01-05T01:30:21Z")

</div>

Probably by making a direct Julia wrapper (at which point it might be seriously faster than using the Python wrapper in Python).

If there is a decent C API, that’s pretty easy to do, otherwise, there are a couple of ways of wrapping C++ code directly, such as [GitHub - JuliaInterop/Cxx.jl: The Julia C++ Interface](https://github.com/Keno/Cxx.jl) and [GitHub - JuliaInterop/CxxWrap.jl: Package to make C++ libraries available in Julia](https://github.com/JuliaInterop/CxxWrap.jl).