The harmonic oscillator


The Harmonic Oscillator

This notebook sets up some basic classical and thermodyanmic simulations with the harmonic oscillator potential.

The Harmonic Oscillator is one of the most important systems in both quantum and classical physics. Here, we'll use it to explore Buckyball's dynamics simulatinos.


In [ ]:
%matplotlib inline
import numpy as np
from matplotlib.pyplot import *

import moldesign as mdt
from moldesign import units as u

First, let's set up the potential. We're using a simple 1-d harmonic oscillator that will be applied to the x coordinate of the first atom in the system. The potential has the form $$V_{HO}(x) = \frac{1}{2} k x^2,$$ where $k$ is the "SPRING_CONSTANT" variable.

As with more realistic chemical systems, we'll need to specify the units of all physical quantities.


In [ ]:
SPRING_CONSTANT = 1.0*u.kcalpermol/(u.angstrom**2)
model = mdt.models.HarmonicOscillator(k=SPRING_CONSTANT)

Next, we'll set up a simple 1-atom system to play with. The atom we use isn't important - we're just going to subject it to the external harmonic potential.


In [ ]:
atom = mdt.Atom(atnum=8, mass=1.0*u.amu)
mol = mdt.Molecule([atom])
mol.set_energy_model(model)

Let's take a look at this model's potential energy and forces. Because this is a 1-dimensional model, we can just scan along the x-coordinate of the first atom.


In [ ]:
scan = mdt.trajectory.Trajectory(mol)
for x in u.angstrom*np.arange(-2.0,2.0,0.05):
    atom.x = x
    mol.calculate()
    scan.new_frame()

In [ ]:
xs = scan.positions[:,0]
subplot(211).set_title('Energy')
plot(xs,scan.potential_energy,label='potential')
ylabel('energy / {}'.format(scan.potential_energy.units))
tick_params(axis='x',labelbottom='off'); grid()
subplot(212).set_title('Force')
plot(xs,scan.forces[:,0],label=u'force')
grid()
xlabel(u'position / Å'); ylabel('force / {}'.format(scan.forces.units))

Modeling

Next, let's actually do something with the model. We'll move the atom far away from the equilibrium position at $x=0$, and then run a basic minimization to find the equilibrium again.


In [ ]:
atom.x = 4.0 * u.angstrom
minimization = mol.minimize(frame_interval=1)

That likely went very quickly - most minimization techniques will perform extremely well with harmonic potentials. Let's take a look at what happened during the minimization:


In [ ]:
plot(minimization.minimization_step, minimization.positions[:,0])
grid(); xlabel('minimization step'); ylabel('atom position')

viewer = minimization.draw3d()
viewer.vdw()
viewer

Dynamics

We can also model the atom's motion under this potential. We'll set up a timestep based "integrator" that solves the Newton's equations of motion to generate a trajectory.


In [ ]:
integrator = mdt.integrators.VelocityVerlet(timestep=0.5*u.fs, frame_interval=30)
mol.set_integrator(integrator)

Next, let's specify the system's initial state - what are the starting positions and velocities?


In [ ]:
mol.positions[0] = 1.0*u.angstrom
mol.time = 0.0*u.fs
mol.integrator.time = 0.0*u.fs
mol.momenta = 0.0 * u.default.momentum

In [ ]:
trajectory = mol.run(1.0*u.ps)

In [ ]:
lone = [atom for atom in mol.atoms if atom.num_bonds == 0]
lone

Running the dynamics should take a few seconds. When it's done, let's start by animating the results:


In [ ]:
trajectory.draw3d()

Next, let's do some deeper analysis - generate some plots to visualize the trajectory. By now, it should be clear why we call this an "oscillator".


In [ ]:
l1 = plot(trajectory.time,trajectory.positions[:,0], label='position')
grid(axis='x', color='black'); grid(color='blue',axis='y')
tick_params(axis='x',labelbottom='off');
ylabel(u'position / Å', color='blue'); yticks(color='blue')
ax2=gca().twinx()
l2 = ax2.plot(trajectory.time,trajectory.momenta[:,0], color='red', label='momentum')
ax2.patch.set_alpha(0.0)
grid(color='red',axis='y'); xlabel('time / fs')
ylabel(u'momentum / {}'.format(u.default.momentum), color='red'); yticks(color='red')
tick_params(axis='x',labelbottom='off')

figure()
pe = trajectory.potential_energy
ke = trajectory.kinetic_energy
plot(trajectory.time, pe,label='potential energy')
plot(trajectory.time, ke,label='kinetic energy')
legend()
grid();xlabel('time / {}'.format(trajectory.time.units));ylabel('energy / {}'.format(pe.units))

In [ ]: