The code.ill.fr has been recreated and upgraded with the latest version this weekend, If you encounter any problem please inform the Helpdesk.

Commit a4fc542a authored by eric pellegrini's avatar eric pellegrini Committed by Remi Perenon
Browse files

Fully rewrite the rbt analysis

parent 76a530a8
......@@ -61,7 +61,6 @@ class RigidBodyTrajectory(IJob):
settings['frames']=('frames', {"dependencies":{'trajectory':'trajectory'}})
settings['atom_selection']=('atom_selection',{"dependencies":{'trajectory':'trajectory'}})
settings['grouping_level']=('grouping_level',{"default" : "atom","dependencies":{'trajectory':'trajectory','atom_selection':'atom_selection'}})
settings['stepwise']=('boolean',{'default':True})
settings['reference']=('integer',{"mini":0})
settings['remove_translation']=('boolean',{'default':False})
settings['output_files']=('output_files', {"formats":["netcdf"]})
......@@ -69,15 +68,12 @@ class RigidBodyTrajectory(IJob):
def initialize(self):
"""
"""
self.numberOfSteps = self.configuration['atom_selection']['selection_length']
self.numberOfSteps = self.configuration['frames']['number']
if (self.configuration['reference']['value'] >= self.configuration['trajectory']['length']):
raise JobError(self,'Invalid reference frame. Must be an integer in [%d,%d[' % (0,self.configuration['trajectory']['length']))
self.nFrames = self.configuration['frames']['number']
self._rbt = {'trajectory':{}}
self._quaternions = numpy.zeros((self.configuration['atom_selection']['selection_length'],self.configuration['frames']['number'], 4),dtype=numpy.float64)
self._coms = numpy.zeros((self.configuration['atom_selection']['selection_length'],self.configuration['frames']['number'], 3),dtype=numpy.float64)
self._fits = numpy.zeros((self.configuration['atom_selection']['selection_length'],self.configuration['frames']['number']),dtype=numpy.float64)
......@@ -89,180 +85,111 @@ class RigidBodyTrajectory(IJob):
for i in range(self.configuration['atom_selection']['selection_length']):
indexes = self.configuration['atom_selection']["indexes"][i]
self._groups.append(Collection([atoms[idx] for idx in indexes]))
self.referenceFrame = self.configuration['reference']['value']
trajectory = self.configuration['trajectory']['instance']
self.referenceConfiguration = trajectory.universe.contiguousObjectConfiguration(conf = trajectory.configuration[self.referenceFrame])
selectedAtoms = Collection()
for indexes in self.configuration['atom_selection']['indexes']:
for idx in indexes:
selectedAtoms.addObject(atoms[idx])
# Create trajectory
self.outputFile = Trajectory(selectedAtoms, self.configuration['output_files']['files'][0], 'w')
# Create the snapshot generator
self.snapshot = SnapshotGenerator(self.configuration['trajectory']['universe'], actions = [TrajectoryOutput(self.outputFile, ["configuration","time"], 0, None, 1)])
def run_step(self, index):
'''
'''
group = self._groups[index]
rbtPerGroup = {}
# Those matrix will store the quaternions and the CMS coming from the RBT trajectory.
rbtPerGroup['quaternions'] = numpy.zeros((self.configuration['frames']['number'], 4), dtype=numpy.float64)
rbtPerGroup['com'] = numpy.zeros((self.configuration['frames']['number'], 3), dtype=numpy.float64)
rbtPerGroup['fit'] = numpy.zeros((self.configuration['frames']['number'],), dtype=numpy.float64)
rbtPerGroup['trajectory'] = {}
# Case of a moving reference.
if self.configuration['stepwise']['value']:
# The reference configuration is always the one of the previous frame excepted for the first frame
# where it is set by definition to the first frame (could we think about a cyclic alternative way ?).
for comp in range(self.configuration['frames']['number']):
if comp == 0:
previousFrame = self.configuration['frames']['value'][0]
else:
previousFrame = self.configuration['frames']['value'][comp-1]
currentFrame = self.configuration['frames']['value'][comp]
self.configuration['trajectory']['instance'].universe.setFromTrajectory(self.configuration['trajectory']['instance'], currentFrame)
refConfig = self.configuration['trajectory']['instance'].universe.contiguousObjectConfiguration()
# refConfig = self.configuration['trajectory']['instance'].configuration[previousFrame]
# The RBT is created just for the current step.
rbt = self.configuration['trajectory']['instance'].readRigidBodyTrajectory(group,
first=currentFrame,
last=currentFrame+1,
skip=1,
reference=refConfig)
# The corresponding quaternions and cms are stored in their corresponding matrix.
rbtPerGroup['quaternions'][comp,:] = rbt.quaternions
rbtPerGroup['com'][comp,:] = rbt.cms
rbtPerGroup['fit'][comp] = rbt.fit
# The simplest case, the reference frame is fixed.
# A unique RBT is performed from first to last skipping skip steps and using refConfig as the reference.
else:
trajectory = self.configuration['trajectory']['instance']
currentFrame = self.configuration['frames']['value'][index]
for group_id, group in enumerate(self._groups):
self.configuration['trajectory']['instance'].universe.setFromTrajectory(self.configuration['trajectory']['instance'], self.configuration['reference']['value'])
refConfig = self.configuration['trajectory']['instance'].universe.contiguousObjectConfiguration()
rbt = trajectory.readRigidBodyTrajectory(group,
first = currentFrame,
last = currentFrame+1,
skip = 1,
reference = self.referenceConfiguration)
centerOfMass = group.centerOfMass(self.referenceConfiguration)
# The rotation matrix corresponding to the selected frame in the RBT.
transfo = Quaternion(rbt.quaternions[0]).asRotation()
if self.configuration['remove_translation']['value']:
# The transformation matrix corresponding to the selected frame in the RBT.
transfo = Translation(centerOfMass)*transfo
# Compose with the CMS translation if the removeTranslation flag is set off.
else:
# The transformation matrix corresponding to the selected frame in the RBT.
transfo = Translation(Vector(self._coms[group_index,frame,:]))*transfo
# Loop over the atoms of the group to set the RBT trajectory.
for atom in group:
# The coordinates of the atoms are centered around the center of mass of the group.
xyz = self.referenceConfiguration[atom] - centerOfMass
atom.setPosition(transfo(Vector(xyz)))
self._quaternions[group_id,index,:] = rbt.quaternions[0]
self._coms[group_id,index,:] = rbt.cms[0]
self._fits[group_id,index] = rbt.fit[0]
self.snapshot(data = {'time' : self.configuration['frames']['time'][currentFrame]})
return index, None
# If a fixed reference has been set. We can already set the reference configuration here.
# refConfig = self.configuration['trajectory']['instance'].configuration[self.configuration['reference']['value']]
# The RBT is created for the whole frame selection
rbt = self.configuration['trajectory']['instance'].readRigidBodyTrajectory(group,
first=self.configuration['frames']['first'],
last=self.configuration['frames']['last']+1,
skip=self.configuration['frames']['step'],
reference=refConfig)
# The corresponding quaternions and cms are stored in their corresponding matrix.
rbtPerGroup['quaternions'] = rbt.quaternions
rbtPerGroup['com'] = rbt.cms
rbtPerGroup['fit'] = rbt.fit
return index, rbtPerGroup
def combine(self, index, x):
"""
"""
self._quaternions[index] = x['quaternions']
self._coms[index] = x['com']
self._fits[index] = x['fit']
pass
def finalize(self):
'''
'''
atoms = sorted(self.configuration['trajectory']['universe'].atomList(), key=operator.attrgetter('index'))
selectedAtoms = Collection()
for indexes in self.configuration['atom_selection']['indexes']:
for idx in indexes:
selectedAtoms.addObject(atoms[idx])
# Create trajectory
outputFile = Trajectory(selectedAtoms, self.configuration['output_files']['files'][0], 'w')
# Create the snapshot generator
snapshot = SnapshotGenerator(self.configuration['trajectory']['universe'], actions = [TrajectoryOutput(outputFile, ["configuration","time"], 0, None, 1)])
# The output is written
for frame in range(self.configuration['frames']['number']):
currentFrame = self.configuration['frames']['value'][frame]
self.configuration['trajectory']['instance'].universe.setFromTrajectory(self.configuration['trajectory']['instance'], currentFrame)
# Case of a moving reference.
if self.configuration['stepwise']['value']:
if frame == 0:
previousFrame = self.configuration['frames']['value'][0]
else:
previousFrame = self.configuration['frames']['value'][frame-1]
refConfig = self.configuration['trajectory']['instance'].configuration[previousFrame]
else:
# If a fixed reference has been set. We can already set the reference configuration here.
refConfig = self.configuration['trajectory']['instance'].configuration[self.configuration['reference']['value']]
for group_index,group in enumerate(self._groups):
centerOfMass = group.centerOfMass(refConfig)
# Loop over the atoms of the group to set the RBT trajectory.
for atom in group:
# The coordinates of the atoms are centered around the center of mass of the group.
xyz = refConfig[atom] - centerOfMass
# The rotation matrix corresponding to the selected frame in the RBT.
transfo = Quaternion(self._quaternions[group_index,frame,:]).asRotation()
if self.configuration['remove_translation']['value']:
# The transformation matrix corresponding to the selected frame in the RBT.
transfo = Translation(centerOfMass)*transfo
# Compose with the CMS translation if the removeTranslation flag is set off.
else:
# The transformation matrix corresponding to the selected frame in the RBT.
transfo = Translation(Vector(self._coms[group_index,frame,:]))*transfo
atom.setPosition(transfo(Vector(xyz)))
snapshot(data = {'time' : self.configuration['frames']['time'][frame]})
outputFile.close()
outputFile = NetCDFFile(self.configuration['output_files']['files'][0], 'a')
outputFile.createDimension('NFRAMES', self.configuration['frames']['number'])
outputFile.createDimension('NGROUPS', self.configuration['atom_selection']['selection_length'])
outputFile.createDimension('NFRAMES', self.configuration['frames']['number'])
outputFile.createDimension('QUATERNIONLENGTH',4)
outputFile.createDimension('XYZ',3)
# The NetCDF variable that stores the quaternions.
QUATERNIONS = outputFile.createVariable('quaternion', numpy.dtype(numpy.float64).char, ('NGROUPS', 'NFRAMES','QUATERNIONLENGTH'))
QUATERNIONS = outputFile.createVariable('quaternions', numpy.dtype(numpy.float64).char, ('NGROUPS','NFRAMES','QUATERNIONLENGTH'))
# The NetCDF variable that stores the centers of mass.
COM = outputFile.createVariable('com', numpy.dtype(numpy.float64).char, ('NGROUPS','NFRAMES','XYZ'))
COM = outputFile.createVariable('coms', numpy.dtype(numpy.float64).char, ('NGROUPS','NFRAMES','XYZ'))
# The NetCDF variable that stores the rigid-body fit.
FIT = outputFile.createVariable('fit', numpy.dtype(numpy.float64).char, ('NGROUPS','NFRAMES'))
FIT = outputFile.createVariable('fits', numpy.dtype(numpy.float64).char, ('NGROUPS','NFRAMES'))
outputFile.info = str(self)
# Loop over the groups.
for comp in range(self.configuration['atom_selection']['selection_length']):
aIndexes = self.configuration['atom_selection']['indexes'][comp]
outputFile.info += 'Group %s: %s\n' % (comp, [index for index in aIndexes])
QUATERNIONS[comp,:,:] = self._quaternions[comp][:,:]
COM[comp,:,:] = self._coms[comp][:,:]
FIT[comp,:] = self._fits[comp][:]
QUATERNIONS.assignValue(self._quaternions[comp,:,:])
COM.assignValue(self._coms[comp,:,:])
FIT.assignValue(self._fits[comp,:])
outputFile.close()
REGISTRY['rbt'] = RigidBodyTrajectory
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment