Python Assignment #1 for MD

# Try out some Python commands. Launch PythonWin, and try typing

print "Hello World"

# You'll find that "cut and paste" doesn't work well: each line needs
# you to put in the carriage return by hand. Uggh.

# Many functions in Python come from "imported" modules

import math
math.sqrt(97)

# This makes Python quite extendable, and keeps the "name space" small
# (Otherwise, names like "vector" and "employee" would already be 
#  defined by someone else)

# The "dir" function gives you lists of allowed command names in a
# module

dir(math)

# If you don't mind having these functions in your name space,
# you can import them directly

from math import *
sqrt(97)

# Let's define a simple "atoms" class
# This is temporary!
# We'll soon compile the C++ atoms class to use instead!

# If you haven't already done so, copy the file
#	F:\Simulations\md\MDPython 
# into your MD directory (probably it's still called MDAssign)

# Use File, Open to look at TAtoms.py
# (you'll have to shift to your MDPython directory). 

# Load the module TAtoms.py into your work space:

import sys
# Substitute your path name here!
sys.path.append("Z:\sethna\Simulations\MDPython")
sys.path					# See what you got!

import TAtoms

# Try looking at some of the features of the temporary Atoms class

dir(TAtoms.Atoms)	# Compare the list with the class definition

TAtoms.Atoms.__doc__	# Documentation lines!

atoms = TAtoms.Atoms(2)	# Make an atoms class with two atoms

# Try a few of the member functions for atoms

atoms.Number()
atoms.positions		# This won't work with our C++ class
atoms.X(1,1)		# but these will
atoms.V(1,1)
atoms.KineticEnergy()

# Now, create an instance of an update strategy
dir(TAtoms.UpdateStrategy)		
# Look familiar? Check in the code: what is mg?

tUpdate = TAtoms.UpdateStrategy(0.1)

tUpdate.Move(atoms)
atoms.positions
atoms.velocities
atoms.KineticEnergy()

# Let's now write an EnergyObserver and a KineticEnergyObserver
# Use "File, Open" to open EnergyObserver.py

# Try loading it before you change it.

import EnergyObserver

dir(EnergyObserver.EnergyObserver)	# ModuleName.ClassName
e = EnergyObserver.EnergyObserver()

# e doesn't do much

e.Update(atoms)
e.Ebar()
e.numE

# Edit EnergyObserver to compute the mean and variance for the energy
# Mean =  = Etot/numE
# Variance = <(E-)^2> = -^2 = E2tot/numE - (Etot/numE)^2 

# When you think you have it right, save the edited file and reload

reload(EnergyObserver)
e = EnergyObserver.EnergyObserver()

e.Update(atoms)
e.numE
e.Ebar()

# You'll likely need to edit, save, and re-load several times.

# Finally, try watching the mean and standard deviation of the energy
# change with time (under the force due to gravity):

for i in xrange(10):
	tUpdate.Move(atoms)
	e.Update(atoms)
	

# Now, we need to set up and run SWIG to get a real
# MD simulation running under Python. To do this, you'll
# need to 
# (1) leave Python
# (2) copy MDSolution into your md (or MDAssign) directory
#	We'll learn the small tricks needed to get your version of the
#	program to SWIG likely next week
# (3) start up a DOS command-prompt window
# (4) cd to C:\Program Files\Microsoft Visual Studio\VC98\bin,
#	and run vcvars32. 
# (5) In that same DOS window, cd back to your MDPython directory, and
#		nmake /f Makefile.msc
#     This should create MD.py and MDc.dll, 
#	which will contain a Python interface to our MD classes
# (6) Restart PythonWin again.

import sys
sys.path.append("Z:\sethna\Simulations\MDPython")	# Your directory here

import MD
dir(MD)
dir(MD.Atoms)		# Doesn't show classes: SWIG mangles names
dir(MD.AtomsPtr)		# End class with Ptr to get class names

atoms = MD.Atoms()
atoms.Number()

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

import EnergyObserver
e = EnergyObserver.EnergyObserver()

e.numE		# Remember?
e.Etot

e.Update(atoms)	# Does it work for real atoms?
e.Ebar()
e.Etot
e.numE

for i in xrange(10):
	verlet.Move(atoms)
	e.Update(atoms)

e.numE
e.Ebar()
e.Variance()

# That's enough for today...