Validating a dihedral potential

In [1]:

    

%matplotlib notebook


    

In [2]:

    

import matplotlib.pyplot as plt


    

In [3]:

    

from lammps import IPyLammps


    

In [4]:

    

L = IPyLammps()


    

    




LAMMPS output is captured by PyLammps wrapper

In [5]:

    

import math

L.units("real")
L.atom_style("molecular")

L.boundary("f f f")
L.neighbor(0.3, "bin")

L.dihedral_style("harmonic")


    

In [6]:

    

L.read_data("data.dihedral")


    

    
Out[6]:





['Reading data file ...',
 '  triclinic box = (-5 -5 -5) to (5 5 5) with tilt (0 0 0)',
 '  1 by 1 by 1 MPI processor grid',
 '  reading atoms ...',
 '  4 atoms',
 '  scanning dihedrals ...',
 '  1 = max dihedrals/atom',
 '  reading dihedrals ...',
 '  1 dihedrals',
 'Finding 1-2 1-3 1-4 neighbors ...',
 ' Special bond factors lj:   0          0          0         ',
 ' Special bond factors coul: 0          0          0         ',
 '  0 = max # of 1-2 neighbors',
 '  0 = max # of 1-3 neighbors',
 '  0 = max # of 1-4 neighbors',
 '  1 = max # of special neighbors']

In [7]:

    

L.pair_style("zero", 5)
L.pair_coeff("*", "*")


    

In [8]:

    

L.mass(1, 1.0)


    

In [9]:

    

L.velocity("all", "set", 0.0, 0.0, 0.0)


    

In [10]:

    

L.run(0);


    

In [11]:

    

L.image(zoom=1.0)


    

    
Out[11]:

In [12]:

    

L.atoms[3].position


    

    
Out[12]:





(1.0000000000000009, 1.0000000000000009, 0.0)

In [13]:

    

L.atoms[3].position = (1.0, 0.0, 1.0)


    

In [14]:

    

L.image(zoom=1.0)


    

    
Out[14]:

In [15]:

    

L.eval("pe")


    

    
Out[15]:





160.0

In [16]:

    

L.atoms[3].position = (1.0, 0.0, -1.0)


    

In [17]:

    

L.run(0);


    

In [18]:

    

phi = [d * math.pi / 180 for d in range(360)]


    

In [19]:

    

pos = [(1.0, math.cos(p), math.sin(p)) for p in phi]


    

In [20]:

    

K = 80.0
d = 1
n = 2
E_analytical = [K * (1 + d * math.cos(n*p)) for p in phi]


    

In [21]:

    

pe = []
for p in pos:
    L.atoms[3].position = p
    L.run(0);
    pe.append(L.eval("pe"))


    

In [22]:

    

plt.plot(range(360), pe, range(360), E_analytical)
plt.xlabel('angle')
plt.ylabel('E')


    

    















    












    
Out[22]:





<matplotlib.text.Text at 0x7f191b1628d0>

In [ ]: