# 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_1 = {
'steps':100,
'skip_steps':1,
'temperature':75.,
'dt': 1 * unit.femtoseconds,
'ensemble':'NVE',
}
parms_2 = {
'steps':500,
'skip_steps':1,
'temperature':95.,
'dt': 1 * unit.femtoseconds,
'ensemble':'NVE',
}
results_1 = simulate_solidwater(parms_1)
results_2 = simulate_solidwater(parms_2)
import matplotlib.pyplot as plt
import numpy as np
print(results_1)
#kinE, nsteps, positions, potE, rOO, time, totalE = *zip(*sorted(results_1.items()))
plt.figure()
plt.subplot(311)
plt.plot(results_1['time'],results_1['kinE'])
plt.title('Kinetic energy vs time')
plt.xlabel('Time')
plt.ylabel('Kinetic energy kJ/mole')
plt.subplot(312)
plt.plot(results_1['time'],results_1['potE'])
plt.title('Potential energy vs time')
plt.xlabel('Time')
plt.ylabel('Potential energy kJ/mole')
plt.subplot(313)
plt.plot(results_1['time'],results_1['totalE'])
plt.title('Total energy vs time')
plt.xlabel('Time')
plt.ylabel('Total energy kJ/mole')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
# plots
# (feel free to separate into more codeblocks)
# (just make sure to enable them when publishing, i.e. double check your output pdf)
Ans >
import statistics
average_kin1 = statistics.mean(results_1['kinE'])
average_kin2 = statistics.mean(results_2['kinE'])
variance_kin1 = statistics.pvariance(results_1['kinE'],average_kin1)
variance_kin2 = statistics.pvariance(results_2['kinE'],average_kin2)
average_pot1 = statistics.mean(results_1['potE'])
average_pot2 = statistics.mean(results_2['potE'])
variance_pot1 = statistics.variance(results_1['potE'],average_pot1)
variance_pot2 = statistics.variance(results_2['potE'],average_pot2)
average_total1 = statistics.mean(results_1['totalE'])
average_total2 = statistics.mean(results_2['totalE'])
variance_total1 = statistics.variance(results_1['totalE'],average_total1)
variance_total2 = statistics.variance(results_2['totalE'],average_total2)
avg = {'Kin1':average_kin1, 'Kin2':average_kin2, 'Pot1':average_pot1,'Pot2':average_pot2, 'Total1':average_total1, 'Total2':average_total2}
energy_avg = list(avg.keys())
values_avg = list(avg.values())
var = {'Kin1':variance_kin1, 'Kin2':variance_kin2, 'Pot1':variance_pot1,'Pot2':variance_pot2, 'Total1':variance_total1, 'Total2':variance_total2}
energy_var = list(var.keys())
values_var = list(var.values())
# 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)
plt.figure()
plt.subplot(311)
plt.plot(results_2['time'],results_2['kinE'])
plt.title('Kinetic energy vs time')
plt.xlabel('Time')
plt.ylabel('Kinetic energy kJ/mole')
plt.subplot(312)
plt.plot(results_2['time'],results_2['potE'])
plt.title('Potential energy vs time')
plt.xlabel('Time')
plt.ylabel('Potential energy kJ/mole')
plt.subplot(313)
plt.plot(results_2['time'],results_2['totalE'])
plt.title('Total energy vs time')
plt.xlabel('Time')
plt.ylabel('Total energy kJ/mole')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
Ans >
#plot distance here
plt.figure()
plt.subplot(211)
plt.plot(results_1['time'],results_1['rOO'])
plt.title('Distance between two oxygens vs time 1')
plt.xlabel('Time')
plt.ylabel('Distance')
plt.subplot(212)
plt.plot(results_2['time'],results_2['rOO'])
plt.title('Distance between two oxygens vs time 2')
plt.xlabel('Time')
plt.ylabel('Distance')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
Ans >
# change your parameters here
parms_1 = {
# place values here
'steps':100,
'skip_steps':1,
'temperature':75.,
'dt': 1. * unit.femtoseconds,
'ensemble':'NVT'
}
parms_2 = {
'steps':100,
'skip_steps':1,
'temperature':95.,
'dt': 1. * unit.femtoseconds,
'ensemble':'NVT',
}
# run the sim here (check usage instructions for help):
results_1 = simulate_solidwater(parms_1)
results_2 = simulate_solidwater(parms_2)
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(311)
plt.plot(results_1['time'],results_1['kinE'])
plt.title('Kinetic energy vs time')
plt.xlabel('Time')
plt.ylabel('Kinetic energy kJ/mole')
plt.subplot(312)
plt.plot(results_1['time'],results_1['potE'])
plt.title('Potential energy vs time')
plt.xlabel('Time')
plt.ylabel('Potential energy kJ/mole')
plt.subplot(313)
plt.plot(results_1['time'],results_1['totalE'])
plt.title('Total energy vs time')
plt.xlabel('Time')
plt.ylabel('Total energy kJ/mole')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
# plots
# (feel free to separate into more codeblocks)
# (just make sure to enable them when publishing, i.e. double check your output pdf)
Ans >
# space for multiple simulation runs
parms_1 = {
# place values here
'steps':300,
'skip_steps':1,
'temperature':75.,
'dt': 1. * unit.femtoseconds,
'ensemble':'NVT'
}
parms_2 = {
'steps': 300,
'skip_steps':1,
'temperature':95.,
'dt': 1. * unit.femtoseconds,
'ensemble':'NVT',
}
results_1 = simulate_solidwater(parms_1)
results_2 = simulate_solidwater(parms_2)
import matplotlib.pyplot as plt
plt.figure()
plt.subplot(311)
plt.plot(results_1['time'],results_1['kinE'])
plt.title('Kinetic energy vs time 1')
plt.xlabel('Time')
plt.ylabel('Kinetic energy kJ/mole')
plt.subplot(312)
plt.plot(results_1['time'],results_1['potE'])
plt.title('Potential energy vs time 1')
plt.xlabel('Time')
plt.ylabel('Potential energy kJ/mole')
plt.subplot(313)
plt.plot(results_1['time'],results_1['totalE'])
plt.title('Total energy vs time 1')
plt.xlabel('Time')
plt.ylabel('Total energy kJ/mole')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
plt.figure()
plt.subplot(311)
plt.plot(results_2['time'],results_2['kinE'])
plt.title('Kinetic energy vs time 2')
plt.xlabel('Time')
plt.ylabel('Kinetic energy kJ/mole')
plt.subplot(312)
plt.plot(results_2['time'],results_2['potE'])
plt.title('Potential energy vs time 2')
plt.xlabel('Time')
plt.ylabel('Potential energy kJ/mole')
plt.subplot(313)
plt.plot(results_2['time'],results_2['totalE'])
plt.title('Total energy vs time 2')
plt.xlabel('Time')
plt.ylabel('Total energy kJ/mole')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
# 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)
import statistics
average_kin1 = statistics.mean(results_1['kinE'])
average_kin2 = statistics.mean(results_2['kinE'])
variance_kin1 = statistics.pvariance(results_1['kinE'],average_kin1)
variance_kin2 = statistics.pvariance(results_2['kinE'],average_kin2)
average_pot1 = statistics.mean(results_1['potE'])
average_pot2 = statistics.mean(results_2['potE'])
variance_pot1 = statistics.variance(results_1['potE'],average_pot1)
variance_pot2 = statistics.variance(results_2['potE'],average_pot2)
average_total1 = statistics.mean(results_1['totalE'])
average_total2 = statistics.mean(results_2['totalE'])
variance_total1 = statistics.variance(results_1['totalE'],average_total1)
variance_total2 = statistics.variance(results_2['totalE'],average_total2)
avg = {'Kin1':average_kin1, 'Kin2':average_kin2, 'Pot1':average_pot1,'Pot2':average_pot2, 'Total1':average_total1, 'Total2':average_total2}
energy_avg = list(avg.keys())
values_avg = list(avg.values())
var = {'Kin1':variance_kin1, 'Kin2':variance_kin2, 'Pot1':variance_pot1,'Pot2':variance_pot2, 'Total1':variance_total1, 'Total2':variance_total2}
energy_var = list(var.keys())
values_var = list(var.values())
#plot distance here
plt.figure()
plt.subplot(211)
plt.plot(results_1['time'],results_1['rOO'])
plt.title('Distance between two oxygens vs time 1')
plt.xlabel('Time')
plt.ylabel('Distance')
plt.subplot(212)
plt.plot(results_2['time'],results_2['rOO'])
plt.title('Distance between two oxygens vs time 2')
plt.xlabel('Time')
plt.ylabel('Distance')
plt.subplots_adjust(wspace=2,
hspace=1)
plt.show()
Ans >
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
results_1 = simulate_argon(parms_1)
from functions import arBox
import numpy as np
results_X = arBox.v_analysis()
standard_error_bars = np.sqrt(results_X['variance']/results_X['nsteps'])
standard_error_bars
import matplotlib.pyplot as plt
plt.figure()
plt.plot(results_X['time'],results_X['potE'])
plt.title('Average potential energy vs time')
plt.xlabel('time')
plt.ylabel('Average potential energy')
plt.show()
Ans >
arBox.gen_pair_dist()
result_file = '/work/Ar_histo'
r, gr = np.loadtxt('Ar_histo', unpack=True)
plt.figure()
plt.plot(r,gr)
plt.title('Radial distribution function vs radius')
plt.xlabel('Radius')
plt.ylabel('Radial distribution function g(r)')
plt.show
Ans > The radial distribution doesn't change when going from N=200 to N=500
n = dr*4.*np.pi*(N/volume)*NN
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':500,
'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'
import matplotlib.pyplot as plt
gOOr = np.loadtxt('OO_histo')
gOHr = np.loadtxt('OH_histo')
plt.figure()
plt.subplot(211)
plt.plot(gOOr)
plt.subplot(212)
plt.plot(gOHr)
plt.show()
Ans >