'''
Normal mode calculation using external apps or libraries.
(c) 2011-2012 Thomas Holder, MPI for Developmental Biology
License: BSD-2-Clause
'''
from io import StringIO
from pymol import cmd, stored, CmdException
[docs]
def normalmodes_pdbmat(selection, cutoff=10.0, force=1.0, mass='COOR',
first=7, last=10, choose='LOWE', substruct='RESI', blocksize=4,
exe='pdbmat', diag_exe='diagrtb',
prefix='mode', states=7, factor=-1, clean=1, quiet=1, async_=-1,
*, _self=cmd, **kwargs):
'''
DESCRIPTION
PDBMAT and DIAGRTB wrapper
Runs "pdbmat" and "diagrtb" and generates perturbed models for modes
"first" to "last". WARNING: May run for a long time.
PDBMAT computes the mass-weighted second derivatives energy matrix,
using Tirion's model, that is, an elastic network model (ENM).
In such models, close particles (atoms) are linked by springs.
http://ecole.modelisation.free.fr/modes.html
NOTES
Only considers ATOM records, if your model contains MSE residues or
ligands that you want to consider, prepare it like this:
mse2met (all) # for MSE residues
alter (hetatm), type='ATOM' # for any hetatm
ARGUMENTS
selection = string: atom selection
cutoff = float: interaction distance cutoff in angstroem {default: 10}
force = float: interaction force constant {default: 1.0}
mass = string: origin of mass values {default: COOR}
first = int: first mode to create perturbed model {default: 7}
last = int: last mode to create perturbed model {default: 10}
choose = string: eigenvalues chosen {default: LOWE}
substruct = string: type of substructuring {default: RESI}
blocksize = int: nb of residues per block {default: 4}
'''
args = [selection, cutoff, force, mass,
first, last, choose, substruct, blocksize,
exe, diag_exe, prefix, states, factor, clean, quiet]
quiet, async_ = int(quiet), int(kwargs.pop('async', async_))
if async_ < 0:
async_ = not quiet
if not async_:
return _normalmodes(*args)
from pymol.wizard.message import Message
wiz = Message(['normalmodes: please wait ...', ''], dismiss=0)
_self.set_wizard(wiz)
import threading
t = threading.Thread(target=_normalmodes, args=args + [wiz])
t.setDaemon(1)
t.start()
def _normalmodes(selection, cutoff, force, mass,
first, last, choose, substruct, blocksize,
exe, diag_exe, prefix, states, factor, clean, quiet, wiz=None,
*, _self=cmd):
import tempfile, subprocess, os, shutil, sys
from chempy import cpv
cutoff, force = float(cutoff), float(force)
first, last, blocksize = int(first), int(last), int(blocksize)
clean, quiet = int(clean), int(quiet)
exe = cmd.exp_path(exe)
diag_exe = cmd.exp_path(diag_exe)
tempdir = tempfile.mkdtemp()
if not quiet:
print(' normalmodes: Temporary directory is', tempdir)
try:
sele_name = _self.get_unused_name('__pdbmat')
except AttributeError:
sele_name = '__pdbmat'
try:
filename = os.path.join(tempdir, 'mobile.pdb')
commandfile = os.path.join(tempdir, 'pdbmat.dat')
_self.select(sele_name, '(%s) and not hetatm' % (selection))
_self.save(filename, sele_name)
f = open(commandfile, 'w')
f.write('''! pdbmat file
Coordinate FILENAME = %s
MATRIx FILENAME = matrix.sdijb
INTERACtion DISTance CUTOF = %.3f
INTERACtion FORCE CONStant = %.3f
Origin of MASS values = %s
Output PRINTing level = 0
MATRix FORMat = BINA
''' % (filename, cutoff, force, mass.upper()))
f.close()
if not quiet:
print(' normalmodes: running', exe, '...')
if wiz is not None:
wiz.message[1] = 'running pdbmat'
_self.refresh_wizard()
process = subprocess.Popen([exe], cwd=tempdir,
universal_newlines=True,
stderr=subprocess.STDOUT, stdout=subprocess.PIPE)
for line in process.stdout:
if not quiet:
sys.stdout.write(line)
try:
natoms = len(open(os.path.join(tempdir, 'pdbmat.xyzm')).readlines())
except Exception as e:
print(e)
natoms = _self.count_atoms(sele_name)
if natoms != _self.count_atoms(sele_name):
raise CmdException('pdbmat did not recognize all atoms')
commandfile = os.path.join(tempdir, 'diagrtb.dat')
f = open(commandfile, 'w')
f.write('''! diagrtb file
MATRIx FILENAME = matrix.sdijb
COORdinates filename = %s
Eigenvector OUTPut filename= diagrtb.eigenfacs
Nb of VECTors required = %d
EigeNVALues chosen = %s
Type of SUBStructuring = %s
Nb of residues per BLOck = %d
Origin of MASS values = %s
Temporary files cleaning = ALL
MATRix FORMat = BINA
Output PRINting level = 0
''' % (filename, last, choose.upper(), substruct.upper(), blocksize, mass.upper()))
f.close()
exe = diag_exe
if not quiet:
print(' normalmodes: running', exe, '...')
if wiz is not None:
wiz.message[1] = 'running diagrtb'
_self.refresh_wizard()
process = subprocess.Popen([exe], cwd=tempdir,
universal_newlines=True,
stderr=subprocess.STDOUT, stdout=subprocess.PIPE)
for line in process.stdout:
if not quiet:
sys.stdout.write(line)
eigenfacs, frequencies = parse_eigenfacs(
os.path.join(tempdir, 'diagrtb.eigenfacs'), last)
if wiz is not None:
wiz.message[1] = 'generating objects'
_self.refresh_wizard()
states = int(states)
factor = float(factor)
if factor < 0:
factor = natoms**0.5
for mode in range(first, last + 1):
name = prefix + '%d' % mode
_self.delete(name)
if not quiet:
print(' normalmodes: object "%s" for mode %d with freq. %.6f' %
(name, mode, frequencies[mode - 1]))
for state in range(1, states + 1):
_self.create(name, sele_name, 1, state, zoom=0)
_self.alter_state(state, name,
'(x,y,z) = cpv.add([x,y,z], cpv.scale(next(myit), myfac))',
space={'cpv': cpv, 'myit': iter(eigenfacs[mode - 1]),
'next': next,
'myfac': factor * (state - (states + 1) / 2.0)})
# if CA only selection, show ribbon trace
if natoms == _self.count_atoms('(%s) and name CA' % sele_name):
_self.set('ribbon_trace_atoms', 1, prefix + '*')
_self.show_as('ribbon', prefix + '*')
# store results
if not hasattr(stored, 'nma_results'):
stored.nma_results = []
stored.nma_results.append({
'facs': eigenfacs,
'freq': frequencies,
'sele': sele_name,
})
except OSError:
print('Cannot execute "%s", please provide full path to executable' % (exe))
except CmdException as e:
print(' normalmodes: failed!', e)
finally:
if clean:
shutil.rmtree(tempdir)
elif not quiet:
print(' normalmodes: Working directory "%s" not removed!' % (tempdir))
# cmd.delete(sele_name)
if wiz is not None:
_self.set_wizard_stack([w for w in _self.get_wizard_stack() if w != wiz])
[docs]
def parse_eigenfacs(filename='diagrtb.eigenfacs', readmax=20):
line_it = iter(open(filename))
eigenfacs = []
values = []
for line in line_it:
a = line.split()
if a[0] == 'VECTOR':
assert a[2] == 'VALUE'
number = int(a[1])
assert number == len(eigenfacs) + 1
if number > readmax:
break
value = float(a[3])
next(line_it)
vector = []
eigenfacs.append(vector)
values.append(value)
else:
a = list(map(float, a))
vector.append(a)
return eigenfacs, values
[docs]
def normalmodes_prody(selection, cutoff=15, first=7, last=10, guide=1,
prefix='prody', states=7, factor=-1, quiet=1, *, _self=cmd):
'''
DESCRIPTION
Anisotropic Network Model (ANM) analysis with ProDy.
Based on:
http://www.csb.pitt.edu/prody/examples/dynamics/enm/anm.html
'''
import prody
first, last, guide = int(first), int(last), int(guide)
states, factor, quiet = int(states), float(factor), int(quiet)
assert first > 6
if guide:
selection = '(%s) and guide and alt A+' % (selection)
tmpsele = _self.get_unused_name('_')
_self.select(tmpsele, selection)
f = StringIO(_self.get_pdbstr(tmpsele))
conf = prody.parsePDBStream(f)
modes = prody.ANM()
modes.buildHessian(conf, float(cutoff))
modes.calcModes(last - first + 1)
if factor < 0:
from math import log
natoms = modes.numAtoms()
factor = log(natoms) * 10
if not quiet:
print(' set factor to %.2f' % (factor))
for mode in range(first, last + 1):
name = prefix + '%d' % mode
_self.delete(name)
if not quiet:
print(' normalmodes: object "%s" for mode %d' % (name, mode))
for state in range(1, states + 1):
xyz_it = iter(modes[mode - 7].getArrayNx3() * (factor *
((state - 1.0) / (states - 1.0) - 0.5)))
_self.create(name, tmpsele, 1, state, zoom=0)
_self.alter_state(state, name, '(x,y,z) = next(xyz_it) + (x,y,z)',
space={'xyz_it': xyz_it, 'next': next})
_self.delete(tmpsele)
if guide:
_self.set('ribbon_trace_atoms', 1, prefix + '*')
_self.show_as('ribbon', prefix + '*')
else:
_self.show_as('lines', prefix + '*')
cmd.extend('normalmodes_pdbmat', normalmodes_pdbmat)
cmd.extend('normalmodes_prody', normalmodes_prody)
cmd.auto_arg[0].update([
('normalmodes_pdbmat', cmd.auto_arg[0]['zoom']),
('normalmodes_prody', cmd.auto_arg[0]['zoom']),
])
# vi: ts=4:sw=4:smarttab:expandtab