"""
This script is intended as a component of the BuildGeometry stage from
desmond's stage.py. The system_builder executable inserts extra ions based on
the charge difference between two FEP ligands, and this script moves those ions
from the ion structures to the correct ligand structure and sets the appropriate
properties such that the ions will be alchemically mutated along with the
ligands.
Copyright Schrodinger, LLC. All rights reserved.
"""
import numpy
from schrodinger import structure
from schrodinger.infra import mm
from . import constants
from .constants import ALCHEMICAL_ION
from .constants import ALCHEMICAL_WATER
# used to indicate an atom is part of any alchemical ion - type structure,
# be it charged or neutral ion or water
# specifically alchemical water
MAP_TO_NEUTRAL_DEFAULT = True
ALCHEMICAL_WATER_DEFAULT = True
ION_ATOMIC_NUMBERS = {-1: 17, 1: 11}
ION_NAMES = {-1: "Cl", 1: "Na"}
[docs]def move_ions_in_file(input_fname,
output_fname,
map_to_neutral_atom=MAP_TO_NEUTRAL_DEFAULT,
use_alchemical_water=ALCHEMICAL_WATER_DEFAULT,
print_fcn=lambda msg: None):
"""
Reads the structures in the input file, moves the alchemical ions to the
base fep ligand block, and writes the new structures to an output file.
Read the structures in the input file, move the alchemical ions to the
base fep ligand block, and write the new structures to an output file.
:param input_fname: the filename of the input file
:type input_fname: str
:param output_fname: the filename of the output file
:type output_fname: str
:param map_to_neutral_atom: see move_alchemical_ions for info
:type map_to_neutral_atom: bool
:param use_alchemical_water: see move_alchemical_ions for info
:type use_alchemical_water: bool
:param print_fcn: function for printing passed to move_alchemical_ions
:type print_fcn: callable
"""
with structure.StructureReader(input_fname) as reader:
all_sts = list(reader)
move_alchemical_ions(all_sts, map_to_neutral_atom, use_alchemical_water,
print_fcn)
with structure.StructureWriter(output_fname) as writer:
writer.extend(all_sts)
def _compute_ligand_translation(ct_ref, ct_mut):
"""
Compute translation vector for the first core atom in both the
reference and mutated structures if reference coordinate properties are defined
for that atom. If reference coordinates are not defined, translation vector
is zero. If all atoms are dummy, the translation vector is zero.
:param ct_ref: first molecule in an FEP calculation
:type ct_ref: stucture.Structure
:param ct_mut: second molecule in an FEP calculation
:type ct_mut: stucture.Structure
:rtype: list[numpy.array]
"""
atoms_list = list(
zip(*[(a.property[constants.FEP_MAPPING], a.index)
for a in ct_mut.atom
if a.property[constants.FEP_MAPPING]]))
if not atoms_list:
return [numpy.array([0.0, 0.0, 0.0])] * 2
ref_core_atoms, mut_core_atoms = atoms_list
trans_vectors = []
for core_atoms, ct in zip((ref_core_atoms, mut_core_atoms),
(ct_ref, ct_mut)):
vector = numpy.array([0.0, 0.0, 0.0])
for a in core_atoms:
atom = ct.atom[a]
# Reference coordinates are defined as atom properties,
# need to check if they exist.
if set(constants.REFERENCE_COORD_PROPERTIES).issubset(
set(atom.property)):
ref_coord = numpy.array([
atom.property[k]
for k in constants.REFERENCE_COORD_PROPERTIES
])
vector = ref_coord - atom.xyz
break
trans_vectors.append(vector)
return trans_vectors
[docs]def move_alchemical_ions(all_sts,
map_to_neutral_atom=MAP_TO_NEUTRAL_DEFAULT,
use_alchemical_water=ALCHEMICAL_WATER_DEFAULT,
print_fcn=lambda msg: None):
"""
Given a list of sts, find the base and mutant ligand and the ion st,
and move the final ions from the ion block to the base ligand such that
the net charge of the ligands are equal. (This allows us to use
system_builder to place the ions in solution, and then set up the
fep-mapping here). Right now this assumes the necessary number of alchemical
ions have been added to the end of the ion block, otherwise it will fail.
:param all_sts: a list of sts, which should contain a base and mutant
ligand and an ion st.
:type all_sts: list of `structure.Structure`
:param map_to_neutral_atom: If True, the alchemical ions added to the
mutant ligand are mapped to a neutral but otherwise equivalent atom
in the reference ligand.
If False (default), they are mapped to dummy atoms.
:type map_to_neutral_atom: bool
:param use_alchemical_water: if True, the solvent block is searched for
tagged water molecules which are moved to the reference ligand. The
oxygen molecules are mapped to ions in the mutant ligand, and the
hydrogens are mapped to dummy atoms.
:type use_alchemical_water: bool
:param print_fcn: function for printing
:type print_fcn: callable
:return: The input list with its structures updated such that alchemical
ions are moved to the base ligand st.
"""
print_fcn("Moving alchemical ions to FEP base ligand")
alchemical_water_st = None
ion_st = None
ligand_0 = None
ligand_1 = None
for idx, st in enumerate(all_sts):
if constants.FEP_FRAGNAME in st.property:
if st.property[constants.FEP_FRAGNAME] == 'none':
ligand_0 = idx
else:
ligand_1 = idx
else:
st_type = st.property.get(constants.CT_TYPE, '')
if use_alchemical_water:
if st_type == 'alchemical_water':
alchemical_water_st = st
elif st_type == 'ion':
ion_st = st
if ligand_0 is None:
raise ValueError("no reference ligand found in input file")
if ligand_1 is None:
raise ValueError("no mutant ligand found in input file")
q0 = all_sts[ligand_0].formal_charge
q1 = all_sts[ligand_1].formal_charge
needed_ions = abs(q1 - q0)
if not needed_ions:
print_fcn("No charge difference in ligands, no alchemical ions needed")
return
if use_alchemical_water:
charged_indices = []
charge = None
if alchemical_water_st is None:
raise ValueError("no solvent structure found in input file")
for atom in alchemical_water_st.atom:
this_charge = atom.property.get(ALCHEMICAL_WATER, None)
if this_charge:
charge = this_charge
charged_indices.append(atom.index)
if not charged_indices:
raise ValueError("Alchemical waters not marked with charge")
charged_alchemical_water_atoms = alchemical_water_st.extract(
charged_indices, copy_props=True)
for atom in charged_alchemical_water_atoms.atom:
atom.atomic_number = ION_ATOMIC_NUMBERS[charge]
atom.property['s_m_pdb_atom_name'] = ION_NAMES[charge]
atom.formal_charge = charge
qi = 0
extracted_ions = alchemical_water_st
extracted_ions_mut = charged_alchemical_water_atoms
else:
if ion_st is None:
# the presence of ion_st should be guaranteed by the fact that
# 1. add_alchemical_ions was called in system_builder_setup
# (should be a prerequisite for this function)
# 2. the delta charge is nonzero, as at least one ion has to be added
# in system builder for the ion structure to be created.
# so reaching this point should be treated as an error
raise ValueError("no ion structure found in input file")
qi = ion_st.formal_charge
total_ions = ion_st.atom_total
# take last needed_ions atoms
indices = range(total_ions, total_ions - needed_ions, -1)
extracted_ions = ion_st.extract(indices, copy_props=True)
extracted_ions_mut = extracted_ions
ion_st.deleteAtoms(indices)
print_fcn("Before move: ")
_print_charges(q0, q1, qi, print_fcn)
# Reference coordinates are set before the simulate box is constructed.
# Molecules might be translated to a different location.
# Need to figure out the translation vector to set reference coordinates
# for alchemical ions. Because snapcore may change either one of the
# ligands, each ligand needs to have its own vector.
trans_vectors = _compute_ligand_translation(all_sts[ligand_0],
all_sts[ligand_1])
move_ions_to_ligands(all_sts, extracted_ions, extracted_ions_mut, ligand_0,
ligand_1, map_to_neutral_atom, trans_vectors,
use_alchemical_water)
q0 = all_sts[ligand_0].formal_charge
q1 = all_sts[ligand_1].formal_charge
qi = ion_st.formal_charge if not use_alchemical_water else 0
print_fcn("After move")
_print_charges(q0, q1, qi, print_fcn)
if q0 != q1:
raise ValueError("ion structure does not contain correct alchemical "
"ions")
# we can safely get rid of this because nothing else can be in this st,
# but do it last so as not to invalidate ligand indices
if use_alchemical_water:
all_sts.remove(alchemical_water_st)
# rebuild fullsystem ct
assert all_sts[0].property[constants.CT_TYPE] == \
constants.CT_TYPE.VAL.FULL_SYSTEM
fsys = all_sts[0]
fsys.deleteAtoms(range(1, fsys.atom_total + 1))
for st in all_sts[1:]:
fsys.extend(st)
return all_sts
[docs]def move_ions_to_ligands(all_sts, extracted_ions, extracted_ions_mut, ligand_0,
ligand_1, map_to_neutral_atom, trans_vectors,
use_alchemical_water):
total_0 = all_sts[ligand_0].atom_total
total_1 = all_sts[ligand_1].atom_total
atoms_per_other = (extracted_ions.atom_total / extracted_ions_mut.atom_total
if use_alchemical_water else 1)
prepare_ions(extracted_ions_mut, total_0, trans_vectors[1], True,
map_to_neutral_atom, use_alchemical_water, atoms_per_other)
all_sts[ligand_1] = structure.Structure(
mm.mmffld_extendCTwithCustomCharges(all_sts[ligand_1].handle,
extracted_ions_mut.handle))
if map_to_neutral_atom or use_alchemical_water:
prepare_ions(extracted_ions, total_1, trans_vectors[0], False,
map_to_neutral_atom, use_alchemical_water,
1 / atoms_per_other)
all_sts[ligand_0] = structure.Structure(
mm.mmffld_extendCTwithCustomCharges(all_sts[ligand_0].handle,
extracted_ions.handle))
[docs]def prepare_ions(extracted_ions,
other_total,
trans_vectors,
is_mutant,
map_to_neutral_atom,
use_alchemical_water,
map_to_every_nth=1):
# use_alchemical_water implies map_to_neutral_atom for mutant ligand
if use_alchemical_water or not is_mutant:
map_to_neutral_atom = True
for i, ion in enumerate(extracted_ions.atom, start=1):
fep_mapping_number = 0
if is_mutant:
qion = ion.formal_charge
if abs(qion) != 1:
# this is currently enforced in system_builder_setup
raise ValueError("Only monovalent ions are allowed")
ion_name = ion.property['s_m_pdb_atom_name']
ion.property[constants.REST_HOTREGION] = 0
ion.property["i_m_residue_number"] = i
ion.property["s_m_pdb_residue_name"] = ion_name
ion.property["s_m_atom_name"] = ion_name
else:
qion = 0
# not is_mutant implies map_to_neutral_atom
# if alchemical_water, only map the atoms marked with alchemical
# charge to the mut ligand ions (these should be the oxygens for
# normal water), everything else is mapped to 0
if (map_to_neutral_atom and (not use_alchemical_water or
ion.property.get(ALCHEMICAL_WATER, 0))):
ion_map_index = int(round((i - 1) * map_to_every_nth)) + 1
fep_mapping_number = other_total + ion_map_index
x, y, z = trans_vectors + ion.xyz
# don't update water charges in non-mutant atom
if not use_alchemical_water or is_mutant:
ion.formal_charge = qion
ion.property["r_m_charge1"] = qion
ion.property["r_m_charge2"] = qion
# since we're adding this to the alchemical structures instead of into
# solvent, we add this property so others can determine if they are in
# fact part of the original alchemical structure rather than
# solvent/ion floating in solution. This is necessary for GCMC, which
# draws a region around the 'ligand'.
ion.property[ALCHEMICAL_ION] = True
ion.property[constants.FEP_MAPPING] = fep_mapping_number
for prop, coord in zip(constants.REFERENCE_COORD_PROPERTIES, [x, y, z]):
ion.property[prop] = coord
if use_alchemical_water:
del ion.property[ALCHEMICAL_WATER]
def _print_charges(q0, q1, qi, print_fcn):
print_fcn('Base ligand charge: {:+}, Mutant ligand charge: {:+}, '
'Ion st charge: {:+}'.format(q0, q1, qi))
if __name__ == '__main__':
import sys
args = sys.argv[1:]
move_ions_in_file(*args)