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 import numpy as np np.random.seed(0)

# 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_1 = { 'steps':100, 'skip_steps':1, 'temperature':75., 'dt': 1 * unit.femtoseconds, 'ensemble':'NVE', } parms_2 = { # second set of parameters 'steps':100, 'skip_steps':1, 'temperature':75., 'dt': 2 * unit.femtoseconds, 'ensemble':'NVE', } parms_3 = { # second set of parameters 'steps':100, 'skip_steps':1, 'temperature':75., 'dt': 2 * unit.femtoseconds, 'ensemble':'NVE', } # run the sim here (check usage instructions for help): results_1 = simulate_solidwater(parms_1) results_2 = simulate_solidwater(parms_2) results_3 = simulate_solidwater(parms_3)

import matplotlib.pyplot as plt kinetic_energy = results_1["kinE"] potential_energy = results_1["potE"] total_energy = results_1["totalE"] time = results_1["time"] # temperature = results_1[""] # plots plt.figure() plt.plot(time, kinetic_energy) plt.title("Kinetic Energy over time") plt.xlabel("time") plt.ylabel("Kinetic Energy") plt.show() plt.plot(time, potential_energy) plt.title("Potential Energy") plt.xlabel("Time") plt.ylabel("Potential Energy") plt.show() plt.plot(time, total_energy) plt.title("Total Energy") plt.xlabel("time") plt.ylabel("Total Energy") plt.show() # (feel free to separate into more codeblocks) # (just make sure to enable them when publishing, i.e. double check your output pdf)

Ans > For the kinetic energy, it decreases and then increases as time progresses. For the

import numpy as np # space for multiple simulation runs kinetic_average = np.mean(kinetic_energy) potential_average = np.mean(potential_energy) total_average = np.mean(total_energy) kinetic_average_2 = np.mean(results_2['kinE'])

# 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)

Ans >

# plot distance here distance = results_1["rOO"] time = results_1["time"] # plotting plt.figure() plt.plot(time, distance) plt.title("Distance between 2 oxygen atoms") plt.xlabel("Time") plt.ylabel("Distance") plt.show() # computing average OO_average = np.mean(distance) # computing variance OO_var = np.var(distance) print(f'The average distance is {OO_average} and the variance is {OO_var}')

Ans > I do not observe any osciallations. The average is 0.276 and the variance is 1.23e-6.

# change your parameters here parms_4 = { # place values here 'steps':100, 'skip_steps':1, 'temperature':75., 'dt': 1 * unit.femtoseconds, 'ensemble':'NVT', } # run the sim here (check usage instructions for help): results_4 = simulate_solidwater(parms_4)

import matplotlib.pyplot as plt kinetic_energy_nvt = results_4['kinE'] potential_energy_nvt = results_4['potE'] total_energy_nvt = results_4['totalE'] time = results_4["time"] # plots # (feel free to separate into more codeblocks) # (just make sure to enable them when publishing, i.e. double check your output pdf) plt.figure() plt.plot(time, kinetic_energy_nvt) plt.title("Kinetic Energy over time") plt.xlabel("Time") plt.ylabel("Kinetic Energy") plt.show() plt.plot(time, potential_energy_nvt) plt.title("Potential Energy") plt.xlabel("Time") plt.ylabel("Potential Energy") plt.show() plt.plot(time, total_energy_nvt) plt.title("Total Energy") plt.xlabel("Time") plt.ylabel("Total Energy") plt.show()

Ans >

# space for multiple simulation runs kinetic_average_nvt = np.mean(kinetic_energy_nvt) potential_average_nvt = np.mean(potential_energy_nvt) total_energy_nvp = np.mean(total_energy_nvt)

# 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)

distance = results_4["rOO"] time = results_4['time'] plt.figure() plt.plot(time, distance) plt.title("Distance between two O atoms") plt.xlabel("Time") plt.ylabel("Distance") plt.show() rOO_avg = np.mean(distance) rOO_var = np.var(distance) print(f"The average distance is {rOO_avg} and {rOO_var}.")

Ans > There are no osciallations. The average distance is 0.279 and the variance is 3.515e-6.

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':150, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':200, 'density':0.9 } parms_2 = { 'steps':150, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':300, 'density':0.9 } parms_3 = { 'steps':200, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':400, 'density':0.9 } parms_4 = { 'steps':150, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':500, 'density':0.9 } # run sim here results_ar_1 = simulate_argon(parms_1)

from functions import arBox results_X_200 = arBox.v_analysis() average_200 = results_X_200["mean"] print(f'The potential energy average for 200 particles is {average_200}.')

results_ar_300 = simulate_argon(parms_2) results_X_300 = arBox.v_analysis() average_300 = results_X_300["mean"] print(f'The potential energy average for 300 particles is {average_300}.')

Ans > As the number of particles increases, the potential energy gets smaller and smaller. The standard error for 200 particles is 1.6508588000892237. The standard error for 300 particles is 2.5645175143446006. The standard error for 400 particles is 2.8188976714367002. The standard error for 500 particles is 1.6906108479074953.

results_ar_1 = simulate_argon(parms_1) r200 = arBox.gen_pair_dist() result_file = '/work/Ar_histo'

[r, g_r] = np.loadtxt("Ar_histo", unpack=True) plt.figure() plt.plot(r,g_r) plt.title("g(r) as a function of distance") plt.ylabel("g(r)") plt.xlabel("x") plt.show()

Ans >

# req = [x for x, y in g_r,r if g_r <= 0.5] use = [(x,y) for x in r for y in g_r if x<=0.5] print(use) sum = 0 dr = 0.01 N = 200 volume = 8.734223 # print(dr) for i in range(len(use)): # 1 is gr # 0 is r distance = (use[i][0]) g = use[i][1] sum += g*(distance)*distance print(dr*4.*np.pi*(N/volume)*sum)

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_200 = { 'steps':100, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':200, 'Box_edge':2.*unit.nanometers } parms_300 = { 'steps':100, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':300, 'Box_edge':2.*unit.nanometers } parms_400 = { 'steps':100, 'skip_steps':10, 'temperature':300., 'equil_steps':100, 'N':400, 'Box_edge':2.*unit.nanometers } # run sim here

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

[r, gr] = np.loadtxt("/work/OO_histo", unpack=True) plt.figure() plt.plot(r, gr) plt.title("g(r) as a function of distance") plt.ylabel("g(r)") plt.xlabel("x") plt.show()

Ans >