Python Assignment #2 for MD

# First, if you haven't SWIGGed your code yet, you'll need to use a
# newer version of mine. Copy my new version of MDSolution into your directory,
# replacing the old one. (I hadn't put in the shifted energy when
# I added the cutoff!) Swig the new code.

# Second, make sure that Exceed is running

import sys
sys.path.append("K:\\Simulations\\xmgr")
import xmgr

graph = xmgr.WindowsXMGRInterface()

graph.StartXMGR()

import math
import Numeric # Efficient numerical arrays package: NumPy

SinX = Numeric.array([[0.0, math.sin(0.0)],
			  [1.5708, math.sin(1.5708),
			  [3.1416, math.sin(3.1416)]])

SinX	# Just an array

graph.Plot(SinX)

X = Numeric.arrayrange(0, 3.1416, 0.31416)
Y = Numeric.sin(X)	# Fast array version

# [X,Y] is transpose of array we need

SinX = Numeric.transpose([X,Y])

sys.path.append("K:\\Simulations\\MDPython")	# Your path here

import MD

atoms = MD.Atoms()
rand = MD.Rand()
nonOverlapRandom = MD.RandomPositionStrategy(1.0,rand)
nonOverlapRandom.Move(atoms)
thermalize = MD.ThermalizeVelocityStrategy(1,rand)
thermalize.Move(atoms)
dt = 0.01
verlet = MD.UpdateStrategy(dt)

# Plot trajectory of atom #3
x3 = Numeric.zeros((100,2), Numeric.Float)	# Allocate array

for i in xrange(Numeric.shape(x3)[0]):
	for j in xrange(10):
		verlet.Move(atoms)
	x3[i][0] = atoms.X(3,0)
	x3[i][1] = atoms.X(3,1)
	
x3

graph.Plot(x3)


# Notice, the periodic boundary conditions make for ugly trajectories

# Plot energy versus time

times = Numeric.arrayrange(0,10,dt)
e = Numeric.zeros(len(times), Numeric.Float)	# Just to allocate space

for i in xrange(len(times)):
	e[i] = atoms.TotalEnergy()
	verlet.Move(atoms)
	
e_vs_t = Numeric.transpose([times,e])

graph.SendCommand("flush")
graph.Plot(e_vs_t)

verlet.SetDt(0.001)
e2 = Numeric.zeros(len(times), Numeric.Float)

for i in xrange(len(times)):
	for j in xrange(10):
		verlet.Move(atoms)
	e2[i] = atoms.TotalEnergy()
	
e2_vs_t = Numeric.transpose([times, e2])
graph.Plot(e_vs_t, e2_vs_t)

# Zoom in with the magnifying glass to see the fluctuations

# Measuring the Specific Heat
# Copy your EnergyObserver.py file into a new file, KineticEnergyObserver.py
# Edit it to calculate the mean and fluctuations in the kinetic energy
# Also, define functions which 
#   calculate the temperature from the kinetic energy
#		Kinetic Energy = nAtoms * DIMENSION * kT/2
#    (assume Boltzmann's constant k = 1)
#   calculate the total specific heat C from the fluctuations in the kinetic energy
#		<(Delta T)^2> = k T^2 [ 2 / (DIMENSION * k * numAtoms) - 1/C]

# Observer the kinetic energy while storing it
import KineticEnergyObserver
keObs = KineticEnergyObserver.KineticEnergyObserver()

# Test your new routine: if it isn't right, edit, reload(KineticEnergyObserver), and repeat
reload(KineticEnergyObserver)
keObs = KineticEnergyObserver.KineticEnergyObserver()
atoms.KineticEnergy()
keObs.Update(atoms)
keObs.numKE
keObs.KEtot
keObs.KEbar()
keObs.KineticTemperature()
verlet.Move(atoms)
keObs.Update(atoms)
keObs.SpecificHeatPerAtom()

# Average over time interval of 100, using every tenth time step

keObs.Reset()
dt = 0.01
verlet.SetDt(dt)
times = Numeric.arrayrange(0,100,10*dt)
ke = Numeric.zeros(len(times), Numeric.Float)

for i in xrange(len(times)):
	for j in xrange(10):
		verlet.Move(atoms)
	ke[i] = atoms.KineticEnergy()
	keObs.Update(atoms)

keObs.KineticTemperature()
keObs.SpecificHeatPerAtom()

ke_vs_t = Numeric.transpose([times, ke])
graph.SendCommand("flush")
graph.Plot(ke_vs_t)

# Let's look at the temperature-dependent specific heat

# Make a few more atoms

atoms = MD.Atoms(50)
nonOverlapRandom.Move(atoms)


# Define a couple of functions for averaging

def Relax(atoms, nRelax, nSkip):
	for i in xrange(nRelax):
		for j in xrange(nSkip):	
			verlet.Move(atoms)
		thermalize.Move(atoms)	# Allows thermalization with heat bath
		
def KEAve(atoms, nAverage, nSkip):
	keObs.Reset()
	for i in xrange(nAverage):
		for j in xrange(nSkip):
			verlet.Move(atoms)
		keObs.Update(atoms)
		
# Specific heat at high temperature

thermalize.SetTemperature(3.0)
Relax(atoms,100,10)
KEAve(atoms,100,10)

keObs.KineticTemperature()
keObs.SpecificHeatPerAtom()	# Should be kT/2 per important DOF: how many DOF per atom?


# Quench slowly to zero temperature	
coolingTemperatures = Numeric.arrayrange(1, 0, -0.1)

conjugateGradient = MD.CGStrategy()

# Before conjugate gradient, quench slowly so atoms
# don't get caught in bad metastable states

def Quench(atoms):
	for i in xrange(len(coolingTemperatures)):
		thermalize.SetTemperature(coolingTemperatures[i])
		thermalize.Move(atoms)
		Relax(atoms,10,10)
		print atoms.TotalEnergy()
	conjugateGradient.Move(atoms)
	thermalize.SetTemperature(0.0)
	thermalize.Move(atoms)
	print atoms.TotalEnergy()

# The best energy I've gotten was -119.4
# You may want to shake it around some if you're energy is much higher

# What's the specific heat per atom at low temperatures?

thermalize.SetTemperature(0.05)
Relax(atoms,100,10)

KEAve(atoms,100,10)
keObs.KineticTemperature()		# Not quite the same?
keObs.SpecificHeatPerAtom()	# Should be 1/2 kT per DOF: how many DOF per atom?


# Now, let's measure the specific heat on heating
heatingTemperatures = Numeric.arrayrange(0.1,1,0.1)
T = Numeric.zeros(len(heatingTemperatures), Numeric.Float)
C = Numeric.zeros(len(T), Numeric.Float)

def Heat(atoms):
	for i in xrange(len(heatingTemperatures)):
		thermalize.SetTemperature(heatingTemperatures[i])
		thermalize.Move(atoms)
		Relax(atoms,10,10)
		KEAve(atoms, 100, 10)
		T[i] = keObs.KineticTemperature()
		C[i] = keObs.SpecificHeatPerAtom()
		print i, T[i], C[i]
	
C_vs_T = Numeric.transpose([T,C])
graph.Plot(C_vs_T)

# Quench and heat a few times. You'll find large fluctuations
# from one run to the next. What do you think might be happening?


# Try higher temperatures

thermalize.SetTemperature(10.0)

Statistical Mechanics: Entropy, Order Parameters, and Complexity, now available at Oxford University Press (USA, Europe).