"""
Classes for creating Grid Generation and Ligand Docking DICE (Impact) input
files from user-friendly keyword/value pairs.
The classes use the mm.mmim* wrappers to convert the keywords and values into
the actual DICE input files.
Job parameters are specified by passing a dictionary of keyword/value pairs to
the class initialization method. In addition, some parameters, such as
constraints, can be specified by calling class methods.
Calling the writeSimplified() method will create the input file that can be
passed to $SCHRODINGER/glide. The value of the JOBNAME keyword will be used
as the base name for the input file.
For a full list of keywords accepted by the Gridgen and Dock classes,
respectively, see the output of ::
$SCHRODINGER/glide -gridgen-keywords
$SCHRODINGER/glide -docking-keywords
Copyright Schrodinger, LLC. All rights reserved.
"""
import io
import json
import os
import re
import sys
import numpy
from schrodinger import structure
from schrodinger.application.glide import utils
from schrodinger.application.inputconfig import InputConfig
from schrodinger.infra import canvas
from schrodinger.infra import mm
from schrodinger.infra import mmim
from schrodinger.infra import mmkv
from schrodinger.job import jobcontrol
from schrodinger.job import util
from schrodinger.structutils import analyze
from schrodinger.structutils import transform
from schrodinger.utils import fileutils
DEBUG = 'SCHRODINGER_PYTHON_DEBUG' in os.environ
#############################################################################
# GLOBAL VARIABLES
#############################################################################
# Types of the constraints (in the receptor); user-friendly strings:
HYDROPHOBIC = 'Hydrophobic'
HBOND_DONOR = 'H-bond donor'
HBOND_ACCEPTOR = 'H-bond acceptor'
METAL = 'Metal'
POSITIONAL = 'Positional'
NOE = 'NOE'
METCOORD = 'Metal coordination'
NA = 'NA'
# Map to the MMIM's int values:
mmim_cons_types = [
NA, HYDROPHOBIC, HBOND_DONOR, HBOND_ACCEPTOR, NA, METAL, POSITIONAL, NOE,
METCOORD, NA
]
# Values for AMIDE_MODE keyword:
FREE = 'free'
FIXED = 'fixed'
PENALIZE = 'penal'
TRANS = 'trans'
# Values for PRECISION keyword:
HTVS = 'HTVS'
SP = 'SP'
XP = 'XP'
# Values for POSE_OUTTYPE keyword:
PV = 'poseviewer'
LIB = 'ligandlib'
PV_SD = 'poseviewer_sd'
LIB_SD = 'ligandlib_sd'
PHASE_SUBSET = 'phase_subset'
# Force field names
OPLS_2005_VER = 14
OPLS_2005 = mm.OPLS_NAME_F14
OPLS3_VER = 16
OPLS3 = mm.OPLS_NAME_F16
# Require all constraints in group:
ALL = -1
# Section names, and old keywords not in the spec which are still allowed:
SECTIONS = {'CONSTRAINT_GROUP', 'FEATURE', 'TORCONS'}
EXTRA_KEYWORDS = {
'READBASE', 'READZIP', 'WRITEZIP', 'GRIDLIG', 'DIELECTRIC', 'ACTXRANGE',
'ACTYRANGE', 'ACTZRANGE', 'REF_LIGAND_SOURCE'
}
######################################################
# Dictionary linking old keyword names to new names:
######################################################
_old_keywords = {} # key: old keyword; value: new keyword
_old_keywords['REF_LIGAND_CT'] = 'REFLIGCT'
_old_keywords['COREATOMS'] = 'CORE_ATOMS'
_old_keywords['CORE_DEF'] = 'CORE_DEFINITION'
_old_keywords['RECVSCALE'] = 'RECEP_VSCALE'
_old_keywords['RECCCUT'] = 'RECEP_CCUT'
_old_keywords['LVDW'] = 'LIG_VSCALE'
_old_keywords['LIGCCUT'] = 'LIG_CCUT'
_old_keywords['OUTTYPE'] = 'POSE_OUTTYPE'
_old_keywords['LIGAND_SOURCE'] = 'LIGSOURCE'
_old_keywords['MODEL_FILE'] = 'RECEP_FILE'
_old_keywords['USE_REFLIGAND'] = 'USE_REF_LIGAND'
_old_keywords['REFLIGSOURCE'] = 'REF_LIGAND_SOURCE'
_old_keywords['POSTDOCKNPOSE'] = 'POSTDOCK_NPOSE'
_old_keywords['_GEOCHECK_FACTOR'] = 'GEOCHECK_FACTOR'
#######################################
# Keywords allowed in special sections
#######################################
torcons_section_allowed_keywords = set(
['ALL', 'SMARTS', 'ATOMS', 'USEVAL', 'TORVAL'])
consgroup_section_allowed_keywords = set(['USE_CONS', 'NREQUIRED_CONS'])
# CLASSES
#############################################################################
[docs]class Grid(object):
"""
Class representing the user selected grid file.
"""
[docs] def __init__(self, grid_file):
self.recep_st_entryid = None
if not grid_file or not os.path.isfile(grid_file):
self.box_center = None
self.inner_box = None
self.outer_box = None
self.recep_st = None
self.constraints = []
self.excluded_volumes = []
else:
# Read the box dimensions from the grid file:
mm.mmim_initialize(mm.error_handler)
abs_grid_file = os.path.abspath(grid_file)
gridbase = os.path.splitext(abs_grid_file)[0]
try:
with fileutils.in_temporary_directory():
self.box_center, self.inner_box, self.outer_box = mm.mmim_get_glide_gridbox(
gridbase)
self.inner_box = [int(i) for i in self.inner_box]
self.constraints = read_constraints(abs_grid_file)
except Exception as err:
msg = "Failed to read grid file: %s" % grid_file
msg += "\n\nException: %s" % str(err)
self.box_center = None
self.inner_box = None
self.outer_box = None
self.constraints = None
raise RuntimeError(msg)
finally:
mm.mmim_terminate()
self._readGridFile(grid_file)
def _readGridFile(self, grid_file):
"""
Read the list of excluded volumes from the specified grid file.
:type grid_file: str
:param grid_file: Path to the grid file (.zip or .grd)
:rtype: List of tuples
:return: Each tuple defines an excluded volume. First item is the name
(str), second item is a tuple of XYZ coordinates, third item is a
radius (float).
"""
grid_file = os.path.abspath(grid_file)
with fileutils.in_temporary_directory():
self.recep_st = utils.get_recep_structure_from_grid(grid_file)
xvol_file = utils.extract_file_from_grid(grid_file, '.gxvol')
self.excluded_volumes = utils.parse_xvol_file(xvol_file)
[docs] def isValid(self):
"""
Whether the specified grid file is valid.
"""
return (self.box_center is not None)
[docs]class Constraint:
""" Class defining a docking constraint in the receptor. """
[docs] def __init__(self, label, constype):
self._label = label
self._type = constype
def __str__(self):
""" Print the label for the constraint. """
return self._label
@property
def name(self):
"""
User-visible constraint name, as it appears in UI and in
Glide input files.
"""
return self._label
[docs] def type(self):
"""
Return the type of the (receptor) constraint as a user-friendly string.
Possible return values:
- 'Hydrophobic'
- 'H-bond donor'
- 'H-bond acceptor'
- 'Metal'
- 'Positional'
- 'NOE'
- 'Metal coordination'
"""
return self._type
class _HBondConstraint(Constraint):
keyword = 'HBOND_CONSTRAINTS'
def __init__(self, label, atoms=None):
super().__init__(label, self._TYPE)
if atoms:
self._atoms = atoms
else:
self._atoms = None
def atoms(self):
return self._atoms
def atom(self):
if self._atoms:
return self._atoms[0]
def sifValue(self):
# Unfortunately, the grid file doesn't store the "nosym" information
# directly. What we can infer however is that if there is more than one
# atom, there was no "nosym". When there is only one atom, either
# "nosym" was used, or the atom just didn't have any symmetrically
# equivalent atoms. If it's the later, "nosym" won't change anything, so
# we can add it safely even if it wasn't used in the original job.
nosym = ' nosym' if len(self._atoms) == 1 else ''
return f'{self._label} {self.atom()}{nosym}'
[docs]class DonorConstraint(_HBondConstraint):
"""
H-Bond constraint that was specified as a *receptor* donor atom.
A ligand must have an ACCEPTOR group that is making an interaction with
this atom to satisfy this constraint.
"""
_TYPE = HBOND_DONOR
[docs]class AcceptorConstraint(_HBondConstraint):
"""
H-Bond constraint that was specified as a *receptor* acceptor atom.
A ligand must have a DONOR group that is making an interaction with
this atom to satisfy this constraint.
"""
_TYPE = HBOND_ACCEPTOR
[docs]class PositionalConstraint(Constraint):
keyword = 'POSIT_CONSTRAINTS'
[docs] def __init__(self, label, x=None, y=None, z=None, radius=None):
Constraint.__init__(self, label, POSITIONAL)
self.x = x
self.y = y
self.z = z
self.radius = radius
[docs] def sifValue(self):
return f'{self._label} {self.x} {self.y} {self.z} {self.radius}'
[docs]class NOEConstraint(Constraint):
keyword = 'NOE_CONSTRAINTS'
[docs] def __init__(self, label, x=None, y=None, z=None, rmin=None, rmax=None):
Constraint.__init__(self, label, NOE)
self.x = x
self.y = y
self.z = z
self.rmin = rmin
self.rmax = rmax
[docs] def sifValue(self):
return (f'{self._label} {self.x} {self.y} {self.z} '
f'{self.rmin} {self.rmax}')
[docs]class MetCoordConstraint(Constraint):
keyword = 'METCOORD_CONSTRAINTS'
sites_keyword = 'METCOORD_SITES'
[docs] def __init__(self, label, xmetc, ymetc, zmetc, rmetc):
Constraint.__init__(self, label, METCOORD)
self.xmetc = xmetc
self.ymetc = ymetc
self.zmetc = zmetc
self.rmetc = rmetc
self.ncoordsites = len(xmetc)
[docs] def sifValue(self):
return f'{self._label} {self.ncoordsites}'
[docs] def siteSifValue(self):
ret = []
for x, y, z, r in zip(self.xmetc, self.ymetc, self.zmetc, self.rmetc):
ret.append(f'{x} {y} {z} {r}')
return ret
[docs]class HydrophobicConstraint(Constraint):
# We no longer support hydrophobic constraints, but we'll leave this
# class here so the read_constraints function will still work with old
# grids.
[docs] def __init__(self, label):
Constraint.__init__(self, label, HYDROPHOBIC)
#############################################################################
#############################################################################
[docs]class FeaturePattern(object):
""" Class defining a constraint feature pattern. """
[docs] def __init__(self, smarts, atoms, exclude=False, nfill=None):
# nfill - number of ligand atoms needed to satisfy hydrophobic
# constraint. Ignored for non-hydrophobic constraints.
# FIXME Isn't <nfill> a contraint-level attribute, not pattern-level?
self.smarts = smarts
self.atoms = atoms
self.exclude = exclude
# num atoms satisfy (hydrophobic constraints only). None=default
self.nfill = nfill
[docs]class Feature(object):
"""
A list of FeaturePattern objects. Each pattern has these properties:
smarts, atoms, exclude.
"""
[docs] def __init__(self, patterns=None):
"""
:param patterns: List of patterns associated with this feature. A
shallow copy of the patterns list will be made for the new Feature
object.
:type patterns: list
"""
self.patterns = []
if patterns:
self.patterns += patterns
[docs]class ConstraintGroup:
""" Class defining a constraint group. """
[docs] def __init__(self, parent, constraints=None, nrequired=ALL):
self.parent = parent
if constraints:
self.constraints = constraints
else:
self.constraints = [
] # labels of constraints to be included in this group
self.nrequired = nrequired # -1 means ALL are required
self.group_handle = None
[docs] def addConstraint(self, label):
if label not in self.constraints:
self.constraints.append(label)
[docs] def setup(self, cons_labels, feature_sets):
# TODO: Update this function to use this module's
# Feature class for feature_sets.
mmim_handle = self.parent.mmim_handle
self.group_handle = mm.mmim_handle_glide_cons_new_group(mmim_handle)
for fulllabel in self.constraints:
# Make sure that all requested constraints are in the *cons file:
if fulllabel in cons_labels:
# Use default feature
conslabel = fulllabel
# TODO: Update this to use the Feature class in this module.
feature_patterns = [] # List of FeaturePattern objects
else:
# Try to use feature name after ":"
s = fulllabel.split(':')
conslabel = ':'.join(s[0:-1])
featurename = s[-1]
if conslabel not in cons_labels:
msg = "glide.py: constraint '%s' is not defined for the grid" % fulllabel
msg += '\nAvailable constraints:' + str(cons_labels)
raise RuntimeError(msg)
if featurename not in feature_sets:
msg = "glide.py: feature '%s' is not defined" % featurename
if feature_sets:
msg += '\nAvailable features: %s' % ', '.join(
feature_sets)
else:
msg += '\nNo features defined in the FEATURE section.'
raise RuntimeError(msg)
feature_patterns = feature_sets[featurename]
# Determine the type of this (receptor) constraint:
type_int = mm.mmim_handle_glide_cons_type_from_label(
mmim_handle, conslabel)
# Map constraint type strings to the MMIM constant value:
constype = mmim_cons_types[type_int]
ignore = 0 # Not used by Glide (Phase only)
headpathandle = -1 # Use default handle
prevpathandle = -1 # No previous pattern yet
smallest_nfill = None # Number of atoms to fill (hydrophobic cons only)
# NOTE: We need to use the smallest number of all available features.
# If none had <nfill> set, we need to use 1.
for pat in feature_patterns:
pathandle = mm.mmim_handle_glide_cons_new_pattern(mmim_handle,\
pat.smarts, len(pat.atoms), pat.atoms, int(pat.exclude), ignore)
if prevpathandle == -1:
# This is the first pattern
headpathandle = pathandle
else:
mm.mmim_handle_glide_cons_link_pattern(
mmim_handle, prevpathandle, pathandle)
prevpathandle = pathandle
if pat.nfill:
pat_nfill = int(pat.nfill)
if pat_nfill < smallest_nfill or smallest_nfill is None:
smallest_nfill = pat_nfill
nfill = 0 # Should be used for non-hydrophobic constraints
if constype == HYDROPHOBIC:
raise ValueError(
"Hydrophobic constraints are no longer supported")
mm.mmim_handle_glide_cons_add_cons_to_group(mmim_handle,
self.group_handle,
conslabel, nfill,
headpathandle)
if self.nrequired == ALL:
# All are required
nreq = len(self.constraints)
else:
nreq = int(self.nrequired)
mm.mmim_handle_glide_cons_set_group_nreq(mmim_handle, self.group_handle,
nreq)
[docs] def cleanup(self):
""" Free the memory associated with this group """
mm.mmim_handle_glide_cons_delete_group(self.parent.mmim_handle,
self.group_handle)
#############################################################################
#############################################################################
[docs]class GlideJob(InputConfig):
"""
Class for writing Glide input files - either Gridgen or Docking. The
keywords passed are used to set MMIM directly with no processing.
"""
[docs] def __init__(self,
jobtype=None,
specs=None,
keywords=None,
createhandle=False):
"""
Specify a dictionary of keywords, where the key is the short keyword
name and the value is what is should be set to. For indexed keywords,
the value needs to be a list.
It is also possible to specify a file path to a simplified input file
that contains the keywords.
"""
self.jobname = None
self._jobtype = jobtype # mm.MMIM_DOCKING_JOB or mm.MMIM_GRIDGEN_JOB; if "None," see below
if createhandle:
self.mmim_handle = mm.mmim_new()
self.mmim_dict = mmim.MMIMDict(self.mmim_handle)
else:
self.mmim_handle = None
self.mmim_dict = None
if isinstance(keywords, dict):
new_keywords = {}
# Make a copy and replace old keywords with new equivalents:
for key, value in keywords.items():
if key in _old_keywords:
new_key = _old_keywords[key]
import warnings
msg = 'glide.py: keyword %s is obsolete. Use %s instead.' % (
key, new_key)
warnings.warn(msg, DeprecationWarning, stacklevel=2)
key = new_key
new_keywords[key] = value
InputConfig.__init__(self, new_keywords, specs)
else: # Input file or no keywords specified
InputConfig.__init__(self, keywords, specs)
self._mungeBeforeValidation()
# Will complain on invalid values:
if specs:
# Validate the InputConfig to specs (convert values to
# expected types and report errors):
self.validateValues(preserve_errors=True, copy=False)
self.defaults_set = set(self.defaults)
if jobtype is None:
tmp = self.get("MMIMJOBTYPE")
if tmp:
self._jobtype = int(tmp)
elif self.get('CG_MODE'):
self._jobtype = mm.MMIM_COMBIGLIDE_JOB
else:
self._jobtype = mm.MMIM_DOCKING_JOB
if DEBUG:
print('**************')
print('GlideJob:', self)
print('**************')
def _mungeBeforeValidation(self):
# hook for hacks needed for backward compatibility and such.
pass
[docs] def validateValues(self, *args, **kwargs):
"""
Validate the keywords against the spec. Unlike the parent
implementation, raise a ValueError if any unknown keywords are present.
"""
super(GlideJob, self).validateValues(*args, **kwargs)
unknown_keywords = []
for k in self.extra_values:
# Keywords that begin with underscore are allowed because that's
# the convention used for unsupported/experimental keywords.
if k not in EXTRA_KEYWORDS and k.split(':')[0] not in SECTIONS \
and not k.startswith('_'):
unknown_keywords.append(k)
if unknown_keywords:
raise ValueError("Unknown SIF keywords: %s" %
(', '.join(unknown_keywords)))
def _applyKeyword(self, keyname, value):
"""
Calls appropriate MMIM APIs to "set" the <keyname> keyword
to the specified value.
Conversion to the correct type is performed if necessary.
"""
if type(value) == type(''):
try:
value = mmim.convert_string(keyname, value)
except ValueError as e:
# Re-raise a bad-key name error:
raise ValueError(e)
if DEBUG:
print('SETTING KEYWORD: "%s" to: %s' % (keyname, value))
try:
self.mmim_dict[keyname] = value
except TypeError as err:
msg = "Keyword %s (value: %s): %s" % (keyname, value, err)
if DEBUG:
print(msg)
raise
else:
raise TypeError(msg)
except mm.MmException as err:
raise Exception(str(err))
def _getAndApplyKeyword(self, username, default=None):
"""
Calls appropriate MMIM APIs to "set" the <username> keyword to the
value that the user specified for this keyword. If the keyword was
not specified, but 'default' is provided, set the value to
'default'.
Returns the value that is set or None if not set.
"""
value = None
try:
if username not in self.defaults_set:
value = self[username]
except KeyError: # Not specified
pass
if value is None:
if default is None:
return None # Did not set the keyword
else:
value = default # Use default
self._applyKeyword(username, value)
return value
[docs] def get(self, key, default=None):
try:
val = self[key]
if val is None:
return default
return val
except KeyError:
return default
[docs] def getAllowedRange(self, key):
"""
Will return a tuple of (min, max) for the allowed range of value. Either
or both of these values may be None.
:return: The allowed min/max range.
:rtype: Tuple of 2 ints/floats.
"""
for spec in read_keywords_file():
if spec['keyword'] == key.lower():
min_value = spec.get('min')
max_value = spec.get('max')
return (min_value, max_value)
raise ValueError('Keyword "%s" is not defined in MMIM' % key)
[docs] def cleanup(self):
""" Clean up constraints objects. """
[docs] def writeSimplified(self, filename=None, yesno=False):
"""
Writes a simplified input file to "<jobname>.inp".
This input file needs to be run via $SCHRODINGER/glide.
"""
self.jobname = self.get('JOBNAME')
self[
'JOBNAME'] = None # So that JOBNAME would not be written to input file
if not filename:
filename = self.jobname + '.in'
else:
if hasattr(filename, 'endswith') \
and not filename.endswith('.inp') \
and not filename.endswith('.in'):
raise ValueError("Glide input file should end with .in or .inp")
self.writeInputFile(filename, ignore_none=True, yesno=yesno)
self['JOBNAME'] = self.jobname
return filename
[docs] def applyContextDependentDefaults(self):
"""
Apply context-dependent defaults for a Glide job. It modifies the job
object, but marks the modified keywords as "default" so they won't
be written by writeSimplified().
"""
defaults_set = set(self.defaults)
for key, val, defaults in self._defaultsTable():
if self[key] == val or (callable(val) and val(self[key])):
for set_key, set_val, override in defaults:
if override or set_key in defaults_set:
self[set_key] = set_val
self.defaults.append(set_key)
[docs] def checkConflicts(self):
"""
Check for invalid keyword combinations. Raises a ValueError if any
are found.
"""
for key1, val1, key2, val2 in self._conflictsTable():
v1, v2 = self[key1], self[key2]
k1 = val1(v1) if callable(val1) else v1 == val1
k2 = val2(v2) if callable(val2) else v2 == val2
if k1 and k2:
raise ValueError('Invalid Glide keyword combination: %s %s '
'is incompatible with %s %s' %
(key1, v1, key2, v2))
[docs] def setupJson(self, jsdict):
"""
Additional setup to do before writing a JSON file. Base implementation
does nothing. This method is expected not to modify 'self', but only
jsdict, which is the dictionary of keywords that ultimately will be
written to the JSON file.
"""
[docs] def checkForceFieldLicense(self, jsdict):
"""
If we are using the default forcefield (OPLS3) and don't have an
OPLS license, fall back to using OPLS_2005.
"""
val = jsdict['FORCEFIELD']
is_opls3 = mm.opls_name_to_version(val) == OPLS3_VER
if ('FORCEFIELD' in self.defaults and is_opls3 and
not mm.mmffld_license_ok(OPLS3_VER)):
jsdict['FORCEFIELD'] = OPLS_2005
[docs] def writeJSON(self, filename=None):
"""
Write the keywords as a JSON representation. Unlike writeSimplified,
this includes the values for all the keywords, not only the
non-default ones.
"""
self.jobname = self.get('JOBNAME')
if not self.jobname:
raise RuntimeError("glide.py: JOBNAME keyword must be specified")
if filename is None:
filename = self['JOBNAME'] + '.js'
else:
if hasattr(filename, 'endswith') and not filename.endswith('.js'):
raise ValueError("Glide input file should end with .js")
self['MMIMJOBTYPE'] = self._jobtype
self.applyContextDependentDefaults()
self.checkConflicts()
jsdict = dict(self)
if int(os.environ.get('SCHRODINGER_GLIDE_SIFTEST', 0)) == 2:
jsdict['_DICELESS_MODE'] = True
self.setupJson(jsdict)
self.checkForceFieldLicense(jsdict)
with open(filename, 'w') as fh:
json.dump(jsdict, fh, indent=1, sort_keys=True)
return filename
[docs]class Gridgen(GlideJob):
"""
Class for writing Glide Grid Generation input files. Currently, the
receptor input file for a Gridgen job (NOT created by this class) must be
named "<jobname>.mae".
"""
[docs] def __init__(self, keywords):
"""
Specify a dictionary of keywords, where the key is the short keyword
name and the value is what is should be set to. For indexed keywords,
the value needs to be a list.
It is also possible to specify a file path to a simplified input file
that contains the keywords.
Constraint sites may be set via the following keywords:
HBOND_CONSTRAINTS
METAL_CONSTRAINTS
POSIT_CONSTRAINTS
NOE_CONSTRAINTS
METCOORD_CONSTRAINTS
METCOORD_SITES
:type keywords: dict or file path
:param keywords: Either a dictionary of keywords or a simplified input file path.
:raise RuntimeError: If no job name specified.
"""
# Setup the configuration specifications:
GlideJob.__init__(self, mm.MMIM_GRIDGEN_JOB,
validation_specs('gridgen'), keywords)
self['JOBTYPE'] = 'GRIDGEN'
self.defaults.append('JOBTYPE')
self._hbond_labels = []
self._metal_labels = []
self._posit_labels = []
self._noe_labels = []
self._metcoord_labels = []
def _defaultsTable(self):
return [
# (key, val, defaults)
(
'PEPTIDE',
True,
[
# (set_key, set_val, override)
('NDIR', mm.GLIDE_PEPTIDE_NDIR, True),
]),
('USECOMPMAE', False, [
('USECOMPMAE', True, False),
]),
]
def _conflictsTable(self):
# Each line describes a potential conflict. If keyword1 has value1
# and keyword2 has value2, there is a conflict. A "value" can also
# be a callable; in that case we use the result of calling it with
# the actual value of the keyword.
return [
# (keyword1, value1, keyword2, value2)
('GLIDE_RECEP_ASLSCALE', True, 'GLIDE_RECEP_MAESCALE', True),
('USEFLEXASL', True, 'USEFLEXMAE', True),
('HBOND_ACCEP_HALO', True, 'HBOND_DONOR_HALO', True),
]
def _gridfile(self):
"""
Return the grid file to use, normally just the value of the GRIDFILE
keyword, but sometimes based on legacy keywords for compatibility.
"""
gridfile = self.get('GRIDFILE') # If None, use <jobname>.zip/.grd
if not gridfile:
if not self.jobname:
raise RuntimeError(
"glide.py: Please specify either GRIDFILE or JOBNAME keyword. Otherwise I don't know where to save the generated grids."
)
# Support for old WRITEZIP keyword:
writezip = self.get('WRITEZIP', True)
# FIXME private function call. This shouldn't be needed, as ideally
# section['WRITEZIP'] should always return a bool type
if mmkv._convert_to_bool(writezip):
gridfile = self.jobname + '.zip'
else:
gridfile = self.jobname + '.grd'
return gridfile
def _recep_file(self):
if self.get('RECEP_FILE'):
recfile = self['RECEP_FILE']
elif self['USECOMPMAE']:
recfile = "%s.maegz" % self['JOBNAME']
else:
recfile = "%s.mae" % self['JOBNAME']
return jobcontrol.get_runtime_path(recfile)
def _grid_center_from_asl(self, grid_center_asl, delete_atoms=False):
recfile = self._recep_file()
st = self._receptorStructure()
atoms = analyze.evaluate_asl(st, grid_center_asl)
xsum = 0.0
ysum = 0.0
zsum = 0.0
for anum in atoms:
xsum += st.atom[anum].x
ysum += st.atom[anum].y
zsum += st.atom[anum].z
natoms = float(len(atoms))
xcent = xsum / natoms
ycent = ysum / natoms
zcent = zsum / natoms
grid_center = (xcent, ycent, zcent)
if delete_atoms:
st.deleteAtoms(atoms)
recfile_basename = fileutils.get_basename(recfile)
recfile_pathame, recfile_ext = fileutils.splitext(recfile)
recfile_noligand = recfile_basename + "_recep" + recfile_ext
st.write(recfile_noligand)
recfile = recfile_noligand
return grid_center, recfile
def _grid_center(self):
"""
Return (grid_center, recep_file). The grid center may come from
any of three places, in order of increasing precedence:
1) via a ligand number reference in RECEP_FILE; it will be deleted.
2) via an ASL evaluated against RECEP_FILE
3) via an explicit GRID_CENTER keyword value (or [XYZ]CENT).
The recep_file normally comes from RECEP_FILE or JOBNAME, but it can
be a munged file if the LIGAND_MOLECULE option is used.
"""
recfile = self._recep_file()
grid_center = self.get('GRID_CENTER')
lig_st = self._ligandStructure()
if lig_st:
grid_center = list(transform.get_centroid(lig_st)[:3])
ligand_molecule = self.get('LIGAND_MOLECULE')
if ligand_molecule:
ligand_asl = "mol.num %d" % ligand_molecule
grid_center, recfile = self._grid_center_from_asl(ligand_asl,
delete_atoms=True)
grid_center_asl = self.get('GRID_CENTER_ASL')
if grid_center_asl:
grid_center, recfile = self._grid_center_from_asl(grid_center_asl)
if 'GRID_CENTER' in self.defaults_set:
if self.get('XCENT') and self.get('YCENT') and self.get('ZCENT'):
grid_center = [self['XCENT'], self['YCENT'], self['ZCENT']]
import warnings
msg = 'glide.py: keywords XCENT, YCENT, and ZCENT are obsolete. Use GRID_CENTER instead.'
warnings.warn(msg, DeprecationWarning, stacklevel=2)
else:
grid_center = self.get('GRID_CENTER')
grid_center = list(map(float, grid_center)) # Convert to floats
return grid_center, recfile
def _inner_box(self):
# If INNERBOX keyword is specified:
inner_box = self.get('INNERBOX')
if inner_box:
if type(inner_box) != type([]): # Not a list
inner_box = [inner_box, inner_box, inner_box]
if len(inner_box) == 1:
# list has one item in it
inner_box *= 3
elif len(inner_box) != 3:
raise ValueError("INNERBOX needs either 1 or 3 values")
# This will convert floats and strings to ints (Ev:71130):
# (works on "int" strings AND "float" strings:
inner_box = list(map(int, list(map(float, inner_box))))
return inner_box
def _outer_box(self):
if 'OUTERBOX' in self.defaults_set and self.get('LIGAND_INDEX'):
st = self._ligandStructure()
_, outer_box = get_default_box_size(st)
else:
outer_box = self.get('OUTERBOX')
if type(outer_box) != type([]): # Not a list
outer_box = [outer_box, outer_box, outer_box]
if len(outer_box) == 1:
# list has one item in it
outer_box *= 3
elif len(outer_box) != 3:
raise ValueError("OUTERBOX needs either 1 or 3 values")
outer_box = list(map(float, outer_box)) # Convert to list of FLOATS
return outer_box
[docs] def setupJson(self, jsdict):
jsdict['GRIDFILE'] = self._gridfile()
jsdict['GRID_CENTER'], jsdict['RECEP_FILE'] = self._grid_center()
jsdict['INNERBOX'] = self._inner_box()
jsdict['OUTERBOX'] = self._outer_box()
jsdict['ACTXRANGE'], jsdict['ACTYRANGE'], jsdict['ACTZRANGE'] = \
jsdict['OUTERBOX']
self._jsonSetupConstraints(jsdict)
def _useSimplifiedConstraints(self):
"""
Return True is using "simplified constraints" keywords such as
HBOND_CONSTRAINTS. Raises a ValueError exception if simplified keywords
are used together with "low-level" keywords such as GLIDECONSATOMS.
"""
low_level_keywords = [
'GLIDE_NTOTALCONS', 'GLIDECONSNAMES', 'GLIDE_NUMPOSITCONS',
'GLIDE_CONS_XPOS', 'GLIDE_CONS_YPOS', 'GLIDE_CONS_ZPOS',
'GLIDE_CONS_RPOS', 'GLIDE_NUMNOECONS', 'GLIDE_CONS_XNOE',
'GLIDE_CONS_YNOE', 'GLIDE_CONS_ZNOE', 'GLIDE_CONS_RNOEMIN',
'GLIDE_CONS_RNOEMAX', 'GLIDE_NUMMETCOORDCONS',
'GLIDE_NUMMETCOORDSITES', 'GLIDE_CONS_XMETCOORD',
'GLIDE_CONS_YMETCOORD', 'GLIDE_CONS_ZMETCOORD',
'GLIDE_CONS_RMETCOORD', 'GLIDECONSATOMS', 'GLIDECONSUSESYMATOMS'
]
high_level_keywords = [
'HBOND_CONSTRAINTS', 'METAL_CONSTRAINTS', 'POSIT_CONSTRAINTS',
'NOE_CONSTRAINTS', 'METCOORD_CONSTRAINTS', 'METCOORD_SITES'
]
used_low_level = [k for k in low_level_keywords if self[k]]
used_high_level = [k for k in high_level_keywords if self[k]]
if used_high_level and used_low_level:
raise ValueError('Invalid combination of simplified and '
'low-level constraint keywords: %s vs %s' %
(used_high_level, used_low_level))
return bool(used_high_level)
def _setupPositConstraints(self, jsdict):
for cons_string in jsdict['POSIT_CONSTRAINTS']:
label, x, y, z, radius = cons_string.split()
jsdict['GLIDE_CONS_XPOS'].append(float(x))
jsdict['GLIDE_CONS_YPOS'].append(float(y))
jsdict['GLIDE_CONS_ZPOS'].append(float(z))
jsdict['GLIDE_CONS_RPOS'].append(float(radius))
jsdict['GLIDECONSNAMES'].append(label)
jsdict['GLIDE_NUMPOSITCONS'] += 1
jsdict['GLIDE_NTOTALCONS'] += 1
def _setupNoeConstraints(self, jsdict):
for cons_string in jsdict['NOE_CONSTRAINTS']:
label, x, y, z, rmin, rmax = cons_string.split()
jsdict['GLIDE_CONS_XNOE'].append(float(x))
jsdict['GLIDE_CONS_YNOE'].append(float(y))
jsdict['GLIDE_CONS_ZNOE'].append(float(z))
jsdict['GLIDE_CONS_RNOEMIN'].append(float(rmin))
jsdict['GLIDE_CONS_RNOEMAX'].append(float(rmax))
jsdict['GLIDECONSNAMES'].append(label)
jsdict['GLIDE_NUMNOECONS'] += 1
jsdict['GLIDE_NTOTALCONS'] += 1
def _setupMetcoordConstraints(self, jsdict):
for cons_string in jsdict['METCOORD_SITES']:
x, y, z, r = cons_string.split()
jsdict['GLIDE_CONS_XMETCOORD'].append(float(x))
jsdict['GLIDE_CONS_YMETCOORD'].append(float(y))
jsdict['GLIDE_CONS_ZMETCOORD'].append(float(z))
jsdict['GLIDE_CONS_RMETCOORD'].append(float(r))
for cons_string in jsdict['METCOORD_CONSTRAINTS']:
label, ncoordsites = cons_string.split()
jsdict['GLIDE_NUMMETCOORDSITES'].append(int(ncoordsites))
jsdict['GLIDECONSNAMES'].append(label)
jsdict['GLIDE_NUMMETCOORDCONS'] += 1
jsdict['GLIDE_NTOTALCONS'] += 1
sumnsites = sum(jsdict['GLIDE_NUMMETCOORDSITES'])
countnsites = len(jsdict['GLIDE_CONS_XMETCOORD'])
if sumnsites != countnsites:
raise ValueError("%i metal-coordination sites specified, "
"%i used by constraints." %
(countnsites, sumnsites))
def _setupReceptorAtomBasedConstraints(self,
jsdict,
keyword,
symmetrize_all=False):
for cons_string in jsdict[keyword]:
label, atomindex, symmetrize = \
self._parseConsString(cons_string)
symmetrize = symmetrize_all or symmetrize
jsdict['GLIDECONSATOMS'].append(int(atomindex))
jsdict['GLIDECONSUSESYMATOMS'].append(symmetrize)
jsdict['GLIDECONSNAMES'].append(label)
jsdict['GLIDE_NTOTALCONS'] += 1
def _jsonSetupConstraints(self, jsdict):
if not self._useSimplifiedConstraints():
return # low-level constraints are in use
# NOTE: The order of the following setup calls is important because the
# backend expects the GLIDECONSNAMES array to have the constraints
# sorted by type, in this order.
jsdict['GLIDECONS'] = True
jsdict['GLIDE_NTOTALCONS'] = 0
jsdict['GLIDECONSNAMES'] = []
self._setupPositConstraints(jsdict)
self._setupNoeConstraints(jsdict)
self._setupMetcoordConstraints(jsdict)
self._setupReceptorAtomBasedConstraints(jsdict, 'HBOND_CONSTRAINTS')
self._setupReceptorAtomBasedConstraints(jsdict,
'METAL_CONSTRAINTS',
symmetrize_all=True)
def _receptorStructure(self):
return structure.Structure.read(self._recep_file())
def _ligandStructure(self):
"""
Return the ligand Structure when using the LIGAND_INDEX keyword, else
None.
"""
index = self.get('LIGAND_INDEX')
if index:
return structure.Structure.read(self._recep_file(), index=index)
else:
return None
def _parseConsString(self, cons_string):
"""
Parse a string with the format "<label> <atom_index> [nosym]" or
"<label> <asl> [nosym]" and return the tuple (label, atomindex,
symmetrize). Raises ValueError if the string is invalid for any reason.
"""
try:
label, rest = cons_string.split(None, 1)
except ValueError:
raise ValueError("Error parsing constraint spec '%s': "
"must be '<label> <atom_index|asl-expr> [nosym]'" %
(cons_string,))
expr, n = re.subn(r'\s+nosym$', '', rest, flags=re.I)
symmetrize = n == 0
try:
atomindex = int(expr)
return label, atomindex, symmetrize
except ValueError:
try:
st = self._receptorStructure()
atoms = analyze.evaluate_asl(st, expr)
except mm.MmException as e:
raise ValueError("Error parsing constraint ASL '%s': %s" %
(cons_string, e.mmlib_message))
if len(atoms) != 1:
raise ValueError(
"ASL expression for constraint "
"must match exactly one atom: '%s' matched %s" %
(expr, atoms))
return label, atoms[0], symmetrize
def _setupHbondConstraints(self):
hbond_cons_list = self.get('HBOND_CONSTRAINTS', [])
if hbond_cons_list:
st = self._receptorStructure()
# Ev:72130:
self._applyKeyword('GLIDECONS', True) # This line is for DICE file
for cons_string in hbond_cons_list:
label, atomindex, symmetrize = \
self._parseConsString(cons_string)
if DEBUG:
print('HBOND CONSTRAINT:', label, atomindex)
# Auto detect if the receptor atom is a donor or an acceptor:
try:
atom = st.atom[atomindex]
except IndexError:
raise ValueError(
"_newHbondConstraint(): Invalid atom index: %i" %
atomindex)
if atom.element == 'H':
constype = mm.MMIM_GLIDECONS_TYPE_DONOR
elif atom.element in ['O', 'N', 'S']:
constype = mm.MMIM_GLIDECONS_TYPE_ACCEP
else:
raise ValueError(
"_newHbondConstraint(): Invalid atom element: %s (atom %i)"
% (atom.element, atomindex))
# FIXME populate this info dynamically from this file:
# mmshare/mmlibs/mmim/mmim_keywords.yaml
if DEBUG:
print("DEBUG Calling mmim_handle_glide_cons_new_hbond()")
print(" ATOM:", atomindex, "TYPE:", constype, "LABEL:",
label)
mm.mmim_handle_glide_cons_new_hbond(self.mmim_handle, atomindex,
constype, label, symmetrize)
self._hbond_labels.append(label)
def _setupConstraints(self):
self._setupHbondConstraints()
for cons_string in self.get('METAL_CONSTRAINTS', []):
label, atomindex = cons_string.split()
mm.mmim_handle_glide_cons_new_metal(self.mmim_handle,
float(atomindex), label)
self._metal_labels.append(label)
# Ev:72130:
self._applyKeyword('GLIDECONS', True) # This line is for DICE file
for cons_string in self.get('POSIT_CONSTRAINTS', []):
label, x, y, z, radius = cons_string.split()
mm.mmim_handle_glide_cons_new_posit(self.mmim_handle, float(x),
float(y), float(z),
float(radius), label)
self._posit_labels.append(label)
# Ev:72130:
self._applyKeyword('GLIDECONS', True) # This line is for DICE file
for cons_string in self.get('NOE_CONSTRAINTS', []):
label, x, y, z, rmin, rmax = cons_string.split()
mm.mmim_handle_glide_cons_new_noe(self.mmim_handle, float(x),
float(y), float(z), float(rmin),
float(rmax), label)
self._noe_labels.append(label)
# Ev:72130:
self._applyKeyword('GLIDECONS', True) # This line is for DICE file
sumnsites = 0
countnsites = 0
xmetc = []
ymetc = []
zmetc = []
rmetc = []
for cons_string in self.get('METCOORD_SITES', []):
x, y, z, r = cons_string.split()
countnsites += 1
xmetc.append(float(x))
ymetc.append(float(y))
zmetc.append(float(z))
rmetc.append(float(r))
for cons_string in self.get('METCOORD_CONSTRAINTS', []):
label, ncoordsites = cons_string.split()
ncoordsites = int(ncoordsites)
if sumnsites + ncoordsites > countnsites:
raise ValueError(
"_setupConstraints(): metal-coordination constraint specifies too many sites: %i (only %i remaining)"
% (ncoordsites, countnsites - sumnsites))
xs = xmetc[sumnsites:sumnsites + ncoordsites]
ys = ymetc[sumnsites:sumnsites + ncoordsites]
zs = zmetc[sumnsites:sumnsites + ncoordsites]
rs = rmetc[sumnsites:sumnsites + ncoordsites]
mm.mmim_handle_glide_cons_new_metcoord(self.mmim_handle,
ncoordsites, xs, ys, zs, rs,
label)
self._metcoord_labels.append(label)
self._applyKeyword('GLIDECONS', True)
sumnsites += ncoordsites
if sumnsites != countnsites:
raise ValueError("_setupConstraints(): %i metal-coordination sites specified, %i used by constraints." % \
(countnsites, sumnsites))
[docs] def cleanup(self):
""" Free the memory associated with constraings """
for label in self._hbond_labels:
mm.mmim_handle_glide_cons_delete_hbond(self.mmim_handle, label)
for label in self._metal_labels:
mm.mmim_handle_glide_cons_delete_metal(self.mmim_handle, label)
for label in self._posit_labels:
mm.mmim_handle_glide_cons_delete_posit(self.mmim_handle, label)
for label in self._noe_labels:
mm.mmim_handle_glide_cons_delete_noe(self.mmim_handle, label)
[docs]class Dock(GlideJob):
"""
Class for generating Glide docking input files.
"""
[docs] def __init__(self, keywords, jobtypes=None):
"""
Specify a dictionary of keywords, where the key is the short keyword
name and the value is what is should be set to. For indexed keywords,
the value needs to be a list.
It is also possible to specify a file path to a simplified input file
that contains the keywords.
:type keywords: dict or file path
:param keywords: Either a dictionary of keywords or a simplified input file path.
"""
# Setup the configuration specifications:
# Pass a "jobtype" of "None" here, so which GlideJob.__init__()
# will replace with keywords['MMIMJOBTYPE'] if that's set, of
# mm.MMIM_DOCKING_JOB otherwise. This allows "special" docking
# job types (e.g. Combiglide, core-hopping) to override the
# default of mm.MMIM_DOCKING_JOB. But we can't get the
# 'MMIMJOBTYPE' setting directly from "keywords" here, because
# that might be a file path (string) rather than a dict.
if jobtypes is None:
jobtypes = ['docking', 'combiglide', 'pcalign', 'chdock']
specs = validation_specs(jobtypes)
GlideJob.__init__(self, None, specs, keywords)
self['JOBTYPE'] = 'DOCKING'
self.defaults.append('JOBTYPE')
self._cons_groups = [
] # List of constraint group objects (populated in setup() and cleared in cleanup() )
# TODO: This should be updated to use the Feature class in this module.
self._feature_sets = {}
# All constraint information is actually stored in the Dock object as
# self['CONSTRAINT_GROUP:X'] sections
# Convert USE_CONS & NREQUIRED_CONS to use groups:
use_cons = self.get('USE_CONS')
if use_cons:
if type(use_cons) == type(''):
use_cons = use_cons.split(' ') # Convert to list
nrequired_str = self['NREQUIRED_CONS']
if nrequired_str == 'ALL' or nrequired_str == 'all' or nrequired_str == 'All':
nrequired = ALL
else:
try:
nrequired = int(nrequired_str)
except:
msg = 'ERROR: Invalid NREQUIRED_CONS value: %s' % nrequired_str
raise RuntimeError(msg)
self['CONSTRAINT_GROUP:1'] = {
'USE_CONS': use_cons,
'NREQUIRED_CONS': nrequired
}
del self['USE_CONS']
del self['NREQUIRED_CONS']
def _mungeBeforeValidation(self):
# GLIDE-3165: be nice if someone uses LIGANDFILE as a list
if isinstance(self.get('LIGANDFILE'), list):
self['LIGANDFILES'] = self.pop('LIGANDFILE')
[docs] def setLigRange(self, start, end):
"""
Shortcut for setting the range of ligands to be read from the
LIGANDFILE.
:param end: If zero, means to go to the end of the file.
"""
self['LIGAND_START'] = start
self['LIGAND_END'] = end
def _getGridProperties(self):
"""
Return the properties from the _recep.mae in the grid file.
If the grid file is not available, return an empty dict.
"""
if self['GRIDFILE']:
gridfile = jobcontrol.get_runtime_path(self['GRIDFILE'])
try:
st = utils.get_recep_structure_from_grid(gridfile)
return dict(st.property)
except IOError:
pass
return {}
def _checkPeptide(self, val):
"""
Return True if the PEPTIDE keyword is true or the grid archive
is available and was generated using peptide mode.
"""
if val:
return True
else:
props = self._getGridProperties()
return props.get(mm.GLIDE_GRID_PEPTIDE_MODE, False)
def _checkFlexOH(self):
"""
Return True if the grid file for a docking job exists and
was generated using rotatable hydroxyl groups.
"""
props = self._getGridProperties()
return props.get(mm.GLIDE_FLEXR_NFLEX, 0) > 0
def _autoLigFormatOneFile(self, filename):
fmt = fileutils.get_structure_file_format(filename)
supported_formats = {
fileutils.MAESTRO,
fileutils.SD,
fileutils.MOL2,
fileutils.SMILES,
fileutils.SMILESCSV,
}
if fmt in supported_formats:
return fmt
elif fmt == fileutils.PDB:
raise ValueError("PDB format is no longer supported by Glide")
else:
return fileutils.MAESTRO
def _autoLigFormat(self, filenames):
formats = {self._autoLigFormatOneFile(f) for f in filenames}
if len(formats) != 1:
raise ValueError("All LIGANDFILES must have the same format; "
"got " + ', '.join(sorted(formats)))
return formats.pop()
def _defaultsTable(self):
return [
# (key, val, [
# (set_key, set_val, override), ...])
('PRECISION', 'Normal', [
('PRECISION', 'SP', True),
]),
('PRECISION', 'Accurate', [
('PRECISION', 'XP', True),
]),
('DOCKING_METHOD', 'refineinput', [
('DOCKING_METHOD', 'mininplace', True),
]),
('DOCKING_METHOD', 'optandscore', [
('DOCKING_METHOD', 'mininplace', True),
]),
('PRECISION', 'HTVS',
[('CORESCALE', 0.05, True), ('MAXSOFT', 50, True),
('MAXHARD', 150, True), ('SAMPLING', 9, True),
('MAXKEEP', mm.GLIDE_MAXKEEP_HTVS_DEFAULT, True),
('MAXREF', mm.GLIDE_MAXREF_HTVS_DEFAULT, True),
('POSTDOCK_ITMAX', 20, True), ('CANONICALIZE', False, False),
('GLIDE_CONS_FINALONLY', True, True),
('FORCEFIELD', OPLS_2005, False)]),
('PRECISION', 'SP', [
('MAXKEEP', mm.GLIDE_MAXKEEP_SP_DEFAULT, False),
('MAXREF', mm.GLIDE_MAXREF_SP_DEFAULT, False),
]),
('PRECISION', 'XP', [
('MAXKEEP', mm.GLIDE_MAXKEEP_XP_DEFAULT, False),
('MAXREF', mm.GLIDE_MAXREF_XP_DEFAULT, False),
('POSTDOCK_NPOSE', mm.GLIDE_POSTDOCK_NPOSES_XP_DEFAULT, False),
]),
('WSCORE', True, [
('MAXKEEP', mm.GLIDE_MAXKEEP_WSCORE_DEFAULT, False),
('MAXREF', mm.GLIDE_MAXREF_WSCORE_DEFAULT, False),
('POSE_DISPLACEMENT', mm.GLIDE_POSE_DISPLACEMENT_WSCORE_DEFAULT,
False),
('POSE_RMSD', mm.GLIDE_POSE_RMSD_WSCORE_DEFAULT, False),
('EXPANDED_SAMPLING', True, False),
('POSES_PER_LIG', mm.GLIDE_POSES_PER_LIG_WSCORE_DEFAULT, True),
('POSTDOCK', False, True),
('PRECISION', 'XP', True),
('CANONICALIZE', False, False),
('WRITE_XP_DESC', True, False),
]),
('MACROCYCLE', True, [
('MAXKEEP', mm.GLIDE_MAXKEEP_MACROCYCLE_DEFAULT, False),
('MAXREF', mm.GLIDE_MAXREF_MACROCYCLE_DEFAULT, False),
('POSTDOCK_NPOSE', mm.GLIDE_POSTDOCK_NPOSES_MACROCYCLE_DEFAULT,
False),
('NPCONF', mm.GLIDE_MACROCYCLE_NPCONF, False),
('MAXCORE', mm.GLIDE_MACROCYCLE_MAXCORE, False),
('BALANCED_STRUCS', True, False),
('RINGCONFCUT', mm.GLIDE_MACROCYCLE_RINGCONFCUT, False),
]),
('PEPTIDE', self._checkPeptide, [
('MAXKEEP', mm.GLIDE_MAXKEEP_PEPTIDE_DEFAULT, False),
('MAXREF', mm.GLIDE_MAXREF_PEPTIDE_DEFAULT, False),
('POSTDOCK_NPOSE', mm.GLIDE_POSTDOCK_NPOSES_PEPTIDE_DEFAULT,
False),
('NDIR', mm.GLIDE_PEPTIDE_NDIR, True),
('NPCONF', mm.GLIDE_PEPTIDE_NPCONF, True),
]),
('DOCKING_METHOD', lambda v: v != 'confgen', [
('MAXKEEP', mm.GLIDE_MAXKEEP_NONCONFGEN_DEFAULT, False),
('MAXREF', mm.GLIDE_MAXREF_NONCONFGEN_DEFAULT, False),
('CANONICALIZE', False, True),
]),
('DOCKING_METHOD', 'inplace', [
('POSTDOCK', False, True),
('NENHANCED_SAMPLING', 1, True),
]),
('DOCKING_METHOD', 'mininplace', [
('GLIDE_CONS_FINALONLY', True, True),
]),
('CORE_POS_MAX_RMSD',
lambda v: v < mm.GLIDE_AUTO_CORE_SNAP_RMSD_CUTOFF, [
('CORE_SNAP', True, False),
]),
('CORE_SNAP', False, [
('CORE_POS_MAX_RMSD',
mm.GLIDE_CORE_POS_MAX_RMSD_NO_SNAP_DEFAULT, False),
]),
('CORE_DEFINITION', 'atomlist', [
('CORE_POS_MAX_RMSD', True, False),
]),
('CORE_DEFINITION', 'mcssmarts', [
('MAXATOMS', canvas.MAX_AGTYPE, False),
]),
('MMIMJOBTYPE', mm.MMIM_CHDOCK_JOB,
[('MAXREF', 5000, True), ('METAL_CUTOFF', 0.0, True),
('MAXATOMS', 300, True), ('MAXROTBONDS', 50, True),
('NREPORT', 10000, True), ('POSE_OUTTYPE', 'ligandlib', False),
('CANONICALIZE', False, False), ('EPIK_PENALTIES', False, False),
('FORCEFIELD', OPLS_2005, False)]),
('MMIMJOBTYPE', mm.MMIM_PCALIGN_JOB,
[('MAXREF', 5000, True), ('METAL_CUTOFF', 0.0, True),
('MAXATOMS', 300, True), ('MAXROTBONDS', 50, True),
('NREPORT', 10000, True), ('POSES_PER_LIG', 400, True),
('POSTDOCK', False, True), ('POSE_OUTTYPE', 'ligandlib', False),
('CANONICALIZE', False, False), ('EPIK_PENALTIES', False, False),
('CHDOCK_FROMTO', False, True),
('FORCEFIELD', OPLS_2005, False)]),
('CG_MODE', bool, [('POSE_OUTTYPE', 'ligandlib', False),
('CANONICALIZE', False, False),
('NREPORT', 150000, True),
('KEEPRAW', True, False),
('COMPRESS_POSES', False, True),
('FORCEFIELD', OPLS_2005, False)]),
('CG_MODE', 'sidechain', [
('POSES_PER_LIG', 10, False),
('NOSORT', True, False),
]),
('USECOMPMAE', False, [
('USECOMPMAE', True, False),
]),
('POSE_OUTTYPE', 'phase_subset', [
('COMPRESS_POSES', False, True),
]),
('REF_LIGAND_FILE', bool, [
('USE_REF_LIGAND', True, False),
]),
('CLIENT_MODULE', bool, [
('NOSORT', True, False),
]),
]
def _conflictsTable(self):
# Each line describes a potential conflict. If keyword1 has value1
# and keyword2 has value2, there is a conflict. A "value" can also
# be a callable; in that case we use the result of calling it with
# the actual value of the keyword.
return [
# (keyword1, value1, keyword2, value2)
('PRECISION', lambda v: v != 'XP', 'WRITE_XP_DESC', True),
('GLIDE_TORCONSFILE', bool, 'DOCKING_METHOD', 'rigid'),
('PRECISION', 'XP', 'SCORE_INPUT_POSE', True),
('PRECISION', 'XP', 'SCORE_MINIMIZED_INPUT_POSE', True),
('PRECISION', 'HTVS', 'PEPTIDE', True),
('PRECISION', 'HTVS', 'DOCKING_METHOD', 'mininplace'),
('WSCORE', True, 'DOCKING_METHOD', 'mininplace'),
('WSCORE', True, 'DOCKING_METHOD', 'rigid'),
('WSCORE', True, 'WS_WMAP_FILE', is_false),
('WSCORE', True, 'PRECISION', lambda v: v != 'XP'),
('HBOND_ACCEP_HALO', True, 'HBOND_DONOR_HALO', True),
('WRITE_XP_DESC', True, 'POSE_OUTTYPE',
lambda v: v.endswith('_sd')),
('GLIDE_ELEMENTS', True, 'PRECISION', 'XP'),
('GLIDE_ELEMENTS', True, 'WSCORE', True),
('GLIDE_ELEMENTS', True, 'CG_MODE', bool),
('GLIDE_ELEMENTS', True, 'MMIMJOBTYPE', mm.MMIM_PCALIGN_JOB),
('GLIDE_ELEMENTS', True, 'MMIMJOBTYPE', mm.MMIM_CHDOCK_JOB),
('GLIDE_ELEMENTS', True, 'FITDEN', True),
('LIGANDFILE', bool, 'LIGANDFILES', bool),
('CORE_RESTRAIN', True, 'USE_REF_LIGAND', False),
('CORE_RESTRAIN', True, 'REF_LIGAND_FILE', is_false),
('DOCKING_METHOD', is_equivalent_to_mininplace, 'CORE_RESTRAIN',
True),
('CORE_DEFINITION', 'mcssmarts', 'MAXATOMS',
lambda v: v > canvas.MAX_AGTYPE),
('LIGPREP', True, 'LIGFORMAT', fileutils.MOL2),
('LIGPREP', True, 'PHASE_DB', bool),
]
[docs] def checkConflicts(self):
super(Dock, self).checkConflicts()
if self._checkFlexOH() and (self['POSE_OUTTYPE'].endswith('_sd') or
self['POSE_OUTTYPE'] == PHASE_SUBSET):
raise ValueError(
"ERROR: Rotatable OH groups in %s are incompatible with "
"POSE_OUTTYPE %s.\n" % (self['GRIDFILE'], self['POSE_OUTTYPE']))
# We'll put the grid incompatibility checks here instead of in
# _conflictsTable to allow for more useful error messages.
props = self._getGridProperties()
grid_type = props.get('s_glide_grid_type')
if self['WSCORE'] and grid_type == 'glide':
raise ValueError("Can't use Glide grids for WScore jobs")
if not self['WSCORE'] and grid_type == 'wscore':
raise ValueError("Can't use WScore grids for Glide jobs")
if self['SHAPE_RESTRAIN'] and not (self['SHAPE_REF_LIGAND_FILE'] or
self['REF_LIGAND_FILE']):
raise ValueError("SHAPE_RESTRAIN requires a reference ligand"
" (SHAPE_REF_LIGAND_FILE or REF_LIGAND_FILE)")
if (self['SHAPE_REFINDEX'] != 1 and not self['SHAPE_REF_LIGAND_FILE']):
raise ValueError("SHAPE_REFINDEX requires SHAPE_REF_LIGAND_FILE")
self._checkConsConflicts()
def _checkConsConflicts(self):
"""
Check that the same constraint does not appear in more than one
constraint group.
"""
groups, _ = self.readConsKeywords()
seen = set()
dups = set()
for group in groups:
names = {cons.rsplit(':', 1)[0] for cons in group.constraints}
dups |= seen & names
seen |= names
if dups:
raise ValueError('Some constraints appear in more than one '
f'group: {sorted(dups)}')
def _setupCombiGlide(self, jsdict):
"""
Copy the values from certain CG_* keywords to the equivalent
"pure Glide" keyword.
"""
cgmap = {
"CG_RINGFLIP": "SAMPLE_RINGS",
"CG_INPUTRING": "INCLUDE_INPUT_RINGS",
"CG_NINVERT": "SAMPLE_N_INVERSIONS",
"CG_INPUTCONF": "INCLUDE_INPUT_CONF",
"CG_POSTDOCK": "POSTDOCK",
"CG_POSTDOCKSTRAIN": "POSTDOCKSTRAIN",
"CG_MAXATOM": "MAXATOMS",
"CG_MAXROTBOND": "MAXROTBONDS",
"CG_COREPOS": "CG_COREPOS",
"CG_POSTDOCKNPOSE": "POSTDOCK_NPOSE",
"CG_LIGVDWSCALE": "LIG_VSCALE",
"CG_LIGVDWCUTOFF": "LIG_CCUT",
"CG_AMIDE": "AMIDE_MODE",
}
for cg_key, g_key in cgmap.items():
jsdict[g_key] = self[cg_key]
[docs] def applyContextDependentDefaults(self):
# Special cases of context-dependent defaults that can't be
# represented by _defaultsTable.
# PRECISION needs to be figured out before calling parent method.
if self['CG_MODE']:
if self['CG_MODE'] == 'enumerate':
if self['CG_ENUM'] == 'combi':
self['PRECISION'] = 'XP'
else:
self['PRECISION'] = self['CG_ENUM'].upper()
self.defaults.append('PRECISION')
elif self['CG_MODE'] == 'molecules' and self['PRECISION'] == 'SP':
self['POSTDOCK'] = False
self.defaults.append('POSTDOCK')
else:
self['PRECISION'] = 'XP'
self.defaults.append('PRECISION')
super(Dock, self).applyContextDependentDefaults()
if self['CORE_SNAP'] and self['CORE_RESTRAIN'] \
and self['CORE_DEFINITION'] == 'all':
self['CORE_DEFINITION'] = 'allheavy'
if self['PRECISION'] == 'XP' \
and self['DOCKING_METHOD'] == 'inplace' \
and 'POSTDOCK' in self.defaults_set:
self['POSTDOCK'] = False
self.defaults.append('CORE_DEFINITION')
if self['PRECISION'] == 'XP' \
and self['DOCKING_METHOD'] == 'mininplace':
self['AMIDE_MODE'] = 'fixed'
self.defaults.append('AMIDE_MODE')
if 'GLIDE_REFLIG_FORMAT' in self.defaults and self['REF_LIGAND_FILE']:
self['GLIDE_REFLIG_FORMAT'] = self._autoLigFormatOneFile(
self['REF_LIGAND_FILE'])
self.defaults.append('GLIDE_REFLIG_FORMAT')
# backend only really knows about LIGANDFILE, so treat
# PCALIGN_PCFILE and CHDOCK_DRESSEDFILE as synonyms.
if self['PCALIGN_PCFILE']:
self['LIGANDFILE'] = self['PCALIGN_PCFILE']
if self['CHDOCK_DRESSEDFILE']:
self['LIGANDFILE'] = self['CHDOCK_DRESSEDFILE']
if 'LIGFORMAT' in self.defaults:
self['LIGFORMAT'] = self._autoLigFormat(self['LIGANDFILES'] or
[self['LIGANDFILE']])
self.defaults.append('LIGFORMAT')
# Glide Elements only supports one-ligand jobs.
if self['GLIDE_ELEMENTS']:
self['LIGAND_END'] = self['LIGAND_START']
# Docking SMILES requires LIGPREP.
if self['LIGFORMAT'] in (fileutils.SMILES, fileutils.SMILESCSV):
self['LIGPREP'] = True
# On-the-fly LigPrep support
if self['LIGPREP']:
options = {
'submodule': self['CLIENT_MODULE'],
'suboptions': self['CLIENT_OPTIONS']
}
if self['CLIENT_MODULE']:
self['NOSORT'] = True
else:
# Undocumented keyword that tells the backend to create a raw
# file despite being in client mode.
self['_WANT_RAW'] = True
self['CLIENT_MODULE'] = \
'schrodinger.application.glide.packages.ligprep'
self['CLIENT_OPTIONS'] = json.dumps(options)
def _setupPhase(self, jsdict):
"""
Make sure that the Phase keywords make sense before writing a JSON file.
"""
phase_db = self.get('PHASE_DB')
if phase_db:
if not self.get('PHASE_SUBSET'):
jsdict['PHASE_SUBSET'] = os.path.join(
phase_db, 'database_master_phase.inp')
if not os.path.isabs(phase_db):
raise RuntimeError('ERROR: PHASE_DB must be an absolute path.')
[docs] def setupJson(self, jsdict):
mm.mmim_initialize(mm.MMERR_DEFAULT_HANDLER)
self.mmim_handle = mm.mmim_new()
self.mmim_dict = mmim.MMIMDict(self.mmim_handle)
self._setupJsonTorcons(jsdict)
jsdict['GRIDFILE'] = self._gridfile()
self.mmim_dict['GRIDFILE'] = jsdict['GRIDFILE']
self._setupConstraints(write_featfile=True)
if self._cons_groups:
jsdict['GLIDECONS'] = True
jsdict['GLIDEUSECONSFEAT'] = True
jsdict['GLIDE_CONS_FEAT_FILE'] = self._featfile()
if self['CG_MODE']:
self._setupCombiGlide(jsdict)
self._setupPhase(jsdict)
self.cleanup()
mm.mmim_delete(self.mmim_handle)
mm.mmim_terminate()
self.mmim_handle = None
self.mmim_dict = None
def _setupJsonTorcons(self, jsdict):
if not any(k.startswith('TORCONS') for k in self.keys()):
return # no torsional constraints blocks
self._setupTorCons()
jsdict["GLIDE_TORCONS_PATTERNS"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_PATTERNS])
jsdict["GLIDE_TORCONS_IATOMS"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_IATOMS])
jsdict["GLIDE_TORCONS_JATOMS"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_JATOMS])
jsdict["GLIDE_TORCONS_KATOMS"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_KATOMS])
jsdict["GLIDE_TORCONS_LATOMS"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_LATOMS])
jsdict["GLIDE_TORCONS_PATTERN_INDEX"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_PATTERN_INDEX])
jsdict["GLIDE_TORCONS_SETVAL"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_SETVAL])
jsdict["GLIDE_TORCONS_VALUES"] = list(
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_VALUES])
torconsfile = "%s.torcons" % self.jobname
jsdict["GLIDE_TORCONSFILE"] = torconsfile
# This mmim call writes the torcons file
mm.mmim_handle_glide_get_torcons(self.mmim_handle, torconsfile)
def _gridfile(self):
"""
Return the grid file to use, normally just the value of the GRIDFILE
keyword, but sometimes based on legacy keywords for compatibility.
"""
if self.get('GRIDFILE'):
return self.get('GRIDFILE')
# Support for old READBASE/READZIP keywords:
gridbase = self.get('READBASE')
readzip = self.get('READZIP', False)
if gridbase:
import warnings
msg = 'glide.py: keyword READBASE is obsolete. " \
"Use GRIDFILE instead.'
warnings.warn(msg, DeprecationWarning, stacklevel=2)
# FIXME private function call. This shouldn't be needed, as
# ideally section['READZIP'] should always return a bool type
if mmkv._convert_to_bool(readzip):
return gridbase + '.zip'
else:
return gridbase + '.grd'
else: # Neither GRIDFILE nor READBASE are specified:
raise RuntimeError("ERROR: grid file was not specified")
def _featfile(self):
return '%s.feat' % self.jobname
[docs] def readConsKeywords(self):
"""
Reads the constraint keywords, parses their values, and returns them
as a tuple of (cons groups, feature sets).
"""
cons_groups = []
feature_sets = {} # key: feature name, value: list of patterns
for key in self.keys(): # Go through all keywords & sections
if key.startswith('CONSTRAINT_GROUP'):
section = self[key]
use_cons = section.get('USE_CONS')
nrequired_str = section.get('NREQUIRED_CONS')
if use_cons is None or nrequired_str is None:
msg = "ERROR: Must define NREQUIRED_CONS and USE_CONS" \
" in section {SECTION}".format(SECTION=key)
raise RuntimeError(msg)
if type(use_cons) == type(''):
use_cons = use_cons.split(' ') # Convert to list
unknown_keywords = set(map(
str.upper, section)) - consgroup_section_allowed_keywords
for bad_key in unknown_keywords:
print('Invalid keyword in section "%s": %s.' %
(key, bad_key),
file=sys.stderr)
if bad_key in self.configspec:
print(
"(%s would have been a valid keyword at the top of the file, before the first subsection.)"
% bad_key,
file=sys.stderr)
if unknown_keywords:
raise RuntimeError('Invalid keywords in section "%s": %s' %
(key, ', '.join(unknown_keywords)))
if nrequired_str == 'ALL':
nrequired = ALL
else:
try:
nrequired = int(nrequired_str)
except:
msg = 'ERROR: Invalid NREQUIRED_CONS value: %s' % nrequired_str
raise RuntimeError(msg)
consgroup = ConstraintGroup(self, use_cons, nrequired)
cons_groups.append(consgroup)
elif key.startswith('FEATURE'):
featurename = key.split(':')[1]
section = self[key]
badfeatkey = False
have_asl = False
have_smarts = False
unknown_keywords = []
for patternkey in list(section):
if patternkey.startswith('PATTERN'):
have_smarts = True
elif patternkey.startswith('ASL'):
have_asl = True
else:
badfeatkey = True
unknown_keywords.append(patternkey)
print('Invalid keyword in section "%s": %s.' %
(key, patternkey),
file=sys.stderr)
if patternkey in self.configspec:
print(
"(%s would have been a valid keyword at the top of the file, before the first subsection.)"
% patternkey,
file=sys.stderr)
if have_smarts and have_asl:
raise RuntimeError("Error: '%s' has both ASL expressions "
"and SMARTS patterns" % key)
if badfeatkey:
raise RuntimeError('Invalid keywords in section "%s": %s' %
(key, ', '.join(unknown_keywords)))
for key, value in section.items():
if key.startswith('ASL'):
# we don't do splitting on spaces here because ASL
# expressions may contain spaces and it would be tricky
# to figure out where to split. We have to rely on
# the input file being quoted correctly using
# InputConfig syntax.
# No atom spec needed for ASL; insert empty atom list
s = [value[0], ''] + value[1:]
atoms = []
else:
# SMARTS-based feature
if type(value) == type([]):
# Convert ConfigObj generated list to string:
value = ' '.join(value)
s = value.split()
atoms = [int(a) for a in s[1].split(',')]
pat = s[0]
if len(s) > 2:
exclude = s[2]
if exclude == 'exclude':
exclude = True
elif exclude == 'include':
exclude = False
else:
exclude = bool(exclude)
else: # Exclude flag is not specified
exclude = False
if len(s) > 3:
nfill = s[3]
else: # Nfill arg is not specified
nfill = None # default
pattern = FeaturePattern(pat, atoms, exclude, nfill)
try:
feature_sets[featurename].append(pattern)
except KeyError:
feature_sets[featurename] = [pattern]
return cons_groups, feature_sets
def _setupConstraints(self, write_featfile=False):
# Calls appropriate MMIM APIs to setup constraints
self._cons_groups, self._feature_sets = self.readConsKeywords()
if self._cons_groups:
# Read receptor constraints from *.cons file:
try:
cons_labels = mm.mmim_handle_glide_read_consfile(
self.mmim_handle)
except mm.MmException as e:
print(str(e))
raise RuntimeError("Failed to read the *.cons file (MMIM)")
for group_obj in self._cons_groups:
group_obj.setup(cons_labels, self._feature_sets)
if write_featfile:
mm.mmim_handle_glide_cons_write_featfile(
self.mmim_handle, self._featfile())
[docs] def cleanup(self):
""" Free up the group handles """
for group_obj in self._cons_groups:
group_obj.cleanup()
self._cons_groups = []
[docs] def addConstraintGroup(self, cons_labels, nrequired):
"""
Add a new constraint group to this Docking instance.
:param cons_labels:
A list of labels for constraints to use.
:param nrequired:
Minimum number of constraints required for match (Can be ALL).
"""
for i in range(1, 100):
groupname = 'CONSTRAINT_GROUP:%i' % i
if not self.get(groupname): # Safe to add this group
self[groupname] = {
'USE_CONS': cons_labels,
'NREQUIRED_CONS': nrequired
}
[docs] def useConstraint(self, label):
"""
Adds specified constraint to the first group.
Default group is #1; constraint type determined automatically.
"""
groupname = 'CONSTRAINT_GROUP:1'
if groupname in self:
self[groupname]['USE_CONS'].append(label)
else:
group = {'USE_CONS': [label], 'NREQUIRED_CONS': ALL}
self[groupname] = group
[docs] def setNRequired(self, nrequired):
"""
Set the minimum number of constraints required for the first group.
A value of ALL means all constraints are required.
"""
groupname = 'CONSTRAINT_GROUP:1'
if groupname not in self:
raise SyntaxError(
"glide.py: setNRequired() called before useConstraint()")
self[groupname]['NREQUIRED_CONS'] = nrequired
def _setupTorCons(self):
# Parses [TORCONS:<n>] sections into mmim array variables, in
# numerical order. The index <n> must range from 000 to 999 in
# order for the sort to give the numerical order we want. (Or
# 00 to 99, or 0000 to 9999, or whatever the number of digits.
# The point is, every <n> must be zero-padded on the left to
# make it the same width as the highest possible <n>.)
self._torcons_keys = []
smarts = []
patindex = []
iatom = []
jatom = []
katom = []
latom = []
flag = []
val = []
for key in self.keys():
if key.startswith('TORCONS'):
self._torcons_keys.append(key)
self._torcons_keys.sort()
npat = 0
ntor = 0
nflag = 0
nval = 0
for key in self._torcons_keys:
section = self[key]
oldntor = ntor
smarts.append(section['SMARTS'])
allbonds = False
# FIXME private function call. This shouldn't be needed, as ideally
# section['ALL'] should always return a bool type
if 'ALL' in section and mmkv._convert_to_bool(section['ALL']):
allbonds = True
else:
atomsets = section['ATOMS']
if type(atomsets) == type(''):
atomsets = atomsets.split(',') # Convert to list
useflags = section['USEVAL']
if type(useflags) == type(''):
useflags = useflags.split(',') # Convert to list
if len(useflags) != len(atomsets):
raise ValueError(
"For pattern %d (%s), n USEVAL flags (%d) != n torsion specs (%d)"
%
(npat + 1, smarts[npat], len(useflags), len(atomsets)))
torvals = section['TORVAL']
if type(torvals) == type(''):
torvals = torvals.split(',') # Convert to list
if len(torvals) != len(atomsets):
raise ValueError(
"For pattern %d (%s), n values (%d) != n torsion specs (%d)"
% (npat + 1, smarts[npat], len(torvals), len(atomsets)))
unknown_keywords = set(map(
str.upper, section)) - torcons_section_allowed_keywords
for bad_key in unknown_keywords:
print('Invalid keyword in section "%s": %s.' % (key, bad_key),
file=sys.stderr)
if bad_key in self.configspec:
print(
"(%s would have been a valid keyword at the top of the file, before the first subsection.)"
% bad_key,
file=sys.stderr)
if unknown_keywords:
raise RuntimeError('Invalid keywords in section "%s": %s' %
(key, ', '.join(unknown_keywords)))
(ct, st) = mm.mmpatty_get_ct_from_pattern(smarts[npat])
if allbonds:
# Iterate through all possible distinct rotatable bonds in the pattern
# ((j,k); don't also take (k,j) *or* any bond specified for a previous
# pattern), selecting one proper dihedral (i bonded to j, l bonded to k)
# per bond.
struc = structure.Structure(ct)
for jstr in struc.atom:
j = jstr.index
for kstr in jstr.bonded_atoms:
k = kstr.index
if (k <= j):
continue
# EV:97482. Skip bonds in rings.
curlevel = mm.mmerr_get_level(mm.MMERR_DEFAULT_HANDLER)
mm.mmerr_level(mm.MMERR_DEFAULT_HANDLER, mm.MMERR_OFF)
isring = mm.mmim_is_ringbond(ct, j, k)
mm.mmerr_level(mm.MMERR_DEFAULT_HANDLER, curlevel)
if (isring):
continue
isdupl = False
for itor in range(oldntor, ntor):
if (j == jatom[itor] and
k == katom[itor]) or (j == katom[itor] and
k == jatom[itor]):
isdupl = True
break
if (isdupl):
continue
foundijkl = False
for istr in jstr.bonded_atoms:
i = istr.index
if (i == k):
continue
for lstr in kstr.bonded_atoms:
l = lstr.index # noqa: E741, bond atom
if (l == i or l == j):
continue
curlevel = mm.mmerr_get_level(
mm.MMERR_DEFAULT_HANDLER)
mm.mmerr_level(mm.MMERR_DEFAULT_HANDLER,
mm.MMERR_OFF)
isdihe = mm.mmim_is_dihedral(ct, i, j, k, l)
mm.mmerr_level(mm.MMERR_DEFAULT_HANDLER,
curlevel)
if (isdihe):
patindex.append(npat)
iatom.append(i)
jatom.append(j)
katom.append(k)
latom.append(l)
ntor += 1
flag.append(0)
nflag += 1
val.append(0.0)
nval += 1
foundijkl = True
break
if (foundijkl
): # only need one i and one l per (j,k)
break
else: # Constrain specified dihedrals only.
curitor = 1
for quartet in atomsets:
strlist = quartet.split('-')
[i, j, k, l] = [int(istr) for istr in strlist]
isdihe = mm.mmim_is_dihedral(ct, i, j, k, l)
if (not isdihe):
raise ValueError(
"Atom indices (%d, %d, %d, %d) do not denote a dihedral in CT from pattern %d (%s)."
% (i, j, k, l, npat + 1, smarts[npat]))
for itor in range(oldntor, ntor):
if (j == jatom[itor] and
k == katom[itor]) or (k == jatom[itor] and
j == katom[itor]):
raise ValueError(
"For pattern %d (%s), bond %d (%d,%d) duplicates bond %d (%d, %d)"
% (npat + 1, smarts[npat], curitor, j, k,
itor + 1, jatom[itor], katom[itor]))
patindex.append(npat)
iatom.append(i)
jatom.append(j)
katom.append(k)
latom.append(l)
ntor += 1
curitor += 1
for curflag in useflags:
flag.append(int(curflag))
nflag += 1
for curval in torvals:
val.append(float(curval))
nval += 1
npat += 1
if ((ct > -1) and (mm.mmct_ct_in_use(ct)) and not allbonds):
mm.mmct_ct_delete(ct)
if (st > -1):
mm.mmstereo_delete(st)
if nflag != ntor:
raise ValueError(
"Total n USEVAL flags (%d) != total n torsion specs (%d)" %
(nflag, ntor))
if nval != ntor:
raise ValueError(
"Total n values (%d) != total n torsion specs (%d)" %
(nval, ntor))
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_PATTERNS] = smarts
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_IATOMS] = iatom
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_JATOMS] = jatom
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_KATOMS] = katom
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_LATOMS] = latom
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_PATTERN_INDEX] = patindex
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_SETVAL] = flag
self.mmim_dict[mm.MMIM_GLIDE_TORCONS_VALUES] = val
#############################################################################
# FUNCTIONS
#############################################################################
[docs]def read_constraints(gridfile):
"""
Reads the receptor constraints from the specified grid file.
The constraints are actually read from the .cons file. The resulting
constraint "type" refers to the atom type in the RECEPTOR, where
applicable. Specifically, H-bond and metal-binding constraints.
The matching ligand atom for an H-bond donor receptor constraint should
be an H-bond acceptor (by default). Likewire the matching ligand ato for
and H-bond acceptor receptor constraint should be a donor (polar H atom).
Other constraint types don't refer to any specific atom type
(Positional, NOE), or require the same atom type in the ligand
as in the receptor (Hydrophobic) so the constraint type is
valid for both molecules (ligand and receptor).
To get the constraint label, simply do str(Constraint).
:param gridfile: The gridfile name; should be either a .grd or .zip file.
:rtype: A list of Constraint objects.
"""
(gridname, ext) = os.path.splitext(gridfile)
(griddir, gridbase) = os.path.split(gridname)
if ext not in ['.grd', '.zip']:
msg = 'glide.py: read_constraints(): invalid gridfile extension:'
msg += gridfile
raise SyntaxError(msg)
# Initialize MMIM:
mm.mmim_initialize(mm.error_handler)
# Create MMIM handle
mmim_handle = mm.mmim_new()
mmim_dict = mmim.MMIMDict(mmim_handle)
# Specify grid file to MMIM
mmim_dict['gridfile'] = gridfile
mmim_dict['jobname'] = gridbase
if DEBUG:
print('Calling mmim_handle_glide_read_consfile(handle)')
# Read the *.cons file. Will return empty list if no *.cons file present:
try:
cons_labels = mm.mmim_handle_glide_read_consfile(mmim_handle)
except mm.MmException as e:
print(str(e))
raise RuntimeError("Failed to read the *.cons file (MMIM)")
if DEBUG:
print(" after mmim_handle_glide_read_consfile() labels:", cons_labels)
constraints = []
for label in cons_labels:
# Determine the constraint type:
type_int = mm.mmim_handle_glide_cons_type_from_label(mmim_handle, label)
# Map constraint type strings to the MMIM constant value:
constype = mmim_cons_types[type_int]
if constype == HBOND_DONOR:
atoms = mm.mmim_handle_glide_cons_get_hbond(mmim_handle, label)
cons = DonorConstraint(label, atoms)
elif constype == HBOND_ACCEPTOR:
atoms = mm.mmim_handle_glide_cons_get_hbond(mmim_handle, label)
cons = AcceptorConstraint(label, atoms)
elif constype == METAL:
atoms = mm.mmim_handle_glide_cons_get_metal(mmim_handle, label)
cons = MetalConstraint(label, atoms)
elif constype == POSITIONAL:
x, y, z, radius = mm.mmim_handle_glide_cons_get_posit(
mmim_handle, label)
cons = PositionalConstraint(label, x, y, z, radius)
elif constype == NOE:
x, y, z, rmin, rmax = mm.mmim_handle_glide_cons_get_noe(
mmim_handle, label)
cons = NOEConstraint(label, x, y, z, rmin, rmax)
elif constype == HYDROPHOBIC:
cons = HydrophobicConstraint(label)
elif constype == METCOORD:
xmetc, ymetc, zmetc, rmetc = mm.mmim_handle_glide_cons_get_metcoord(
mmim_handle, label)
cons = MetCoordConstraint(label, xmetc, ymetc, zmetc, rmetc)
else:
raise ValueError("Invalid constraint type: %i" % type_int)
constraints.append(cons)
# Delete the MMIM handle:
mm.mmim_delete(mmim_handle)
# Terminate MMIM:
mm.mmim_terminate()
return constraints
[docs]def read_yaml():
"""
Read mmim_keywords.json and return a list of objects. The read_yaml()
name was retained for compatibility.
"""
return read_keywords_file()
_keywords_data = None
_keywords_dict = None
[docs]def read_keywords_file():
"""
Read the mmim_keywords.json file from the mmshare data directory
and return a list of objects.
"""
# The new function name is format-agnostic, just in case we change file
# formats again. :-)
global _keywords_data
if _keywords_data is None:
mmshare_data = util.hunt('mmshare', dir='data')
fname = os.path.join(mmshare_data, 'mmim_keywords.json')
with open(fname) as fh:
_keywords_data = json.load(fh)
return _keywords_data
[docs]def read_keywords_file_as_dict():
"""
Read the mmim_keywords.json file from the mmshare data directory
and return a dict of dicts, where the key is the keyword name and
the value is the spec for the keyword (see get_keyword for more details).
"""
global _keywords_dict
if _keywords_dict is None:
_keywords_dict = {k['keyword'].upper(): k for k in read_keywords_file()}
return _keywords_dict
[docs]def get_keyword(keyword):
"""
Return the specification of a Glide keyword as a dictionary, with keys such
as "keyword", "type", "default", "allowed", "min", "max", and "comment". For
the full list and description for each field, see the comments in
mmim_keywords.yaml.
"""
return read_keywords_file_as_dict()[keyword]
mmim_to_inputconfig_type_map = {
'boolean': 'boolean',
'boolean_array': 'bool_list',
'integer': 'integer',
'integer_array': 'int_list',
'float': 'float',
'float_array': 'float_list',
'string': 'string',
'string_array': 'string_list',
}
[docs]def get_keywords_dicts(jobtypes, include_undocumented=False):
"""
Return the raw mmim_keywords.yaml data filtered to include only
the keywords for the requested job types.
"""
jobtypes_set = set(jobtypes) | set(['all'])
keywords = [
k for k in read_keywords_file() if set(k['jobtypes']) &
jobtypes_set and (include_undocumented or not k.get('undocumented'))
]
return keywords
[docs]def keyword_spec(jobtypes, include_undocumented=False):
"""
Return a string with the specification of the mmim keywords that are valid
for the jobtypes given. Jobtypes should be an iterable with items such as
'docking', 'multi_docking', or 'gridgen'. Keywords for jobtype 'all' are
always included.
By default, undocumented keywords are not included in the return value, but
that behavior can be overridden by setting 'include_undocumented' to True.
"""
keywords_dicts = get_keywords_dicts(jobtypes, include_undocumented)
keywords_dicts.sort(key=lambda x: x['keyword'])
maxlen = max(len(k['keyword']) for k in keywords_dicts)
fh = io.StringIO()
for k in keywords_dicts:
default = k.get('default')
if k['type'] == 'string' and default is not None:
default = "'%s'" % default
elif isinstance(default, list):
default = "list(%s)" % (", ".join(str(d) for d in default))
comment = ""
if 'comment' in k and k['comment']:
comment = " # %s" % k['comment']
typename = mmim_to_inputconfig_type_map[k['type']]
if include_undocumented:
allowed = k.get('allowed')
else:
allowed = k.get('show_allowed') or k.get('allowed')
allowed_str = ""
if allowed:
quoted_allowed = ["'%s'" % value for value in allowed]
# trailing comma needed because allowed_str will always be followed
# by default value specification
allowed_str = ", ".join(quoted_allowed) + ", "
typename = 'option'
in_paren_str = allowed_str + "default=%s" % default
min_value = k.get('min')
max_value = k.get('max')
if min_value is not None:
in_paren_str += ', min=%s' % min_value
if max_value is not None:
in_paren_str += ', max=%s' % max_value
line = "%-*s = %s(%s)%s" % (maxlen, k['keyword'].upper(), typename,
in_paren_str, comment)
print(line, file=fh)
return fh.getvalue()
[docs]def validation_specs(jobtype):
"""
Return a list of validation specs for a given jobtype (e.g., 'docking' or
'gridgen') or list of jobtypes. This includes not only the keywords
returned by keyword_spec(jobtype), but also some extra keywords that are
supported by glide.py but not by mmim.
"""
if type(jobtype) == str:
jobtypes = [jobtype]
else:
jobtypes = jobtype
jobtype = jobtypes[0]
specs_str = keyword_spec(jobtypes, include_undocumented=True)
specs = specs_str.split('\n')
specs.pop() # rm traling blank due to split of EOL-terminated specs_str
return specs
[docs]def get_keywords_set(jobtypes):
"""Return the set of all keywords valid for the given job types."""
keywords_dicts = get_keywords_dicts(jobtypes, include_undocumented=True)
return {k['keyword'].upper() for k in keywords_dicts}
[docs]def is_gridgen_keyword(keyword):
"""Return True is a keyword is valid for gridgen jobs."""
return keyword in get_keywords_set({'gridgen'})
[docs]def is_docking_keyword(keyword):
"""Return True is a keyword is valid for docking jobs."""
return keyword in get_keywords_set({'docking'})
[docs]def calc_ligand_size(coords):
"""
Given a list of (x, y, z) coordinates for all atoms in a ligand,
calculate the ligand size. This value should be added to the
innerbox to get the default outerbox.
:type coords: List of 3-item tuples.
:param coords: List of coordinates for all ligand atoms.
:rtype: float
:return: The size of the ligand.
"""
ligand_coords = numpy.array(coords)
natoms = ligand_coords.shape[0]
# Get the length of the ligand by finding the longest distance
# between two atoms
dist_matrix = (numpy.tile(ligand_coords,
(1, natoms)) - numpy.tile(ligand_coords.flatten(),
(natoms, 1)))
dist_matrix = dist_matrix.reshape((natoms, natoms, 3))
ligand_length = numpy.sqrt(numpy.amax(numpy.sum(dist_matrix**2, 2)))
# The ligand size is 4 + 80% of the ligand length:
return 4.0 + 0.8 * ligand_length
[docs]def get_default_box_size(struct, ligand_asl='all'):
"""
Gets the default box size for a grid. This returns the outer and inner
box lengths.
:param struct: The structure containing the reference ligand
:type struct: structure.StructureReader
:param ligand_asl: The ASL for the ligand. Default: "all"
:type ligand_asl: str
:rtype: tuple of lists, e.g. ([10, 10, 10], [-5.9863, -1.887 , -0.2703])
"""
# Set the default inner box but overwrite it if it is in the
# mmim_keywords file.
inner_box = [10.0, 10.0, 10.0]
for spec in read_keywords_file():
if spec['keyword'] == 'innerbox' and spec['default']:
innerbox = spec['default']
# Get an array of the ligand coordinates (make sure to cast as numpy.array)
ligand_coords = []
for atom in analyze.get_atoms_from_asl(struct, ligand_asl):
ligand_coords.append(atom.xyz)
lig_size = calc_ligand_size(ligand_coords)
# Pad the inner box by adding the ligand size to each
# of the box's dimensions:
outer_box = numpy.array(inner_box) + lig_size
return (inner_box, outer_box.tolist())
[docs]def is_gridgen_job(keywords):
"""
Return True if the keywords suggest that this is a gridgen job.
:type keywords: iterable (typically a dict)
:rtype: bool
"""
# All gridgen jobs should have at least one of these keywords.
gridgen_keywords = {
'XCENT', 'GRID_CENTER', 'GRID_CENTER_ASL', 'LIGAND_MOLECULE',
'LIGAND_INDEX'
}
return bool(gridgen_keywords.intersection(keywords))
[docs]def get_glide_job(input_file):
"""
Return a Gridgen or Dock object given an SIF input filename or keyword
dictionary. If it is a filename, use it to set JOBNAME (e.g. for
"dir/myjob.in", set it to "myjob").
"""
keywords = InputConfig(input_file)
if is_gridgen_job(keywords):
jobtype = Gridgen
else:
jobtype = Dock
glide_job = jobtype(keywords)
if isinstance(input_file, str):
glide_job['JOBNAME'], ext = os.path.splitext(
os.path.basename(input_file))
return glide_job
[docs]def is_false(val):
return not val
[docs]def is_equivalent_to_mininplace(val):
""" Checks if val is equivalent to mininplace"""
return val in ["mininplace", "optandscore", "refineinput"]
# Used by GUIs:
MAX_ATOMS_SUPPORTED = get_keyword('MAXATOMS')['max']
#EOF