Outputting the electrostatic potential field in LAMMPS


The spatial distribution of the electrostatic potential can be computed from a LAMMPS simulation relatively easily using the atc user package in LAMMPS, but this does not appear to be well documented, so here I will explain how I got it to work. This method can also be used to output other continuum properties such as charge density, temperature, etc. and thus could be more broadly useful as well.

(Update: For 1D cases, you may also want to consider the alternative method described in this blog post.)

First, it is necessary to install the user-atc package. See https://lammps.sandia.gov/doc/Build_extras.html#user-atc. In the src/ directory of my LAMMPS source code directory used to build LAMMPS, I executed:

make libatc args=“-m mpi”

Since I had already manually installed BLAS and LAPACK for other purposes, I then modified lib/atc/Makefile.lammps by changing the following variables (you should modify the paths to your LAPACK and BLAS directories appropriately):

user-atc_SYSLIB = -llapack -lblas -lgfortran
user-atc_SYSPATH =  -L/home/alta/software/lapack-3.8.0

Alternatively, in the src/ directory, execute:

make liblinalg args=“-m mpi”

and don’t change lib/atc/Makefile.lammps, which should get atc to use the built-in linear algebra libraries.

Then, in the src/ directory, execute:

make yes-manybody
make yes-user-atc
make mpi

This produces the executable lmp_mpi which is used to run LAMMPS.

Next, edit your LAMMPS input script. For further references, see the fix atc documentation and examples in the examples/USER/atc/ directory of the LAMMPS source code. For example, to compute a 2D spatial map of electrostatic potential in the xz plane:

units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls

pair_style lj/cut/coul/long 12.0 12.0
pair_modify mix geometric
kspace_style pppm 1.0e-4

read_data data.initial

variable NMESHX equal 100
variable NMESHY equal 1
variable NMESHZ equal 50
lattice sc 1 # Necessary to define lattice for ATC
fix ATC all atc field
fix_modify ATC mesh create ${NMESHX} ${NMESHY} ${NMESHZ} BOX p p p
fix_modify ATC atom_element_map eulerian 1
fix_modify ATC fields add electric_potential
fix_modify ATC output out.atc_potential 1 text

rerun dump.lammpstrj dump x y z q

This will output the electric potential on a uniform 100 x 1 x 50 grid on your simulation box every 1 timestep. This can be used with either the run command or rerun command (for example, rerun on a dump file of a completed simulation). Here, the dimension of 1 in the y direction means values are averaged over the y direction, allowing a 2D field to be generated. You can also set for example NMESHX to be 1, which then generates a 1D profile averaged over the y and z directions.

After running this, there will be 3 output text files: out.atc_potential.DATA, out.atc_potential.MAP, and out.atc_potential.XYZ. The main results are in the .DATA file, where the last columns contain the values of the outputted fields. These files can then be parsed using for example a python script, for analysis and plotting.

Miscellaneous tips:

  • Make sure the mesh create boundary conditions match your box boundary conditions! If for example your box is p p f, make sure mesh also uses p p f.
  • I have only been able to get this to work properly for the rerun command and not with run (i.e. during a simulation while it’s running, only as post-processing).
  • As always, double check that your results make sense (periodic systems give periodic results, etc.). For 1D systems, you can check your results against manually integrating the charge density profile twice (which may actually be a preferable strategy to this fix atc method, for 1D at least).
  • Be sure to average your data over a sufficient amount of time because electrostatic potential is a value that can fluctuate quite considerably.


  1. I was looking for a way to do this for a while now. Thanks very much! Also, I believe we met a number of years ago when I was an REU student at Princeton. Hope you are doing well.


    • Hi Jon, glad to hear that this was helpful! Yes, I remember meeting you at Princeton! You did some really excellent work with Minnie Liu and Craig Arnold. Hope you are doing well too!


      • Do you mind describing how you produced the first plot? I am trying to make a similar planar averaged potential in a capacitor cell. What ATC settings did you use?


  2. Hi Jon, Actually, for the 1D potential profiles such as that plot, I typically did not use fix atc but integrated the charge density profile twice. I’ll send you an email so I can attach a Python script, and perhaps later I’ll update this blog post with more information as well, since someone else has also asked me a similar question after seeing this blog post.


  3. hey mate, thanks- this is easy to understand compared to some other docs. Do you have some insight into how to connect boundary conditions to this too (i.e Dirichlet)??


  4. Nice guide, thank you! Do you know of any software to visualize the volumetric data produced by this fix? Or to convert the output to known formats, e.g., CUBE.


    • Thanks, glad you liked it!

      1. I cannot remember the exact version I used any more. I did this work more than 2 years ago, so probably a version from around 2018 or 2019. Although electric_potential may not be currently listed in that documentation link, it should still be available in newer versions. To check what fields are available through “fix atc”, you can look at the source code on github: https://github.com/lammps/lammps/blob/master/lib/atc/ATC_TypeDefs.h . This can be found by searching for “electric_potential” in the LAMMPS github, which also can be a useful way to find examples: https://github.com/lammps/lammps/search?q=electric_potential.
      2. Since I did this work a few years ago, I can’t remember off the top of my head, but my guess is that it would depend on what units you are using for your LAMMPS simulation overall (https://docs.lammps.org/units.html). It should also be possible to double check the units by calculating the electrostatic potential in a different way and comparing/figuring it out from that (for example, see https://altafang.com/2020/10/13/calculating-1d-electrostatic-potential-profile-from-md-simulations/). I can look back at my old files this weekend to see if I can give a better answer to this question for the “units real” case.


      • Thanks a lot. Although I will ask LAMMPS developers to add the description for the electric_potential field, your information for the unit of the electrostatic potential of this case is also helpful for me and might be also helpful for someone else later.


      • I just looked back at my old files, and I think the units of the electric_potential in out.atc_potential.DATA are just Volts (V), when using “units real”.


    • Yes, I have tried using the CONP module before. It has now been a while since I did that so I don’t remember all the details any more, but I think what you need to do is modify your LAMMPS Makefile, for example src/MAKE/Makefile.mpi, and change the line that says “LIB =” to something like the following: “LIB = -llapack -lblas”. You may also need to make sure that lapack and blas are either on your LD_LIBRARY_PATH or add to that LIB = line something like -L/pathToYourLapackInstallation (replace with the actual path to your lapack installation) similarly to what I suggest in my blog post above for lib/atc/Makefile.lammps. Best of luck!


  5. Hi Alta,
    This article helps me a lot! I really appreciate your sharing. And through this method, I have successfully output the 1D potential and nextstep I am going to output a 2D potential. However, there still is a problem bothers me. How do you determine the zero potential? And is there a way to achieve zero potential at infinity?
    Thanks again for your kindly help!


    • Glad this was helpful! From a theoretical perspective, you are free to shift the electrostatic potential everywhere by a constant. Thus, you can just define the potential at infinity to be zero and then shift all of the potential values by the initial value at infinity, so that the value at infinity becomes zero.


      • Thank you very much for your help. But how do you obtain the initial value of the potential at infinity? Sometimes it is difficult to determine the potential at infinity in the system under study especially when it is a narrow system. So is there a corresponding command to automatically realize zero potential at infinity?


  6. You’re welcome! I do not know of a command to automatically make the potential at infinity be zero. I am not sure what your system looks like, but I imagine that you would need to make sure your simulation box size is large enough so that the potential value approaches a constant value as you approach infinity. To better understand any finite size effects/artifacts arising from having a finite simulation box but trying to simulate something that goes out to infinity, you could even try running multiple simulations with different box sizes/setups to see if your results converge.


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s