Source code for schrodinger.application.desmond.stage.app.fragment_linking.stage
from typing import List
from typing import Optional
from schrodinger import structure
from schrodinger.application.desmond import cmj
from schrodinger.application.desmond import struc
from schrodinger.application.desmond.constants import FEP_ABSOLUTE_ENERGY
from schrodinger.application.desmond.constants import FEP_ABSOLUTE_LIGAND
from schrodinger.application.desmond.constants import FEP_FRAGNAME
from schrodinger.application.desmond.constants import FEP_MAPPING
from schrodinger.application.desmond.constants import FEP_RESTRAIN
from schrodinger.application.desmond.constants import FEP_RESTRAINT_CORRECTION
from schrodinger.application.desmond.constants import \
FRAGMENT_LINKING_HYDRATION_JOBS
from schrodinger.application.desmond.constants import FRAGMENT_LINKING_JOBS
from schrodinger.application.desmond.constants import FRAGMENT_LINKING_SIGN
from schrodinger.application.desmond.stage.prepare import forcefield
from schrodinger.utils import sea
__all__ = ['FragmentLinkingPrimer']
[docs]class FragmentLinkingPrimer(cmj.StructureStageBase):
NAME = "fragment_linking_primer"
PARAM = cmj._create_param_when_needed([
"""
DATA = {
job = %s
restraint = {
force_constants = [1.0 78.0]
temperature = 300.0
lambda_schedule_name = "alchemical_restraint"
}
}
VALIDATE = {
job = {type = enum range = [%s]}
restraint = {
force_constants = {type = list size = 2 elem = {type = float+}}
temperature = {type = float+}
lambda_schedule_name = {type = str1}
}
}
""" % (FRAGMENT_LINKING_JOBS.VAL.SOLVENT, FRAGMENT_LINKING_JOBS.VAL)
])
[docs] def run(self, jobname: str, mae_fname: str) -> Optional[str]:
cts = list(structure.StructureReader(mae_fname))
# Determine which structure contains the fragments
ref_ict = None
mut_ict = None
for ict, ct in enumerate(cts):
if ct.mol_total == 1:
ref_ict = ict
elif ct.mol_total == 2:
mut_ict = ict
else:
ref_ict = None
mut_ict = None
break
is_direction_fragments_to_linked = ref_ict > mut_ict
# Check for invalid input structures
if ref_ict is None or mut_ict is None or len(cts) != 2:
error_msg = (
f"ERROR: {self.NAME} failed: Invalid structures."
f" Must have one ct containing the whole molecule and one ct containing two fragments."
)
return None
if 0 not in _charges_per_mol(cts[mut_ict]):
error_msg = (
f"ERROR: {self.NAME} failed: Invalid structures."
f" At least one fragment must be neutral, both have formal charges."
)
return None
# Clear existing restraint property
forcefield.clear_restraints(cts[ref_ict])
forcefield.clear_restraints(cts[mut_ict])
# Set restraints if used for this job
alchemical = False
restrain = False
lambda_state = None
if self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.SOLVENT:
restrain = True
alchemical = True
lambda_state = 0 if is_direction_fragments_to_linked else 1
elif self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.RESTRAINED_FRAGMENT_HYDRATION:
restrain = True
elif self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.FRAGMENT_HYDRATION:
restrain = False
if restrain:
fc = self.param.restraint.force_constants.val
sigma = [0, 0] # Fixed for now, not present in correction term
schedule_name = self.param.restraint.lambda_schedule_name.val
new_restraints = _get_fragment_restraints(cts[mut_ict], fc, sigma,
schedule_name,
lambda_state)
self._log("Multisim settings for fragment linking restraints:\n"
f"{new_restraints}")
forcefield.encode_restraints(cts[mut_ict], new_restraints)
if self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.RESTRAINED_FRAGMENT_HYDRATION:
cts[mut_ict].property[
FEP_RESTRAINT_CORRECTION] = forcefield.calculate_restraint_correction_term(
new_restraints, self.param.restraint.temperature.val)
# Determine the sign for the hydration dG based on the edge direction
if self.param.job.val in FRAGMENT_LINKING_HYDRATION_JOBS:
sign = 1 if is_direction_fragments_to_linked else -1
cts[mut_ict].property[FRAGMENT_LINKING_SIGN] = sign
# Prepare the output structures
if self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.SOLVENT:
out_cts = cts
elif self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.RESTRAINED_FRAGMENT_HYDRATION:
# Clear relative FEP property
mut_ct = cts[mut_ict]
struc.delete_structure_properties(mut_ct, FEP_FRAGNAME)
# Set up absolute FEP
mut_ct_abs = cts[mut_ict].copy()
mol_for_hydration = _mol_for_hydration(mut_ct_abs)
# Must be set explicitly on all atoms
for atom in mut_ct_abs.atom:
atom.property[FEP_ABSOLUTE_ENERGY] = 0
atom.property[FEP_ABSOLUTE_LIGAND] = 0
for atom in mol_for_hydration.atom:
atom.property[FEP_ABSOLUTE_ENERGY] = 1
atom.property[FEP_ABSOLUTE_LIGAND] = 1
out_cts = [mut_ct_abs]
elif self.param.job.val == FRAGMENT_LINKING_JOBS.VAL.FRAGMENT_HYDRATION:
# Extract the smallest fragment
mol_for_hydration = _mol_for_hydration(cts[mut_ict])
frag_ct = mol_for_hydration.extractStructure(copy_props=True)
# Clear relative FEP property
struc.delete_structure_properties(frag_ct, FEP_FRAGNAME)
# Set up absolute FEP
frag_ct_abs = frag_ct.copy()
for atom in frag_ct_abs.atom:
atom.property[FEP_ABSOLUTE_ENERGY] = 1
atom.property[FEP_ABSOLUTE_LIGAND] = 1
out_cts = [frag_ct_abs]
for ct in out_cts:
ct.property[FRAGMENT_LINKING_JOBS] = self.param.job.val
out_fname = f"{jobname}-out.mae"
struc.write_structures(out_cts, out_fname)
return out_fname
def _charges_per_mol(ct: structure.Structure) -> List[int]:
"""
Given a structure, return the formal charge for each
molecule in that structure.
"""
return [sum(a.formal_charge for a in m.atom) for m in ct.molecule]
def _mol_for_hydration(ct: structure.Structure) -> structure._Molecule:
"""
Given a structure containing two fragments, return the molecule
that has the smallest number of atoms if both fragments are
neutral. If one fragment is neutral and the other is not,
return the neutral one. Both fragments must not be charged.
:param ct: Input structure.
"""
charges = _charges_per_mol(ct)
if sum(charges) == 0:
return min(ct.molecule, key=lambda mol: len(mol))
else:
# Molecule is 1-indexed
return ct.molecule[charges.index(0) + 1]
def _get_fragment_restraints(
ct: structure.Structure, fc: List[float], sigma: List[float],
schedule_name: str,
lambda_state: Optional[int]) -> forcefield.Restraints:
"""
Given a structure containing two fragments, find a optimal distance
and angle restraints.
:param ct: Input structure, modified in place.
:param fc: List of force constants. The first element is for the stretch term,
and the second is for the angle terms.
:param sigma: List of sigmas.
:param lambda_state: If not None, the lambda that corresponds to the fully
restrained state. If None, use a fixed restraint.
"""
# These modules have dependency on the Desmond product.
from schrodinger.application.desmond.packages import restraint
# Calculate restraints
# Only consider core heavy atoms for each fragment
mol_nums = sorted(set(atom.molecule_number for atom in ct.atom))
frag0_asl = f"(m.n {mol_nums[0]}) and (not atom.att 1) and (not atom.{FEP_MAPPING} 0)"
frag1_asl = f"(m.n {mol_nums[1]}) and (not atom.att 1) and (not atom.{FEP_MAPPING} 0)"
flr = restraint.FragmentLinkingRestraint.findRestraint(
ct, fragment0_asl=frag0_asl, fragment1_asl=frag1_asl)
# Mark the restrained atoms
ct.atom[flr.A].property[FEP_RESTRAIN] = 1
ct.atom[flr.a].property[FEP_RESTRAIN] = 2
ct.atom[flr.B].property[FEP_RESTRAIN] = 3
ct.atom[flr.b].property[FEP_RESTRAIN] = 4
atom_asl_dict = {
'A': f"a.{FEP_RESTRAIN} 1",
'a': f"a.{FEP_RESTRAIN} 2",
'B': f"a.{FEP_RESTRAIN} 3",
'b': f"a.{FEP_RESTRAIN} 4",
}
restraint_str = flr.asMsjSetting(fc, sigma, schedule_name, lambda_state,
atom_asl_dict)
return sea.Map(f"restraint = [{restraint_str}]").restraint