My previous posts Calculating 1D electrostatic potential profile from MD simulations and Outputting the electrostatic potential field in LAMMPS are continuing to generate interest and questions, and in response to one recent question, I went back to look at my Python scripts I used for my paper “Simulation Study of the Capacitance and Charging Mechanisms of Ionic Liquid Mixtures near Carbon Electrodes”.

Although in my previous blog post, I provided LAMMPS commands to calculate the 1D charge density profile, in fact for at least some of my paper, I used a different method to calculate the charge density directly from the LAMMPS dump file. One reason I did this was to more easily experiment with things like different bin widths and discarding some initial time snapshots. Below is Python code to do this:

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

import numpy as np

def get_charge_density(filename, n_atom_line, z_col, q_col, bin_width, area, z_min, z_max): 
    """Directly calculate charge density profile from a LAMMPS dump file

    Inputs:
        filename : String. Filename of the LAMMPS dump file
        n_atom_line : Integer. Number of values on a line that contains the full atom information
        z_col : Integer. The index of the column with the value of z
        q_col : Integer. The index of the column with the value of q
        bin_width : Float. Width of bins along z
        area : Float. Surface area (width in x multiplied by width in y)
        z_min : Float. Minimum z value to consider
        z_max : Float. Maximum z value to consider
    
    Outputs:
        bin_centers : Numpy array of bin centers along z
        charge_dens_profile : Numpy array of charge density (charge per volume)
    """  
    nbins = int((z_max-z_min)/bin_width)
    slcharges = []
    with open(filename, 'r') as f:
        for line in f:
            if line == "ITEM: TIMESTEP\n":
                slcharges.append(np.zeros(nbins))
            splitline = line.split()
            if len(splitline) == n_atom_line:
                z = float(splitline[z_col])
                q = float(splitline[q_col])
                i = int(np.floor((z-z_min)/(z_max-z_min)*nbins))
                slcharges[-1][i] += q
    bins = np.linspace(z_min, z_max, nbins+1)
    bin_centers = 0.5*(bins[1:] + bins[:-1]) # Return bin centers
    slcharges = np.array(slcharges) 
    avg_charge_dens = np.average(slcharges, axis=0) # Here you can change it to slcharges[n_skip:, :] to skip some initial frames
    charge_dens_profile = avg_charge_dens/(area*(bin_centers[1]-bin_centers[0]))
    return bin_centers, charge_dens_profile

This can then be used as input to the Python script in this post to calculate the electrostatic potential profile. Below is an example of what the LAMMPS dump file I used with the above code looked like. For example, I used n_atom_line = 10, z_col = 5, and q_col = 9 for the LAMMPS input command:

dump TRAJ all custom ${ndump} dump.lammpstrj id mol type x y z vx vy vz q