[MMTK] the universe adding object can not be viewed

dpy99 dpy99 at mails.tsinghua.edu.cn
Fri Jan 20 01:37:36 CET 2006


I want to use MMTK and its molecular dynamics to simulate the simplest assembly of
atoms. I make two kinds of atoms named as A and B as the database described, and I
make a molecule AB two. My aim is to view the trajectory varying with time, as the
atom A and atom B could be used to build a Molecule AB, two atoms are removed and
a new molecule is added. The universe.addObjective(newobjects) is included in the
script. But I can not view the universe again. the data base of atom A atom B and
molecule AB is presented as follows:
##A in database
name= 'AtomA'
symbol='A'
mass=20
color='red'
vdW_radius=0.2

LJ_radius = 3.4*Ang
LJ_energy = 120*K*k_B

b_coherent = 1.909*fm
b_incoherent = 0.*fm
b_coherent_deut = 1.909*fm
b_incoherent_deut = 0.*fm
######
##B in database
 name= 'AtomB'
symbol='B'
mass=20
color='green'
vdW_radius=0.2

LJ_radius = 3.5*Ang
LJ_energy = 120*K*k_B

b_coherent = 1.909*fm
b_incoherent = 0.*fm
b_coherent_deut = 1.909*fm
b_incoherent_deut = 0.*fm
##########
##AB in database
name = 'AB'

structure = \
   "  A\n" + \
   " | \n" + \
   "B\n"

A  = Atom('A')
B = Atom('B')
bonds = [Bond(A, B)]

pdbmap = [('AB', {'A': A, 'B': B})]
configurations = {
    'default': Cartesian({A: (-0.85, 0., -0.0),
                          B: (0.85, 0., -0.0)
                         })
    }
#######
next is the python script: 
#AB dynamics with full molecular dynamics
from MMTK import *
from MMTK.Collection import *
from MMTK.ForceFields import LennardJonesForceField
from MMTK.Environment import NoseThermostat, AndersenBarostat
from MMTK.Trajectory import Trajectory, TrajectoryOutput, LogOutput
from MMTK.Dynamics import VelocityVerletIntegrator, VelocityScaler, \
                          TranslationRemover, BarostatReset
from Scientific.Geometry.Objects3D import SCLattice
from MMTK.Visualization import view
from Scientific.Geometry import Vector
import sys 
# Define parameters of the lattice
edge_length = 0.8*Units.nm
lattice_size = (5, 5, 5)
# Construct the universe
universe = OrthorhombicPeriodicUniverse((edge_length*lattice_size[0],
                                         edge_length*lattice_size[1],
                                         edge_length*lattice_size[2]),
                                        LennardJonesForceField(15.*Units.Ang))

for point in SCLattice(edge_length, lattice_size):
    universe.addObject(Atom('A', position = point))
    universe.addObject(Atom('B', position = point+Vector(0.4,0.0,0.)))

# Visualization
view(universe)
# Set the initial items list
OneItemList=[]
AtomAList=[]
AtomBList=[]
cycles=0
temperature = 300.*Units.K
pressure = 1.*Units.atm
universe.thermostat = NoseThermostat(temperature)
universe.barostat = AndersenBarostat(pressure)
while cycles<10: 
# add molecular dynamics
# Initialize velocities.
       universe.initializeVelocitiesToTemperature(temperature)

# Create trajectory and integrator.
       trajectory = Trajectory(universe, "AB_test.nc", "w", "AB NPT test")
       integrator = VelocityVerletIntegrator(universe, delta_t=10*Units.fs)
# Periodical actions for trajectory output and text log output.
       output_actions = [TrajectoryOutput(trajectory,
                                   ('configuration', 'energy', 'thermodynamic',
                                    'time', 'auxiliary'), 0, None, 20),
                  LogOutput("AB_test.log", ('time', 'energy'), 0, None, 100)]

# Do some equilibration steps, rescaling velocities and resetting the
# barostat in regular intervals.
       integrator(steps = 200,
           actions = [TranslationRemover(0, None, 100),
                      VelocityScaler(temperature, 0.1*temperature,
                                     0, None, 100),
                      BarostatReset(100)] + output_actions)

# Do some "production" steps.
       integrator(steps = 200,
           actions = [TranslationRemover(0, None, 100)] + output_actions)


# Construct the lists
       universe.acquireReadStateLock
       for objects in universe:
           print objects.name
           if objects.numberOfAtoms()==1:
                if objects.name is 'AtomA':
                    AtomAList.append(objects)
                else : AtomBList.append(objects)
       print OneItemList
       print AtomAList
       print AtomBList
       universe.releaseReadStateLock
# Assembly procedue
       for objects1 in AtomAList: 
           for objects2 in AtomBList:
               if universe.distance(objects1.position,objects2.position)> 1.72:
                    print 'break out'
                    break
               elif universe.distance(objects1.position,objects2.position)<1.70:
                    print  'break out again'
                    break
               else:
                    v1=objects1.centerOfMass()
                    v2=objects2.centerOfMass()
                    print (v1+v2)/2
                    newobject=Molecule('AB', position=(v1+v2)/2)
                    print newobject.position()
                    print universe
                    universe.acquireWriteStateLock()
                    universe.acquireConfigurationChangeLock()
                    universe.removeObject(objects1)
                    universe.removeObject(objects2)
                    print universe
                    universe.addObject(newobject)
                    print universe
                    universe.view()
                    universe.releaseConfigurationChangeLock()
                    universe.releaseWriteStateLock()
                    AtomBList.remove(objects2)
                    AtomAList.remove(objects1)
                    OneItamList.append(newobject)
                    print OneItamList
                    print AtomAList
                    print AtomBList
                    break
       universe.acquireReadStateLock
       universe.releaseReadStateLock
       universe.setBondConstraints( )
       print cycles
       cycles=cycles+1
# Close the trajectory.
trajectory.close()

The result gives error, and I don¡¯t know how to fix it could you give some clue.

[2.5384242048822721, -0.33740634135007008, -1.1721006703959536]
[2.5384242048822721, -0.33740634135007008, -1.1721006703959536]
OrthorhombicPeriodicUniverse  containing 250 objects.
OrthorhombicPeriodicUniverse  containing 248 objects.
OrthorhombicPeriodicUniverse  containing 249 objects.
Traceback (most recent call last):
  File
"D:\MOLECU~1\PYTHON~1.3\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py
", line 310, in RunScript
    exec codeObject in __main__.__dict__
  File "D:\Molecular dynamics\temp\toMMTKmaillist.py", line 101, in ?
    universe.view()
  File "D:\MOLECU~1\PYTHON~1.3\Lib\site-packages\MMTK\Collection.py", line 365, in
view
    configuration)
  File "D:\MOLECU~1\PYTHON~1.3\Lib\site-packages\MMTK\Universe.py", line 753, in
contiguousObjectConfiguration
    offset = self.contiguousObjectOffset(objects, conf)
  File "D:\MOLECU~1\PYTHON~1.3\Lib\site-packages\MMTK\Universe.py", line 1099, in
contiguousObjectOffset
    pairs.append((pairs[-1][1], mpairs[0][0]))
IndexError: list index out of range


Great appreciation given advanced. 





More information about the mmtk mailing list