[MMTK] Trying to write a Monte Carlo code -- problems with coordinate changes

[Georg Ganzenmueller]

Hi there,

I'm new to this mailing list because I just started tinkering a little 
bit with MMTK.
I'm trying to write a small Monte Carlo NVT code and I have problems 
with the MMTK methods of translating / rotating a molecule.
The code (please see below) runs initially O.K. but after a few 
iterations in isweep, the energy rises to very high numbers. Inspection 
of the final configuration of the molecules with rasmol shows molecules 
torn apart (that is, hydrogens at one place, oygens at some other 
place). This does not happen, if I use ONLY rotation or ONLY translation 
attempts.

Could someone please point me out to the mistake I have in my code.

Thanks a lot!
Georg

%--------------------------------------------------------------------------------%
from MMTK import *
from Scientific.Geometry.Objects3D import SCLattice
from MMTK.Visualization import view
from MMTK.Random import *
from MMTK.ForceFields import Amber94ForceField
from MMTK.Units import *

from random import Random, random
from Numeric import exp
import timing
import time

# Define parameters of the lattice
edge_length = 0.5*Units.nm
lattice_size = (3, 3, 3)

# Construct the universe
universe = CubicPeriodicUniverse(edge_length*lattice_size[0],
                                         Amber94ForceField())


# Add the water molecules
for point in SCLattice(edge_length, lattice_size):
    universe.addObject(Molecule('water', position = point))

print 'Initial energy = ', universe.energy()

nmolecules = len(universe)
print 'number of chemical objects in universe = ', nmolecules

allmolecules = universe.objectList()

# initialize random number generator with seed
rv = Random(0)


# start Markov chain
timing.start()
for isweep in range(10):
    for i in range(nmolecules):
        molecule = allmolecules[i] # pick a molecule

        # rotation
        angle = 5.0*rv.random()
        zaxis  = randomDirection()
        vold = universe.energy()    # old pot. energy
        molecule.rotateAroundCenter(zaxis,angle) # rotate
        vnew = universe.energy(small_change=0) # new pot. energy
        deltaE = (vnew - vold) / (298.15 * k_B) # difference in energy, 
scaled by beta
        u_random = rv.random() # get a random number on [0,1)
        accept = 0
        if deltaE < 75.0:
            if deltaE < 0.0:
                accept = 1
            elif exp(-deltaE) > u_random:
                accept = 1
        #
        if accept == 0:
            molecule.rotateAroundCenter(zaxis,-angle) # restore old state

        #
        #
        #
        # translation
        molecule = allmolecules[i]
        trn = randomDirection()
        vold = universe.energy()
        molecule.translateBy(trn)
        vnew = universe.energy()
        deltaE = (vnew - vold) / (298.15 * k_B)
        u_random = rv.random()
        accept = 0
        if deltaE < 75.0:
            if deltaE < 0.0:
                accept = 1
            elif exp(-deltaE) > u_random:
                accept = 1
        #
        if accept == 0:
            molecule.translateBy(-trn)
   
    print 'sweep # ', isweep, ' energy = ', universe.energy()

#
timing.finish()
print "seconds:", timing.seconds()
print "milliseconds:", timing.milli()
#
view(universe)
#
%-----------------------------------------------------------------------%