[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