# We now create a polyisoprene molecule
polyisoprene = pubchem_atoms_search(smiles="CC=C(C)CCC=C(C)CCC=C(C)CCC=C(C)CCC=C(C)CCC=C(C)C")
view(polyisoprene, viewer='x3d')
polyisoprene.calc = ani.calculator
polyisoprene.get_potential_energy()
# This initialises the velocities by drawing from the appropriate distribution at T=300K
MaxwellBoltzmannDistribution(polyisoprene, temperature_K=500)
# We now create a molecular dynamics object that can be used to run Langevin dynamics
# The parameters after the structure are the time step, the temperature, and the friction constant (in units of picoseconds)
dynamics = Langevin(polyisoprene, temperature_K=500, timestep=0.5*units.fs, friction=0.01)
# While the dynamics is running, we want to collect some data! This is achieved by creating a function to the dynamics object
# which gets called after some number of steps, and it can report to us what is happenning to the molecule, record its trajectory, etc.
distance = []
time = []
xyzfile = open('polyisoprene2.xyz', 'w') # the file we are going to record the structures to, the visualiser application "Ovito" can read such XYZ files.
def report():
print("Time: {:.3f} fs | Distance 8-13: {d} ".format(dynamics.get_time()/units.fs, d = np.linalg.norm((polyisoprene.get_positions()[8,:]-polyisoprene.get_positions()[13,:])) ))
ase.io.write(xyzfile, polyisoprene, format="extxyz")
def savearray():
distance.append(np.linalg.norm((polyisoprene.get_positions()[8,:]-polyisoprene.get_positions()[13,:])))
time.append(dynamics.get_time()/units.fs)
dynamics.attach(report, interval=10) # the interval argument specifies how many steps of dynamics are run between each call to record the trajectory
dynamics.attach(savearray, interval= 10)
dynamics.run(500)
xyzfile.close()
del dynamics
import matplotlib.pyplot as plt
import statsmodels.tsa.stattools as sm
import statsmodels.graphics.tsaplots as smg
plt.plot(time, distance)
plt.ylabel('Distance Between 8-13 Atoms')
plt.xlabel('Time')
plt.figure()
smg.plot_acf(np.array(distance)) # Calculates and plots ACF
plt.xlabel("Tau")
plt.ylabel("ACF")
# we will need the ability to run molecular dynamics with constraints: the positions of two atoms at the ends of the molecule should be fixed
# the following module enables this
from ase.constraints import FixAtoms
# suppose we have identified that the atoms with indices 27 and 29 are the ones we are going to held fixed, the following
# command achieves this:
polyisoprene.set_constraint(FixAtoms(indices=[27,29]) )
# We now create a molecular dynamics object that can be used to run Langevin dynamics
# the parameters after the structure are the time step, the temperature, and the friction constant (in units of picoseconds)
dynamics2 = Langevin(polyisoprene, temperature_K=500, timestep=0.5*units.fs, friction=0.01)
# While the dynamics is running, we want to collect some data! This is achieved by creating a function to the dynamics object
# which gets called after some number of steps, and it can report to us what is happenning to the molecule, record its trajectory, etc.
distance = []
time = []
xyzfile = open('polyisoprenefixed.xyz', 'w') # the file we are going to record the structures to, the visualiser application "Ovito" can read such XYZ files.
def report():
print("Fixed Positions 27: {twoseven} | Fixed Positions 29: {twonine}".format(twoseven = polyisoprene.get_positions()[27,:], twonine = polyisoprene.get_positions()[29,:]))
ase.io.write(xyzfile, polyisoprene, format="extxyz")
def savearray():
distance.append(np.linalg.norm((polyisoprene.get_positions()[8,:]-polyisoprene.get_positions()[13,:])))
time.append(dynamics2.get_time()/units.fs)
dynamics2.attach(report, interval=10) # the interval argument specifies how many steps of dynamics are run between each call to record the trajectory
dynamics2.attach(savearray, interval= 10)
dynamics2.run(500)
xyzfile.close()
del dynamics2
# We now create a molecular dynamics object that can be used to run Langevin dynamics
# the parameters after the structure are the time step, the temperature, and the friction constant (in units of picoseconds)
dynamics3 = Langevin(polyisoprene, temperature_K=300, timestep=0.5*units.fs, friction=0.01)
# While the dynamics is running, we want to collect some data! This is achieved by creating a function to the dynamics object
# which gets called after some number of steps, and it can report to us what is happenning to the molecule, record its trajectory, etc.
distance = []
time = []
xyzfile = open('polyisoprenestretch.xyz', 'w') # the file we are going to record the structures to, the visualiser application "Ovito" can read such XYZ files.
def report():
print("Distance: ", np.linalg.norm((polyisoprene.get_positions()[27,:]-polyisoprene.get_positions()[29,:])))
ase.io.write(xyzfile, polyisoprene, format="extxyz")
def savearray():
distance.append(np.linalg.norm((polyisoprene.get_positions()[27,:]-polyisoprene.get_positions()[29,:])))
time.append(dynamics3.get_time()/units.fs)
def stretch():
distance = np.linalg.norm((polyisoprene.get_positions()[27,:]-polyisoprene.get_positions()[29,:]))
newdis = 1.02*distance # Stretches by 2%
polyisoprene.set_distance(27, 29, newdis, fix = 0) #
dynamics3.attach(report, interval= 10) # the interval argument specifies how many steps of dynamics are run between each call to record the trajectory
dynamics3.attach(stretch, interval = 10)
dynamics3.run(100)
xyzfile.close()
del dynamics3
# Just Library Imports:
# import basic atomistic simulation modules
import numpy as np
import ase
from ase import Atoms
from ase.build import bulk
from ase.data.pubchem import pubchem_atoms_search
from ase.visualize import view
# import an energy model
import sys
sys.path.insert(0, "ANI")
import ani
# now import the modules we need to run molecular dynamics and fix atoms
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md.verlet import VelocityVerlet
from ase import units
from ase.md.langevin import Langevin
from ase.io.trajectory import Trajectory
from ase.constraints import FixAtoms
# Plotting and ACF
import matplotlib.pyplot as plt
import statsmodels.tsa.stattools as sm
import statsmodels.graphics.tsaplots as smg
# This is for supplement created files
# This opens a handle to your file, in 'r' read mode
file_handle = open('/work/lmp170.forces', 'r')
# Read in all the lines of your file into a list of lines
lines_list = file_handle.readlines()
#print(lines_list)
fhead = [] # Initializes vector we save force measurements to
ftail = []
forcebetween = []
for i in range(1, len(lines_list)):
f1, f2, f3, f4, f5, f6, f7, f8 = (float(val) for val in lines_list[i].split())
fhead.append(f1)
ftail.append(f4)
forcebetween.append(f1-f4)
forcebetweennp = np.array(forcebetween)
# Graph
plt.figure()
plt.plot(forcebetween) # Just the mess of forces
plt.title("Force Over Iterations")
plt.ylabel("Force")
# ACF
x = np.arange(37)
plt.figure()
acf_f = sm.acf(forcebetween) # Points of ACF
smg.plot_acf(forcebetweennp) # Plot of ACF w/ line things
plt.plot(acf_f, 'ro') # Just points of ACF
tauexp = 11 # Best fit envelope variable
envelope = np.exp(-x/tauexp) # Eq for best fit envelope
plt.plot(x,envelope,'b') # Plot best fit envelope
plt.title("ACF Restoring Force")
plt.ylabel("Correlation")
plt.xlabel("Tau")
print("Tauexp: ", tauexp)
print("Tauint: ", 1/2+sum(acf_f))
# This is for supplement created files
count = 0
tauint = []
errorbars = []
avgforce = np.zeros(7)
# This opens a handle to your file, in 'r' read mode
logs = ["150", "170", "190", "210", "230", "250", "270"]
samplingrates = [45, 45, 25, 30, 25, 50, 50]
for j in logs:
file_handle = open('/work/lmp'+j+'.forces', 'r')
# Read in all the lines of your file into a list of lines
lines_list = file_handle.readlines()
print(lines_list)
print(j)
# Extract dimensions from line. Put forces in intialized f array.
fhead = [] # Initializes vector we save force measurement to
ftail = []
forcebetween = []
for i in range(1, len(lines_list)):
f1, f2, f3, f4, f5, f6, f7, f8 = (float(val) for val in lines_list[i].split())
fhead.append(f1)
ftail.append(f4)
forcebetween.append(f1-f4)
avgforce[count] = sum(forcebetween)/len(forcebetween)
acf_int = sm.acf(forcebetween)
tauint.append((1/2+sum(acf_int))*samplingrates[count])
errorbars.append(np.std(forcebetween)/(np.sqrt(len(forcebetween)/(2*tauint[count]))))
count = count + 1
plt.figure()
stretch = [150, 170, 190, 210, 230, 250, 270]
plt.plot(stretch, avgforce, 'o')
plt.plot(stretch, avgforce, 'r')
plt.errorbar(stretch, avgforce, yerr = errorbars, ecolor='b')
plt.xlabel("Stretch")
plt.ylabel("Average Restoring Force")
logs = ["150", "170", "190", "210", "230", "250", "270", "290", "310"]
controlopen = open('/work/lmp150.log', 'r')
controlline = controlopen.readlines()
controlstart = controlline[624]
controlend = controlline[2849].split()[0]
allpoteng = []
avgpoteng = []
for j in logs:
file_handle = open('/work/lmp'+j+'.log', 'r')
# Read in all the lines of your file into a list of lines
lines_list = file_handle.readlines()
for i in range(1, len(lines_list)):
if lines_list[i] == controlstart:
#print(i)
callstartpoint = i
if len(lines_list[i].split()) > 0:
if lines_list[i].split()[0] == controlend:
callendpoint = i
poteng = []
for k in range(callstartpoint+1, callendpoint):
poteng.append(float(lines_list[k].split()[2]))
allpoteng.append(poteng) # List of lists
avgpoteng.append(sum(poteng)/len(poteng))
plt.figure()
plt.plot(allpoteng[0])
plt.plot(allpoteng[1])
plt.plot(allpoteng[2])
plt.plot(allpoteng[3])
plt.plot(allpoteng[4])
plt.plot(allpoteng[5])
plt.ylabel("Potential Energy")
plt.title("Potential Energy Over Iterations: Enthropic")
plt.legend(logs[0:6])
plt.figure()
plt.plot(logs, avgpoteng, 'o')
plt.plot(logs, avgpoteng)
plt.xlabel("Stretch Length")
plt.ylabel("Average Potential Energy")
plt.title("Average Potential Energy Over Stretch Length")
plt.figure()
plt.plot(allpoteng[6])
plt.plot(allpoteng[7])
plt.plot(allpoteng[8])
plt.ylabel("Potential Energy")
plt.title("Potential Energy Over Iterations: Enthalpic")
plt.legend(logs[6:9])