"""
Deprecated: `structutils.minimize` has been deprecated. Use
`schrodinger.forcefield.minimizer` instead.
A deprecated module for force field energy evaluation and minimization.
The `minimize.Minimizer` class in this module is stateful and can cause
problematic interactions with existing instances of mmffld.
`schrodinger.forcefield.minimizer` meets most needs, though is not a drop-in
replacement. If you need a feature that it lacks, please file a SHARED ticket.
"""
import contextlib
import warnings
from schrodinger import structure
from schrodinger.forcefield import OPLSVersion
from schrodinger.forcefield import common
from schrodinger.infra import mm
from schrodinger.infra import mmerr
# NOTE: Different constant names are used here for backwards compatability
OPLS_2005 = OPLSVersion.F14
OPLS3 = OPLSVersion.F16
MinBFGS = mm.MMFfldMinBFGS
MinAUTO = mm.MMFfldMinAUTO
MinCG = mm.MMFfldMinCG
#
# All output from the minimizer is stored as mmct properties with rather
# long names. The following functions provide for easier access to these
# properties given a Structure and the force field version used.
#
[docs]def minimize_structure(struct: structure.Structure, **kwargs):
"""
Deprecated: `structutils.minimize.minimize_structure` has been
deprecated. Use `schrodinger.forcefield.minimizer.minimize_structure`
instead.
"""
with minimizer_context(struct=struct, **kwargs) as minimizer:
return minimizer.minimize()
#
# Functions to create properties in the Minimizer class.
#
def _minimizer_int_property(mmffld_option, doc=None, readonly=True):
"""
A function to generate an int property for minimizer settings.
"""
def prop_get(self):
int_value, float_value, string_value = \
mm.mmffld_getOption(self.handle, mmffld_option)
return int_value
def prop_set(self, value):
if readonly and self._struct:
# If used by anything but the constructor (Ev:81067)
raise RuntimeError("This Minimizer property is read-only")
int_value = int(value)
float_value = 0.0
string_value = ""
mm.mmffld_setOption(self.handle, mmffld_option, int_value, float_value,
string_value)
return property(prop_get, prop_set, doc=doc)
def _minimizer_float_property(mmffld_option, doc=None, readonly=True):
"""
A function to generate an float property for minimizer settings.
"""
def prop_get(self):
int_value, float_value, string_value = \
mm.mmffld_getOption(self.handle, mmffld_option)
return float_value
def prop_set(self, value):
if readonly and self._struct:
# If used by anything but the constructor (Ev:81067)
raise RuntimeError("This Minimizer property is read-only")
int_value = 0
float_value = value
string_value = ""
mm.mmffld_setOption(self.handle, mmffld_option, int_value, float_value,
string_value)
if mmffld_option == mm.MMFfldOption_SYS_CUTOFF:
# PBC may cause non-bonded cutoff value to change, so back up
# user's initial value. Doing the back up ensures that we can
# re-calculate the cutoff for each structure.
self._users_nonbonded_cutoff = value
return property(prop_get, prop_set, doc=doc)
[docs]class Minimizer(object):
"""
deprecated: `structutils.minimize.Minimizer` has been deprecated.
Use `schrodinger.forcefield.minimizer` instead.
A class to provide force field energy and minimization.
For simple cases the minimize_structure function can be used. Direct use
of this class is mostly useful for cases where multiple conformers are
to be minimized, as the updateCoordinates() method allows one to avoid
running atom-typing for every conformer.
"""
[docs] def __init__(self,
ffld_version=None,
struct=None,
cleanup=True,
honor_pbc=True,
**kwargs):
"""
Initialize.
Additional keyword arguments will be used to set properties of the
instance. See the __init__() method for supported property names.
:type ffld_version: integer module constant
:param ffld_version: Use one of the valid force fields. Default is to
use the force field that is selected in Maestro preferences.
:type struct: schrodinger.structure.Structure
:param struct:
The structure to be minimized.
:type cleanup: bool
:param cleanup:
If True, attempts to automatically clean up the structure will
be made. (This uses the C function `mmlewis_apply()`.) Note
that this can modify the atom types of the passed in structure.
:type min_method: enum
:param min_method: The minimizer method. Valid values are MinBFGS,
MinAUTO and MinCG. Default is MinAUTO, which uses BFGS if number of
atoms is less than 500, CG otherwise.
:type max_steps: int
:param max_steps: The maximum number of steps to run the mimization
for. Default is 1500.
:type energy_criterion: float
:param energy_criterion: The energy convergence criterion.
Default is 5.0e-09.
:type gradient_criterion: float
:param gradient_criterion: The gradient convergence criterion.
Default is 0.05.
:type verbose: bool
:param verbose: Printing verbosity. Default is False.
:type energy_no_electrostatics: bool
:param energy_no_electrostatics: Whether to turn off electrostatics.
Default is False, i.e. to use electrostatics.
NOTE: Can not be modified after the instance is created
:type energy_lj_repulsive_only: bool
:param energy_lj_repulsive_only: Whether to use only the repulsive
portion of the Lennard-Jones term. Default is False.
NOTE: Can not be modified after the instance is created
:type nonbonded_cutoff: float
:param nonbonded_cutoff: The cutoff for non-bonded interactions. It
will be automatically scaled down for small BPC boxes.
Default is 14.0.
:type maintain_planarity: bool
:param maintain_planarity: Enable scaling of forces to artificially
maintain planarity of certain sub-structures. This sets associated
scale factors to their default values. Default is False.
NOTE: Can not be modified after the instance is created
:type dielectric_constant: float
:param dielectric_constant: The strength of the constant dielectric
field used in the energy calculation. The default is 1.0 (vacuum).
:type no_zob_restrain: bool
:param no_zob_restrain: The default (False) is to restrain zero-order
metal bonds to their input value. Setting this to True causes the
forcefield parameters to be taken from lookup results if possible.
:type charges_from_ct: bool
:param charges_from_ct: The default (False) is to use the
charges from the force field. Setting this to True causes the
charges from the ct to be used.
:param honor_pbc: Honor Periodic Boundary Conditions, if defined as
properties in the structure. Default is True.
:type honor_pbc: bool
"""
self._exit_stack = contextlib.ExitStack()
self._exit_stack.enter_context(common.mmffld_environment())
if ffld_version is None:
ffld_version = OPLS3 if mm.mmffld_license_ok(OPLS3) else OPLS_2005
self._ffld_version = ffld_version
# Default value, which will get over-written if the user specified
# nonbonded_cutoff option in kwargs:
self._users_nonbonded_cutoff = 14.0
self._honor_pbc = honor_pbc
self.handle = mm.mmffld_new(self._ffld_version)
self._struct = None # forces setStructure() to reset the CT
self.energy_no_electrostatics = False
self.cleanup = cleanup
# Set attributes of Minimizer based on specified keywords:
# Must do this before we enter the structure via setStructure()
self.setOptions(**kwargs)
# Will call mmffld_enterMol(), after which options can no longer
# be changed (Ev:81067):
self.setStructure(struct)
def __del__(self):
self.cleanup_resources()
[docs] def cleanup_resources(self):
# Run this code only if constructor of the class was successful:
if hasattr(self, 'handle') and self.handle is not None:
if hasattr(self, '_struct'):
if self._struct is not None:
mm.mmffld_deleteMol(self.handle, self._struct)
self._struct = None
mm.mmffld_delete(self.handle)
self.handle = None
self._exit_stack.close()
[docs] def getStructure(self):
"""
Return the current structure.
"""
return self._struct
[docs] def updateCoordinates(self, struct):
"""
Update the coordinates of the current structure with the values from
the provided struct.
This allows an additional conformer to be minimized without
re-running atomtyping. It is the caller's responsibility to make
sure the molecules are equivalent and have the same atomtypes (i.e.
same charges, connectivity).
"""
if self._struct != struct:
for ix in range(1, len(struct.atom) + 1):
self._struct.atom[ix].xyz = struct.atom[ix].xyz
[docs] def minimize(self):
"""
Minimize the provided Structure.
Note that this method will require a valid product license.
Currently, MacroModel, GLIDE, Impact, or PLOP will suffice.
:return: minimization results
:rtype: mm.MinimizeResults
:raise RuntimeError: If minimization fails
"""
return mm.mmffld_minimize(self.handle)
[docs] def getTotalEnergy(self):
"""
get the total energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldTOT_ENE:
return ene_list[mm.MMFfldTOT_ENE]
else:
msg = 'getTotalEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getStretchEnergy(self):
"""
get the stretch energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldSTR_ENE:
return ene_list[mm.MMFfldSTR_ENE]
else:
msg = 'getStretchEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getBendingEnergy(self):
"""
get the beinding energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldBND_ENE:
return ene_list[mm.MMFfldBND_ENE]
else:
msg = 'getBendingEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getTorsionEnergy(self):
"""
get the torsion energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldTOR_ENE:
return ene_list[mm.MMFfldTOR_ENE]
else:
msg = 'getTorsionEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getImpTorsionEnergy(self):
"""
get the imporoper torsion energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldIMP_ENE:
return ene_list[mm.MMFfldIMP_ENE]
else:
msg = 'getImpTorsionEnergy: requested energy not available '
raise ValueError(msg)
[docs] def get14LJEnergy(self):
"""
get the 1,4-Lennard Jones energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldLJ14_ENE:
return ene_list[mm.MMFfldLJ14_ENE]
else:
msg = 'get14LJEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getLJEnergy(self):
"""
get the Lennard Jones energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldLJ_ENE:
return ene_list[mm.MMFfldLJ_ENE]
else:
msg = 'getLJEnergy: requested energy not available '
raise ValueError(msg)
[docs] def get14EleEnergy(self):
"""
get the 1,4-electrostatic energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldELE14_ENE:
return ene_list[mm.MMFfldELE14_ENE]
else:
msg = 'get14EleEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getEleEnergy(self):
"""
get the electrostatic energy of the system
"""
ene_list = mm.mmffld_getTotalEnergyTerms(self.handle)
if len(ene_list) > mm.MMFfldELE_ENE:
return ene_list[mm.MMFfldELE_ENE]
else:
msg = 'getEleEnergy: requested energy not available '
raise ValueError(msg)
[docs] def getBondedEnergyComponents(self, stretch_item_list, bending_item_list,
torsion_item_list, imp_torsion_item_list):
"""
Get energy components for bonded terms given in an list of
stretches/bend/torsions/improper torsions (given as integer numbers)
in _item_list. The _energy_list arguments are ignored.
Returns the corresponding energies as a tuple of double lists:
(stretch_energy_list, bending_energy_list, torsion_energy_list,
imp_torsion_energy_list).
Atoms are counted starting at zero.
NOTE: The API of this function may change in the next release.
"""
return mm.mmffld_getBondedEnergyComponents(self.handle,
stretch_item_list,
bending_item_list,
torsion_item_list,
imp_torsion_item_list)
[docs] def getNBEnergyForTwoAtomLists(self, atom_list_i, atom_list_j):
"""
Get the non-bonded energy components for two arbitrary list of atoms.
The total Coulomb and LJ interaction energies are returned.
"""
return mm.mmffld_getNBEnergyForTwoAtomLists(self.handle, atom_list_i,
atom_list_j)
[docs] def getTotalAtom(self):
""" get the total number of atoms """
return mm.mmffld_getParameterTotal(self.handle, mm.MMFfldNATOM)
[docs] def getTotalStretch(self):
""" get the total number of stretch """
return mm.mmffld_getParameterTotal(self.handle, mm.MMFfldNSTR)
[docs] def getStretch(self, istr):
""" get the stretch type information and parameters """
(i, j, fc, r0, icb, quality,
comment) = mm.mmffld_getStretch_reveal_ok(self.handle, istr)
return (i, j, fc, r0, quality, comment)
[docs] def getTotalBending(self):
""" get the total number of bending """
return mm.mmffld_getParameterTotal(self.handle, mm.MMFfldNBND)
[docs] def getBending(self, ibnd):
""" get the bending type information and parameters """
(i, j, k, fc, theta0, ict, quality,
comment) = mm.mmffld_getBending_reveal_ok(self.handle, ibnd)
return (i, j, k, fc, theta0, quality, comment)
[docs] def getTotalTorsion(self):
""" get the total number of torsion """
return mm.mmffld_getParameterTotal(self.handle, mm.MMFfldNTOR)
[docs] def getTorsion(self, itor):
""" get the torsion type information and parameters """
(i, j, k, l, vs, type, quality,
comment) = mm.mmffld_getTorsion_reveal_ok(self.handle, itor)
return (i, j, k, l, vs, quality, comment)
[docs] def getTotalImpTorsion(self):
""" get the total number of improper torsion """
return mm.mmffld_getParameterTotal(self.handle, mm.MMFfldNIMPTOR)
[docs] def getImpTorsion(self, imp_tor):
""" get the improper torsion type information and parameters """
(i, j, k, l, v2, type, quality,
comment) = mm.mmffld_getImpTorsion_reveal_ok(self.handle, imp_tor)
return (i, j, k, l, v2, quality, comment)
[docs] def getTotal14Pair(self):
""" get the total number of 1,4-pairs """
return mm.mmffld_getParameterTotal(self.handle, mm.MMFfldN14PAIR)
[docs] def getNB14Pairs(self, i14pair):
""" get the 1,4-pair information """
(i, j) = mm.mmffld_getNB14Pairs(self.handle, i14pair)
return (i, j)
[docs] def getExcludedAtomList(self, iatom):
""" get the excluded atoms """
try:
return mm.mmffld_getExcludedAtomList(self.handle, iatom)
except:
msg = 'getExcludedAtomList cannot obtain exlusions for atom %d ' % iatom
raise ValueError(msg)
[docs] def setStructure(self, struct):
"""
Change the structure to be minimized.
This method generally won't need to be called by the end user. To
minimize different molecules, it is preferable to create separate
Minimizer instances.
If only the coordinates of the structure changed, it is much faster
to call updateCoordinates() instead of setStructure().
"""
if self._struct == struct:
return
# Change the structure only if a different CT handle is specified
if self._struct is not None:
# Delete the previous structure (if present):
mm.mmffld_deleteMol(self.handle, self._struct)
self._struct = None
if self.cleanup:
ifo, r, c = mm.mmlewis_get_option(mm.MMLEWIS_OPTION_IFO)
del_dummies, r, c = mm.mmlewis_get_option(
mm.MMLEWIS_OPTION_MMLEWIS_APPLY_DELETE_DUMMIES)
# replicating the mmlewis treatment that all backends do
# except for deleting dummy atoms.
with mmerr.Capture(handle=mm.mmlewis_get_error_handler()):
try:
mm.mmlewis_set_option(mm.MMLEWIS_OPTION_IFO, 5, 0.0, "")
mm.mmlewis_set_option(
mm.MMLEWIS_OPTION_MMLEWIS_APPLY_DELETE_DUMMIES, 0, 0.0,
"")
mm.mmlewis_apply(struct)
except mm.MmException as e:
if e.rc != mm.MMLEWIS_ERR:
raise
raise ValueError(
"Lewis structure check failed on the "
"specified structure.\nCheck formal charges and bond "
f"orders.\n{e}")
finally:
mm.mmlewis_set_option(mm.MMLEWIS_OPTION_IFO, ifo, 0.0, "")
mm.mmlewis_set_option(
mm.MMLEWIS_OPTION_MMLEWIS_APPLY_DELETE_DUMMIES,
del_dummies, 0.0, "")
if self._honor_pbc:
nbc = self._users_nonbonded_cutoff
# Import xtal here to avoid having circular imports
from schrodinger.application.matsci.nano import xtal
try:
avec, bvec, cvec, alpha, beta, gamma = xtal.get_params_from_chorus(
xtal.get_chorus_properties(struct))
except KeyError:
# PBC properties are missing
pass
else:
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_BX, 0,
avec, "")
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_BY, 0,
bvec, "")
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_BZ, 0,
cvec, "")
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_ALPHA, 0,
alpha, "")
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_BETA, 0,
beta, "")
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_GAMMA, 0,
gamma, "")
# Limit the non-bonded intersction cutoff to 45% of the smallest box
# dimension. This prevents the situation where atom pairs have 2 sets
# of NB interactions (one within the box, and one across a PBC).
minvec = min([avec, bvec, cvec])
nbc = min(nbc, minvec * .45)
# Set cutoff value, overwriting one set for the previous structure,
# without overwriting self._users_nonbonded_cutoff:
mm.mmffld_setOption(self.handle, mm.MMFfldOption_SYS_CUTOFF, 0, nbc,
"")
try:
mm.mmffld_enterMol(self.handle, struct)
except mm.MmException as e:
if e.rc == mm.MMFFLD_ERR:
msg = "Force field failed on specified structure."
raise ValueError(msg)
self._struct = struct
return
[docs] def getEnergy(self, recalc=None):
"""
Run a zero-step minimization to calculate the energy, and return it.
Units are kcal/mol.
:type recalc: bool
:param recalc: Deprecated, and will propably go away.
"""
if recalc is not None:
warnings.warn("Minimizer.getEnergy(): recalc option is obsolete.",
stacklevel=2)
# Temporarily set max_steps to zero to get an energy evaluation at
# the current structure.
orig_max_steps = self.max_steps
energy = None
try:
self.max_steps = 0
min_res = self.minimize()
energy = min_res.potential_energy
finally:
self.max_steps = orig_max_steps
return energy
[docs] def setOptions(self, **kwargs):
"""
Set Minimizer options by keyword.
To be used only internally by the constructor.
"""
for attr, value in kwargs.items():
setattr(self, attr, value)
[docs] def getOptions(self):
"""
Return the values for all properties that control Minimizer
behavior.
"""
# Find all public class attributes that have an fset method that is
# not None; these are properties with set methods. This ignores
# read-only properties.
options = {}
for attr in dir(Minimizer):
if not attr.startswith(("_", "next")):
if getattr(Minimizer.__dict__[attr], "fset", None):
options[attr] = getattr(self, attr)
return options
[docs] def setLigandCleanupOptions(self):
"""
Obsolete.
"""
msg = "Minimizer.setLigandCleanupOptions() can no longer be used."
msg += "\nPlease set energy_no_electrostatics and energy_lj_repulsive_only to True when creating the Minimizer class"
raise RuntimeError(msg)
[docs] def addTorsionRestraint(
self,
i,
j,
k,
l, # noqa: E741, bond atom
force_constant,
target,
flat_bottom=0.0):
"""
Define a torsional restraint.
:type i: int
:type j: int
:type k: int
:type l: int
:param i, j, k, l: atoms defining the restraint
:type force_constant: float
:param force_constant: force constant
:type target: float
:param target: target value
:type flat_bottom: float
:param flat_bottom: flat-bottom half width (default 0.0)
"""
mm.mmffld_enterRestraint(self.handle, mm.MMFfldTorsionRest,
mm.MMCT_INVALID_CT, i, j, k, l, force_constant,
target, flat_bottom)
[docs] def addAngleRestraint(self,
i,
j,
k,
force_constant,
target,
flat_bottom=0.0):
"""
Define an angle restraint.
:type i: int
:type j: int
:type k: int
:param i, j, k: atoms defining the restraint
:type force_constant: float
:param force_constant: force constant
:type target: float
:param target: target value
:type flat_bottom: float
:param flat_bottom: flat-bottom half width (default 0.0)
"""
mm.mmffld_enterRestraint(self.handle, mm.MMFfldAngleRest,
mm.MMCT_INVALID_CT, i, j, k, 0, force_constant,
target, flat_bottom)
[docs] def addDistanceRestraint(self,
i,
j,
force_constant,
target,
flat_bottom=0.0):
"""
Define a distance restraint.
:type i: int
:type j: int
:param i, j: atoms defining the restraint
:type force_constant: float
:param force_constant: force constant
:type target: float
:param target: target value
:type flat_bottom: float
:param flat_bottom: flat-bottom half width (default 0.0)
"""
mm.mmffld_enterRestraint(self.handle, mm.MMFfldDistRest,
mm.MMCT_INVALID_CT, i, j, 0, 0, force_constant,
target, flat_bottom)
[docs] def addPosRestraint(self, i, force_constant, flat_bottom=0.0):
"""
Define a positional restraint.
:type i: int
:param i: atom defining the restraint
:type force_constant: float
:param force_constant: force constant
:type flat_bottom: float
:param flat_bottom: flat-bottom half width (default 0.0)
"""
mm.mmffld_enterRestraint(self.handle, mm.MMFfldPosRest,
mm.MMCT_INVALID_CT, i, 0, 0, 0, force_constant,
0.0, flat_bottom)
[docs] def addPosFrozen(self, i):
"""
Define a frozen atom.
:type i: int
:param i: atom to be frozen
"""
mm.mmffld_enterRestraint(self.handle, mm.MMFfldPosFrozen,
mm.MMCT_INVALID_CT, i, 0, 0, 0, 0.0, 0.0, 0.0)
[docs] def deleteAllRestraints(self, rest_type=mm.MMFfldAllRestType):
"""
Delete all the restraints for a given type (default all types).
:type rest_type: mm Constant
:param rest_type: One of:
- mm.MMFfldTorsionRest
- mm.MMFfldPosRest
- mm.MMFfldPosFrozen
- mm.MMFfldAllRestType
"""
mm.mmffld_deleteAllRestraints(self.handle, rest_type,
mm.MMCT_INVALID_CT)
[docs] def printParameters(self):
"""
Print all parameters for the entered molecules.
"""
mm.mmffld_printParameters(self.handle)
# Settings control how the minimizer behaves and are exposed here as
# properties generated by the _minimizer_property function.
min_method = _minimizer_int_property(
mm.MMFfldOption_MIN_METHOD,
doc="See the constructor method for documentation.",
readonly=False # Can be set after calling enterMol()
)
max_steps = _minimizer_int_property(
mm.MMFfldOption_MIN_MAX_STEP,
doc="See the constructor method for documentation.",
readonly=False # Can be set after calling enterMol()
)
energy_criterion = _minimizer_float_property(
mm.MMFfldOption_MIN_CONV_ENER,
doc="See the constructor method for documentation.",
readonly=False # Can be set after calling enterMol()
)
gradient_criterion = _minimizer_float_property(
mm.MMFfldOption_MIN_CONV_GRAD,
doc="See the constructor method for documentation.",
readonly=False # Can be set after calling enterMol()
)
verbose = _minimizer_int_property(
mm.MMFfldOption_MIN_VERBOSE,
doc="See the constructor method for documentation.",
readonly=False # Can be set after calling enterMol()
)
energy_no_electrostatics = _minimizer_int_property(
mm.MMFfldOption_ENE_NO_ELE,
doc="See the constructor method for documentation.",
readonly=True # Can NOT be set after enterMol() was called
)
energy_lj_repulsive_only = _minimizer_int_property(
mm.MMFfldOption_ENE_LJ_REPULSIVE_ONLY,
doc="See the constructor method for documentation.",
readonly=True # Can NOT be set after enterMol() was called
)
nonbonded_cutoff = _minimizer_float_property(
mm.MMFfldOption_SYS_CUTOFF,
doc="See the constructor method for documentation.",
readonly=False # Can be set after calling enterMol()
)
maintain_planarity = _minimizer_int_property(
mm.MMFfldOption_TYPER_PLANAR_SC_ENABLE,
doc="See the constructor method for documentation.",
readonly=True # Can NOT be set after enterMol() was called
)
no_zob_restrain = _minimizer_int_property(
mm.MMFfldOption_TYPER_NO_RESTRAIN_ZOB,
doc="See the constructor method for documentation.",
readonly=True # Can NOT be set after enterMol() was called
)
charges_from_ct = _minimizer_int_property(
mm.MMFfldOption_ENE_CHARGES_FROM_CT,
doc="See the constructor method for documentation.",
readonly=True # Can NOT be set after enterMol() was called
)
# (EV 126347, 126368)
dielectric_constant = _minimizer_float_property(
mm.MMFfldOption_ENE_DIELECTRIC,
doc="See the constructor method for documentation.",
readonly=True)
# Results are stored in the MMCT, and are made available here as
# Minimizer properties.
[docs] def getBondedEnergies(self, atom_subset):
"""
Return bonded energies for bonded interactions within the given
atom subset.
:type atom_subset: List of atom indices.
:param atom_subset: Atoms to calculate bonded energies for.
:return: (stretch-energy, bend-energy, torsion-energy, imp-tor-energy)
:rtype: tuple
"""
if self._ffld_version == OPLS3:
raise ValueError("getBondedEnergy() is not supported by OPLS3")
# We need a zero-indexed list here:
subset_atoms_array = [ia - 1 for ia in atom_subset]
total_estr, total_ebnd, total_etor, total_eimptor = mm.mmffld_getBondedEnergiesForSubset(
self.handle, subset_atoms_array)
return total_estr, total_ebnd, total_etor, total_eimptor
[docs] def getBondedInteractionEnergy(self,
atomset1,
atomset2=None,
verbose=False):
"""
Return the bonded component of the interaction energy between
the two given atom sets. If the second subset is not specified,
instead use all other atoms from the CT.
:param atomset1: Atoms from the first subset.
:type atomset1: List of atom indices.
:param atomset2: Atoms from the second subset.
:type atomset2: List of atom indices.
:return: The bonded interaction energy.
:rtype: float
"""
def consider_bond(i, j):
if atomset2 is None:
return (i in atomset1 and
j not in atomset1) or (j in atomset1 and
i not in atomset1)
else:
return (i in atomset1 and j in atomset2) or (i in atomset2 and
j in atomset1)
def consider_bend(i, j, k):
# Consider bends only if 1 of the atoms is in one set, and the others are in the other:
if atomset2 is None:
if i in atomset1 and k not in atomset1:
if j in atomset1 or j not in atomset1:
return True
elif k in atomset1 and i not in atomset1:
if j in atomset1 or j not in atomset1:
return True
else:
if i in atomset1 and k in atomset2:
if j in atomset1 or j in atomset2:
return True
elif i in atomset2 and k in atomset1:
if j in atomset1 or j in atomset2:
return True
return False
if verbose:
print("getBondedInteractionEnergy()")
print(" Number of atoms in subset1:", len(atomset1))
if atomset2 is None:
print(" Number of atoms in subset2 (rest of the CT):",
self._struct.atom_total - len(atomset1))
else:
print(" Number of atoms in subset2:", len(atomset2))
# Make a list of stretches that we want to consider.
# Each item is an index to the Minimizer's internal stretch list
stretch_item_list = []
# Add only those stretches that we want to consider:
for istr in range(self.getTotalStretch()):
(i, j, fc, r0, quality, comment) = self.getStretch(istr)
if consider_bond(i + 1, j + 1):
#debug
# print 'getBondedInteractionEnergy istr:', istr, ': ', i+1, j+1
stretch_item_list.append(istr)
# Make a list of bending that we want to consider:
# Each item is an index to the Minimizer's internal bending list
bending_item_list = []
for ibnd in range(self.getTotalBending()):
(i, j, k, fc, theta0, quality, comment) = self.getBending(ibnd)
# Consider bends only if 1 of the atoms is in one set, and the others are in the other:
if consider_bend(i + 1, j + 1, k + 1):
#debug
# print 'getBondedInteractionEnergy ibnd:', ibnd, ': ', i+1, j+1, k+1
bending_item_list.append(ibnd)
# Make a list of torsions that we want to consider:
# Each item is an index to the Minimizer's internal torsions list
torsion_item_list = []
for itor in range(self.getTotalTorsion()):
(i, j, k, l, vs, quality, comment) = self.getTorsion(itor)
# The two atoms involved in the torsion bond need to be in the separate subsets:
if consider_bond(j + 1, k + 1):
#debug
# print 'getBondedInteractionEnergy itor:', itor, ': ', i+1, j+1, k+1, l+1
torsion_item_list.append(itor)
# Make a list of improper torsions that we want to consider:
# Each item is an index to the Minimizer's internal improper torsions list
imp_torsion_item_list = []
for iitor in range(self.getTotalImpTorsion()):
(i, j, k, l, v2, quality, comment) = self.getImpTorsion(iitor)
# The two atoms involved in the torsion bond need to be in the separate subsets:
if consider_bond(j + 1, k + 1):
imp_torsion_item_list.append(iitor)
# Empty lists (vectors) that the minimizer will populate:
# get the bonded energy components for the selected bonded types
ret = self.getBondedEnergyComponents(stretch_item_list,
bending_item_list,
torsion_item_list,
imp_torsion_item_list)
stretch_energy_list, bending_energy_list, torsion_energy_list, imp_torsion_energy_list = ret
# Sum up all the relevant stretch components:
total_estr = sum(stretch_energy_list)
# Sum up all the relevant bending components:
total_ebnd = sum(bending_energy_list)
# Sum up all the relevant torsion components:
total_etor = sum(torsion_energy_list)
# Sum up all the relevant improper-torsion components:
total_eimptor = sum(imp_torsion_energy_list)
return total_estr + total_ebnd + total_etor + total_eimptor
[docs] def getNonBondedEnergies(self, atom_subset, verbose=False):
"""
Given a list of atoms representing a substructure (e.g. a residue),
return a tuple of:
1. the sum of non-bonded energies WITHIN the set.
2. the sum of non-bonded energies between the set atoms and everything
else.
:type atom_subset: list(int)
:param atom_subset: Atoms indicies to calculate non-bonded energies for
:return: (internal non-bonded energy in kcal/mol, non-bonded
interaction energy in kcal/mol)
:rtype: tuple
"""
if verbose:
print("enter getNonBondedEnergies()")
print('number of CT atoms:', len(self._struct.atom))
print('considering num atoms:', len(atom_subset))
if verbose:
print('build lists...')
# Make a list of atoms that we need to consider:
# (We need a zero-indexed list here)
subset_atoms_array = [ia - 1 for ia in atom_subset]
if verbose:
print('Calculating non-bonded energy components...')
intra_nb_energy = 0.0
inter_nb_energy = 0.0
intra_nb_energy, inter_nb_energy = mm.mmffld_getNonBondedEnergiesForSubset(
self.handle, subset_atoms_array)
# Ev:131438 must cast to float, to avoid returning a numby float:
return float(intra_nb_energy), float(inter_nb_energy)
[docs] def getSelfEnergy(self, atom_subset, include_intra_nb=True):
"""
Return the internal energy of the the given atoms. This is the sum
of all bonded interactions between the given atoms.
:type atom_subset: List of atom indices.
:param atom_subset: Atoms to calculate self-energy for.
:type include_intra_nb: bool
:param include_intra_nb:
If True, include non-bonded interaction energy within the given
subset. If False, include only bonded interactions.
:return: Total self-energy
:rtype: float
"""
if self._ffld_version == OPLS3:
raise ValueError("getSelfEnergy() does not support OPLS3")
# Calculate the bonded energies WITHIN the subset:
total_estr, total_ebnd, total_etor, total_eimptor = self.getBondedEnergies(
atom_subset)
total_energy = total_estr + total_ebnd + total_etor + total_eimptor
if include_intra_nb:
intra_nb_energy, inter_nb_energy = self.getNonBondedEnergies(
atom_subset, verbose=False)
total_energy += intra_nb_energy
return total_energy
def _getInteractionEnergyBetweenSubsets(self,
atom_subset1,
atom_subset2,
include_bonded=False,
verbose=False):
"""
Returns the interaction energy between the two atom groups.
Private method, use getInteractionEnergy() instead.
:type atom_subset1: List of atom indices.
:param atom_subset1: First atom set.
:type atom_subset2: List of atom indices.
:param atom_subset2: Second atom set.
:type include_bonded: Boolean.
:param include_bonded: Whether to also include bonded terms. By
default, only non-bonded interactions are included, even if the
atom sets are covalently bonded, and any bonded interactions are
excluded.
:return: Interaction energy betwen the given atom sets.
(key: j atom; value: interaction energy)).
:rtype: float
"""
st = self._struct
if verbose:
print("enter _getInteractionEnergyBetweenSubsets()")
print('total CT atoms:', len(st.atom))
print('atoms in set 1:', len(atom_subset1))
print('atoms in set 2:', len(atom_subset2))
print('include_bonded:', include_bonded)
print('build lists...')
# Make lists of atoms that we need to consider:
# (We need a zero-indexed list here)
subset1_atoms_array = [ia - 1 for ia in atom_subset1]
subset2_atoms_array = [ia - 1 for ia in atom_subset2]
if verbose:
print('Calculating non-bonded energy components...')
E_LJ, E_C = self.getNBEnergyForTwoAtomLists(subset1_atoms_array,
subset2_atoms_array)
tot_inter_energy = E_LJ + E_C
if verbose:
print(" Total considered NB interaction energy: %.4f" %
tot_inter_energy)
if include_bonded:
bonded_inter_energy = self.getBondedInteractionEnergy(
atom_subset1, atom_subset2)
tot_inter_energy += bonded_inter_energy
if verbose:
print(" Total considered bonded interaction energy: %.4f" %
bonded_inter_energy)
print(" Total considered interaction energy: %.4f" %
tot_inter_energy)
return tot_inter_energy
[docs] def getInteractionEnergy(self,
subset_atoms,
consider_atoms=None,
include_bonded=False,
verbose=False):
"""
Calculate the energy of interaction between the <subset_atoms> and
the <consider_atoms>. If <consider_atoms> is not specified, then
returns a sum of all interactions between the <subset_atoms> and the
rest of the atoms.
Note that in the default case where <consider_atoms> is None, atoms
outside of the cutoff are not considered.
:type subset_atoms: List of atom indices.
:param subset_atoms: First atom set.
:type consider_atoms: List of atom indices.
:param consider_atoms: Second atom set (optional)
:type include_bonded: Boolean.
:param include_bonded: Whether to also include bonded terms. By default,
only non-bonded interactions are included, even if the atom sets
are covalently bonded, and any bonded interactions are excluded.
:return: Interaction energy betwen the given atom sets (or between the
given atom sets and the rest of the structure).
:rtype: float
"""
if consider_atoms is not None:
return self._getInteractionEnergyBetweenSubsets(
subset_atoms, consider_atoms, include_bonded, verbose)
if verbose:
print("enter getInteractionEnergy()")
print(" Number of subset_atoms:", len(subset_atoms))
print(" include_bonded:", include_bonded)
intra_nb_energy, inter_nb_energy = self.getNonBondedEnergies(
subset_atoms, verbose=verbose)
tot_inter_energy = inter_nb_energy
if verbose:
print(" Total considered NB interaction energy: %.4f" %
tot_inter_energy)
if include_bonded:
bonded_inter_energy = self.getBondedInteractionEnergy(subset_atoms)
tot_inter_energy += bonded_inter_energy
if verbose:
print(" Total considered bonded interaction energy: %.4f" %
bonded_inter_energy)
print(" Total considered interaction energy: %.4f" %
tot_inter_energy)
return tot_inter_energy
[docs]@contextlib.contextmanager
def minimizer_context(**kwargs):
"""
Context manager to patch the Minimizer class. ::
with minimizer_context(**kwargs) as min:
min.setStructure(st)
min.minimize()
"""
minimizer = Minimizer(**kwargs)
try:
yield minimizer
finally:
minimizer.cleanup_resources()
[docs]def compute_energy(struct, **kwargs):
"""
Compute the energy in kcal/mol, store it on the structure, and return it.
:type struct: `schrodinger.structure.Structure`
:param struct: the structure for which to compute the energy
:rtype: float
:return: the energy in kcal/mol
"""
kwargs['struct'] = struct
with minimizer_context(**kwargs) as mizer:
return mizer.getEnergy()