In [ ]:
%matplotlib notebook
In [ ]:
import matplotlib.pyplot as plt
In [ ]:
from lammps import IPyLammps
In [ ]:
L = IPyLammps()
In [ ]:
import math
L.units("real")
L.atom_style("molecular")
L.boundary("f f f")
L.neighbor(0.3, "bin")
L.dihedral_style("harmonic")
In [ ]:
L.read_data("data.dihedral")
In [ ]:
L.pair_style("zero", 5)
L.pair_coeff("*", "*")
In [ ]:
L.mass(1, 1.0)
In [ ]:
L.velocity("all", "set", 0.0, 0.0, 0.0)
In [ ]:
L.run(0);
In [ ]:
L.image(zoom=1.0)
In [ ]:
L.atoms[3].position
In [ ]:
L.atoms[3].position = (1.0, 0.0, 1.0)
In [ ]:
L.image(zoom=1.0)
In [ ]:
L.eval("pe")
In [ ]:
L.atoms[3].position = (1.0, 0.0, -1.0)
In [ ]:
L.run(0);
In [ ]:
phi = [d * math.pi / 180 for d in range(360)]
In [ ]:
pos = [(1.0, math.cos(p), math.sin(p)) for p in phi]
In [ ]:
K = 80.0
d = 1
n = 2
E_analytical = [K * (1 + d * math.cos(n*p)) for p in phi]
In [ ]:
pe = []
for p in pos:
L.atoms[3].position = p
L.run(0);
pe.append(L.eval("pe"))
In [ ]:
plt.plot(range(360), pe, range(360), E_analytical)
plt.xlabel('angle')
plt.ylabel('E')
In [ ]: