This PyMol script is capable of generating a set of a molecule’s conformations, which can be further used to calculate their potential energy or any other properties. The script enables one to vary several internal degrees of freedom (DOFs) synchronously as well as scanning over a multidimensional grid in a space of DOFs.
To use the script save the code into a file (say «scan.py»), and issue the following command from PyMOL:
PyMOL> run /path/to/scan.py
Now load the molecule into PyMOL:
PyMOL> load mymol.pdb, mymol
The molecule is loaded into the object «mymol»
You need to know atoms IDs in order to define the DOFs of interest. Label atoms with their IDs using the following command:
PyMOL> label mymol, ID
You also need to define a set of atoms to move. Select the atoms of interest using a mouse or any other method. Assume the selection name is «sele».
The following command generates a set of conformations by displacing a group of the atoms along the line passing through the atoms with IDs 3 and 7:
PyMOL> scan mymol, lin(3 7) sele -1.0 3.0 5
As a result 5 new states of the object «mymol» are created. The first state of the object contains the original conformation. The second state of the «mymol» contains the conformation with «sele» group of atoms displaced along the line in backward direction by -1.0 angtsrom. Further states contain the conformations with «sele» atoms shifted gradually to +3.0 position along the reference line.
The following command rotates the «sele» atoms around the line passing through atoms with IDs 3 and 7 by 360 degrees in 36 steps (torsion angle scan):
PyMOL> scan mymol, dih(3 7) sele 0 360 37
The following command scans the covalent angle defined by atoms 1, 3 and 7. The angle varies from -30 deg to 30 deg from original position in 30 steps:
PyMOL> scan mymol, ang(1 3 7) sele -30 30 31
The following command varies two former DOFs simultaneously in 100 steps:
PyMOL> scan mymol, lin(3 7) sele -1.0 3.0 & dih(3 7) sele 0 360 100
The following command scans over a 2-dimensional grid in the space of two DOFs by displacing «sele» atoms along vector 3->7 and rotating them around vector 3->7:
PyMOL> scan mymol, lin(3 7) sele -1.0 3.0 5 & dih(3 7) sele 0 360 100
The general invocation pattern of the «scan» command is as follows:
Usage: scan obj, scan_str
obj: object to modify
scan_str: coord_type(atom1 atom2 ...) sele_name start stop [& coord_type(atom1 atom2 ...) ... [& ...]] nsteps [| coord_type ...]
Here is the scan.py code:
from pymol import cmd, stored
import re
import numpy as np
#from scandim import *
class scandim:
"Scan dimension class"
def __init__(self, obj, type, atlist, sele, start, stop, npt, friend=None):
self.obj = obj
self.type = type
self.atlist = atlist
self.sele = sele
self.vals = []
self.i = 0
for i in range(0,int(npt)):
self.vals.append(float(start) + i*(float(stop) - float(start))/(int(npt)-1))
self.friend = friend # recursive list of dependent coords
def M(self, S=1):
t = self.type
stored.m = {'H':1.01,
'C':12.01,
'N':14.01,
'O':16.,
'S':32.06,
'Fe':55.84}
stored.atm = []
for id in self.atlist:
cmd.iterate_state(S, self.obj + " & id " + str(id), "stored.atm.append(np.array([x, y, z]))")
stored.M = 0.0
if t == 'lin':
cmd.iterate_state(S, self.sele, 'stored.M += stored.m[elem]')
elif t in ['ang', 'dih']:
if t == 'ang':
v1 = stored.atm[1] - stored.atm[0]
v2 = stored.atm[2] - stored.atm[1]
v = np.cross(v1,v2)
v /= np.sqrt((v**2).sum())
else:
v = stored.atm[1] - stored.atm[0]
v /= np.sqrt((v**2).sum())
stored.v = v
stored.o = stored.atm[1]
cmd.iterate_state(S, self.sele, 'a = np.array([x,y,z])-stored.o; l = np.cross(a, stored.v); stored.M += stored.m[elem] * np.dot(l,l)')
else:
raise NameError, "Unknown coord type '%s'" % self.type
return t,stored.M
def inc(self):
if self.friend: self.friend.inc()
self.i += 1
if self.i >= len(self.vals):
self.i = 0
return True
else:
return False
def apply(self, state):
obj=self.obj
stored.atm = []
for id in self.atlist:
cmd.iterate_state(state, obj + " & id " + str(id), "stored.atm.append(np.array([x, y, z]))")
if self.type in ['lin','dih']:
o = stored.atm[1]
v = stored.atm[1] - stored.atm[0]
v /= np.sqrt((v**2).sum())
if self.type == 'lin':
v *= self.vals[self.i]
cmd.translate([v[0],v[1],v[2]], selection=self.sele, state=state, camera=0)
else:
cmd.rotate([v[0],v[1],v[2]], self.vals[self.i], selection=self.sele, origin=[o[0],o[1],o[2]], state=state, camera=0)
elif self.type == 'ang':
o = stored.atm[1]
v1 = stored.atm[1] - stored.atm[0]
v2 = stored.atm[2] - stored.atm[1]
v = np.cross(v1,v2)
v /= np.sqrt((v**2).sum())
cmd.rotate([v[0],v[1],v[2]], self.vals[self.i], selection=self.sele, origin=[o[0],o[1],o[2]], state=state, camera=0)
else:
raise NameError, "Unknown coord type '%s'" % self.type
if self.friend: self.friend.apply(state)
return
def scan(obj, scanstr, save="" ):
"""
Usage: scan obj, scan_str
obj: object to modify
scan_str: coord_type(atom1 atom2 ...) sele_name start stop [& coord_type(atom1 atom2 ...) ... [& ...]] nsteps [| coord_type ...] """
# parse input and generate list of independent coordinates
dims = []
scanpat = "^\s*(?P<typ>\w+)\s*\((?P<atoms>[\d\s]+)\)\s+(?P<mov>[\w\d_-]+)\s+(?P<start>-?\d+(\.\d+)?)\s+(?P<stop>-?\d+(\.\d+)?)\s+(?P<nstep>-?\d+)\s*$"
friendpat = "^\s*(?P<typ>\w+)\s*\((?P<atoms>[\d\s]+)\)\s+(?P<mov>[\w\d_-]+)\s+(?P<start>-?\d+(\.\d+)?)\s+(?P<stop>-?\d+(\.\d+)?)\s*$"
# split and iterate independent coords
for s in scanstr.split('|'):
# deal with dependent coords
ss = s.split('&')
# last of dependent coordinates ends with numpoints (parse using 'scanpat' pattern)
i = len(ss)-1
m = re.match(scanpat, ss[i])
if (m):
nstep = m.groupdict()['nstep']
friend = scandim(obj, m.groupdict()['typ'], m.groupdict()['atoms'].split(), m.groupdict()['mov'], m.groupdict()['start'], m.groupdict()['stop'], nstep)
else:
raise NameError, "Unrecognized scan string '%s'" % ss[i]
# form recursive list of dependent coords (using 'friendpat' pattern)
while (i):
i -= 1
m = re.match(friendpat, ss[i])
if m:
friend = scandim(obj, m.groupdict()['typ'], m.groupdict()['atoms'].split(), m.groupdict()['mov'], m.groupdict()['start'], m.groupdict()['stop'], nstep, friend)
else:
raise NameError, "Unrecognized scan string '%s'" % ss[i]
# add first dependent coord (last iterated) to the list of independent coords
dims.append(friend)
# calculate total number of points
NP = 1
for i in range(len(dims)):
NP *= len(dims[i].vals)
print "Total number of points: %d" % NP
# iterate through mesh points
for i in range(NP):
cmd.create(obj, obj, 1, i+2)
# apply transformations
for dim in dims:
dim.apply(i+2)
# save conformation if requested
if save:
log = open(save + ".txt", 'w')
print >>log, scanstr
log.close()
fname=save
for dim in dims:
fname += "_" + ("%02d" % dim.i)
fname += ".pdb"
cmd.save(fname, obj, i+2, "pdb")
# increment indices
d = len(dims)-1
while (dims[d].inc()):
d -= 1
if d < 0:
d = len(dims)-1
cmd.extend("scan", scan)