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 lib–atc 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
-L/home/alta/software/BLAS-3.8.0
Alternatively, in the src/ directory, execute:
make lib–linalg 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 patom_style full
bond_style harmonic
angle_style harmonic
dihedral_style oplspair_style lj/cut/coul/long 12.0 12.0
pair_modify mix geometric
kspace_style pppm 1.0e-4read_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
region BOX block EDGE EDGE EDGE EDGE EDGE EDGE
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 textrerun 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.
Alta Fang 
this article helped me – thank you!
LikeLike
You’re welcome, glad it was helpful!
LikeLike
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.
LikeLike
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!
LikeLike
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?
LikeLike
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.
LikeLike
[…] people have reached out to me to ask me questions after seeing my previous blog post, “Outputting the electrostatic potential field in LAMMPS“, leading to some good […]
LikeLike
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)??
LikeLike
Glad to hear this has been helpful! According to the fix atc documentation, the boundary conditions available for this method are only either periodic or fixed: https://lammps.sandia.gov/doc/atc_mesh_create.html, which should match the overall LAMMPS boundary conditions set using the boundary command (https://lammps.sandia.gov/doc/boundary.html). If you want to specify fixed electrostatic potential at certain positions, that is actually quite challenging. To do that, you could consider the constant potential method: https://github.com/zhenxingwang/lammps-conp. However, that can be considerably more computationally expensive, and I recommend thinking carefully about the physics of your system — for example, specifying a constant potential would only make sense on the surface of an electrically conductive material.
LikeLike
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.
LikeLike
Sure, glad it is helpful! In the past I have used VisIt to visualize 3D fields: https://visit-dav.github.io/visit-website/ , which is open-source. My colleagues have used Tecplot 360: https://www.tecplot.com/products/tecplot-360/ although that is commercial software.
LikeLike
Hello Alta,
Thank you for a nice post!
I have two questions:
1. I cannot see the field electric_potential in the LAMMPS manual (https://docs.lammps.org/atc_hardy_fields.html), what was the version you used?
2. What is the unit of the electrostatic potential in the out.atc_potential.DATA?
LikeLike
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.
LikeLike
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.
LikeLike
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”.
LikeLike
Hi, I am quite excited to see this. I do wonder if you have had success implementing the CONP module (https://github.com/zhenxingwang/lammps-conp) in LAMMPS. I hit a roadblock at the step 3 of 4 during the process i.e. “adding library files of LAPACK and BLAS to the link flag in LAMMPS Makefile”
LikeLike
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!
LikeLike
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!
LikeLike
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.
LikeLike
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?
LikeLike
Thank you so much for your help. It inspired me a lot. I hope everything goes well with you!
LikeLike
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.
LikeLike