Running Lammps simulation

Using constant potential method to study electric double layers

Posted by Vasilii Triandafilidi on October 11, 2018

Running a constant potential simulation and analyzing the results via VMD

  1. Running the simulation:
~/bin/lmp_mpi -in example_input

Here lammps outputs the relevant information

dump 2 all dcd 500 lmp.dcd

This is the potential defined in lammps

#fix [ID] all conp [Nevery] [η] [Molecule-ID 1] [Molecule-ID 2] [Potential 1] [Potential 2] [Method] [Log] [Matrix]
#Potential 1 = -0.5 V
#Potential 2 = 0.5 V
fix e all conp 1 1.979 81 82 -0.5 0.5 inv iter
  1. Analyzing the potential via VMD

After that lets read the dcd file and analyze the results of the simulation. To do that lets create a viz.vmd file:

topo readlammpsdata example_datafile full waitfor all
animate read dcd lmp.dcd waitfor all
mol modstyle 0 0 VDW 0.600000 12.000000

Then lets open the file:

vmd -e viz.vmd

There we will see this:

Screen Shot 2018-10-11 at 3.08.20 PM

To analyze the potential distribution, we head to the Analysis/Volmap tool where we see this window. We update the selection, choose coulombmsm, select the molecule 0 and compute for all frames

Screen Shot 2018-10-11 at 3.08.28 PM

This simulation produces a volmap_out.dx file that we will read using Python module gridData (make sure to install it beforehand

  1. Plotting the results

Analyzing results of the simulation

import pandas as pd
import numpy as np
import string
from gridData import Grid
import matplotlib.pyplot as plt'fivethirtyeight')

def analyze_dx(_f):
    """analyzes dx _file

        _f (string): filename with extension of the dx file that we want to analyze

        pandas: pandas dataframe columns=['coord_phi','phiz']
    g = Grid(_f)
    nx, ny, nz = g.grid.shape
    Lx = g.edges[0][0]
    Ly = g.edges[1][0]
    Lzl = g.edges[2][0]
    Lzr = g.edges[2][-1]

    g_all = np.zeros(nz)

    count = 0

    for i in range(nx):
        for j in range(ny):
            g_all += g.grid[i][j]
            count += 1
    r_array = np.linspace(Lzl, Lzr, nz)

    yy = g_all/count

    A = np.vstack([r_array, yy])
    df = pd.DataFrame(A.T, columns=['coord_phi','phiz'])
    return df

filename = 'volmap_out'
df = analyze_dx(filename + '.dx')
plt.plot(df['coord_phi'], df['phiz'],'o--')

This script will produce this plot:



  1. Analyzing the plotted data

The data is very noisy, the terminal values are but it is important to notice that terminal values are of opposite sign and are equal. After proper equilibration the potential will converge to a proper values of -0.5-0.5

  1. Automating the calculation:
set a [atomselect top all]
volmap coulombmsm $a -res 1.0 -allframes -combine avg -o volmap_out.dx
volmap coulombmsm $a -res 1.0 -allframes -combine stdev -o volmap_out_stdv.dx

filename = 'volmap_out'
df = parse_file_info.analyze_dx(filename + '.dx')
df_std = parse_file_info.analyze_dx(filename + '_stdv.dx')

plt.errorbar(df['coord_phi'], df['phiz'], yerr=df_std['phiz'], fmt='o')

we are getting this:


Now we can see how “un-equilibrated” the simulation is