# Electrostatics in Julia

**URL:** <https://discourse.julialang.org/t/electrostatics-in-julia/41633>\
**Category:** Modelling & Simulations\
**Created:** [June 17, 2020, 10:14pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633 "2020-06-17T22:14:19Z")\
**Posts on this page:** 13\
**Page:** 2

<div class="post-metadata">

**Author:** ![jgreener64](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgreener64/32/2483_2.png) [@jgreener64](https://discourse.julialang.org/u/jgreener64)\
**Post date:** [June 21, 2020, 11:13am UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/21 "2020-06-21T11:13:23Z")

</div>

> [@cortner](#):
>
> Off topic: Maybe molly would like to join JuliaMolSim? I’m still hoping we can make Molly and JuLIP compatible.

Yes, that seems its natural home, let’s transfer it there.

I think making me a member of the org will let me make the transfer, could someone do that? Thanks.

---

<div class="post-metadata">

**Author:** ![cortner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cortner/32/204_2.png) [@cortner](https://discourse.julialang.org/u/cortner)\
**Post date:** [June 21, 2020, 1:40pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/22 "2020-06-21T13:40:15Z")

</div>

I’ve just given you “owner” privileges.

Note - the purpose of this organisation is two-fold. (1) just collect packages relevant to molecular simulation; (2) manage our own registry [MolSim](https://github.com/JuliaMolSim/MolSim) where we register our packages essentially without and constraints and quality-checks other than our own, so that we can have version control for our research projects. (This doesn’t prevent us from registering polished version with General as well.)

---

<div class="post-metadata">

**Author:** ![jgreener64](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgreener64/32/2483_2.png) [@jgreener64](https://discourse.julialang.org/u/jgreener64)\
**Post date:** [June 21, 2020, 3:34pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/23 "2020-06-21T15:34:18Z")

</div>

Great, I’ve made the switch.

---

<div class="post-metadata">

**Author:** ![Cory\_Simon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cory_simon/32/3950_2.png) [@Cory\_Simon](https://discourse.julialang.org/u/Cory_Simon)\
**Post date:** [November 28, 2020, 3:28am UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/24 "2020-11-28T03:28:05Z")

</div>

I wrote an Ewald sums code here:  
[https://github.com/SimonEnsemble/PorousMaterials.jl/blob/master/src/electrostatics.jl](https://github.com/SimonEnsemble/PorousMaterials.jl/blob/master/src/electrostatics.jl)  
despite my best efforts to speed it up (resulting in very ugly code), it is too slow for production. also the default choices for the number of replications in k-space are not optimal. I tested it with some NIST data (see [https://github.com/SimonEnsemble/PorousMaterials.jl/blob/master/test/electrostatics.jl](https://github.com/SimonEnsemble/PorousMaterials.jl/blob/master/test/electrostatics.jl)). feel free to copy any of my code.  
I would love to both use an EWalds.jl package and contribute to one.  
we have a package for crystal structures here:  
[https://github.com/SimonEnsemble/Xtals.jl](https://github.com/SimonEnsemble/Xtals.jl)  
and we simulate gas adsorption in MOFs (GCMC sims) with PorousMaterials.jl for our research:  
[https://github.com/SimonEnsemble/PorousMaterials.jl](https://github.com/SimonEnsemble/PorousMaterials.jl)  
(only handles LJ potentials and point charges).

Cory

---

<div class="post-metadata">

**Author:** ![mfh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfh/32/13787_2.png) [@mfh](https://discourse.julialang.org/u/mfh)\
**Post date:** [November 28, 2020, 8:45am UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/25 "2020-11-28T08:45:25Z")

</div>

A bit off topic, but regarding Xtals.jl: I’ve been thinking about something similar as a light-weight ASE to provide more convenience for starting calculations in DFTK. Is that the intention you have in mind with the package?

---

<div class="post-metadata">

**Author:** ![Cory\_Simon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cory_simon/32/3950_2.png) [@Cory\_Simon](https://discourse.julialang.org/u/Cory_Simon)\
**Post date:** [November 28, 2020, 6:52pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/26 "2020-11-28T18:52:31Z")

</div>

we built the package over the last three years, out of pure selfishness-- whatever we need to get our research done in Julia!  
that would be awesome if Xtals.jl turns into a more general purpose package for crystals, PRs welcome!  
we have the capability to infer bonding graphs and work with those (another package on top of Xtals.jl coming soon). the only qualification is that there are speed considerations for how to represent atoms and charges for molecular simulations, in PorousMaterials.jl. so I wouldn’t want those to change, unless they were faster/as fast.  
plz try out the package and let us know what you think!

---

<div class="post-metadata">

**Author:** ![mfh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfh/32/13787_2.png) [@mfh](https://discourse.julialang.org/u/mfh)\
**Post date:** [November 29, 2020, 7:38am UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/27 "2020-11-29T07:38:07Z")

</div>

Cool. Thanks. Will do 😉

---

<div class="post-metadata">

**Author:** ![cortner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cortner/32/204_2.png) [@cortner](https://discourse.julialang.org/u/cortner)\
**Post date:** [November 30, 2020, 7:02pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/28 "2020-11-30T19:02:45Z")

</div>

I’ve got some functionality for that in JuLIP and wouldn’t mind at all getting rid of it and using another oackage

---

<div class="post-metadata">

**Author:** ![pxshen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pxshen/32/38776_2.png) [@pxshen](https://discourse.julialang.org/u/pxshen)\
**Post date:** [March 2, 2021, 6:19pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/29 "2021-03-02T18:19:39Z")

</div>

Bumping this thread. Any updates?  
I might implement a barebones particle mesh since it’s just couple lines of fourier transform. Probably similar to Michael’s DFTK.jl ewald.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 2, 2021, 7:59pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/30 "2021-03-02T19:59:21Z")

</div>

The situation hasn’t changed from my end. If you or somebody else want to register a package and maintain it we’d be happy to have DFTK depend on it, and I would imagine a number of people are in the same situation.

---

<div class="post-metadata">

**Author:** ![pxshen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pxshen/32/38776_2.png) [@pxshen](https://discourse.julialang.org/u/pxshen)\
**Post date:** [March 4, 2021, 4:23pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/31 "2021-03-04T16:23:13Z")

</div>

Thanks I’ll take a look.

---

<div class="post-metadata">

**Author:** ![Cory\_Simon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cory_simon/32/3950_2.png) [@Cory\_Simon](https://discourse.julialang.org/u/Cory_Simon)\
**Post date:** [March 12, 2021, 4:11am UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/32 "2021-03-12T04:11:34Z")

</div>

getting serious about this. we do molecular simulations of gases adsorbing in porous crystals such as MOFs with [PorousMaterials.jl](https://github.com/SimonEnsemble/PorousMaterials.jl). need to compute electrostatic potential inside the MOF _very_ quickly to enable simulation of molecular models with point charges on them.

to start a detailed discussion:

# goal

a simple function to compute electrostatic potential at a point given in fractional coords `xf`, arising from `charges` in a periodic `box`.

```julia-auto
phi = electrostatic_potential(xf::Array{Float64, 1}, charges::Charges, box::Box, ewald_settings::Settings)

```

is this what you are all looking for too?

but we need to call the above function hundreds of thousands of times, involving many loops over ~100 atoms—same `charges` different `xf`’s. for this function to be of reasonable speed for us, we need to pre-compute (and re-use) the k-vectors and pre-allocate memory for arrays that are used over and over in the for loop over the charges to compute e^{i n k ⋅ r}.  
my code for this is [here](https://github.com/SimonEnsemble/PorousMaterials.jl/blob/master/src/electrostatics.jl)– very ugly. also has several tricks that take advantage of (a) symmetry to halve the computations and (b) recurrence relations to speed up e^{i n k ⋅ r}.

anyway, I propose the following for building blocks for this.

# how to store the point charges

this is based on our `Xtals.jl` package [here](https://github.com/SimonEnsemble/Xtals.jl). can we use `Xtals.jl` or is there a more flexible way to represent charges for the community?

- a `Box` to tell us about the periodicity of the system. see [here](https://github.com/SimonEnsemble/Xtals.jl/blob/master/src/box.jl#L1-L37).

- a `Charges` structure to tell us where the point charges are in the box and the magnitude of the charges. we use [this](https://github.com/SimonEnsemble/Xtals.jl/blob/master/src/matter.jl#L174-L178) because we distinguish between fractional and cartesian coordinates in our code to avoid bugs (treat them as separate types). maybe not needed for ya’ll? if not, we can do more general:

```julia-auto
struct Charges
    q::Array{Float64, 1} # values of the point charges
    xf::Array{Float64, 2} # fractional coordinates of charges in the box
end

```

and maybe put info about the box in there too.

in practice I found the above much faster than storing charges as `charges::Array{Charge, 1}` with

```julia-auto
struct Charge
   q::Float64
   xf::Array{Float64, 1}
end

```

since the former allows for quick array operations through `broadcasting` (memory locality thing?).

for an `Ewald.jl` would we start over in defining these basic building blocks or use `Xtals.jl`?

# Offsetarrays

I found OffsetArrays.jl to be very helpful for making the code match the math since it allows `for k = -2:2` and to index arrays as `x[k]`. would this give a performance hit, though?

# resources

Allen & Tildesley give a [sample Ewald code](https://github.com/glennklockwood/allen-tildesley/blob/master/ewald.f) that I found helpful.

these [notes](http://micro.stanford.edu/mediawiki/images/4/46/Ewald_notes.pdf) for understanding Ewald.

---

<div class="post-metadata">

**Author:** ![jgreener64](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jgreener64/32/2483_2.png) [@jgreener64](https://discourse.julialang.org/u/jgreener64)\
**Post date:** [November 3, 2025, 5:38pm UTC](https://discourse.julialang.org/t/electrostatics-in-julia/41633/33 "2025-11-03T17:38:40Z")

</div>

[Molly.jl](https://github.com/JuliaMolSim/Molly.jl) now has PME code, compatible with CPU and GPU backends via [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl): [Molly.jl/src/interactions/ewald.jl at master · JuliaMolSim/Molly.jl · GitHub](https://github.com/JuliaMolSim/Molly.jl/blob/master/src/interactions/ewald.jl).

The implementation is based on OpenMM and seems to be decently performant, taking \< 1 ms per step for 15k atoms on a RTX 4090 GPU. It is differentiable with Enzyme.

It is rather tied into the Molly infrastructure since the neighbour list is used to calculate the short-range contributions. Here is an example:

```julia
using Molly # Requires v0.23.0
using GLMakie

n_atoms = 100
atoms = [Atom(mass=10.0u"g/mol", charge=(i % 2 == 0 ? 0.3 : -0.3))
         for i in 1:n_atoms]

boundary = CubicBoundary(2.0u"nm")
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm")

dist_cutoff = 0.9u"nm"
coul = CoulombEwald(dist_cutoff=dist_cutoff)
pme = PME(dist_cutoff, atoms, boundary)

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    pairwise_inters=(coul,),
    general_inters=(pme,),
    loggers=(coords=CoordinatesLogger(10),),
)

temp = 1.0u"K"
simulator = VelocityVerlet(
    dt=0.001u"ps",
    coupling=AndersenThermostat(temp, 1.0u"ps"),
)

random_velocities!(sys, temp)
simulate!(sys, simulator, 2_000)

visualize(sys.loggers.coords, boundary, "sim_pme.mp4";
          color=[(i % 2 == 0 ? :red : :blue) for i in 1:n_atoms])

```

[Previous page](https://discourse.julialang.org/t/electrostatics-in-julia/41633.md?page=1)