Multiple 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 discussions.

To these people, I have actually recommended using a different approach to calculate the electrostatic potential, namely using a Python script to integrate the 1D charge density profile twice, rather than using the “fix atc” feature of LAMMPS. This is the approach I actually used to produce the first plot in my previous post and also the approach used in the publication “Simulation Study of the Capacitance and Charging Mechanisms of Ionic Liquid Mixtures near Carbon Electrodes”.

Here, as a follow up to my earlier post, I will share the excerpt of LAMMPS commands and the Python script I sent to those who contacted me.

An example set of LAMMPS commands to output charge density are shown below, where some of the values are variables that need to be substituted in (such as q_bin_width), and the assumption is that the charge density profile is calculated along the x direction (so adjust the commands below if needed):

compute CHARGE all property/atom q
compute CHUNK_XBIN_ALL all chunk/atom bin/1d x lower ${q_bin_width}
fix CHARGE_DENSITY_FIX all ave/chunk ${Nevery} 1 ${Nfreq} CHUNK_XBIN_ALL c_CHARGE file out.charge_density

For further information, please see the LAMMPS documentation, particularly the pages on “fix ave/chunk” and “compute chunk/atom”. Note that you need to divide by surface area to get charge density per volume for the next step of integrating twice to get the electrostatic potential.

An alternative method for calculating the 1D charge density profile is to use the LAMMPS “dump” command to output the position and charge of each particle in the simulation, at various times, and then write a Python script to parse the dump file and calculate the charge density profile by histogramming/binning the data. I actually tended to use this latter method in my publications. For code to do this, see this post.

Below is the Python script to calculate the electric potential from the charge density profile.

"""Alta Fang
alta.y.fang@gmail.com
The below code snippet has the 0BSD license
(https://opensource.org/licenses/0BSD)"""

import scipy.integrate

def integrate_poisson_1D(r, charge_dens_profile, periodic=True):
    """Integrate Poisson's equation twice to get electrostatic potential from charge density.

    Inputs:
    r : 1D numpy array. Values of coordinates, in Angstroms.
    charge_dens_profile : 1D numpy array. Values of charge density, 
        in e*Angstrom^-3.
    periodic : Boolean. Optional keyword. If True, adds linear function to 
        ensure values at boundaries are periodic.
    Outputs:
    phi : Numpy array of same length as r. Electrostatic potential in V. 
        The value of phi at the first r value is set to be 0.
    """

    eps_0_factor = 8.854e-12/1.602e-19*1e-10 # Note: be careful about units!
    int1 = scipy.integrate.cumtrapz(charge_dens_profile, r, initial=0)
    int2 = -scipy.integrate.cumtrapz(int1, r, initial=0)/eps_0_factor
    if periodic:
        # Ensure periodic boundary conditions by adding a linear function such
        # that the values at the boundaries are equal
        phi = int2 - ((int2[-1]-int2[0])/(r[-1]-r[0])*(r-r[0]) + int2[0])
    else:
        phi = int2
    return phi

Note that you should double check the units to make sure they correspond to the units your data has! Also, you should modify the “periodic” flag as needed.

The above code corresponds to the following equations:

As can be seen from the equations above, there are two constants that are determined by boundary conditions. Thus, when the system is periodic along the x direction, the code adds a linear function so that the values of the potential are equal at both ends of the box.

Another thing to note is that the the input charge density is the charge per volume (not the 1D charge per length) because the integration (dx” and dx’) contributes the Angstrom² dimensions.