"""
Utilities for working with coarse grain structures
Copyright Schrodinger, LLC. All rights reserved.
"""
import base64
import itertools
import math
import random
import re
from collections import Counter
from collections import OrderedDict
from collections import namedtuple
from past.utils import old_div
import numpy
from schrodinger import structure
from schrodinger.application.desmond import cms
from schrodinger.application.matsci import msprops
from schrodinger.application.matsci import msutils
from schrodinger.infra import mm
from schrodinger.structutils import analyze
from schrodinger.structutils import color
from schrodinger.structutils import rings
from schrodinger.structutils import transform
from schrodinger.utils import fileutils
NB_TYPE_KEY = 's_matsci_cg_non-bond_type'
ANGLE_TYPE_KEY = 's_matsci_cg_angle_type'
FF_CUTOFF_PROPERTY = 'r_matsci_cg_ff_cutoff_(Ang)'
TIME_STEP_PROPERTY = 'r_matsci_cg_time_step_(fs)'
GENERAL_SITE_TYPE = 'General'
MARTINI_SITE_TYPE = 'Martini'
SITE_TYPES = [GENERAL_SITE_TYPE, MARTINI_SITE_TYPE]
HARMONIC = 'Harmonic'
ANGLES_HARMONIC_KEY = 'harm'
TRIGONOMETRIC = 'Trigonometric'
ANGLES_TRIGONOMETRIC_KEY = 'mrtn'
ANGLE_POTENTIALS = [HARMONIC, TRIGONOMETRIC]
ANGLES_KEYS = [ANGLES_HARMONIC_KEY, ANGLES_TRIGONOMETRIC_KEY]
OPLS = 'OPLS'
DIHEDRAL_POTENTIALS = [OPLS, TRIGONOMETRIC]
LENNARD_JONES = 'Lennard-Jones'
DISSIPATIVE_PARTICLE = 'Repulsive harmonic'
SHIFTED_LENNARD_JONES = 'Shifted Lennard-Jones'
NONBOND_POTENTIALS = [
LENNARD_JONES, DISSIPATIVE_PARTICLE, SHIFTED_LENNARD_JONES
]
CG_UNIQUE_BOND_PARTICLE_1_LABEL_KEY = 's_matsci_cg_particle_1_label'
CG_UNIQUE_BOND_PARTICLE_2_LABEL_KEY = 's_matsci_cg_particle_2_label'
CG_UNIQUE_PARTICLE_LABEL_KEY = 's_matsci_cg_particle_label'
# see MATSCI-3442 - use backend hidden properties
CG_CAPPING_PARTICLE_KEY = 's' + mm.M2IO_BACKEND_HIDDEN_PROP_PATTERN + 'cg_capping_particle_%s'
CG_PARTICLE_KEY = 's' + mm.M2IO_BACKEND_HIDDEN_PROP_PATTERN + 'cg_particle_%s'
CG_PARTICLE_KEY_KEY = 's_matsci_cg_particle_key'
CG_PARTICLE_PARENT_INDICES_KEY = 's_matsci_cg_particle_parent_indices'
CG_PARTICLE_PARENT_BONDS_KEY = 's_matsci_cg_particle_parent_bonds'
CG_ELEMENT_SYMBOL = 'C'
CG_ATOMIC_NUMBER = 6
DEFAULT_MMOD_ATOM_TYPE = 400
DEFAULT_RADIUS = 3. # Angstrom
DEFAULT_MASS = 12.
DEFAULT_MARTINI_MASS = 72.
DEFAULT_RGB = (100, 100, 100)
# note that the key s_matsci_cg_particle_color_name is hardcoded in
# mmshare/mmlibs/colorscheme/schemes/coarse_grain.sch
CG_PARTICLE_COLOR_NAME_KEY = 's_matsci_cg_particle_color_name'
# Header text for Desmond .msj files running CG simulations
CUTOFF = 'cutoff'
R_CLONE = 'r_clone'
FAR_TYPE = 'far_type'
MARGIN = 'margin'
KAPPA = 'kappa'
FAR_TYPE_NONE = 'none'
FAR_TYPE_PME = 'pme'
MSJ_LJ_COARSE_GRAIN_HEADER = """
#
# Universal settings for Lennard-Jones coarse grain systems
#
glue = none
coulomb_method = pme
cutoff_radius = {cutoff}
backend = {{
force.nonbonded.far.type = {far_type}
force.nonbonded.near.type = polynomial
force.nonbonded.near.average_dispersion = 0.0
force.nonbonded.near.correct_average_dispersion = false
ensemble.barostat.tau = 6.0
ensemble.thermostat.tau = 3.0
global_cell.margin = {margin}
global_cell.r_clone = {r_clone}
}}
#
# End settings for coarse grain systems
#
"""
MSJ_RH_COARSE_GRAIN_HEADER = """
#
# Universal settings for repulsive harmonic coarse grain systems
#
glue = none
coulomb_method = pme
cutoff_radius = {cutoff}
backend = {{
force.nonbonded.far.type = {far_type}
force.nonbonded.far.order = [4 4 4]
force.nonbonded.far.n_k = [144 144 144]
force.nonbonded.far.kappa = [{kappa} {kappa} {kappa}]
force.nonbonded.near.type = polynomial
force.nonbonded.near.average_dispersion = 0.0
force.nonbonded.near.correct_average_dispersion = false
global_cell.margin = {margin}
global_cell.r_clone = {r_clone}
}}
#
# End settings for coarse grain systems
#
"""
# as of MATSCI-4712 this is the same as the LJ settings except coulomb_method
MSJ_SLJ_COARSE_GRAIN_HEADER = """
#
# Universal settings for shifted Lennard-Jones coarse grain systems
#
glue = none
coulomb_method = useries
cutoff_radius = {cutoff}
backend = {{
force.nonbonded.far.type = {far_type}
force.nonbonded.near.type = polynomial
force.nonbonded.near.average_dispersion = 0.0
force.nonbonded.near.correct_average_dispersion = false
ensemble.barostat.tau = 6.0
ensemble.thermostat.tau = 3.0
global_cell.margin = {margin}
global_cell.r_clone = {r_clone}
}}
#
# End settings for coarse grain systems
#
"""
# see mmlibs/mmct/mmct_priv.h: MMCT_MAXCHARGEF, MMCT_MINCHARGEF
MAX_INFRA_FF_CHARGE = 18
MIN_INFRA_FF_CHARGE = -18
# Exclusion types for forcefield nonbonds
EXCLUDE_NONE = 'None'
EXCLUDE_BONDS = 'Bonds'
EXCLUDE_BONDS_ANGLES = 'Bonds and 1-3 Angles'
EXCLUDE_BONDS_ANGLES_DIHEDRALS = 'Bonds, 1-3 Angles, and 1-4 Dihedrals'
EXCLUDE_RINGS = 'Exclude pairs within rings and fused-ring systems'
EXCLUDE_SEPARATOR = ';'
# these markers must not overlap with valid particle name punctuation
# (see smartsutils.VALID_NAME_PUNCTUATION)
SHADOW_START = '{'
SHADOW_END = '}'
SHADOW_PATTERN = SHADOW_START + r'\d+' + SHADOW_END
RH_CUTOFF_FACTOR_DIHEDRALS = 4.0
RH_CUTOFF_FACTOR_ANGLES = 2.75
RH_CUTOFF_FACTOR_BONDS = 1.5
ParticleInfo = namedtuple('ParticleInfo',
['key', 'name', 'rgb', 'radius', 'mass'])
DEFAULT_EPSILON = 1.000
DEFAULT_A_PARAMETER = 0.1500
# Martini 2.0 data taken from Marrink, et. al., JPCB 111 7812 2007
# (see the Martini homepage for more information: http://cgmartini.nl/)
# Martini charge utils
MARTINI_CHARGE_SUBTYPE_START = 'Q'
CHIRAL_SEPARATOR = '@'
# CGFF Builder JSON keys for parameters
FORCE_CONSTANT = 'force_constant'
BOND_LENGTH = 'bond_length'
MEAN_ANGLE = 'mean_angle'
[docs]def get_base64_str_from_structure(st):
"""
Return a base64 string from the given structure.
:type st: schrodinger.structure.Structure
:param st: the structure
:rtype: str
:return: the base64 string
"""
ascii_str = st.writeToString(fileutils.MAESTRO)
ascii_bytes = ascii_str.encode()
base64_bytes = base64.b64encode(ascii_bytes)
base64_str = base64_bytes.decode()
return base64_str
[docs]def get_structure_from_base64_str(base64_str):
"""
Return a structure from the given base64 string.
:type astr: str
:param astr: the base64 string
:rtype: schrodinger.structure.Structure
:return: the structure
"""
base64_bytes = base64_str.encode()
ascii_bytes = base64.b64decode(base64_bytes)
ascii_str = ascii_bytes.decode()
return next(structure.StructureReader.fromString(ascii_str))
[docs]def is_martini_charge_type(atype):
"""
Return True if the given Martini type is a charge type.
:type atype: str
:param atype: the Martini type
:rtype: bool
:return: True if a Martini charge type
"""
return atype.startswith(MARTINI_CHARGE_SUBTYPE_START)
# Martini ring utils
MARTINI_RING_SUBTYPE_START = 'S'
[docs]def is_martini_ring_type(atype):
"""
Return True if the given Martini type is a ring type.
:type atype: str
:param atype: the Martini type
:rtype: bool
:return: True if a Martini ring type
"""
is_std_ring_type = atype.startswith(MARTINI_RING_SUBTYPE_START)
is_extra_ring_type = atype in MARTINI_EXTRA_RING_SUBTYPES
return is_std_ring_type or is_extra_ring_type
[docs]def get_martini_ring_type_from_type(atype):
"""
Return the Martini ring type from the type.
:type atype: str
:param atype: the Martini type
:rtype: str
:return: the Martini ring type
"""
if not is_martini_ring_type(atype):
return MARTINI_RING_SUBTYPE_START + atype
else:
return atype
# Martini antifreeze utils
MARTINI_ANTIFREEZE_SUBTYPE_START = 'B'
[docs]def is_martini_antifreeze_type(atype):
"""
Return True if the given Martini type is an antifreeze
type.
:type atype: str
:param atype: the Martini type
:rtype: bool
:return: True if a Martini antifreeze type
"""
return atype.startswith(MARTINI_ANTIFREEZE_SUBTYPE_START)
[docs]def get_martini_antifreeze_type_from_type(atype):
"""
Return the Martini antifreeze type from the type.
:type atype: str
:param atype: the Martini type
:rtype: str
:return: the Martini antifreeze type
"""
if not is_martini_antifreeze_type(atype):
return MARTINI_ANTIFREEZE_SUBTYPE_START + atype
else:
return atype
# Martini general utils
[docs]def get_martini_type_from_super_type(atype):
"""
Return the Martini type from the super type, i.e.
ring or antifreeze type.
:type atype: str
:param atype: the Martini super type
:rtype: str
:return: the Martini type
"""
if atype not in MARTINI_EXTRA_RING_SUBTYPES and \
(is_martini_ring_type(atype) or is_martini_antifreeze_type(atype)):
return atype[1:]
else:
return atype
[docs]def get_bead_count(struct, excluded_atom_ids=None):
"""
Gets a dictionary containing the counts of the different "beads" in the
coarse-grained structure.
:type struct: `schrodinger.structure.Structure`
:param struct: the structure to get formula from
:type excluded_atom_ids: None or a list of atom ids
:param excluded_atom_ids: exclude these atoms when computing the formula
:rtype: dict
:return: A dictionary whose keys are the "bead"/atom names and whose values
are the number of beads/atoms
"""
excluded_atom_ids = set() if excluded_atom_ids is None else set(
excluded_atom_ids)
counter = Counter(atom.name
for atom in struct.atom
if atom.index not in excluded_atom_ids)
bead_counter = dict(counter)
return bead_counter
# Martini basic types
MARTINI_QDA_SUBTYPE = 'Qda'
MARTINI_QD_SUBTYPE = 'Qd'
MARTINI_QA_SUBTYPE = 'Qa'
MARTINI_Q0_SUBTYPE = 'Q0'
MARTINI_P5_SUBTYPE = 'P5'
MARTINI_P4_SUBTYPE = 'P4'
MARTINI_P3_SUBTYPE = 'P3'
MARTINI_P2_SUBTYPE = 'P2'
MARTINI_P1_SUBTYPE = 'P1'
MARTINI_NDA_SUBTYPE = 'Nda'
MARTINI_ND_SUBTYPE = 'Nd'
MARTINI_NA_SUBTYPE = 'Na'
MARTINI_N0_SUBTYPE = 'N0'
MARTINI_C5_SUBTYPE = 'C5'
MARTINI_C4_SUBTYPE = 'C4'
MARTINI_C3_SUBTYPE = 'C3'
MARTINI_C2_SUBTYPE = 'C2'
MARTINI_C1_SUBTYPE = 'C1'
MARTINI_BASIC_SUBTYPES = [
MARTINI_QDA_SUBTYPE, MARTINI_QD_SUBTYPE, MARTINI_QA_SUBTYPE,
MARTINI_Q0_SUBTYPE, MARTINI_P5_SUBTYPE, MARTINI_P4_SUBTYPE,
MARTINI_P3_SUBTYPE, MARTINI_P2_SUBTYPE, MARTINI_P1_SUBTYPE,
MARTINI_NDA_SUBTYPE, MARTINI_ND_SUBTYPE, MARTINI_NA_SUBTYPE,
MARTINI_N0_SUBTYPE, MARTINI_C5_SUBTYPE, MARTINI_C4_SUBTYPE,
MARTINI_C3_SUBTYPE, MARTINI_C2_SUBTYPE, MARTINI_C1_SUBTYPE
]
# Martini extra ring types
MARTINI_CNP_SUBTYPE = 'CNP'
MARTINI_EXTRA_RING_SUBTYPES = [MARTINI_CNP_SUBTYPE]
# Martini ring types
MARTINI_RING_SUBTYPES = list(
map(get_martini_ring_type_from_type, MARTINI_BASIC_SUBTYPES))
MARTINI_SQDA_SUBTYPE = MARTINI_RING_SUBTYPES[0]
MARTINI_SQD_SUBTYPE = MARTINI_RING_SUBTYPES[1]
MARTINI_SQA_SUBTYPE = MARTINI_RING_SUBTYPES[2]
MARTINI_SQ0_SUBTYPE = MARTINI_RING_SUBTYPES[3]
MARTINI_SP5_SUBTYPE = MARTINI_RING_SUBTYPES[4]
MARTINI_SP4_SUBTYPE = MARTINI_RING_SUBTYPES[5]
MARTINI_SP3_SUBTYPE = MARTINI_RING_SUBTYPES[6]
MARTINI_SP2_SUBTYPE = MARTINI_RING_SUBTYPES[7]
MARTINI_SP1_SUBTYPE = MARTINI_RING_SUBTYPES[8]
MARTINI_SNDA_SUBTYPE = MARTINI_RING_SUBTYPES[9]
MARTINI_SND_SUBTYPE = MARTINI_RING_SUBTYPES[10]
MARTINI_SNA_SUBTYPE = MARTINI_RING_SUBTYPES[11]
MARTINI_SN0_SUBTYPE = MARTINI_RING_SUBTYPES[12]
MARTINI_SC5_SUBTYPE = MARTINI_RING_SUBTYPES[13]
MARTINI_SC4_SUBTYPE = MARTINI_RING_SUBTYPES[14]
MARTINI_SC3_SUBTYPE = MARTINI_RING_SUBTYPES[15]
MARTINI_SC2_SUBTYPE = MARTINI_RING_SUBTYPES[16]
MARTINI_SC1_SUBTYPE = MARTINI_RING_SUBTYPES[17]
# Martini extra types
MARTINI_BP4_SUBTYPE = get_martini_antifreeze_type_from_type(MARTINI_P4_SUBTYPE)
MARTINI_CUSTOM_SUBTYPE = 'Custom'
DEFAULT_MARTINI_SUBTYPE = MARTINI_C1_SUBTYPE
# Martini all types
MARTINI_SUBTYPES = MARTINI_BASIC_SUBTYPES + MARTINI_RING_SUBTYPES + \
[MARTINI_BP4_SUBTYPE, MARTINI_CNP_SUBTYPE, MARTINI_CUSTOM_SUBTYPE]
# Martini epsilon parameters in kcal/mol
MARTINI_O_EPSILON = 1.338
MARTINI_I_EPSILON = 1.195
MARTINI_II_EPSILON = 1.076
MARTINI_III_EPSILON = 0.956
MARTINI_IV_EPSILON = 0.837
MARTINI_V_EPSILON = 0.741
MARTINI_VI_EPSILON = 0.645
MARTINI_VII_EPSILON = 0.550
MARTINI_VIII_EPSILON = 0.478
MARTINI_IX_EPSILON = 0.478
# yapf: disable
MARTINI_EPSILONS = {
(MARTINI_QDA_SUBTYPE, MARTINI_QDA_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_QD_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_QA_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_Q0_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_P5_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_P4_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_QDA_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_QD_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_QA_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_Q0_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_P5_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_P4_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_QA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_Q0_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_P5_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_P4_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_Q0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_P5_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_P4_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IX_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_P5_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_P4_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_O_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VIII_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_P4_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VIII_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_P3_SUBTYPE): MARTINI_I_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_P2_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VII_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_P1_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P1_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_NDA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_NDA_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_ND_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_II_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_ND_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_NA_SUBTYPE): MARTINI_III_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_NA_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_N0_SUBTYPE, MARTINI_N0_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_N0_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_N0_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_N0_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_N0_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_N0_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_C5_SUBTYPE, MARTINI_C5_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C5_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C5_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C5_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_C5_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_C4_SUBTYPE, MARTINI_C4_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C4_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C4_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_C4_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_C3_SUBTYPE, MARTINI_C3_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C3_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C3_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C2_SUBTYPE, MARTINI_C2_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C2_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C1_SUBTYPE, MARTINI_C1_SUBTYPE): MARTINI_IV_EPSILON
}
# these are mostly one-offs not worth classifying, doesn't include
# all possible pairs
MARTINI_EXTRA_EPSILONS = {
(MARTINI_CNP_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P5_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P4_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_P3_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_V_EPSILON,
(MARTINI_P2_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.777,
(MARTINI_P1_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.777,
(MARTINI_NDA_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.789,
(MARTINI_ND_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.777,
(MARTINI_NA_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.777,
(MARTINI_N0_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.813,
(MARTINI_C5_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.813,
(MARTINI_SC5_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.627,
(MARTINI_C4_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_SC4_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.645,
(MARTINI_C3_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_IV_EPSILON,
(MARTINI_C2_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.789,
(MARTINI_C1_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.753,
(MARTINI_SC1_SUBTYPE, MARTINI_CNP_SUBTYPE): 0.627,
(MARTINI_QDA_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_QD_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_QA_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON,
(MARTINI_Q0_SUBTYPE, MARTINI_CNP_SUBTYPE): MARTINI_VI_EPSILON
}
MARTINI_EPSILONS.update(MARTINI_EXTRA_EPSILONS)
# yapf: enable
MARTINI_EPSILONS = dict([
(tuple(sorted(k)), v) for k, v in MARTINI_EPSILONS.items()
])
MARTINI_RING_EPSILON_FACTOR = 0.75
# Martini sigma parameters in Ang.
# for the Martini sigmas:
# (1) is a general Martini sigma
# (2) is a Qx-C1,C2 x = da, a, d, or 0 Martini sigma
# (3) is a ring-ring (aromatic or not) Martini sigma
# (4) is an antifreeze BP4-P4 Martini sigma
# (5) is a Polystyrene Martini sigma (modeled with
# slightly smaller than normal particles) (see MATSCI-4889)
MARTINI_1_SIGMA = 4.7
MARTINI_2_SIGMA = 6.2
MARTINI_3_SIGMA = 4.3
MARTINI_4_SIGMA = 5.7
MARTINI_5_SIGMA = 4.1
MARTINI_6_SIGMA = 4.5
MARTINI_SIGMAS = [
MARTINI_5_SIGMA, MARTINI_3_SIGMA, MARTINI_6_SIGMA, MARTINI_1_SIGMA,
MARTINI_4_SIGMA, MARTINI_2_SIGMA
]
# in units of Angstrom
MARTINI_DEFAULT_NONRING_BOND_EQ_LENGTH = 4.7
MARTINI_DEFAULT_RING_BOND_EQ_LENGTH = 4.3
# in units of (kcal/mol)/Ang.^2
MARTINI_DEFAULT_NONRING_BOND_FORCE_CONSTANT = 2.988
MARTINI_DEFAULT_RING_BOND_FORCE_CONSTANT = 11.950
InternalInfo = namedtuple('InternalInfo', ['names', 'indices', 'data'])
[docs]def split_shadow_bond_tag(name):
"""
Split the shadow bond tag from the given name.
:type name: str
:param name: the particle name
:rtype: str, str
:return: the given name less the shadow bond tag, the shadow
bond tag (empty string if not present)
"""
match = re.search(SHADOW_PATTERN, name)
if match:
shadow_tag = match.group(0)
else:
shadow_tag = ''
return name.replace(shadow_tag, ''), shadow_tag
[docs]def get_martini_subtypes(type_tuple, martini_subtypes):
"""
Return a list of Martini subtypes from the given
type tuple.
:type type_tuple: tuple
:param type_tuple: a type tuple, for example bond type
tuple, containing particle names
:type martini_subtypes: dict
:param martini_subtypes: keys are site type tuples and values
are Martini subtypes
:rtype: list
:return: contains Martini subtypes
"""
# MATSCI-5191 - strip shadow bond tags
f_martini_subtypes = []
for atype in type_tuple:
base, tag = split_shadow_bond_tag(atype)
try:
martini_subtype = martini_subtypes[(base,)]
except KeyError:
# see MATSCI-5763, this either means that the given type_tuple
# was uniquified in the CG particulate GUI or a real exception,
# if the former then because the uniquifying labels are
# user specified and only stored in the CT bonds loop over
# site names as it is inexpensive, do it this way now as
# opposed to more directly earlier in the code so that potential
# future uses with angles, etc. are covered
for atuple in martini_subtypes:
site_type = atuple[0]
if site_type in base:
base = site_type
break
else:
raise
martini_subtype = martini_subtypes[(base,)]
f_martini_subtypes.append(martini_subtype)
return f_martini_subtypes
[docs]def get_martini_eq_length(pair_type, martini_subtypes):
"""
Return the Martini equilibrium length parameter.
:type pair_type: tuple
:param pair_type: the bond pair type
:type martini_subtypes: dict
:param martini_subtypes: keys are site type tuples and values
are Martini subtypes
:rtype: float
:return: the equilibrium length parameter in Angstrom
"""
types = get_martini_subtypes(pair_type, martini_subtypes)
both_are_ring_types = all(map(is_martini_ring_type, types))
if both_are_ring_types:
return MARTINI_DEFAULT_RING_BOND_EQ_LENGTH
else:
return MARTINI_DEFAULT_NONRING_BOND_EQ_LENGTH
[docs]def get_martini_bond_force_constant(pair_type, martini_subtypes):
"""
Return the Martini bond force constant parameter.
:type pair_type: tuple
:param pair_type: the bond pair type
:type martini_subtypes: dict
:param martini_subtypes: keys are site type tuples and values
are Martini subtypes
:rtype: float
:return: the bond force constant parameter in (kcal/mol)/Ang.^2
"""
types = get_martini_subtypes(pair_type, martini_subtypes)
both_are_ring_types = all(map(is_martini_ring_type, types))
if both_are_ring_types:
return MARTINI_DEFAULT_RING_BOND_FORCE_CONSTANT
else:
return MARTINI_DEFAULT_NONRING_BOND_FORCE_CONSTANT
[docs]def get_martini_epsilon(pair_type, martini_subtypes):
"""
Return the Martini epsilon parameter.
:type pair_type: tuple
:param pair_type: the nonbond pair type
:type martini_subtypes: dict
:param martini_subtypes: keys are site type tuples and values
are Martini subtypes
:rtype: float
:return: the epsilon parameter in kcal/mol
"""
types = tuple(sorted(martini_subtypes[(i,)] for i in pair_type))
if MARTINI_CUSTOM_SUBTYPE in types:
return DEFAULT_EPSILON
elif MARTINI_BP4_SUBTYPE in types and MARTINI_P4_SUBTYPE in types:
return MARTINI_O_EPSILON
epsilon = MARTINI_EPSILONS.get(types)
if epsilon is not None:
return epsilon
both_are_ring_types = all(map(is_martini_ring_type, types))
types = tuple(sorted(map(get_martini_type_from_super_type, types)))
epsilon = MARTINI_EPSILONS.get(types, DEFAULT_EPSILON)
if both_are_ring_types:
epsilon *= MARTINI_RING_EPSILON_FACTOR
return epsilon
[docs]def get_martini_sigma(pair_type, martini_subtypes):
"""
Return the Martini sigma parameter.
:type pair_type: tuple
:param pair_type: the nonbond pair type
:type martini_subtypes: dict
:param martini_subtypes: keys are site type tuples and values
are Martini subtypes
:rtype: float
:return: the sigma parameter in Ang.
"""
types = [martini_subtypes[(i,)] for i in pair_type]
if MARTINI_CUSTOM_SUBTYPE in types:
return MARTINI_1_SIGMA
elif MARTINI_BP4_SUBTYPE in types and MARTINI_P4_SUBTYPE in types:
return MARTINI_4_SIGMA
elif all(map(is_martini_ring_type, types)):
return MARTINI_3_SIGMA
elif any(map(is_martini_charge_type, types)) and \
(MARTINI_C1_SUBTYPE in types or MARTINI_C2_SUBTYPE in types):
return MARTINI_2_SIGMA
else:
return MARTINI_1_SIGMA
[docs]def get_combined_exclusion_str(bad_str, exclude_rings=False):
"""
Return a combined exclusion string.
:type bad_str: str
:param bad_str: a bond, angle, or dihedral string, i.e. one
of EXCLUDE_BONDS, EXCLUDE_BONDS_ANGLES, or
EXCLUDE_BONDS_ANGLES_DIHEDRALS
:type exclude_rings: bool
:param exclude_rings: whether pairs in rings and fused-ring
systems are also excluded
:rtype: str
:return: the combined exclusion string
"""
if exclude_rings:
return EXCLUDE_SEPARATOR.join([bad_str, EXCLUDE_RINGS])
else:
return bad_str
[docs]def parse_combined_exclusion_str(combined_str):
"""
Parse a combined exclusion string.
:type combined_str: str
:param combined_str: a combined exclusion string obtained
from get_combined_exclusion_str
:rtype: str, bool
:return: a bond, angle, or dihedral string, i.e. one
of EXCLUDE_BONDS, EXCLUDE_BONDS_ANGLES, or
EXCLUDE_BONDS_ANGLES_DIHEDRALS, and whether pairs in
rings and fused-ring systems are also excluded
"""
try:
bad_str, ring_str = combined_str.split(EXCLUDE_SEPARATOR)
except ValueError:
bad_str, ring_str = combined_str, ''
return bad_str, bool(ring_str)
[docs]def get_type_dict_from_st(st, type_dict=None):
"""
Return a sorted type dict for the given structure.
:type st: `schrodinger.structure.Structure`
:param st: the structure from which to obtain the type dict
:type type_dict: OrderedDict or None
:param type_dict: specify a type dict to include in the
returned type dict (note that this type dict takes precedence over
the type dict data for the given structure) or None if there isn't one
:rtype: OrderedDict
:return: the sorted type dict, keys are particle names, values are
ParticleInfo
"""
if type_dict is None:
f_type_dict = OrderedDict()
else:
f_type_dict = type_dict.copy()
for atom in st.atom:
if atom.name not in f_type_dict:
key = atom.name
name = atom.name
rgb = atom.color.rgb
radius = mm.mmct_atom_get_display_radius(st, atom.index)
mass = mm.mmct_atom_get_coarse_grain_mass(st, atom.index)
info = ParticleInfo(key=key,
name=name,
rgb=rgb,
radius=radius,
mass=mass)
f_type_dict[name] = info
f_type_dict = OrderedDict(sorted(f_type_dict.items()))
return f_type_dict
def _get_rh_cg_cutoff_factor(model):
"""
Get the scale factor for the cutoff when computing rclone in repulsive
harmonic CG models. The scale factor is based on whether there are FF terms
for dihedrals or angles.
:type model: `schrodinger.application.desmond.cms.Cms`
:param model: The CMS structure
:rtype: float
:return: The scale factor for the cutoff term
"""
try:
structs = model.comp_ct[:]
except AttributeError:
# model is just a Structure, probably, but could be an FFIOStructure
structs = [model]
structs = [x for x in structs if hasattr(x, 'ffio')]
for struct in structs:
if len(struct.ffio.dihedral) > 0:
return RH_CUTOFF_FACTOR_DIHEDRALS
elif len(struct.ffio.angle) > 0:
return RH_CUTOFF_FACTOR_ANGLES
return RH_CUTOFF_FACTOR_BONDS
[docs]def system_has_charged_particles(model):
"""
Check if any particle in the given system has a charge based on the sites
FFIO block
:type model: `schrodinger.application.desmond.cms.Cms` or
`schrodinger.application.desmond.ffiostructure.FFIOStructure`
:param model: The system to check the FFIO site block for charges
:rtype: bool
:return: True if an FFIO sites block exists in the given system and contains
at least one charged particle. False if there is either no such block or the
block contains only neutral particles
"""
if isinstance(model, cms.Cms):
structs = model.comp_ct[:]
else:
structs = [model]
for struct in structs:
try:
sites = struct.ffio.site
except AttributeError:
# Not an FFIOStructure object with a sites block
continue
for site in sites:
if site.charge != 0:
return True
return False
[docs]def get_cg_cutoff_data(model):
"""
Return the coarse grain cutoff data.
:type model: `schrodinger.application.desmond.cms.Cms`
:param model: The model structure
:raise ValueError: if model is incorrect
:rtype: dict
:return: the coarse grain cutoff data, each key in the dictionary
is a format specifier from the MSJ_*_COARSE_GRAIN_HEADER string
and the value is the value for that format specifier
"""
if not msutils.is_coarse_grain(model):
msg = 'Model must be a coarse-grained model.'
raise ValueError(msg)
cutoff = model.property.get(FF_CUTOFF_PROPERTY)
if cutoff is None:
msg = 'Cutoff must be defined as a structure property.'
raise ValueError(msg)
# as of MATSCI-4712 for shifted LJ use margin,
# r_clone, and r_clone_dist_modifier from LJ
far_type = FAR_TYPE_NONE
if is_repulsive_harmonic(model):
margin = cutoff
r_clone = cutoff * _get_rh_cg_cutoff_factor(model)
r_clone_dist_modifier = 0.5 * cutoff
if system_has_charged_particles(model):
far_type = FAR_TYPE_PME
elif is_shifted_lennard_jones(model):
margin = 2.00
r_clone = 3.00 + old_div(float(cutoff), 2)
r_clone_dist_modifier = 0.5
else:
margin = 2.00
r_clone = 3.00 + old_div(float(cutoff), 2)
r_clone_dist_modifier = 0.5
if system_has_charged_particles(model):
far_type = FAR_TYPE_PME
max_dist = compute_max_internal_distance(model)
if r_clone <= max_dist:
r_clone = max_dist + r_clone_dist_modifier
data = {}
data[CUTOFF] = cutoff
data[R_CLONE] = r_clone
data[FAR_TYPE] = far_type
data[MARGIN] = margin
if is_repulsive_harmonic(model):
epsilon = 1E-9
toler = math.sqrt(abs(math.log(epsilon * cutoff)))
sigma_ewald = 2 * cutoff / math.sqrt(
abs(math.log(epsilon * cutoff * toler)))
kappa = 1 / (math.sqrt(6) * sigma_ewald)
data[KAPPA] = kappa
return data
[docs]def compute_max_internal_distance(struct):
"""
Compute the largest distance between any pair of atoms that are a) bonded,
b) 1-3 in an angle or c) 1-4 in a dihedral
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to compute the distance over
:rtype: float
:return: The largest distance between any pair of atoms involved in a bond,
angle or torsion
"""
graph = analyze.create_nx_graph(struct)
pairs = set()
for genner in (analyze.bond_iterator, analyze.angle_iterator,
analyze.torsion_iterator):
for indexes in genner(struct, nx_graph=graph):
pairs.add((indexes[0], indexes[-1]))
if not pairs:
return 0.0
else:
return max([struct.measure(*x) for x in pairs])
[docs]def validate_no_cg(struct=None, structs=None, rows=None):
"""
Fail a validation check if any of the structures is a coarse grain structure
One and only one of struct, structs or rows should be given
:type struct: `schrodinger.structure.Structure`
:param struct: The structure, if only a single structure is to be tested
:type structs: list
:param structs: A list of structures to test
:type rows: iterator of `schrodinger.project.ProjectRow`
:param rows: Project Table rows to check
:rtype: (False, msg) or True
:return: False and error message if the structure is coarse grain, True if
everything is OK. Return value is acceptable for an AF2 valiator.
"""
if sum([1 for x in [struct, structs, rows] if x]) != 1:
raise RuntimeError('One and only one of struct, structs and rows must '
'be given')
msg = '%s is a coarse-grained structure and cannot be used with this panel'
if struct:
structs = [struct]
elif rows:
structs = (x.getStructure() for x in rows)
for astruct in structs:
# We do by_atom=True because we don't know the origin of the structure
if msutils.is_coarse_grain(astruct, by_atom=True):
title = astruct.title or 'Input'
return (False, msg % title)
return True
[docs]def set_as_coarse_grain(struct):
"""
Mark struct as a coarse grain structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to mark
"""
struct.property[msprops.REP_TYPE_KEY] = msprops.CG_REP_TYPE
[docs]def set_nonbond_potential_type(struct, nbtype):
"""
Mark the structure with the non-bond pontential type
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to mark
:type nbtype: str
:param nbtype: Should be a module constant in the NONBOND_POTENTIALS list
:raise ValueError: if nbtype is not in NONBOND_POTENTIALS
:raise TypeError: if struct is not a coarse grain system
"""
if nbtype not in NONBOND_POTENTIALS:
raise ValueError('Non-bond potential type must be one of the types in '
'the NONBOND_POTENTIALS list')
if not msutils.is_coarse_grain(struct):
raise TypeError('Non-bond potential type can only be set on coarse '
'grain structures')
struct.property[NB_TYPE_KEY] = nbtype
[docs]def set_angle_potential_type(struct, atype):
"""
Mark the structure with the angle pontential type
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to mark
:type atype: str
:param atype: Should be a module constant in the ANGLES_KEYS list
:raise ValueError: if atype is not in ANGLES_KEYS
:raise TypeError: if struct is not a coarse grain system
"""
if atype not in ANGLES_KEYS:
raise ValueError('Angle potential type must be one of the types in '
'the ANGLES_KEYS list')
if not msutils.is_coarse_grain(struct):
raise TypeError('Angle potential type can only be set on coarse '
'grain structures')
struct.property[ANGLE_TYPE_KEY] = atype
[docs]def get_nonbond_potential_type(struct):
"""
Return the type of non-bond potential for this structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: str or None
:return: The value of the NB_TYPE_KEY property on the structure, or None if
the structure is not a coarse grain structure or does not have this property
set
"""
if not msutils.is_coarse_grain(struct):
return None
return struct.property.get(NB_TYPE_KEY)
[docs]def get_angle_potential_type(struct):
"""
Return the type of angle potential for this structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: str or None
:return: The value of the ANGLE_TYPE_KEY property on the structure, or None if
the structure is not a coarse grain structure or does not have this property
set
"""
if not msutils.is_coarse_grain(struct):
return None
return struct.property.get(ANGLE_TYPE_KEY)
[docs]def is_lennard_jones(struct):
"""
Check if this is a Lennard-Jones coarse grain structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: bool
:return: True if this is a coarse grain Lennard-Jones structure, False if
not
"""
return get_nonbond_potential_type(struct) == LENNARD_JONES
[docs]def is_repulsive_harmonic(struct):
"""
Check if this is a repulsive harmonic coarse grain structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: bool
:return: True if this is a coarse grain repulsive harmoic structure, False
if not
"""
return get_nonbond_potential_type(struct) == DISSIPATIVE_PARTICLE
[docs]def is_shifted_lennard_jones(struct):
"""
Check if this is a shifted Lennard-Jones coarse grain structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: bool
:return: True if this is a coarse grain shifted Lennard-Jones structure, False if
not
"""
return get_nonbond_potential_type(struct) == SHIFTED_LENNARD_JONES
[docs]def is_harmonic_angle(struct):
"""
Check if this is a harmonic angle coarse grain structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: bool
:return: True if this is a coarse grain harmonic angle structure, False if
not
"""
return get_angle_potential_type(struct) == ANGLES_HARMONIC_KEY
[docs]def is_trigonometric_angle(struct):
"""
Check if this is a trigonometric angle coarse grain structure
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:rtype: bool
:return: True if this is a coarse grain trigonometric angle structure, False if
not
"""
return get_angle_potential_type(struct) == ANGLES_TRIGONOMETRIC_KEY
[docs]def set_atom_coarse_grain_properties(struct,
atom,
name,
rgb=DEFAULT_RGB,
atom_type=DEFAULT_MMOD_ATOM_TYPE,
formal_charge=0.0,
partial_charge=0.0,
radius=DEFAULT_RADIUS,
mass=DEFAULT_MASS):
"""
Set the properties required for a coarse grain particle atom
:type structure: `schrodinger.structure.Structure`
:param structure: The structure containing the atom
:type atom: `schrodinger.structure._StructureAtom`
:param atom: The atom to set properties on
:type name: str
:param name: The name of the coarse grain particle.
:type rgb: tuple
:param rgb: The 0-255 tuple of red, green and blue that defines the atom
color
:type atom_type: int
:param atom_type: The mmod atom type
:type formal_charge: int
:param formal_charge: The formal charge
:type partial_charge: float
:param partial_charge: The partial charge
:type radius: float
:param radius: The particle radius, in Angstrom
:type mass: float
:param mass: The particle mass
"""
atom.element = CG_ELEMENT_SYMBOL
atcolor = color.Color(rgb)
atom.color = atcolor
atom.property[CG_PARTICLE_COLOR_NAME_KEY] = atcolor.name
atom.name = name
atom.atomic_number = CG_ATOMIC_NUMBER
# Set atom type after atomic number, because atomic number retypes
atom.atom_type = atom_type
atom.formal_charge = max(min(formal_charge, MAX_INFRA_FF_CHARGE),
MIN_INFRA_FF_CHARGE)
atom.partial_charge = partial_charge
index = atom.index
mm.mmct_atom_set_display_radius(struct, index, radius)
# note that mmct_atom_set_coarse_grain_particle must be called
# after setting the atom.atom_type or warnings regarding
# the calling of mmct_atom_set_isotope and mmct_atom_set_atomic_number
# on atoms marked as coarse grain particles will be raised, also
# mmct_atom_set_coarse_grain_particle must be called before
# mmct_atom_set_coarse_grain_mass
mm.mmct_atom_set_coarse_grain_particle(struct, index, True)
mm.mmct_atom_set_coarse_grain_mass(struct, index, mass)
[docs]def get_cg_molecule_hash(mol):
"""
Return a hash for the given CG molecule.
:type mol: schrodinger.structure._Molecule
:param mol: the molecule
:rtype: str
:return: the hash
"""
# the disordered system builder marks molecules by marking
# their atoms with the following pdb property, yet MATSCI-5775
# reveals that this alone is insufficient and so prepend
# an ordered name hash
pdb_res_names = set()
atom_name_hash = ''
for atom in mol.atom:
pdb_res_names.add(atom.property.get('s_m_pdb_residue_name', '').strip())
atom_name_hash += atom.name
# We need to account for bonding of particles to ensure we don't have an
# isomer. MATSCI-11290
atom_name_hash += ".".join(
str(x.number_by_molecule) for x in atom.bonded_atoms)
pdb_res_names = ''.join(sorted(pdb_res_names))
return atom_name_hash + '_' + pdb_res_names
[docs]def create_cg_molecule_name(mol):
"""
Create a name for the passed molecule by getting the count of each type of
bead in it
The format is "A:4, B:5, ..." where A and B are bead names and 4 and 5 are
the counts of each in the molecule
:param `structure._Molecule` mol: The molecule to get the name for
:rtype: str
:return: The molecule name
"""
bead_names = [atom.name for atom in mol.atom]
counts_dict = Counter(bead_names)
parts = []
for key in sorted(counts_dict):
parts.append(f"{key}: {counts_dict[key]}")
return ', '.join(parts)
[docs]def collect_unique_structures(st, hash_getter=None):
"""
Using the given structure and hasher collect unique instances
of molecule structures as the keys of a dictionary with their
values being lists of tuples of atom indices for all molecule
instances.
:type st: schrodinger.structure.Structure
:param st: the structure containing the molecules
:type hash_getter: function
:param hash_getter: the function used to obtain a unique
hash for a given molecule, takes schrodinger.structure._Molecule,
should return None if a hash can not be created
:rtype: dict
:return: keys are schrodinger.structure.Structure, values
are lists of atom index tuples
"""
if hash_getter is None:
hash_getter = get_cg_molecule_hash
# first collect unique structures by hash
hash_dict = {}
for mol in st.molecule:
ahash = hash_getter(mol)
if ahash is None:
ahash = str(mol.number)
idxs = mol.getAtomIndices()
hash_dict.setdefault(ahash, []).append(idxs)
# instead of storing hashes in keys store a unique
# instance of the structure object
st_dict = {}
for idxs in hash_dict.values():
first_mol_first_idx = idxs[0][0]
first_mol = st.atom[first_mol_first_idx].getMolecule()
st_dict[first_mol.extractStructure()] = idxs
return st_dict
[docs]def get_internal_info(atoms, data_getter=None, uniqueify=False):
"""
Return information about the given internal.
:type atoms: list
:param atoms: contains schrodinger.structure._StructureAtom
that define the internal
:type data_getter: function or None
:param data_getter: if given used to obtain data for the given internal,
takes a list of schrodinger.struture._StructureAtom and returns a
data tuple
:type uniqueify: bool
:param uniqueify: if True then relevant internal
coordinates will be uniqueified by type
:rtype: InternalInfo
:return: the internal info
"""
if data_getter is None:
data_getter = lambda x: None
atoms = order_internal_by_name(atoms)
names = get_names(atoms, uniqueify=uniqueify)
indices = tuple(atom.index for atom in atoms)
data = data_getter(atoms)
return InternalInfo(names=names, indices=indices, data=data)
[docs]def internals_iterator(sts_idxs_dict,
idx_getter,
data_getter=None,
uniqueify=False):
"""
Iterator that yields InternalInfo for each internal in the system
defined by the given dictionary of unique molecular structures and
atom indices of all molecule instances. Internal indices and data
are obtained with the given getters.
:type sts_idxs_dict: dict
:param sts_idxs_dict: keys are schrodinger.structure.Structure,
values are lists of atom index tuples, together they define
the system
:type idx_getter: function
:param idx_getter: gets indices of internals, takes a
schrodinger.structure.Structure as kwarg 'struct' and
yields index tuples
:type data_getter: function or None
:param data_getter: gets data of internals, takes a list of
schrodinger.struture._StructureAtom and returns a data tuple
:type uniqueify: bool
:param uniqueify: if True then relevant internal
coordinates will be uniqueified by type
:rtype: InternalInfo
:return: each iteration yields an InternalInfo
"""
if data_getter is None:
data_getter = lambda x: None
for st, idxs in sts_idxs_dict.items():
# define internal info for the given unique structure
t_idx_getter = idx_getter(struct=st)
if not t_idx_getter:
continue
for internal_idxs in t_idx_getter:
internal_atoms = [st.atom[i] for i in internal_idxs]
info = get_internal_info(internal_atoms,
data_getter=data_getter,
uniqueify=uniqueify)
# copy internal info to all structure instances, the following
# maps indices of the unique molecule (keys) to indices of all of
# its instances (values)
for jdxs in idxs:
amap = dict(zip(list(range(1, len(jdxs) + 1)), jdxs))
kdxs = tuple(amap[i] for i in info.indices)
yield InternalInfo(names=info.names,
indices=kdxs,
data=info.data)
return ()
[docs]def order_internal_by_name(atoms):
"""
Order the atoms of the given internal coordinate by name.
:type atoms: list
:param atoms: contains structure._StructureAtom
of the internal coordinate to be ordered
:rtype: list
:return: atoms sorted by name
"""
def wrong_order(at1, at2):
# Atoms should be ordered alphabetically by type or numerically by
# index if the type is identical
return (at1.name > at2.name or
(at1.name == at2.name and at1.index > at2.index))
if wrong_order(atoms[0], atoms[-1]):
if atoms[0].name == atoms[-1].name and len(atoms) == 4:
# Order dihedrals with same end types by the middle atoms
if wrong_order(atoms[1], atoms[2]):
atoms.reverse()
else:
atoms.reverse()
return atoms
[docs]def get_names(atoms, uniqueify=False):
"""
Return the atom names for the given list of atoms.
:type atoms: list
:param atoms: contains structure._StructureAtom
:type uniqueify: bool
:param uniqueify: if True then relevant internal
coordinates will be uniqueified by type
:rtype: tuple
:return: contains atom names
"""
# collect bonding pairs and initialize labels (creating two extra
# labels eases the final labels build), the number of pairs is
# one less than the number of atoms and the number of labels is
# two more than the number of pairs (the number of labels is one
# more than the number of atoms)
pairs = [(atoms[i], atoms[i + 1]) for i in range(len(atoms) - 1)]
labels = [('', '')] * (len(pairs) + 2)
# update labels if we are uniqueifying, the bond properties are
# sorted in the same way as order_internal_by_name does by
# particle name so the ordering may need to be reversed
if uniqueify:
# Labels tatic bonds set in the particulate gui
bond_1_key = CG_UNIQUE_BOND_PARTICLE_1_LABEL_KEY
bond_2_key = CG_UNIQUE_BOND_PARTICLE_2_LABEL_KEY
# Labels atoms that have a "phantom" bond assign via the FF gui
atom_key = CG_UNIQUE_PARTICLE_LABEL_KEY
for idx, pair in enumerate(pairs):
st = pair[0].structure
if st.areBound(*pair):
bond = st.getBond(*pair)
i_label = bond.property.get(bond_1_key, '')
j_label = bond.property.get(bond_2_key, '')
else:
i_label = pair[0].property.get(atom_key, '')
j_label = pair[1].property.get(atom_key, '')
if i_label and j_label:
alist = [i_label, j_label]
if list(pair) != order_internal_by_name(list(pair)):
alist.reverse()
labels[idx + 1] = alist
# do the final labels build, convert these bond labels to atom labels,
# the number of labels is one more than the number of atoms
labels = [
(labels[idx][1], labels[idx + 1][0]) for idx in range(len(labels) - 1)
]
# create final names
get_name = lambda x, y: y[0] + x.name + y[1]
names = tuple(get_name(*pair) for pair in zip(atoms, labels))
return names
[docs]def site_iterator(struct=None):
"""
Create a site iterator from the given structure.
:type struct: schrodinger.structure.Structure
:param struct: the structure
:rtype: tuple
:return: each iteration yields a tuple containing a single index
"""
if not struct:
raise ValueError('structure must be given')
for atom in struct.atom:
yield (atom.index,)
[docs]def fused_ring_nonbond_iterator(struct=None):
"""
Create a fused ring nonbond iterator from the given structure.
:type struct: schrodinger.structure.Structure
:param struct: the structure
:rtype: tuple
:return: each iteration yields a nonbond index pair tuple
"""
if not struct:
raise ValueError('structure must be given')
fused_ring_idxs = set()
for fused_ring in FusedRing.getFusedRings(struct):
fused_ring_idxs.add(tuple(sorted(fused_ring.getAtomIndices())))
for idxs in fused_ring_idxs:
return itertools.combinations(idxs, 2)
[docs]class Internals(object):
"""
Manage the internal coordinates of a structure.
"""
[docs] def __init__(self, astructure, uniqueify=False):
"""
Create an instance.
:type astructure: structure.Structure
:param astructure: the structure for which the internal
coordinates are needed
:type uniqueify: bool
:param uniqueify: if True then relevant internal
coordinates will be uniqueified by type
"""
self.astructure = astructure
self.uniqueify = uniqueify
self.sts_idxs_dict = collect_unique_structures(self.astructure)
self.sites = OrderedDict()
self.bonds = OrderedDict()
self.angles = None
self.dihedrals = None
self.impropers = None
self.nonbonds = OrderedDict()
self.bond_nonbonds = None
self.angle_nonbonds = None
self.dihedral_nonbonds = None
self.fused_ring_nonbonds = None
# these only contain data for fake bonds, angles, and dihedrals
# resulting from the user adding bonds in the CG FF assignment gui
# DefineNewBondDialog
self.newbond_nonbonds = {}
self.newangle_nonbonds = {}
self.newdihedral_nonbonds = {}
# This points to the current included nonbonds data
self.included_nonbonds = self.nonbonds
# This points to all non-bonds that are currently excluded
self.excluded_nonbonds = None
self._removeUniqueParticleLabels()
[docs] def avgInternals(self, internals):
"""
Average the values of the internals in the given
dictionary of internals.
:type internals: dict
:param internals: keys are tuples of names of the constitutive
atoms, values is a list of pair tuples containing (1) a tuple
of atom indices and (2) a tuple of internal values
:rtype: dict
:return: keys are tuples of names of the constitutive
atoms, values are pair tuples containing (1) tuple of tuples
containing atom indices for all of the internals
with the given key and (2) a tuple of average internal values
"""
for key, value in internals.items():
idxs, all_parameters = list(zip(*value))
avg_all_parameters = []
for parameters in zip(*all_parameters):
if isinstance(parameters[0], list):
avg = []
for comp in zip(*parameters):
avg_comp = old_div(float(sum(comp)), len(comp))
avg.append(int(round(avg_comp)))
elif isinstance(parameters[0], float):
avg = old_div(sum(parameters), len(parameters))
avg_all_parameters.append(avg)
avg_all_parameters = tuple(avg_all_parameters)
internals[key] = (idxs, avg_all_parameters)
return internals
[docs] def recursivelySort(self, internals):
"""
Return the given dictionary of internal coordinates recursively
sorted by the names of the constitutive atoms.
:type internals: dict
:param internals: keys are tuples of names of the constitutive
atoms, values are pair tuples containing (1) tuple of tuples
containing atom indices for all of the internals
with the given key and (2) a tuple of average internal values
:rtype: OrderedDict
:return: internals recursively sorted
"""
names = list(internals)
afunc = lambda x: tuple([x[i] for i in range(len(names[0]))])
names = sorted(names, key=afunc)
internals_sorted = OrderedDict()
for name in names:
internals_sorted[name] = internals[name]
return internals_sorted
[docs] def setInternals(self, internals, info, internals_excluded=None):
"""
Set internals.
:type internals: dict
:param internals: keys are tuples of names of the constitutive
atoms, values is a list of pair tuples containing (1) a tuple
of atom indices and (2) a tuple of internal values
:type info: InternalInfo
:param info: information for this internal
:type internals_excluded: dict or None
:param internals_excluded: contains internals that should be excluded
from being set. keys are tuples of names of the constitutive atoms,
values is a list of pair tuples containing (1) a tuple of atom indices
and (2) a tuple of internal values or None if there isn't one
:rtype: bool
:return: True if the internal was set, or False if it was not set due to
already appearing in the internals_excluded dict
"""
if internals_excluded is None:
internals_excluded = {}
names = info.names
idxs = info.indices
data = info.data
if idxs in internals_excluded.get(names, ([], None))[0]:
return False
if data is None:
data = (self.astructure.measure(*idxs),)
atuple = (idxs, data)
value = internals.setdefault(names, [])
if atuple not in value:
value.append(atuple)
return True
[docs] def averageAndSort(self, internals):
"""
Average and sort the given internals.
:type internals: dict
:param internals: keys are tuples of names of the constitutive
atoms, values is a list of pair tuples containing (1) a tuple
of atom indices and (2) a tuple of internal values
:rtype: dict
:return: sorted keys where keys are tuples of names of the constitutive
atoms, values are pair tuples containing (1) tuple of tuples
containing atom indices for all of the internals
with the given key and (2) a tuple of average internal values
"""
internals = self.avgInternals(internals)
internals = self.recursivelySort(internals)
return internals
[docs] def setSites(self):
"""
Set the sites dictionary.
"""
def get_base_data(atoms):
mass = atoms[0].atomic_weight
charge = atoms[0].partial_charge
if not charge and atoms[0].formal_charge:
charge = float(atoms[0].formal_charge)
radius = atoms[0].radius
color = list(atoms[0].color.rgb)
return (mass, charge, radius, color)
sites = {}
for info in internals_iterator(self.sts_idxs_dict,
site_iterator,
data_getter=get_base_data):
self.setInternals(sites, info)
self.sites = self.averageAndSort(sites)
[docs] def setBonds(self):
"""
Set the bonds dictionary.
"""
bonds = {}
for info in internals_iterator(self.sts_idxs_dict,
analyze.bond_iterator,
uniqueify=self.uniqueify):
self.setInternals(bonds, info)
self.bonds = self.averageAndSort(bonds)
[docs] def setAngles(self):
"""
Set angles dictionary.
"""
angles = {}
for info in internals_iterator(self.sts_idxs_dict,
analyze.angle_iterator,
uniqueify=self.uniqueify):
self.setInternals(angles, info)
self.angles = self.averageAndSort(angles)
[docs] def setDihedrals(self):
"""
Set dihedrals dictionary.
"""
dihedrals = {}
for info in internals_iterator(self.sts_idxs_dict,
analyze.torsion_iterator,
uniqueify=self.uniqueify):
self.setInternals(dihedrals, info)
self.dihedrals = self.averageAndSort(dihedrals)
[docs] def setImpropers(self):
"""
Set impropers dictionary.
"""
impropers = {}
for info in internals_iterator(self.sts_idxs_dict,
analyze.improper_dihedral_iterator,
uniqueify=self.uniqueify):
self.setInternals(impropers, info)
self.impropers = self.averageAndSort(impropers)
[docs] def getNewAngles(self, new_bonds):
"""
Return the new angles.
:type new_bonds: list
:param new_bonds: contains schrodinger.structure._StructureAtom pair
tuples for new bonds
:rtype: list
:return: contains schrodinger.structure._StructureAtom triples tuples
for new angles
"""
# we do not consider angles involving multiple new bonds
new_angles = []
for aatom, batom in new_bonds:
for catom in aatom.bonded_atoms:
new_angles.append((catom, aatom, batom))
for catom in batom.bonded_atoms:
new_angles.append((aatom, batom, catom))
return new_angles
[docs] def getNewDihedrals(self, new_bonds):
"""
Return the new dihedrals.
:type new_bonds: list
:param new_bonds: contains schrodinger.structure._StructureAtom pair
tuples for new bonds
:rtype: list
:return: contains schrodinger.structure._StructureAtom quadruples tuples
for new dihedrals
"""
# we do not consider dihedrals involving multiple new bonds
new_dihedrals = []
for catom, datom in new_bonds:
for batom in catom.bonded_atoms:
for aatom in batom.bonded_atoms:
if aatom.index != catom.index and aatom.index != datom.index:
new_dihedrals.append((aatom, batom, catom, datom))
for aatom in catom.bonded_atoms:
for batom in datom.bonded_atoms:
if aatom.index != batom.index:
new_dihedrals.append((aatom, catom, datom, batom))
for aatom in datom.bonded_atoms:
for batom in aatom.bonded_atoms:
if batom.index != datom.index and batom.index != catom.index:
new_dihedrals.append((catom, datom, aatom, batom))
return new_dihedrals
[docs] def getNewImpropers(self, new_bonds):
"""
Return the new impropers.
:type new_bonds: list
:param new_bonds: contains schrodinger.structure._StructureAtom pair
tuples for new bonds
:rtype: set
:return: contains schrodinger.structure._StructureAtom quadruples tuples
for new impropers
"""
# we do not consider impropers involving multiple new bonds
new_impropers = set()
for batom, catom in new_bonds:
atoms = [batom, catom]
atoms.extend(batom.bonded_atoms)
atoms.extend(catom.bonded_atoms)
sub_graph = analyze.create_nx_graph(self.astructure, atoms=atoms)
sub_graph.add_edge(batom.index, catom.index)
for idx_quadruple in analyze.improper_dihedral_iterator(
nx_graph=sub_graph):
if batom.index in idx_quadruple and catom.index in idx_quadruple:
atom_quadruple = tuple(
self.astructure.atom[idx] for idx in idx_quadruple)
new_impropers.add(atom_quadruple)
return new_impropers
[docs] def removeInternalsWithNewBonds(self):
"""
Remove internals with new bonds.
"""
def remove(internals):
"""
Do the actual work of removing any of the given internals
that involve particle types marked with the new bond tag.
:type internals: OrderedDict or None
:param internals: the internals with keys being tuples
of particle types
"""
if not internals:
return
for key in list(internals):
for atype in key:
base, tag = split_shadow_bond_tag(atype)
if tag:
internals.pop(key)
break
remove(self.bonds)
remove(self.angles)
remove(self.dihedrals)
remove(self.impropers)
self.newbond_nonbonds = {}
self.newangle_nonbonds = {}
self.newdihedral_nonbonds = {}
self._removeUniqueParticleLabels()
def _removeUniqueParticleLabels(self):
"""
Remove the unique particle labels from the structure.
"""
for atom in self.astructure.atom:
atom.property.pop(CG_UNIQUE_PARTICLE_LABEL_KEY, None)
[docs] def findNewBonds(self,
name_1,
name_2,
bonds_deep,
exclude_what=EXCLUDE_BONDS):
"""
Return the new bond pairs.
:type name_1: str
:param name_1: the name of the first particle
:type name_2: str
:param name_2: the name of the second particle
:type bonds_deep: int
:param bonds_deep: the number of bonds separating the particles
:type exclude_what: str
:param exclude_what: What nonbonds to exclude when regenerating the
lists of included and excluded nonbonds after adding these bonds -
should be one of the EXCLUDE_* constants, potentially combined
for rings with get_combined_exclusion_str
:rtype: list
:return: contains schrodinger.structure._StructureAtom pair
tuples for new bonds
"""
atom_key = CG_UNIQUE_PARTICLE_LABEL_KEY
new_bonds = set()
asl = ('(withinbonds {deep} atom.n {ind}) and '
'(beyondbonds {deepm1} atom.n {ind}) and '
'atom.name {name2}')
for aatom in self.astructure.atom:
if aatom.name == name_1:
this_asl = asl.format(deep=bonds_deep,
ind=aatom.index,
deepm1=bonds_deep - 1,
name2=name_2)
idxs = analyze.evaluate_asl(self.astructure, this_asl)
for idx in idxs:
batom = self.astructure.atom[idx]
pair = tuple(sorted([aatom, batom], key=lambda x: x.index))
if not self.astructure.areBound(*pair):
for atom in pair:
tag = SHADOW_START + str(bonds_deep) + SHADOW_END
atom.property[atom_key] = tag
new_bonds.add(pair)
new_bonds = list(new_bonds)
if new_bonds:
self.updateInternalsForNewBonds(new_bonds,
exclude_what=exclude_what)
return new_bonds
[docs] def updateFusedRingIdxsForNewBonds(self, new_bonds_tuples):
"""
Update the fused ring idxs when a new set of bonds has been defined.
:type new_bonds_tuples: list
:param new_bonds_tuples: contains schrodinger.structure._StructureAtom
pair tuples for new bonds
"""
# we do not consider fused rings involving multiple new bonds
st = self.astructure.copy()
for aatom, batom in new_bonds_tuples:
st.addBond(aatom, batom, 1)
sts_idxs_dict = {}
for mol in st.molecule:
sts_idxs_dict[mol.extractStructure()] = [mol.getAtomIndices()]
self.fused_ring_nonbonds = None
self.setFusedRingNonBonds(sts_idxs_dict=sts_idxs_dict)
def _getNewNonbonds(self, adict):
"""
Return new 1-2, 1-3, or 1-4 nonbonds formed from the given internal
dict of new bonds, angles, or dihedrals, respectively.
:type adict: dict
:param adict: sorted keys where keys are tuples of names of the constitutive
atoms possibly marked with the shadow bond tag, values are pair tuples
containing (1) tuple of tuples containing atom indices for all of the
internals with the given key and (2) a tuple of average internal
values
:rtype: dict
:return: keys are pair tuples of names of the constitutive nonbonding atoms
void of any shadow bond tags, values are pair tuples containing (1)
tuple of tuples containing atom indices for all of the nonbonds
with the given key and (2) a tuple of average internal values
(which for nonbonds is empty but kept for consistency with other
datastructures)
"""
# by "new" it is meant that shadow bonds are present, atom
# names in the incoming internals dict may be marked with
# the shadow bond tag which is necessary to distinguish them
# from their real bonded instances, shadow bond tags are not
# needed in the nonbond dict
all_new_nonbonds = {}
for key, value in adict.items():
# for each new bond, angle, or dihedral type, form
# the relevant nonbond type
new_key = []
for name in (key[0], key[-1]):
base, tag = split_shadow_bond_tag(name)
new_key.append(base)
new_key = tuple(new_key)
# grab the nonbond instances
new_idxs = tuple((idxs[0], idxs[-1]) for idxs in value[0])
# form a dict entry using an empty average internal value tuple
# (kept for consistency with other data structures)
new_nonbonds = {new_key: (new_idxs, ())}
# merge this entry with all entries to prevent potential
# duplication, etc.
all_new_nonbonds = self.mergeNonbonds(all_new_nonbonds,
new_nonbonds)
return all_new_nonbonds
[docs] def updateInternalsForNewBonds(self,
new_bonds_tuples,
exclude_what=EXCLUDE_BONDS):
"""
Update the internals when a new set of bonds has been defined
:type new_bonds_tuples: list
:param new_bonds_tuples: contains schrodinger.structure._StructureAtom
pair tuples for new bonds
:type exclude_what: str
:param exclude_what: What nonbonds to exclude when regenerating the
lists of included and excluded nonbonds after adding these bonds -
should be one of the EXCLUDE_* constants, potentially combined for
rings with get_combined_exclusion_str
"""
# the behavior of uniqueify when defining internals here is different
# because the new bonds are not actually bound, because of this
# the new bonds dictionary will only contain a single key-value pair
new_angles_tuples = self.getNewAngles(new_bonds_tuples)
new_dihedrals_tuples = self.getNewDihedrals(new_bonds_tuples)
new_impropers_tuples = self.getNewImpropers(new_bonds_tuples)
dicts = []
for tuples in [
new_bonds_tuples, new_angles_tuples, new_dihedrals_tuples,
new_impropers_tuples
]:
internals = {}
for atuple in tuples:
info = get_internal_info(list(atuple), uniqueify=self.uniqueify)
self.setInternals(internals, info)
dicts.append(self.averageAndSort(internals))
self.bonds.update(dicts[0])
attrs = ('angles', 'dihedrals', 'impropers')
for attr, ndict in zip(attrs, dicts[1:]):
if not ndict:
continue
odict = getattr(self, attr)
if odict is None:
setattr(self, attr, ndict)
else:
odict.update(ndict)
# the dicts now contain types and instances of new bonds, angles,
# and dihedrals (new meaning containing shadow bonds) formed on this
# function call, where types are marked with a shadow bond tag to
# distinguish them from their real bonded instances, now determine
# the 1-2 bond, 1-3 angle, and 1-4 dihedral nonbonds resulting
# from the new internals and update the stored attributes from
# previous function calls
newbond_nonbonds, newangle_nonbonds, newdihedral_nonbonds = map(
self._getNewNonbonds, dicts[:-1])
self.newbond_nonbonds = self.mergeNonbonds(self.newbond_nonbonds,
newbond_nonbonds)
self.newangle_nonbonds = self.mergeNonbonds(self.newangle_nonbonds,
newangle_nonbonds)
self.newdihedral_nonbonds = self.mergeNonbonds(
self.newdihedral_nonbonds, newdihedral_nonbonds)
self.updateFusedRingIdxsForNewBonds(new_bonds_tuples)
self.setIncludedExcludedNonBonds(exclude_what=exclude_what)
def _findNonBondPairs(self, exclude=None):
"""
This finds all pairs of particle types that have at least one non-bond
interaction. Unlike other types of internals (sites, bonds, angles, etc)
we do not want to enumerate all instances. Therefore, for each pair of
types, we just check pairs of those atoms until we find a single
non-bonded example. Once we have that one example, that's all we need to
add that pair to the non-bond pair type list.
:type exclude: dict or None
:param exclude: contains internals that should not be considered as
example pairs, probably due to the 1-3 angle nonbond setting.
keys are tuples of names of the constitutive atoms, values is a list of
pair tuples containing (1) a tuple of atom indices and (2) a tuple of
internal values or None if there isn't one
:rtype: dict
:return: keys are tuples of names of the constitutive atoms, values is a
list of pair tuples containing (1) a tuple of atom indices and (2) a
tuple of internal values
"""
# Don't enumerate all instances of nonbonds due to time (MATSCI-3342)
nonbonds = {}
# we want to find one example of each (type i, type j) pair of
# non-bonded particles. we don't need to enumerate all pairs of
# (type i, type j) non-bonds, we only need to know that at least one
# example exists.
types = list(self.sites)
# loop i over all types
for itypes_index, itypes in enumerate(types):
# loop j over all types starting at i
for jtypes in types[itypes_index:]:
pair_found = False
# loop over all atoms of type i
for idata in self.sites[itypes][0]:
# data is stored as (index, ), this gets the i atom index
atom_idex = idata[0]
atom_i = self.astructure.atom[atom_idex]
# loop over all atoms of type j
for jdata in self.sites[jtypes][0]:
# get the j atom index
atom_jdex = jdata[0]
atom_j = self.astructure.atom[atom_jdex]
if atom_idex == atom_jdex:
continue
# if atom i and atom j are not bound, we have the one
# example of this (type i, type j) pair we need
if not self.astructure.areBound(atom_i, atom_j):
info = get_internal_info([atom_i, atom_j],
data_getter=lambda x: ())
added = self.setInternals(
nonbonds, info, internals_excluded=exclude)
if added:
pair_found = True
# we no longer need to loop over j atoms
break
if pair_found:
# we no longer need to loop over i atoms
break
return nonbonds
[docs] def setNonBonds(self, exclude_newbonds=True):
"""
Set nonbonds dictionary.
:type exclude_newbonds: bool
:param exclude_newbonds: If True, exclude new bonds added by the user.
If False, only exclude pairs that are actually bonded.
"""
if not self.sites:
self.setSites()
# This call is technically not necessary but keeps the excluded_nonbonds
# property up to date and it is fast. The ForceField class depends on
# it.
self.setBondNonBonds()
if exclude_newbonds:
exclude = self.newbond_nonbonds
else:
exclude = {}
self.excluded_nonbonds = self.mergeNonbonds(self.bond_nonbonds, exclude)
included_nonbonds = self._findNonBondPairs(
exclude=self.excluded_nonbonds)
self.included_nonbonds = self.nonbonds = self.averageAndSort(
included_nonbonds)
[docs] def setBasic(self, exclude_what=EXCLUDE_BONDS):
"""
Set basic dictionaries.
:type exclude_what: str
:param exclude_what: What nonbonds to exclude - should be one of the
EXCLUDE_* constants, potentially combined for rings with
get_combined_exclusion_str
"""
self.setSites()
self.setBonds()
self.setAngles()
self.setDihedrals()
self.setImpropers()
self.setIncludedExcludedNonBonds(exclude_what=exclude_what)
[docs] def setFusedRingNonBonds(self, sts_idxs_dict=None):
"""
Set fused ring nonbonds dictionary.
:type sts_idxs_dict: dict
:param sts_idxs_dict: keys are schrodinger.structure.Structure,
values are lists of atom index tuples, together they define
the system
"""
if self.fused_ring_nonbonds is not None:
return
if sts_idxs_dict is None:
sts_idxs_dict = self.sts_idxs_dict
fused_ring_nonbonds = {}
data_getter = lambda x: ()
for info in internals_iterator(sts_idxs_dict,
fused_ring_nonbond_iterator,
data_getter=data_getter):
self.setInternals(fused_ring_nonbonds, info)
self.fused_ring_nonbonds = self.avgInternals(fused_ring_nonbonds)
[docs] def handleRingExclusions(self):
"""
Handle ring exclusions.
"""
if not self.sites:
self.setSites()
self.setFusedRingNonBonds()
self.excluded_nonbonds = self.mergeNonbonds(self.excluded_nonbonds,
self.fused_ring_nonbonds)
included_nonbonds = self._findNonBondPairs(
exclude=self.excluded_nonbonds)
self.included_nonbonds = self.averageAndSort(included_nonbonds)
[docs] def setIncludedExcludedNonBonds(self, exclude_what=EXCLUDE_BONDS):
"""
Define which non-bond interactions are included and excluded
:type exclude_what: str
:param exclude_what: What to exclude - should be one of the EXCLUDE_*
constants, potentially combined for rings with
get_combined_exclusion_str
"""
if exclude_what == EXCLUDE_NONE:
self.setAllPairs()
else:
bad_str, exclude_rings = parse_combined_exclusion_str(exclude_what)
if bad_str == EXCLUDE_BONDS_ANGLES_DIHEDRALS:
self.setNonAngleNonDihedralNonBonds()
elif bad_str == EXCLUDE_BONDS_ANGLES:
self.setNonAngleNonBonds()
elif bad_str == EXCLUDE_BONDS:
self.setNonBonds()
if exclude_rings:
self.handleRingExclusions()
[docs] def setAngleNonBonds(self):
"""
Set angle nonbonds dictionary.
"""
if self.angle_nonbonds is not None:
return
angle_nonbonds = {}
for info in internals_iterator(self.sts_idxs_dict,
analyze.angle_iterator):
atom_double = [
self.astructure.atom[i]
for i in (info.indices[0], info.indices[-1])
]
info = get_internal_info(atom_double, data_getter=lambda x: ())
self.setInternals(angle_nonbonds, info)
self.angle_nonbonds = self.avgInternals(angle_nonbonds)
[docs] def setNonAngleNonBonds(self):
"""
Set the nonangle nonbonds dictionary.
"""
if not self.sites:
self.setSites()
# exclude all bonding pairs and 1-3 angle pairs, first
# collect those pairs from new bonds and new angles, new
# meaning from shadow bonds, then update the collection
# with pairs from real bonds and angles
excluded_nonbonds = self.mergeNonbonds(self.newbond_nonbonds,
self.newangle_nonbonds)
self.setBondNonBonds()
excluded_nonbonds = self.mergeNonbonds(excluded_nonbonds,
self.bond_nonbonds)
self.setAngleNonBonds()
self.excluded_nonbonds = self.mergeNonbonds(excluded_nonbonds,
self.angle_nonbonds)
included_nonbonds = self._findNonBondPairs(
exclude=self.excluded_nonbonds)
self.included_nonbonds = self.averageAndSort(included_nonbonds)
[docs] def setDihedralNonBonds(self):
"""
Set dihedral nonbonds dictionary.
"""
# improper dihedrals are covered via regular angles
# and therefore do not require special treatment
if self.dihedral_nonbonds is not None:
return
dihedral_nonbonds = {}
for info in internals_iterator(self.sts_idxs_dict,
analyze.torsion_iterator):
atom_double = [
self.astructure.atom[i]
for i in (info.indices[0], info.indices[-1])
]
info = get_internal_info(atom_double, data_getter=lambda x: ())
self.setInternals(dihedral_nonbonds, info)
self.dihedral_nonbonds = self.avgInternals(dihedral_nonbonds)
[docs] def setNonAngleNonDihedralNonBonds(self):
"""
Set the nonangle nondihedral nonbonds dictionary.
"""
if not self.sites:
self.setSites()
# exclude all bonding pairs, 1-3 angle pairs, and 1-4 dihedral
# pairs, first collect those pairs from new bonds, new angles,
# and new dihedrals, new meaning from shadow bonds, then update
# the collection with pairs from real bonds, angles, and dihedrals
excluded_nonbonds = self.mergeNonbonds(self.newbond_nonbonds,
self.newangle_nonbonds)
excluded_nonbonds = self.mergeNonbonds(excluded_nonbonds,
self.newdihedral_nonbonds)
self.setBondNonBonds()
excluded_nonbonds = self.mergeNonbonds(excluded_nonbonds,
self.bond_nonbonds)
self.setAngleNonBonds()
excluded_nonbonds = self.mergeNonbonds(excluded_nonbonds,
self.angle_nonbonds)
self.setDihedralNonBonds()
self.excluded_nonbonds = self.mergeNonbonds(excluded_nonbonds,
self.dihedral_nonbonds)
included_nonbonds = self._findNonBondPairs(
exclude=self.excluded_nonbonds)
self.included_nonbonds = self.averageAndSort(included_nonbonds)
[docs] def setBondNonBonds(self):
"""
Set the bond nonbonds dictionary.
"""
if self.bond_nonbonds is not None:
return
bond_nonbonds = {}
data_getter = lambda x: ()
for info in internals_iterator(self.sts_idxs_dict,
analyze.bond_iterator,
data_getter=data_getter):
self.setInternals(bond_nonbonds, info)
self.bond_nonbonds = self.avgInternals(bond_nonbonds)
[docs] def setAllPairs(self):
"""
Set the all pairs dictionary.
"""
# Regenerate the set of nonbonds, ensuring that any NewBonds added by
# the user are ignored
self.setNonBonds(exclude_newbonds=False)
self.setBondNonBonds()
included_nonbonds = self.mergeNonbonds(self.nonbonds,
self.bond_nonbonds)
self.included_nonbonds = self.recursivelySort(included_nonbonds)
self.excluded_nonbonds = {}
[docs] def mergeNonbonds(self, nonbonds, new_nonbonds):
"""
Return a copy of the first nonbonds dict merged with the second.
:type nonbonds: dict
:param nonbonds: keys are tuples of names of the constitutive atoms, values
is a list of pair tuples containing (1) a tuple of atom indices and (2)
a tuple of internal values or None if there isn't one
:type new_nonbonds: dict
:param new_nonbonds: keys are tuples of names of the constitutive atoms, values
is a list of pair tuples containing (1) a tuple of atom indices and (2)
a tuple of internal values or None if there isn't one
:rtype: dict
:return: merged nonbonds
"""
tmp_nonbonds = nonbonds.copy()
for key, value in new_nonbonds.items():
ovalue = tmp_nonbonds.get(key)
if ovalue:
tmp_value = tuple([i for i in value[0] if i not in ovalue[0]])
tmp_nonbonds[key] = (ovalue[0] + tmp_value, ())
else:
tmp_nonbonds[key] = value
return tmp_nonbonds
[docs]class FusedRing(object):
"""
Manage a fused ring.
"""
[docs] def __init__(self, rings):
"""
Create an instance.
:type rings: list
:param rings: contains schrodinger.structure._Ring
"""
self.rings = rings
self.edge = self._getEdges()
def _getEdges(self):
"""
Return a list of edges.
:rtype: list
:return: contains schrodinger.structure._StructureBond
"""
edges = []
for ring in self.rings:
for edge in ring.edge:
if edge not in edges:
edges.append(edge)
return edges
[docs] def getAtomIndices(self):
"""
Return a list of atom indices.
:rtype: list
:return: atom indices
"""
idxs = []
for ring in self.rings:
for idx in ring.getAtomIndices():
if idx not in idxs:
idxs.append(idx)
return idxs
@staticmethod
def _getIdxsRingsDict(st):
"""
Return an ordered dictionary where keys are tuples of atom
indices and values are lists of ring objects for the
given atoms.
:type st: structure.Structure
:param st: the structure
:rtype: OrderedDict
:return: keys are tuples of atom indices, values are
lists of schrodinger.structure._Ring
"""
# algorithm is:
# (1) loop over rings and consider each as a fusable candidate
# (2) in a macro-loop iteratively iterate through a collection
# of previously categorized rings looking for rings with which
# to fuse the current candidate and create a larger candidate
# (2a) if the current candidate can not be fused then put a new
# entry in the collection of categorized rings
# (2b) if the current candidate can be fused then replace the
# data in the collection of categorized rings with the fused
# equivalent and continue the macro-loop in step (2) until
# the larger candidate can not be fused any longer and step
# (2a) is reached
fix_it = lambda x: tuple(sorted(list(x)))
idxs_rings_dict = OrderedDict()
for ring in st.ring:
fusable_idxs = set(ring.getAtomIndices())
fusable_rings = [ring]
while True:
for idxs in list(idxs_rings_dict):
jdxs = set(idxs)
common = fusable_idxs.intersection(jdxs)
uncommon = fusable_idxs.symmetric_difference(jdxs)
if common and uncommon:
fusable_idxs = fusable_idxs.union(jdxs)
fusable_rings.extend(idxs_rings_dict.pop(idxs))
break
else:
idxs_rings_dict[fix_it(fusable_idxs)] = fusable_rings
break
return idxs_rings_dict
[docs] @staticmethod
def getFusedRings(st):
"""
Return a list of fused ring objects for
the given structure.
:type st: structure.Structure
:param st: the structure
:rtype: list
:return: contains FusedRing
"""
idxs_rings_dict = FusedRing._getIdxsRingsDict(st)
return [FusedRing(rings) for rings in idxs_rings_dict.values()]
[docs]def setCGBondLengths(st):
"""
Set the bond lengths of the given CG structure
according to the particle radii.
:type st: structure.Structure
:param st: the CG structure
:rtype: structure.Structure
:return: the CG structure with the new bond lengths
"""
# this function uses VDW radii rather than covalent radii
get_len = lambda x, y: sum([i.radius for i in (x, y)])
new_st = st.copy()
# first handle bonds in rings and fused ring systems by just
# uniformly scaling them using their average bond length
# scale factor and moving the attached non-ring atoms
# appropriately, then handle the remaining bonds by skipping
# bonds in rings for which .adjust will silently fail
for fring in FusedRing.getFusedRings(new_st):
factors = []
for bond in fring.edge:
length = get_len(bond.atom1, bond.atom2)
factors.append(old_div(length, bond.length))
factor = numpy.mean(factors)
fring_idxs = fring.getAtomIndices()
center = transform.get_centroid(new_st, atom_list=fring_idxs)
center = center[:3]
for idx in fring_idxs:
atom = new_st.atom[idx]
all_moving_atoms = set([atom.index])
for oatom in atom.bonded_atoms:
if oatom.index not in fring_idxs:
all_moving_atoms.update(new_st.getMovingAtoms(atom, oatom))
x, y, z = center + factor * (numpy.array(atom.xyz) - center) - \
numpy.array(atom.xyz)
transform.translate_structure(new_st,
x=x,
y=y,
z=z,
atom_index_list=all_moving_atoms)
ring_bonds = rings.find_ring_bonds(new_st)
for bond in new_st.bond:
atom1, atom2 = bond.atom1, bond.atom2
if {atom1.index, atom2.index} not in ring_bonds:
length = get_len(atom1, atom2)
new_st.adjust(length, atom1, atom2)
return new_st
[docs]def fix_linear_angles(struct, internals=None, angles=None):
"""
Return a copy of the given structure with linear angles
fixed.
:type struct: structure.Structure
:param struct: the structure whose angles will be fixed
:type internals: Internals or None
:param internals: the internals of the structure, if None
it will be defined
:type angles: list or None
:param angles: contains tuples of atom names of the angles
to be fixed, if None then all will be fixed
:rtype: structure.Structure
:return: the structure with angles fixed
"""
# prevent linearity and allow a bending force to be
# defined, do this by randomly perturbing the vertex
# particle rather than via .adjust (see MATSCI-4340),
# note that this function just perturbs the structure
# meaning that it doesn't redefine the PBC nor update
# the internals object
linear_thresh = 1.e-3
is_linear = lambda x: abs(180. - x) <= linear_thresh
st = struct.copy()
if internals is None:
internals = Internals(st)
internals.setAngles()
if angles is None:
angles = list(internals.angles)
to_fixs = []
for angle in angles:
for idxs in internals.angles[angle][0]:
value = st.measure(*idxs)
if is_linear(value):
to_fixs.append(idxs)
# don't share random state and force seed
myrandom = random.Random()
myrandom.seed(0)
noise_lb = 1.e-3
noise_ub = 1.e-2
get_random_noise = lambda: myrandom.choice([-1, 1]) * myrandom.uniform(
noise_lb, noise_ub)
for to_fix in to_fixs:
while is_linear(st.measure(*to_fix)):
atom = st.atom[to_fix[1]]
atom.x += get_random_noise()
atom.y += get_random_noise()
atom.z += get_random_noise()
return st
[docs]def get_unique_bonds(st):
"""
Return a list of unique bonds for the given coarse grain
structure.
:type st: `schrodinger.structure.Structure`
:param st: the coarse grain structure
:rtype: list[`schrodinger.structure._StructureBond`]
:return: unique bonds
"""
assert msutils.is_coarse_grain(st)
key_1 = CG_UNIQUE_BOND_PARTICLE_1_LABEL_KEY
key_2 = CG_UNIQUE_BOND_PARTICLE_2_LABEL_KEY
unique_bonds = []
for bond in st.bond:
for prop_key in [key_1, key_2]:
if bond.property.get(prop_key):
unique_bonds.append(bond)
break
return unique_bonds
[docs]def get_mass(st):
"""
Return the mass of the given coarse grain structure.
:type st: `schrodinger.structure.Structure`
:param st: the coarse grain structure
:rtype: float
:return: the mass in g/mol
"""
assert msutils.is_coarse_grain(st)
return sum(
mm.mmct_atom_get_coarse_grain_mass(st, atom.index) for atom in st.atom)