Dear developpers,
I am trying to make an extension that allows normal mode calculations with MMTK into chimera. I have several problems beginning with the following.
I have first used the MMTKinter.py and added a new function called NM_ffm that calculates normal modes with a full atom force field. This works perfectly (although I am not yet able to make an animation in chimera that follows the vectors calculated...that will be a future question for you). Then I added another function (NM_mcm) that allows me to do NM calculations with the mechanical constraint (see http://dirac.cnrs-orleans.fr/MMTK/using-mmtk/mmtk-example-scripts/normal-modes/normal-modes-with-mechanical-constraints). For this I need to define a subspace that requires a universe.protein object (the functions are given at the end of the mail). Unfortunately, I can not pass the molecular object from the chimera space as a universe.protein object. We thought in putting self.universe.protein call into the _makeUniverse function in MMTKinter.py but that did not work. We tried other combinations but we now face the fact that we are blocked. May be we need another function in order to define subspaces but yet we still don't know how to pass chimera infos to MMTK.
Could you give me a hand so that I can define correctly subspaces on chimera molecules?
All the best and thank in advance for any help,
Victor
*************
Computational Biotechnological Chemistry @transmet
####Function for normal modes calculations in full atom ff
def NM_ffm(self, nsteps, stepsize=0.02,
interval=None, action=None, **kw):
from chimera import replyobj
timestamp("_minimize")
from MMTK import Units
from MMTK.ForceFields.Amber import AmberData
from MMTK.Minimization import SteepestDescentMinimizer
from MMTK.Trajectory import LogOutput
# from JD to make NM calculations
from MMTK.NormalModes import NormalModes
from chimera import runCommand as run
from chimera import selection as sele
from chimera import Point, Vector
import os, sys
# JD'end
import sys
if not interval:
actions = []
else:
actions = [ LogOutput(sys.stdout, ["energy"],
interval, None, interval) ]
kw["step_size"] = stepsize * Units.Ang
minimizer = SteepestDescentMinimizer(self.universe,
actions=actions, **kw)
if action is None or not interval:
interval = None
msg = "Initial energy: %f" % self.universe.energy()
replyobj.status(msg)
replyobj.info(msg)
saveNormalizeName = AmberData._normalizeName
AmberData._normalizeName = simpleNormalizeName
remaining = nsteps
while remaining > 0:
timestamp(" minimize interval")
if interval is None:
realSteps = remaining
else:
realSteps = min(remaining, interval)
minimizer(steps=realSteps)
remaining -= realSteps
if action is not None:
action(self)
timestamp(" finished %d steps" % realSteps)
msg = "Finished %d of %d minimization steps" % (
nsteps - remaining, nsteps)
replyobj.status(msg)
replyobj.info(msg)
replyobj.info("\n")
#Calculate normal modes JD's trucho
def despl(vec):
run("sel #0")
atoms=sele.currentAtoms()
for a in atoms:
a.setCoord(a.xformCoord()-10*(vec))
modes = NormalModes(self.universe)
for mode in modes:
print mode
print mode.array
L1=[]
L2=[]
L3=[]
for mode in modes:
y=6
u=len(mode.array)
while y<u:
L1.append(mode.array[y][0])
L2.append(mode.array[y][1])
L3.append(mode.array[y][2])
y+=1
y=0
run("sel:all")
atoms=sele.currentAtoms()
run("~sel")
u=len(L1)
a=0
o=len(atoms)
while y<u:
vec=Vector(float(L1[y]),float(L2[y]),float(L3[y]))
atoms[a].setCoord(atoms[a].coord()-(vec))
a+=1
if a==o:
a=0
raw_input()
else:
pass
y+=1
#JD's end
return modes
AmberData._normalizeName = saveNormalizeName
timestamp("end _minimize")
####Function for mechanical constraint
def NM_mcm(self, nsteps, stepsize=0.02,
interval=None, action=None, **kw):
from MMTK import *
from MMTK.Proteins import Protein
from MMTK.ForceFields import Amber94ForceField
from MMTK.NormalModes import NormalModes, SubspaceNormalModes
from MMTK.Subspace import RigidMotionSubspace
from MMTK.Minimization import ConjugateGradientMinimizer
from MMTK.Trajectory import StandardLogOutput
import numpy
# Victor Modifications
from chimera import replyobj
timestamp("_minimize")
from MMTK import Units
from MMTK.ForceFields.Amber import AmberData
from MMTK.Minimization import SteepestDescentMinimizer
from MMTK.Trajectory import LogOutput
from MMTK.NormalModes import NormalModes
from chimera import runCommand as run
from chimera import selection as sele
from chimera import Point, Vector
import os, sys
#/end
if not interval:
actions = []
else:
actions = [ LogOutput(sys.stdout, ["energy"],
interval, None, interval) ]
kw["step_size"] = stepsize * Units.Ang
minimizer = SteepestDescentMinimizer(self.universe,
actions=actions, **kw)
if action is None or not interval:
interval = None
msg = "Initial energy: %f" % self.universe.energy()
replyobj.status(msg)
replyobj.info(msg)
saveNormalizeName = AmberData._normalizeName
AmberData._normalizeName = simpleNormalizeName
remaining = nsteps
while remaining > 0:
timestamp(" minimize interval")
if interval is None:
realSteps = remaining
else:
realSteps = min(remaining, interval)
minimizer(steps=realSteps)
remaining -= realSteps
if action is not None:
action(self)
timestamp(" finished %d steps" % realSteps)
msg = "Finished %d of %d minimization steps" % (
nsteps - remaining, nsteps)
replyobj.status(msg)
replyobj.info(msg)
replyobj.info("\n")
# Construct system
#self.universe = InfiniteUniverse(Amber94ForceField())
#self.universe.protein = self.universe
# Set up the subspace: rigid-body translation and rotation for each residue
self.subspace = RigidMotionSubspace(self.universe, self.universe.protein.residues())