"""
Collection of functions to repair structures to be used for desmond simulation
Copyright Schrodinger, LLC. All rights reserved.
"""
import os
import numpy
from typing import Optional
from typing import List
from typing import Tuple
from schrodinger.application.desmond import smarts
from schrodinger.infra import mm
from schrodinger import structure
from schrodinger.structutils.analyze import evaluate_asl
from schrodinger.application.desmond import constants
LIPIDS_PDB_TYPE = [
"DPPC", "DPPE", "DPPS", "POPC", "POPE", "POPS", "DOPC", "DOPE", "DOPS",
"DMPC", "PIP"
]
[docs]def create_cms_from_mae(input_fname: str, output_fname: str, membrane_asl: str,
solvent_asl: str, cosolvent_asl: str):
"""
This function reads in a structure file which can be either a single
full-system CT or a list of component CTs. If the input is a full-
system CT, then we woud pass it directly to `decompose_to_comp_ct()`.
However, if the inputs are a list of component CTs, we would generate a
full-system CT and then and then pass it to `decomopse_to_comp_ct()`.
"""
if os.path.exists(output_fname):
os.remove(output_fname)
streader = structure.StructureReader(input_fname)
full_st = None
st_list = []
for st in streader:
if st.property.get(
constants.CT_TYPE) == constants.CT_TYPE.VAL.FULL_SYSTEM:
full_st = st
break
else:
st_list.append(st)
if full_st:
if constants.NUM_COMPONENT in full_st.property:
del full_st.property[constants.NUM_COMPONENT]
else:
full_st = structure.create_new_structure()
for st in st_list:
full_st.extend(st)
for st in decompose_to_comp_ct(full_st, membrane_asl, solvent_asl,
cosolvent_asl):
st.append(output_fname)
[docs]def create_full_system_from_comp(input_fname, output_fname):
streader = structure.StructureReader(input_fname)
st_list = []
for st in streader:
if (constants.CT_TYPE in st.property and st.property[constants.CT_TYPE]
== constants.CT_TYPE.VAL.FULL_SYSTEM):
continue
else:
st_list.append(st)
full_st = structure.create_new_structure()
for st in st_list:
full_st.extend(st)
full_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.FULL_SYSTEM
full_st.write(output_fname)
for st in st_list:
st.append(output_fname)
[docs]def decompose_to_comp_ct(st: structure.Structure,
membrane_asl: str = "",
solvent_asl: str = "",
cosolvent_asl: str = "") -> List[structure.Structure]:
"""
Decompose st to component cts, solute, membrane, cosolvent,
solvent (if exists any).
:param st: Structure object representing a full-system CT
:return: list of non-empty structures that would comprize a CMS:
[full_system, membrane_st, ..., water_st]
"""
repair_box_vector(st)
membrane_atoms, membrane_st = _extract_membrane(st, membrane_asl)
ion_atoms, ion_st = _extract_ion(st)
water_atoms, water_st = _extract_water(st, solvent_asl)
cosolvent_atoms, cosolvent_st = _extract_cosolvent(st, cosolvent_asl)
atoms_to_delete = membrane_atoms + ion_atoms + water_atoms + cosolvent_atoms
if len(atoms_to_delete) != len(set(atoms_to_delete)):
raise ValueError("Selection ASLs matched the same atoms. Please "
"make sure that all selections matches are unique.")
solute_st = st.copy()
solute_st.deleteAtoms(atoms_to_delete)
solute_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.SOLUTE
# add atom group for biasing potential on the waters if membrane exists
if membrane_st and water_st:
for a in water_st.atom:
a.property[constants.FEP_ABSOLUTE_ENERGY] = 1
full_st = solute_st.copy() if solute_st else structure.Structure()
st_list = [solute_st]
for comp_st in [membrane_st, ion_st, cosolvent_st, water_st]:
if comp_st:
full_st.extend(comp_st)
st_list.append(comp_st)
full_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.FULL_SYSTEM
return [full_st] + st_list
[docs]def repair_box_vector(st: structure.Structure):
"""
repair box vectors when they are corrupted.
create box vectors based on system size
"""
for name in constants.SIM_BOX:
if name not in st.property:
print(f"WARNING: Property {name} is missing. Generating box "
"vectors based on system dimension.")
break
else:
return
for name in constants.SIM_BOX:
st.property[name] = 0.0
atom_xyz_array = st.getXYZ(copy=False)
xmin, ymin, zmin = numpy.min(atom_xyz_array, 0)
xmax, ymax, zmax = numpy.max(atom_xyz_array, 0)
vals = (xmax - xmin, ymax - ymin, zmax - zmin)
for prop, val in zip(constants.SIM_BOX_DIAGONAL, vals):
st.property[prop] = val
def _extract_water(
st: structure.Structure,
asl: str = "") -> Tuple[List[int], Optional[structure.Structure]]:
"""
Extract water molecules from full-system structure.
:return: a tuple of atom indices corresponding to full-system CT, followed
by an extracted Structure object.
"""
atoms = evaluate_asl(st, asl) if asl else \
smarts.evaluate_net_smarts(st, constants.WATER_SMARTS,
constants.ZOB_WATER_SMARTS)
if not atoms:
return [], None
tmp_st = st.extract(atoms, True)
# reorder atoms such that each all atoms in the same molecule
# are together
atom_order = [0] + [a.index for m in tmp_st.molecule for a in m.atom]
new_st = tmp_st.copy()
mm.mmct_ct_reorder(new_st, tmp_st, atom_order)
res_name = new_st.atom[1].pdbres
new_st.title = res_name + " water box"
new_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.SOLVENT
new_st.property[constants.NUM_COMPONENT] = 1
# overwrite residue # with molecule #
# viparr fails when two or more molecules have same residue #
for a in new_st.atom:
a.resnum = a.molecule_number
return atoms, new_st
def _extract_membrane(
st: structure.Structure,
asl: str = "") -> Tuple[List[int], Optional[structure.Structure]]:
"""
Extract Membrane lipids from full-system st.
assign reference custom charges for lipid molecules found at
'system_builder/data/lipid_charge.mae'
:return: a tuple of atom indices corresponding to full-system CT, followed
by an extracted Structure object.
"""
if not asl:
asl = 'res.ptype '
for r in LIPIDS_PDB_TYPE:
asl += f'\"{r}\" '
atoms = evaluate_asl(st, asl)
if not atoms:
return [], None
membrane_st = st.extract(atoms, True)
resname = membrane_st.atom[1].pdbres
membrane_st.title = resname + " bilayer"
membrane_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.MEMBRANE
homogeneous = set(a.pdbres for a in membrane_st.atom)
if len(homogeneous) == 1:
membrane_st.property[constants.NUM_COMPONENT] = 1
return atoms, membrane_st
def _extract_ion(
st: structure.Structure
) -> Tuple[List[int], Optional[structure.Structure]]:
"""
Extract monatomic ions from full-system st
:return: a tuple of atom indices corresponding to full-system CT, followed
by an extracted Structure object.
"""
atoms = evaluate_asl(st, constants.MONATOMIC_IONS_ASL)
if not atoms:
return [], None
new_st = st.extract(atoms, True)
res_name = new_st.atom[1].pdbres
new_st.title = res_name
new_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.ION
# overwrite residue # with molecule #
# viparr fails when two or more molecules have same residue #
for a in new_st.atom:
a.resnum = a.molecule_number
return atoms, new_st
def _extract_cosolvent(
st: structure.Structure,
asl: str = "") -> Tuple[List[int], Optional[structure.Structure]]:
"""
Extract cosolvent atoms from full-system st
:return: a tuple of atom indices corresponding to full-system CT, followed
by an extracted Structure object.
"""
if not asl:
return [], None
atoms = evaluate_asl(st, asl)
if not atoms:
return [], None
cosolvent_st = st.extract(atoms, True)
cosolvent_st.title = cosolvent_st.atom[1].pdbres.strip()
cosolvent_st.property[constants.CT_TYPE] = constants.CT_TYPE.VAL.COSOLVENT
# overwrite residue # with molecule #
# viparr fails when two or more molecules have same residue #
for a in cosolvent_st.atom:
a.resnum = a.molecule_number
return atoms, cosolvent_st