Assignment 3 - Duplicate

# environment setup - packages # NOTE: takes a couple of minutes to run !conda install -c conda-forge openmm !conda install -c conda-forge openmmtools=0.20.3 !conda install -c conda-forge mdtraj

# environment setup - scripts import sys sys.path.append('./work/') from openmm import unit

# make sure to run this codeblock :) def simulate_solidwater(parms): from functions import water2 parm_error = water2.check_parms(parms) if not parm_error: water2.prepare_system(parms) water2.minimize() water2.prepare_sim(parms) water2.run_sim(parms) results = water2.analyse() return(results) else: results = None

# change your parameters here parms_water_nve = { "parms_0" : { # original simulation parameters 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVE', }, "parms_1" : { # same parameters as parms_1, but modified number of steps 'steps':500, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVE', }, "parms_2" : { # same parameters as parms_1, but modified number of steps 'steps':1000, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVE', }, "parms_3" : { # same parameters as parms_1, but modified number of steps 'steps':1500, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVE', }, "parms_4" : { # same parameters as parms_1, but modified number of steps 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 2. * unit.femtoseconds, 'ensemble':'NVE', }, "parms_5" : { # same parameters as parms_1, but modified number of steps 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 5. * unit.femtoseconds, 'ensemble':'NVE', }, "parms_6" : { # same parameters as parms_1, but modified number of steps 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 10. * unit.femtoseconds, 'ensemble':'NVE', } }

import matplotlib.pyplot as plt from cycler import cycler plt.style.use('seaborn-v0_8-whitegrid') # plots # (feel free to separate into more codeblocks) # (just make sure to enable them when publishing, i.e. double check your output pdf) plt.plot(results_1["time"],results_1["potE"],label="Potential Energy") plt.plot(results_1["time"],results_1["kinE"],label="Kinetic Energy") plt.plot(results_1["time"],results_1["totalE"],label="Total Energy") plt.legend(loc='best') plt.xlabel("Time (s)"),plt.ylabel("Energy (kJ/mol)") plt.title("Energy of system over time interval") plt.show()

Ans > The total energy is the sum of potential and kinetic energy. Total energy is constant, which aligns with our expectation that energy is conserved in an isolated system. Kinetic energy and potential energy vary at the same rate but with opposite signs, so that KE + PE is constant.

# DEFINE VARIABLES FOR STORING ALL RESULTS results_water_nve = [] results_water_nve_var,results_water_nve_avg = [], [] for parm_name in list(parms_water_nve.keys()): # RUN THE SIMULATION FOR VARIOUS TIME LENGTHS, HOLDING THE TIMESTEP CONSTANT results_i = simulate_solidwater(parms_water_nve[parm_name]); # COMPUTE AVERAGE AND VARIANCE FOR ALL RESULTS var_i = variances(results_i); avg_i = averages(results_i); results_water_nve.append(results_i); results_water_nve_var.append(var_i); results_water_nve_avg.append(avg_i);

# plots for average and variance # (feel free to separate into more codeblocks) # (just make sure to enable them when publishing, i.e. double check your output pdf) # number of plots to generate N = len(results_water_nve) for i in range(N): # this loop iterates over every simulation result obtained plt.rc('axes', prop_cycle=(cycler('color','ccmmyy'))) time = results_water_nve[i]["time"] for field in ["potE","kinE","totalE"]: result = results_water_nve[i][field] avg = results_water_nve_avg[i][field]*np.ones_like(time) var = results_water_nve_var[i][field]*np.ones_like(time) plt.plot(time,result,alpha=0.75,linewidth=0.5,label=field) #plt.errorbar(time,y=avg,yerr=var,fmt='none') plt.fill_between(time,avg-var,avg+var,alpha=0.17,color='black',label="variance "+field) plt.plot(time,avg,alpha=1,linestyle='dotted',label="average "+field) plt.legend(bbox_to_anchor=(1.1, 1.05)) plt.xlabel("Time (s)"),plt.ylabel("Energy (kJ/mol)") plt.title("Energy of system over time interval, {} steps, dt={:.2e}".format(len(time),time[-1]/len(time))) plt.show()

Ans > Vary the length of the simulation for a constant time step (dt), what do you observe? I observe that the variance and average are highly dependent on the simulation itself. The variance is not obviously approaching zero with a large step number (at least for the step sizes that I'm looking at). This is somewhat counterintuitive, as the central limit theorem suggests that the variance would approach zero as the sum of data approaches the mean. For a fixed number of time steps, what happens when you vary the distance between them (dt)? When we vary the dt, we are able to plot over a larger timespan, at the cost of precision of each calculated value for a given t=t'. If dt is too large, the simulation might not converge. Which quantity do you expect will be conserved? We expect the the Total Energy to be conserved for an NVE system since the system is considered isolated.

# plot distance here for i in range(N): # this loop iterates over every simulation result obtained #plt.rc('axes', prop_cycle=(cycler('color','ccmmyy'))) time = results_water_nve[i]["time"] plt.plot(time, results_water_nve[i]['rOO']) #plt.legend(bbox_to_anchor=(1.1, 1.05)) plt.xlabel("Time (s)"),plt.ylabel("Distance [nm]") plt.title("Energy of system over time interval, {} steps, dt={:.2e}".format(len(time),time[-1]/len(time))) plt.show() print("Average: ",np.mean(results_water_nve[i]["rOO"])) print("Variance: ",np.var(results_water_nve[i]["rOO"]))

Ans > yes we observe oscillations.

parms_water_nvt = { "parms_0" : { # original simulation parameters 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVT', }, "parms_1" : { # same parameters as parms_1, but modified number of steps 'steps':500, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVT', }, "parms_2" : { # same parameters as parms_1, but modified number of steps 'steps':1000, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVT', }, "parms_3" : { # same parameters as parms_1, but modified number of steps 'steps':1500, 'skip_steps':1, 'temperature':275., 'dt': 1. * unit.femtoseconds, 'ensemble':'NVT', }, "parms_4" : { # same parameters as parms_1, but modified number of steps 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 2. * unit.femtoseconds, 'ensemble':'NVT', }, "parms_5" : { # same parameters as parms_1, but modified number of steps 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 5. * unit.femtoseconds, 'ensemble':'NVT', }, "parms_6" : { # same parameters as parms_1, but modified number of steps 'steps':100, 'skip_steps':1, 'temperature':275., 'dt': 10. * unit.femtoseconds, 'ensemble':'NVT', }, } # DEFINE VARIABLES FOR STORING ALL RESULTS results_water_nvt = [] results_water_nvt_var,results_water_nvt_avg = [], [] for parm_name in list(parms_water_nvt.keys()): # RUN THE SIMULATION FOR VARIOUS TIME LENGTHS, HOLDING THE TIMESTEP CONSTANT results_i = simulate_solidwater(parms_water_nvt[parm_name]); # COMPUTE AVERAGE AND VARIANCE FOR ALL RESULTS var_i = variances(results_i); avg_i = averages(results_i); results_water_nvt.append(results_i); results_water_nvt_var.append(var_i); results_water_nvt_avg.append(avg_i);

import matplotlib.pyplot as plt # plots # (feel free to separate into more codeblocks) # (just make sure to enable them when publishing, i.e. double check your output pdf) # number of plots to generate N = len(results_water_nvt) for i in range(N): # this loop iterates over every simulation result obtained plt.rc('axes', prop_cycle=(cycler('color','cmy'))) time = results_water_nvt[i]["time"] for field in ["potE","kinE","totalE"]: result = results_water_nvt[i][field] avg = results_water_nvt_avg[i][field]*np.ones_like(time) var = results_water_nvt_var[i][field]*np.ones_like(time) plt.plot(time,result,alpha=0.75,linewidth=0.5,label=field) #plt.errorbar(time,y=avg,yerr=var,fmt='none') #plt.fill_between(time,avg-var,avg+var,alpha=0.17,color='black',label="variance "+field) #plt.plot(time,avg,alpha=1,linestyle='dotted',label="average "+field) plt.legend(bbox_to_anchor=(1.5, 1.05)) plt.xlabel("Time (s)"),plt.ylabel("Energy (kJ/mol)") plt.title("Energy of system over time interval, {} steps, dt ={:.2e}".format(len(time),time[-1]/len(time))) plt.show()

Ans > The total energy is no longer conserved since an NVE system is not isolated. However, we still have the total energy as the sum of kinetic and potnetial.

# DEFINE METHODS FOR COMPUTING AVERAGES AND VARIANCES OF KE, PE, AND TOTE FOR A SIMULATION RESULTS OUTPUT import numpy as np def variances(results): variances = { "potE" : np.var(results["potE"]), "kinE" : np.var(results["kinE"]), "totalE" : np.var(results["totalE"]) } return variances def averages(results): averages = { "potE" : np.average(results["potE"]), "kinE" : np.average(results["kinE"]), "totalE" : np.average(results["totalE"]) } return averages

for i in range(N): # this loop iterates over every simulation result obtained #plt.rc('axes', prop_cycle=(cycler('color','ccmmyy'))) time = results_water_nvt[i]["time"] plt.plot(time, results_water_nve[i]['rOO']) #plt.legend(bbox_to_anchor=(1.1, 1.05)) plt.xlabel("Time (s)"),plt.ylabel("Distance [nm]") plt.title("Energy of system over time interval, {} steps, dt={:.2e}".format(len(time),time[-1]/len(time))) plt.show() print("Average: ",np.mean(results_water_nvt[i]["rOO"])) print("Variance: ",np.var(results_water_nvt[i]["rOO"]))

Ans > yes we observe oscillations.

def simulate_argon(parms): from functions import arBox parm_error = arBox.check_parms(parms) if not parm_error: arBox.prepare_system(parms) arBox.minimize() arBox.equilibrate(parms) arBox.prepare_sim(parms) arBox.run_sim(parms) results = arBox.v_analysis() return(results) else: results = None

parms_1 = { 'steps':100, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':200, 'density':0.9 } # run sim here

from functions import arBox results_X = arBox.v_analysis()

Ans >

arBox.gen_pair_dist() result_file = '/work/Ar_histo'

Ans >

def simulate_liquidwater(parms): from functions import waterBox parm_error = waterBox.check_parms(parms) if not parm_error: waterBox.prepare_system(parms) waterBox.minimize() waterBox.equilibrate(parms) waterBox.prepare_sim(parms) waterBox.run_sim(parms)

parms_X = { 'steps':100, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':200, 'Box_edge':2.*unit.nanometers } # run sim here simulate_liquidwater(parms_X)

from functions import waterBox waterBox.gen_pair_dist() results_OO = '/work/OO_histo' results_OH = '/work/OH_histo'

Ans >