<jupyter_start><jupyter_text>This is an introductory tutorial for the Harmonic Oscilator with Openmm  The objectives are:1. Set up a system consisting of harmonic oscilator2. Change the initial position of the particle and calculate the potential energy and force3. Minimize the potential energy4. Perform molecular dynamics simulations at different timesteps to monitor integration errors Homework:1. Create two alternative harmonic oscilators with doubled mass (or doubled force constants) 2. Repeat point 2 and 3 with the two harmonic oscillators: How does the potential energy and force change3. Choose the appropiate time step for the new harmonic oscilator4. Integrate equation of motion with verlocity verlet<jupyter_code>from simtk import openmm, unit
from simtk.openmm import app
import numpy as np
from openmmtools import testsystems
import matplotlib.pyplot as plt
import scipy as scipy
from scipy.stats import linregress
%matplotlib inline
Ho_1 = testsystems.HarmonicOscillator()
help(Ho_1)
position = Ho_1.positions
print(position)<jupyter_output>[[0. 0. 0.]] A<jupyter_text>Now we change the position (point 2)<jupyter_code>new_position = [[1 , 0, 0]]
Ho_1.positions = new_position
print(Ho_1.positions)<jupyter_output>[[1, 0, 0]]<jupyter_text>It is better to use the unit in the position. Standard unit in the harmonic oscilator are kcal/mol and kcal/mol/AA^2<jupyter_code>new_position = [[1 * unit.angstrom , 0, 0]]
Ho_1.positions = new_position
print(Ho_1.positions)<jupyter_output>[[Quantity(value=1, unit=angstrom), 0, 0]]<jupyter_text>Now we define an integrator class with the verlet intergrator and a timestep o 1 fs.<jupyter_code>timestep = 1.0 * unit.femtosecond
integrator = openmm.VerletIntegrator(timestep)<jupyter_output><empty_output><jupyter_text>With the integrator we define a context class which stores the current state of the simulation.<jupyter_code>context = openmm.Context(Ho_1.system,integrator)<jupyter_output><empty_output><jupyter_text>Before doing anything with context class we have to define the initial position<jupyter_code>context.setPositions(Ho_1.positions)<jupyter_output><empty_output><jupyter_text>With this initial position we want to know some properties of the harmonic oscilator in the current state as:1. Potential energy2. ForceThis is done with the getState function of context class and the output is a copy of the current state which contains  all variables which were set to True. This variables can be obtained with the getXXX function of the new class<jupyter_code>context.getState?
state = context.getState(getEnergy=True, getForces=True)
potential_energy = state.getPotentialEnergy()
print(potential_energy)<jupyter_output>209.20001220703125 kJ/mol<jupyter_text>Question: Is this value correct?<jupyter_code>print(state.getForces())<jupyter_output>[Vec3(x=-4184.0, y=0.0, z=0.0)] kJ/(nm mol)<jupyter_text>Let´s reduce the position in x to only 0.1 angstrom<jupyter_code>position = [[0.1 * unit.angstrom, 0., 0.]]
context.setPositions(position)
state = context.getState(getEnergy=True, getForces=True)
print(state.getPotentialEnergy(), state.getForces())<jupyter_output>2.0920000076293945 kJ/mol [Vec3(x=-418.3999938964844, y=0.0, z=0.0)] kJ/(nm mol)<jupyter_text>Question: Why is the energy 100 times smaller when we divided the position by 10?  Now let´s minimize the potential energy of the harmonic oscilator system with the LocalEnergyMinimizer class (point 3)<jupyter_code>openmm.LocalEnergyMinimizer.minimize(context)
state = context.getState(getPositions=True,getEnergy=True)
state.getPositions(), state.getPotentialEnergy()<jupyter_output><empty_output><jupyter_text>We are at the bottom of the harmonic oscilator again.  Let´s go back to the original positions and do some molecular dynamics at constant energy (microcanonical ensemble) (point 4) The total simulation time will be 3 times the period of the oscilator which is given as:<jupyter_code>T = 2 * np.pi * np.sqrt(Ho_1.mass / Ho_1.K)
print(T)<jupyter_output>0.19414738754194238 ps<jupyter_text>First start a new simulation with a new context class using a timestep the factor times smaller than the period<jupyter_code>del context, integrator
factor = 10
timestep = T / factor
integrator = openmm.VerletIntegrator(timestep)
position = [[0.1 * unit.nanometer, 0., 0.]]
context = openmm.Context(Ho_1.system,integrator)
context.setPositions(position)
total_iterations = factor 
for iteration in range(total_iterations + 1):
    state = context.getState(getEnergy=True,getPositions=True)   
    print('%8.3f ps : position %12.6f nm' % (state.getTime() / unit.picoseconds, state.getPositions()[0][0] / unit.nanometer))
    integrator.step(1)<jupyter_output>0.000 ps : position     0.100000 nm
   0.019 ps : position     0.060522 nm
   0.039 ps : position    -0.002850 nm
   0.058 ps : position    -0.065096 nm
   0.078 ps : position    -0.101644 nm
   0.097 ps : position    -0.098064 nm
   0.116 ps : position    -0.055770 nm
   0.136 ps : position     0.008541 nm
   0.155 ps : position     0.069480 nm
   0.175 ps : position     0.102989 nm
   0.194 ps : position     0.095840 nm<jupyter_text>Questions: 1.) Why is the final position not exactly the one at the beginning after one period?2.) What can you do to make them equal each other? Let´s look at the potential enerygy after one periods of oscilation<jupyter_code>del context, integrator
factor = 10
timestep = T / factor
integrator = openmm.VerletIntegrator(timestep)
position = [[0.1 * unit.nanometer, 0., 0.]]
context = openmm.Context(Ho_1.system,integrator)
context.setPositions(position)
total_iterations = factor  + 1
for iteration in range(total_iterations):
    state = context.getState(getEnergy=True, getPositions=True)   
    print('%8.3f ps : energy %12.6f kJ/mol' % (state.getTime() / unit.picoseconds, state.getPotentialEnergy() / unit.kilojoule_per_mole))
    integrator.step(1)<jupyter_output>0.000 ps : energy   209.200012 kJ/mol
   0.019 ps : energy    76.627052 kJ/mol
   0.039 ps : energy     0.169900 kJ/mol
   0.058 ps : energy    88.648636 kJ/mol
   0.078 ps : energy   216.133057 kJ/mol
   0.097 ps : energy   201.176819 kJ/mol
   0.116 ps : energy    65.066917 kJ/mol
   0.136 ps : energy     1.526111 kJ/mol
   0.155 ps : energy   100.990929 kJ/mol
   0.175 ps : energy   221.894882 kJ/mol
   0.194 ps : energy   192.157486 kJ/mol<jupyter_text>Let´s look at the force after one period of oscilation<jupyter_code>del context, integrator
factor = 10
timestep = T / factor
integrator = openmm.VerletIntegrator(timestep)
position = [[0.1 * unit.nanometer, 0., 0.]]
context = openmm.Context(Ho_1.system,integrator)
context.setPositions(position)
total_iterations = factor
for iteration in range(total_iterations + 1):
    state = context.getState(getEnergy=True,getForces=True)   
    print('%8.3f ps : position %12.6f kJ/mol/nm' % (state.getTime() / unit.picoseconds, state.getForces()[0][0] / unit.kilojoule_per_mole *unit.nanometer))
    integrator.step(1)<jupyter_output>0.194 ps : position -4009.954834 kJ/mol/nm
   0.214 ps : position -2127.763916 kJ/mol/nm
   0.233 ps : position   594.434692 kJ/mol/nm
   0.252 ps : position  3081.959717 kJ/mol/nm
   0.272 ps : position  4352.775879 kJ/mol/nm
   0.291 ps : position  3905.184570 kJ/mol/nm
   0.311 ps : position  1915.888062 kJ/mol/nm
   0.330 ps : position  -829.770447 kJ/mol/nm
   0.349 ps : position -3247.848633 kJ/mol/nm
   0.369 ps : position -4383.727539 kJ/mol/nm
   0.388 ps : position -3788.979980 kJ/mol/nm<jupyter_text>Store the position, forces and energy in arrays to analyze them<jupyter_code>del context, integrator
factor = 10
timestep = T / factor
integrator = openmm.VerletIntegrator(timestep)
position = [[0.1 * unit.nanometer, 0., 0.]]
context = openmm.Context(Ho_1.system,integrator)
context.setPositions(position)
time = []
positions = []
energy = []
forces = []
total_iterations = factor
for iteration in range(total_iterations):
    state = context.getState(getEnergy=True,getPositions=True, getForces=True)
    time.append(state.getTime() / unit.picoseconds)
    position = (state.getPositions(asNumpy=True) / unit.nanometer)
    positions.append(position[0][0])
    force = (state.getForces(asNumpy=True) / unit.kilojoule_per_mole * unit.nanometer)
    forces.append(force[0][0])
    energy.append(state.getPotentialEnergy() / unit.kilojoule_per_mole)
    integrator.step(1)<jupyter_output><empty_output><jupyter_text>Let´s plot the position as function of time<jupyter_code>plt.plot(time,positions,'-o')<jupyter_output><empty_output><jupyter_text>Now let´s plot the potential energy<jupyter_code>plt.plot(time,energy,'-o')<jupyter_output><empty_output><jupyter_text>Plot the force as fucntion of time<jupyter_code>plt.plot(time,forces,'o-')<jupyter_output><empty_output><jupyter_text>After one period the forces should be zero if the integration is done correctly.<jupyter_code>sum(forces), forces[0]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09890232974910394, 'fraction_numerical': 0.06843637992831542, 'mean_word_length': 4.976572958500669, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '565774', 'license': 'MIT', 'n_tokens_mistral': 3483, 'n_tokens_neox': 3049, 'n_words': 840, 'path': 'harmonic_oscilator.ipynb', 'repo_name': 'QCMM/Advanced_md', 'size': '8928'}