Source code for schrodinger.application.matsci.cgforcefield

"""
Constants, functions and classes for assigning coarse-grained force fields to
structures

Copyright Schrodinger, LLC. All rights reserved.
"""

import argparse
import glob
import json
import math
import os
from collections import OrderedDict
from collections import namedtuple
from past.utils import old_div

from schrodinger import structure
from schrodinger.application.desmond import cms
from schrodinger.application.desmond import ffiostructure
from schrodinger.application.matsci import coarsegrain
from schrodinger.application.matsci import codeutils
from schrodinger.application.matsci import enc
from schrodinger.application.matsci import jobutils
from schrodinger.application.matsci import msutils
from schrodinger.application.matsci import parserutils
from schrodinger.application.matsci.nano import xtal
from schrodinger.infra import mm
from schrodinger.job.util import hunt
from schrodinger.utils import fileutils

# Database parameter keys
SITE_MASS_KEY = 'mass/(g/mol)'
SITE_CHARGE_KEY = 'charge'
SITE_MARTINI_SUBTYPE_KEY = 'martini_type'
SITE_RADIUS_KEY = 'radius/Ang.'
SITE_COLOR_KEY = 'rgb_color'
SITE_SMARTS_KEY = 'atomistic_smarts'
BOND_EQ_LENGTH_KEY = 'eq_length/Ang.'
BOND_FORCE_CONSTANT_KEY = 'force_constant/(kcal/mol)/Ang.^2'
ANGLE_EQ_ANGLE_KEY = 'eq_angle/deg.'
ANGLE_HARMONIC_FORCE_CONSTANT_KEY = 'force_constant/(kcal/mol)/rad.^2'
ANGLE_TRIGONOMETRIC_FORCE_CONSTANT_KEY = 'force_constant/(kcal/mol)'
DIHEDRAL_V1_KEY = 'v1/(kcal/mol)'
DIHEDRAL_V2_KEY = 'v2/(kcal/mol)'
DIHEDRAL_V3_KEY = 'v3/(kcal/mol)'
DIHEDRAL_V4_KEY = 'v4/(kcal/mol)'
DIHEDRAL_EQ_ANGLE_1_KEY = 'eq_angle_1/deg.'
DIHEDRAL_EQ_ANGLE_2_KEY = 'eq_angle_2/deg.'
DIHEDRAL_EQ_ANGLE_3_KEY = 'eq_angle_3/deg.'
DIHEDRAL_EQ_ANGLE_4_KEY = 'eq_angle_4/deg.'
DIHEDRAL_FORCE_CONSTANT_1_KEY = 'force_constant_1/(kcal/mol)'
DIHEDRAL_FORCE_CONSTANT_2_KEY = 'force_constant_2/(kcal/mol)'
DIHEDRAL_FORCE_CONSTANT_3_KEY = 'force_constant_3/(kcal/mol)'
DIHEDRAL_FORCE_CONSTANT_4_KEY = 'force_constant_4/(kcal/mol)'
DIHEDRAL_MULT_1_KEY = 'n_1'
DIHEDRAL_MULT_2_KEY = 'n_2'
DIHEDRAL_MULT_3_KEY = 'n_3'
DIHEDRAL_MULT_4_KEY = 'n_4'
IMPROPER_EQ_ANGLE_KEY = 'eq_angle/deg.'
IMPROPER_FORCE_CONSTANT_KEY = 'force_constant/(kcal/mol)/rad.^2'
VDW_M_PARAMETER_KEY = 'm'
VDW_N_PARAMETER_KEY = 'n'
VDW_EPSILON_KEY = 'epsilon/(kcal/mol)'
VDW_SIGMA_KEY = 'sigma/Ang.'
VDW_A_PARAMETER_KEY = 'a/(kcal/mol)/Ang.^2'

# shortcuts
DIHEDRAL_TRIG_1_TRIPLE = (DIHEDRAL_EQ_ANGLE_1_KEY,
                          DIHEDRAL_FORCE_CONSTANT_1_KEY, DIHEDRAL_MULT_1_KEY)
DIHEDRAL_TRIG_2_TRIPLE = (DIHEDRAL_EQ_ANGLE_2_KEY,
                          DIHEDRAL_FORCE_CONSTANT_2_KEY, DIHEDRAL_MULT_2_KEY)
DIHEDRAL_TRIG_3_TRIPLE = (DIHEDRAL_EQ_ANGLE_3_KEY,
                          DIHEDRAL_FORCE_CONSTANT_3_KEY, DIHEDRAL_MULT_3_KEY)
DIHEDRAL_TRIG_4_TRIPLE = (DIHEDRAL_EQ_ANGLE_4_KEY,
                          DIHEDRAL_FORCE_CONSTANT_4_KEY, DIHEDRAL_MULT_4_KEY)
DIHEDRAL_TRIG_TRIPLES = (DIHEDRAL_TRIG_1_TRIPLE, DIHEDRAL_TRIG_2_TRIPLE,
                         DIHEDRAL_TRIG_3_TRIPLE, DIHEDRAL_TRIG_4_TRIPLE)

# dummy keys are used as placeholders when having to add
# multiple instances of a given internal to the ffio block
DIHEDRAL_EQ_ANGLE_DUMMY_KEY = 'eq_angle_dummy/deg.'
DIHEDRAL_FORCE_CONSTANT_DUMMY_KEY = 'force_constant_dummy/(kcal/mol)'
DIHEDRAL_MULT_DUMMY_KEY = 'n_dummy'
DIHEDRAL_TRIG_DUMMY_TRIPLE = (DIHEDRAL_EQ_ANGLE_DUMMY_KEY,
                              DIHEDRAL_FORCE_CONSTANT_DUMMY_KEY,
                              DIHEDRAL_MULT_DUMMY_KEY)

# Misc. database keys
SITES_INTERNAL_KEY = 'sites'
SITES_GENERAL_ATOM_KEY = 'atom'
SITES_MARTINI_ATOM_KEY = 'martini_atom'
BONDS_INTERNAL_KEY = 'bonds'
BONDS_HARMONIC_KEY = 'harm'
ANGLES_INTERNAL_KEY = 'angles'
ANGLES_HARMONIC_KEY = coarsegrain.ANGLES_HARMONIC_KEY
ANGLES_TRIGONOMETRIC_KEY = coarsegrain.ANGLES_TRIGONOMETRIC_KEY
DIHEDRALS_INTERNAL_KEY = 'dihedrals'
DIHEDRALS_OPLS_PROPER_KEY = 'opls_proper'
DIHEDRALS_TRIGONOMETRIC_PROPER_KEY = 'proper_trig'
IMPROPERS_INTERNAL_KEY = 'impropers'
IMPROPERS_HARMONIC_KEY = 'improper_harm'
NONBONDS_INTERNAL_KEY = 'nonbonds'
NONBONDS_LENNARD_JONES_KEY = coarsegrain.LENNARD_JONES
NONBONDS_DISSIPATIVE_PARTICLE_KEY = coarsegrain.DISSIPATIVE_PARTICLE
NONBONDS_SHIFTED_LENNARD_JONES_KEY = coarsegrain.SHIFTED_LENNARD_JONES
EXCLUSIONS_INTERNAL_KEY = 'exclude'
DIELECTRIC_INTERNAL_KEY = 'dielectric_constant'
CUTOFF_INTERNAL_KEY = 'cutoff'
REDUCED_DENSITY_INTERNAL_KEY = 'reduced_density'
DESCRIPTION_INTERNAL_KEY = 'description'
SHADOWBOND_TYPE_INTERNAL_KEY = 'shadow_bonds'
GAMMA_DISSIPATIVE_PARTICLE_INTERNAL_KEY = 'dpd_reduced_gamma'
FF_TYPE_INTERNAL_KEY = 'force_field_type'

# Separator of forcefield data type
# example string :- PE,PE
FORCEFIELD_DATA_TYPE_SEPARATOR = ','

# Separator of forcefield data type used in hdf key
# example string :- PE_PE
FORCEFIELD_DATA_TYPE_KEY_SEPARATOR = '_'

# type-to-text maps
get_rev_dict = lambda x: dict([(j, i) for i, j in x.items()])
SITES_TYPE_TEXT_MAP = {
    SITES_GENERAL_ATOM_KEY: coarsegrain.GENERAL_SITE_TYPE,
    SITES_MARTINI_ATOM_KEY: coarsegrain.MARTINI_SITE_TYPE
}
SITES_TEXT_TYPE_MAP = get_rev_dict(SITES_TYPE_TEXT_MAP)
ANGLES_TYPE_TEXT_MAP = {
    ANGLES_HARMONIC_KEY: coarsegrain.HARMONIC,
    ANGLES_TRIGONOMETRIC_KEY: coarsegrain.TRIGONOMETRIC
}
ANGLES_TEXT_TYPE_MAP = get_rev_dict(ANGLES_TYPE_TEXT_MAP)
DIHEDRALS_TYPE_TEXT_MAP = {
    DIHEDRALS_OPLS_PROPER_KEY: coarsegrain.OPLS,
    DIHEDRALS_TRIGONOMETRIC_PROPER_KEY: coarsegrain.TRIGONOMETRIC
}
DIHEDRALS_TEXT_TYPE_MAP = get_rev_dict(DIHEDRALS_TYPE_TEXT_MAP)

# ffio parameters
FFIO_NAME = 'general'
FFIO_COMBINING_RULE = 'geometric'
FFIO_SITE_TYPE = 'atom'
FFIO_SITE_CHARGE = 0
FFIO_VDW_FUNCT = 'polynomial_cij'
FFIO_VDW_NTERMS = 16
FFIO_VDW_TYPE_NAME_BASE = '%s_vdw_%s'
FFIO_VDW_COMBINED_FUNCT = 'polynomial_cij'
FFIO_VDW_COMBINED_T1 = 1
FFIO_PROPERTY_KEY = 'r_ffio_c%s'
FFIO_CT_TYPE_KEY = 's_ffio_ct_type'
FFIO_FULL_SYSTEM = 'full_system'
FFIO_SOLUTE = 'solute'

TYPE_SEPARATOR = ','

# Ewald parameters
# the following polynomial_cij data is obtained by
# fitting to the Ewald real space Coulomb contribution,
# i.e. k_e*erfc(width*r/sqrt(2))/r, where k_e = 1/(4*pi*epsilon_0)
# = 332.0634283 is the Coulomb constant converted from atomic units
# of 1 Hartree*Bohr/(elementary_charge)**2 to
# (kcal/mol)*Ang./(elementary_charge)**2, and using a width obtained
# from some reference cut off via NonbondRow.getWidthParameter and an
# epsilon of EWALD_EPSILON, the data are (index, exponent, prefactor)
# tuples where index is the r_ffio_c<index>, exponent is the r**exponent,
# and prefactor is the constant multiplying r**exponent
EWALD_EPSILON = 1.00E-9
EWALD_REF_CUTOFF_1 = 12.00
EWALD_FIT = {
    EWALD_REF_CUTOFF_1: [
        (1, 0, -60.29447747264311630033), (3, -2, 756.48469519180366660294),
        (5, -4, -1668.18483484391572346794), (7, -6, 3347.27774485723011821392),
        (9, -8, -4185.97700151554909098195), (11, -10,
                                              2823.83598188037967702257),
        (13, -12, -780.26168456680193230568), (14, 1, 10.95682830288469666868),
        (15, 2, -0.74647138663679613035), (16, 3, 0.01801261680487299330)
    ]
} # yapf: disable

# shifted Lennard-Jones parameters
# the following polynomial_cij data is obtained by
# fitting to the shifted Lennard-Jones, found here:
# http://www.cgmartini.nl/index.php/faq (about half way down
# in "Which shift function is used?"), with the following quadruples
# of (m_exponent, n_exponent, r_shift, r_cutoff) and various sigma,
# the Martini FF uses the data in SHIFTED_LJ_REF_KEY_1 with sigma
# SHIFTED_LJ_REF_SIGMA_X X = 1, 2, 3, 4, and 5,
# the data are (index, exponent, prefactor) tuples where index
# is the r_ffio_c<index>, exponent is the r**exponent,
# and prefactor is the constant multiplying r**exponent
# 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)
SHIFTED_LJ_REF_M_1 = 12
SHIFTED_LJ_REF_N_1 = 6
SHIFTED_LJ_REF_SHIFT_1 = 9.
SHIFTED_LJ_REF_CUTOFF_1 = 12.
SHIFTED_LJ_REF_KEY_1 = (SHIFTED_LJ_REF_M_1, SHIFTED_LJ_REF_N_1,
                        SHIFTED_LJ_REF_SHIFT_1, SHIFTED_LJ_REF_CUTOFF_1)
SHIFTED_LJ_REF_SIGMA_1 = coarsegrain.MARTINI_1_SIGMA
SHIFTED_LJ_REF_SIGMA_2 = coarsegrain.MARTINI_2_SIGMA
SHIFTED_LJ_REF_SIGMA_3 = coarsegrain.MARTINI_3_SIGMA
SHIFTED_LJ_REF_SIGMA_4 = coarsegrain.MARTINI_4_SIGMA
SHIFTED_LJ_REF_SIGMA_5 = coarsegrain.MARTINI_5_SIGMA
SHIFTED_LJ_REF_SIGMA_6 = coarsegrain.MARTINI_6_SIGMA
# yapf: disable
SHIFTED_LJ_FIT = {
    SHIFTED_LJ_REF_SIGMA_1: {SHIFTED_LJ_REF_KEY_1:
        [(1, 0, -73.39828345035263623686),
         (2, -1, 3078.25438778890338653582),
         (3, -2, -58741.17745161147468024865),
         (4, -3, 648196.54037688521202653646),
         (5, -4, -4416992.42427262477576732635),
         (6, -5, 18288900.95356648042798042297),
         (7, -6, -39646057.39946962147951126099),
         (8, -7, 8260493.41963949706405401230),
         (9, -8, 152938526.81131097674369812012),
         (10, -9, -249799717.60283806920051574707),
         (13, -12, 874556930.72175145149230957031),
         (16, 3, 0.00112237709062471187)]},
    SHIFTED_LJ_REF_SIGMA_2: {SHIFTED_LJ_REF_KEY_1:
        [(1, 0, -383.42256089551301556639),
         (2, -1, 15675.40687285328385769390),
         (3, -2, -288178.32800226396648213267),
         (4, -3, 3010910.53871261747553944588),
         (5, -4, -18921120.42772997543215751648),
         (6, -5, 69204176.69233772158622741699),
         (7, -6, -121703422.64002545177936553955),
         (9, -8, 235458516.45386704802513122559),
         (13, -12, 9910545882.02453231811523437500),
         (16, 3, 0.00608905594939570041)]},
    SHIFTED_LJ_REF_SIGMA_3: {SHIFTED_LJ_REF_KEY_1:
        [(1, 0, -34.06475554138776828950),
         (2, -1, 1429.81047397730799275450),
         (3, -2, -27539.30168704309471650049),
         (4, -3, 311750.58101170684676617384),
         (5, -4, -2249550.70831884769722819328),
         (6, -5, 10573170.05618538893759250641),
         (7, -6, -31736713.10633733496069908142),
         (8, -7, 56120868.41938952356576919556),
         (9, -8, -46705251.43325433880090713501),
         (10, -9, 3999798.13562615634873509407),
         (13, -12, 206041770.29576241970062255859),
         (16, 3, 0.00052782307087695509)]},
    SHIFTED_LJ_REF_SIGMA_4: {SHIFTED_LJ_REF_KEY_1:
        [(1, 0, -204.48163128283857759016),
         (2, -1, 8180.27961600910202832893),
         (3, -2, -146211.95554758931393735111),
         (4, -3, 1471162.00692251301370561123),
         (5, -4, -8765775.50950740277767181396),
         (6, -5, 29530680.27191623672842979431),
         (7, -6, -44840328.59478991478681564331),
         (9, -8, -18626101.44420809298753738403),
         (10, -9, 275352338.58366030454635620117),
         (13, -12, 2361213317.45370149612426757812),
         (16, 3, 0.00338370618131226706)]},
    SHIFTED_LJ_REF_SIGMA_5: {SHIFTED_LJ_REF_KEY_1:
        [(1, 0, -39.03073697306906097992),
         (2, -1, 1742.57189917887126284768),
         (3, -2, -36259.12575324094359530136),
         (4, -3, 452327.86035708559211343527),
         (5, -4, -3694956.60914124595001339912),
         (6, -5, 20434373.27112276107072830200),
         (7, -6, -76611144.68549674749374389648),
         (8, -7, 188609050.98408910632133483887),
         (9, -8, -280014654.98826020956039428711),
         (10, -9, 195377734.85506445169448852539),
         (13, -12, 2488094.45196075504645705223),
         (16, 3, 0.00053697348759807427)]},
    SHIFTED_LJ_REF_SIGMA_6: {SHIFTED_LJ_REF_KEY_1:
        [(1, 0, -51.19725759172219881066),
         (2, -1, 2152.30007194916152002406),
         (3, -2, -41383.06854359569842927158),
         (4, -3, 464678.72201346338260918856),
         (5, -4, -3286472.17592053301632404327),
         (6, -5, 14784164.58295695297420024872),
         (7, -6, -40208412.44410600513219833374),
         (8, -7, 54114540.17910003662109375000),
         (9, -8, -1161886.85101260547526180744),
         (10, -9, -64705322.71410297602415084839),
         (13, -12, 428826909.41746771335601806641),
         (16, 3, 0.00078565416984644604)]}
}
# yapf: enable

# shifted Coulomb parameters
# the following polynomial_cij data is obtained by
# fitting to the shifted Coulomb, found here:
# http://www.cgmartini.nl/index.php/faq (about half way down
# in "Which shift function is used?"), with alpha = 1 for
# Coulomb and with the following pairs of (r_shift, r_cutoff),
# the Martini FF uses the data in SHIFTED_COULOMB_REF_KEY_1,
# the data are (index, exponent, prefactor) tuples where index
# is the r_ffio_c<index>, exponent is the r**exponent,
# and prefactor is the constant multiplying r**exponent
SHIFTED_COULOMB_REF_SHIFT_1 = 0.
SHIFTED_COULOMB_REF_CUTOFF_1 = 12.
SHIFTED_COULOMB_REF_KEY_1 = (SHIFTED_COULOMB_REF_SHIFT_1,
                             SHIFTED_COULOMB_REF_CUTOFF_1)
SHIFTED_COULOMB_FIT = {
    SHIFTED_COULOMB_REF_KEY_1: [
        (1, 0, -180.42332686399348062878), (2, -1, 9672.26238257898330630269),
        (3, -2, -210386.88374458774342201650), (4, -3,
                                                2389202.76275781635195016861),
        (5, -4,
         -14976514.80119523406028747559), (6, -5,
                                           46201526.30586227774620056152),
        (8, -7,
         -501674356.84652894735336303711), (9, -8,
                                            1585493173.73096442222595214844),
        (10, -9,
         -1750270815.15244007110595703125), (13, -12,
                                             1772552232.45128130912780761719)
    ]
} # yapf: disable

# handle the coarse grain force field parameter files
FF_PARAMETERS_DIR = 'coarse_grain_force_field_parameters'
FF_PARAMETERS_INSTALLED_PATH = os.path.join(hunt('mmshare', dir='data'),
                                            FF_PARAMETERS_DIR)
ORIG_FF_LOCAL_PATH = os.path.join(msutils.get_matsci_user_data_dir(),
                                  FF_PARAMETERS_DIR)
FF_PARAMETERS_LOCAL_PATH = ORIG_FF_LOCAL_PATH
if not os.path.exists(FF_PARAMETERS_LOCAL_PATH):
    try:
        os.makedirs(FF_PARAMETERS_LOCAL_PATH, exist_ok=True)
    except OSError:
        # MATSCI-6198
        raise
        FF_PARAMETERS_LOCAL_PATH = None
FF_PARAMETERS_PATHS = [FF_PARAMETERS_INSTALLED_PATH, FF_PARAMETERS_LOCAL_PATH]

FF_LOCATION_TYPES_TO_PARAMETERS_PATHS_DICT = dict(
    zip(parserutils.CG_FF_LOCATION_TYPES, FF_PARAMETERS_PATHS))

FORCE_FIELD_FILE_EXT = '_cgff.json'

DESCRIPTION_KEY = 's_matsci_ff_description'
FF_TYPE_KEY = 's_matsci_ff_type'

# gamma in reduced units
DPD_REDUCED_GAMMA = 4.5
DPD_REDUCED_GAMMA_PROP_KEY = 'r_matsci_dpd_gamma'

# FF types and maps
GENERAL_TYPE = 'general'
MARTINI_TYPE = 'martini'
SITE_TEXT_TO_FF_TYPE_MAP = {
    coarsegrain.GENERAL_SITE_TYPE: GENERAL_TYPE,
    coarsegrain.MARTINI_SITE_TYPE: MARTINI_TYPE
}
FF_TYPE_TO_SITE_TEXT_MAP = get_rev_dict(SITE_TEXT_TO_FF_TYPE_MAP)

MARTINI_FACTOR = 0.5

# TODO:  These FF constants were moved to parserutils and then renamed. To keep
# these pointers working, we made pointers between the original variable names
# and the new variable names in parserutils (e.g.,
# parserutils.INSTALLED_FF_LOCATION_TYPE points to
# parserutils.INSTALLED_CG_FF_LOCATION_TYPE). When the pointers here are
# removed, we should also remove the pointers in parserutils. See MATSCI-11895.
MOVED_VARIABLES = (  # module, remove_release, variables
    ('parserutils', '22-3', {
        'INSTALLED_FF_LOCATION_TYPE', 'LOCAL_FF_LOCATION_TYPE',
        'FF_LOCATION_TYPES'
    }),)


def __getattr__(name):
    """
    If a variable doesn't exist in the module, check the moved variables

    :param str name: The variable name

    :raise AttributeError: If `name` is not a moved variable

    :rtype: Any
    :return: The moved variable
    """
    try:
        return codeutils.check_moved_variables(name, MOVED_VARIABLES)
    except AttributeError:
        raise AttributeError(f"module '{__name__}' has no attribute '{name}'")


[docs]def get_density_from_scale_factor(density, scale_factor): """ Get the coarse grain density in g/cm^3 after applying scale factor to it :type density: float :param density: current density :type scale_factor: float :param scale_factor: scaling factor to scale the density by :rtype: float :return: scaled density in g/cm^3 """ return density / scale_factor**3
[docs]def get_scale_factor_from_density(old_density, new_density): """ Get the scale factor for change from old density to the new density :type old_density: float :param old_density: old density of the system in g/cm^3 :type new_density: float :param new_density: new (desired) density of the system in g/cm^3 :rtype: float :return: scale factor to change the old density to new density """ return (old_density / new_density)**(1 / 3)
[docs]def get_reduced_density_from_cutoff(nparticles, volume, cutoff): """ Get reduced density (number) from DPD cutoff in Ang. :type nparticles: int :param nparticles: number of particles in the system :type volume: float :param volume: volume of the system :type cutoff: float :param cutoff: DPD potential cutoff :rtype: float :return: reduced density of the DPD simulation """ return (nparticles / volume) * cutoff**3
[docs]def get_cutoff_from_reduced_density(nparticles, volume, reduced_density): """ Get the DPD cutoff in Ang. from the reduced density :type nparticles: int :param nparticles: number of particles in the system :type volume: float :param volume: volume of the system :type reduced_density: float :param reduced_density: reduced density of the DPD simulation :rtype: float :return: DPD potential cutoff """ return (reduced_density * volume / nparticles)**(1 / 3)
[docs]def get_scaled_volume(initial_volume, scale_factor): """ Return the new volume for coarse grain system after applying scale factor :type initial_volume: float :param initial_volume: volume before applying the scale factor :type scale_factor: float :param scale_factor: scale factor to apply volume :rtype: float :return: new volume """ return initial_volume * scale_factor**3
[docs]def write_force_field_to_json(ff_data, ff_file_path): """ Write force field data to a json file :type ff_data: dict :param ff_data: The dictionary that defines the force field :type ff_file_path: str :param ff_file_path: filename path for the json file """ with open(ff_file_path, 'w') as afile: json.dump(ff_data, afile, sort_keys=True, indent=4, separators=(',', ':'))
[docs]def scale_cell(astructure, scale_factor): """ Scale the cell of the given structure according to the given scale factor. :type astructure: structure.Structure :param astructure: the structure whose cell will be scaled :type scale_factor: float :param scale_factor: the uniform scale factor """ vecs = xtal.get_lattice_vectors( *xtal.get_lattice_param_properties(astructure)) fracs = xtal.trans_cart_to_frac_from_vecs(astructure.getXYZ(), *vecs) # scale parameters a, b, c, alpha, beta, gamma = xtal.get_lattice_param_properties(astructure) a *= scale_factor b *= scale_factor c *= scale_factor pdb_params = [a, b, c, alpha, beta, gamma] chorus_params = xtal.get_chorus_from_params(pdb_params) xtal.set_pbc_properties(astructure, chorus_params) scaled_vecs = xtal.get_vectors_from_chorus(astructure) astructure.setXYZ(xtal.trans_frac_to_cart_from_vecs(fracs, *scaled_vecs)) # clear ASU and fractional properties xtal.clear_asu_and_fractional_properties(astructure)
[docs]def is_like_lennard_jones(nonbond_type): """ Return True if the given type is like Lennard-Jones. :type nonbond_type: str :param nonbond_type: the nonbond type :rtype: bool :return: True if like Lennard-Jones """ return nonbond_type == NONBONDS_LENNARD_JONES_KEY or \ nonbond_type == NONBONDS_SHIFTED_LENNARD_JONES_KEY
[docs]def get_type_dict_from_ff_file(ff_file): """ Return a sorted type dict for the given force field file. :type ff_file: str :param ff_file: path to the force field file from which to obtain the type dict :rtype: OrderedDict :return: the sorted type dict, keys are particle names, values are ParticleInfo """ type_dict = OrderedDict() ff_dict = load_force_field_parameters(ff_file) site_type = ForceField.getSiteType(ff_dict) sites_dict = ff_dict[SITES_INTERNAL_KEY] for name, value in sites_dict.items(): if name not in type_dict: key = name # see MATSCI-4345 - include site radius and color in FF files # but provide backward compatibility rgb = tuple(value[site_type].get(SITE_COLOR_KEY, [255, 0, 0])) radius = value[site_type].get(SITE_RADIUS_KEY, coarsegrain.DEFAULT_RADIUS) mass = value[site_type][SITE_MASS_KEY] info = coarsegrain.ParticleInfo(key=key, name=name, rgb=rgb, radius=radius, mass=mass) type_dict[name] = info type_dict = OrderedDict(sorted(type_dict.items())) return type_dict
[docs]def get_all_force_field_paths(allow_enc=False): """ Get a dictionary of force field names and force field file paths for all defined CG force fields, keyed at the top level by location, INSTALLED_CG_FF_LOCATION_TYPE for FF_PARAMETERS_INSTALLED_PATH and LOCAL_CG_FF_LOCATION_TYPE for FF_PARAMETERS_LOCAL_PATH. :type allow_enc: bool :param allow_enc: whether to allow encrypted force field files :rtype: dict :return: top level keys are either INSTALLED_CG_FF_LOCATION_TYPE or LOCAL_CG_FF_LOCATION_TYPE while inner OrderedDicts have force field names as keys and absolute paths to force field files as values, inner dictionaries are sorted """ all_paths = {} for atype, location in FF_LOCATION_TYPES_TO_PARAMETERS_PATHS_DICT.items(): if location is not None: unsorted_paths = {} pattern = os.path.join(location, '*' + FORCE_FIELD_FILE_EXT) glob_paths = glob.glob(pattern) if allow_enc: pattern += enc.ENC_FILE_EXT glob_paths += glob.glob(pattern) for path in glob_paths: unsorted_paths[get_force_field_name(path)] = path paths = OrderedDict() for key in sorted(list(unsorted_paths)): paths[key] = unsorted_paths[key] if paths: all_paths[atype] = paths return all_paths
[docs]def add_local_type_force_field_to_job_builder(builder, ffname, local_type_includes_cwd=True): """ Set up the Launch API Job Builder to use the given local type force field. :type builder: `schrodinger.job.launchapi.JobSpecificationArgsBuilder` :param builder: The job builder to use for setting the FF as an input file :type ffname: str :param ffname: The name of the force field to use (should be the user-facing name, not the force field file name) :type local_type_includes_cwd: bool :param local_type_includes_cwd: whether the LOCAL_CG_FF_LOCATION_TYPE also includes the CWD """ # by local type it is meant that this force field is # not in the standard Schrodinger installation and # hence needs to be explicity added to Job Control path = get_force_field_file_path( ffname, location_type=parserutils.LOCAL_CG_FF_LOCATION_TYPE, local_type_includes_cwd=local_type_includes_cwd, check_existence=True) # Launch API behavior is to copy files from absolute paths right into the # job folder, which is what we want - we'll be able to access the file in # the working directory via just the file name. builder.setInputFile(path)
[docs]def force_field_name_to_file_name(ffname): """ Convert the user-facing force field name to a file name - without the full path :type ffname: str :param ffname: The name of the force field :rtype: str :return: The name of the force field file """ return ffname + FORCE_FIELD_FILE_EXT
[docs]def get_force_field_file_path( force_field_name, location_type=parserutils.LOCAL_CG_FF_LOCATION_TYPE, local_type_includes_cwd=False, check_existence=False): """ Return the force field file path given the force field name and location type. :type force_field_name: str :param force_field_name: the force field name :type location_type: str :param location_type: the location type, either INSTALLED_CG_FF_LOCATION_TYPE or LOCAL_CG_FF_LOCATION_TYPE :type local_type_includes_cwd: bool :param local_type_includes_cwd: whether the LOCAL_CG_FF_LOCATION_TYPE also includes the CWD :type check_existence: bool :param check_existence: check the existence of the path :raise: ValueError: If the directory the file should be found in does not exist and cannot be made, or if an existence check is being performed and the path doesn't exist :rtype: str :return: the force field file path """ # if local_type_includes_cwd is True and a force field file # with the given name exists in the CWD as well as in the # standard FF_PARAMETERS_LOCAL_PATH location then the former # takes precedence file_name = force_field_name_to_file_name(force_field_name) path = None if location_type == parserutils.LOCAL_CG_FF_LOCATION_TYPE and \ local_type_includes_cwd and os.path.exists(file_name): path = '' if path is None: path = FF_LOCATION_TYPES_TO_PARAMETERS_PATHS_DICT[location_type] if path is None: raise ValueError('Unable to locate or create the local ' 'coarse-grained force field directory, which ' f'should be at {ORIG_FF_LOCAL_PATH}') path = os.path.join(path, file_name) enc_path = path + enc.ENC_FILE_EXT if not os.path.exists(path) and os.path.exists(enc_path): path = enc_path if check_existence and not os.path.exists(path): local_locations = FF_PARAMETERS_LOCAL_PATH if local_type_includes_cwd: cwd = os.getcwd() local_locations = ('current working directory of ' '{cwd} and {std}').format(cwd=cwd, std=local_locations) msg = ('The force field {ffld} can not be found at the ' '\'{atype}\' location. For location type ' '{install_type} a standard directory in the ' 'Schrodinger installation is searched while for ' 'location type {local_type} the following are ' 'searched: {local_locations}. The Coarse Grain ' 'Force Field Assignment GUI may be used to create ' 'force field files in the {std} location.') msg = msg.format(ffld=force_field_name, atype=location_type, install_type=parserutils.INSTALLED_CG_FF_LOCATION_TYPE, local_type=parserutils.LOCAL_CG_FF_LOCATION_TYPE, local_locations=local_locations, std=FF_PARAMETERS_LOCAL_PATH) raise ValueError(msg) return path
[docs]def get_force_field_name(force_field_file_path): """ Return the force field name given the force field file path. :type force_field_file_path: str :param force_field_file_path: the force field file path :rtype: str :return: the force field name """ name = os.path.split(force_field_file_path)[1] ext = FORCE_FIELD_FILE_EXT if name.endswith(enc.ENC_FILE_EXT): ext += enc.ENC_FILE_EXT return name.replace(ext, '')
[docs]def load_force_field_parameters(path): """ Load force field parameters. :type path: str :param path: The path to the force field file to read :rtype: dict or None :return: the force field parameters dictionary or None if there isn't one """ try: text = enc.get_unencrypted_text(path, enc.CG_KEY_FILE_NAME) except ValueError: # this means the given path is not encrypted ff_parameters_dict = {} if os.path.exists(path): with open(path, 'r') as afile: ff_parameters_dict = json.load(afile) else: ff_parameters_dict = json.loads(text) return ff_parameters_dict
[docs]def get_type_name(types): """ Get combined name for a parameter based on the given particle names :type types: tuple :param types: A tuple of particle names that should be joined together for form the type name :rtype: str :return: The particle names in types joined by the TYPE_SEPARATOR character """ return TYPE_SEPARATOR.join(types)
[docs]def split_type_name(name): """ Split a type name up into particle types based on the TYPE_SEPARATOR delimiter :type name: str :param name: A string with the particle type names joined by the TYPE_SEPARATOR delimiter :rtype: list :return: The list of particle type names that formed name """ return name.split(TYPE_SEPARATOR)
[docs]def get_vdw_type_name(ff_type, name): """ Get the VDW type name based on the given FF type and site name. :type ff_type: str :param ff_type: The FF type :type name: str :param name: The name of the site for this vdw type :rtype: str :return: The VDW type name for this particle name """ return FFIO_VDW_TYPE_NAME_BASE % (ff_type, name)
[docs]class CGFFIOStructure(object): """ Manage building a coarse grain structure. """
[docs] def __init__(self, st): """ Create an instance. :type st: structure.Structure :param st: the structure for which the coarse grain structure is needed """ self.st = st.copy() self.ffio_block = None self.createFFIOBlock()
[docs] def createFFIOBlock(self): """ Create a ffio block in the structure. """ mm.mmffio_initialize(mm.error_handler) ffio_block = mm.mmffio_ff_new() mm.mmffio_ff_mmct_put(ffio_block, self.st) mm.mmffio_terminate() st_str = self.st.writeToString(fileutils.MAESTRO) st_reader = structure.StructureReader.fromString( st_str, format=fileutils.MAESTRO) self.st = next(st_reader) self.ffio_block = ffiostructure.FFIOStructure(self.st)
[docs] def setName(self, name=FFIO_NAME): """ Set the name. :type name: str :param name: the name """ self.ffio_block.ffio.name = name
[docs] def setCombiningRule(self, combining_rule=FFIO_COMBINING_RULE): """ Set the combining rule. :type combining_rule: str :param combining_rule: the combining rule """ self.ffio_block.ffio.combining_rule = combining_rule
[docs] def setDescription(self, description=None): """ Set the description. :type description: str or None :param description: the description or None if there isn't one """ if description is not None: self.ffio_block.ffio.property[DESCRIPTION_KEY] = description
[docs] def setFFType(self, ff_type=GENERAL_TYPE): """ Set the FF type. :type ff_type: str or None :param ff_type: the FF type, if None then GENERAL_TYPE is used """ if not ff_type: ff_type = GENERAL_TYPE self.ffio_block.ffio.property[FF_TYPE_KEY] = ff_type
[docs] def writeCMS(self, st, file_name): """ Write `*cms` file. :type st: structure.Structure :param st: the full system structure :type file_name: str :param file_name: the file name """ # structure must be passed in because the .st attribute # necessarily has an empty ffio block st.property[FFIO_CT_TYPE_KEY] = FFIO_FULL_SYSTEM st.write(file_name) self.ffio_block.property[FFIO_CT_TYPE_KEY] = FFIO_SOLUTE ffiostructure.write_cms(self.ffio_block, file_name, mode=mm.M2IO_APPEND)
[docs]class FFIOParam(object): """ The base class for FFIO parameters such as bonds, angles, VDW parameters, exclusions, etc. """ # In base classes, this list will be the set of dictionary keys in the inner # dictionary of the FF dictionary in the order the values should be # converted to the FFIO cx parameters PARAM_ORDER = [] # The name of the FFIO parameters in the order they are given to the # addToFFST method. These are *generally* site-specific data and not part of # the force field (such as atom indexes). NAMED_PROPERTIES = ['ai', 'aj', 'ak', 'al'] # This dictionary should be filled with keys that are FFIO parameter names # and values that are constant values that do not vary based on the involved # particle types CONSTANT_PROPERTIES = {} # The value of the funct FFIO parameter for classes that do not take this # value from the database dictionary FUNCTION = "" # The starting index of the FFIO cx parameters CSTART = 1
[docs] def __init__(self, data, ff_type=GENERAL_TYPE): """ Create an instance :type data: dict or None :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. These values in most classes are converted to FFIO cx parameters using the PARAM_ORDER list. :type ff_type: str :param ff_type: the force field type, either GENERAL_TYPE or MARTINI_TYPE """ # Set the default value for the .funct parameter self.function = self.FUNCTION self.ff_type = ff_type self.extractData(data)
[docs] def extractData(self, data): """ Extract the data from the database parameters :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. These values in most classes are converted to FFIO cx parameters using the PARAM_ORDER list. """ if data: self.function, pdict = list(data.items())[0] self.params = [pdict[x] for x in self.PARAM_ORDER] else: self.params = []
[docs] def getObject(self, fst): """ Reimplement in subclasses to return the proper FFIO block object for the class :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOSubBlock` :return: The block object for this parameter """ raise NotImplementedError
[docs] def addParams(self, fobject): """ Add the parameters taken from the database data to the FFIO block. These parameters generally have generic cx (x=integer) names. :type fobject: ffiostructure._FFIOSubBlock :param fobject: The object to add the parameters to """ for index, param in enumerate(self.params, self.CSTART): fobject.property[FFIO_PROPERTY_KEY % index] = param
[docs] def addNamedProperties(self, fobject, named_props): """ Add the parameters passed to the addToFST method to the FFIO block. These parameters generally have custom names. :type fobject: ffiostructure._FFIOSubBlock :param fobject: The object to add the parameters to :type named_props: tuple :param named_props: Property values in the same order as the NAMED_PROPERTIES list """ for name, prop in zip(self.NAMED_PROPERTIES, named_props): setattr(fobject, name, prop) if self.function: fobject.funct = self.function
[docs] def addConstantProperties(self, fobject): """ Add any properties that are constant for all instances of this class :type fobject: ffiostructure._FFIOSubBlock :param fobject: The object to add the parameters to """ if self.CONSTANT_PROPERTIES: for key, value in self.CONSTANT_PROPERTIES.items(): setattr(fobject, key, value)
[docs] def addToFFST(self, fst, *named_props): """ Add a block with associated properties to the CGFFIO structure :type fst: `CGFFIOStructure` :param fst: The ffio structure to add the block to @keywrd named_props: method arguments. There should be as many arguments and in the same order as the NAMED_PROPERTIES list :rtype: `schrodinger.application.desmond.ffiostructure._FFIOSubBlock` :return: The FFIO object with all parameters added to it. """ fobject = self.getObject(fst) self.addNamedProperties(fobject, named_props) self.addConstantProperties(fobject) self.addParams(fobject) return fobject
[docs]class FFIOBond(FFIOParam): """ The class that handles the FFIO bond block """ PARAM_ORDER = [BOND_EQ_LENGTH_KEY, BOND_FORCE_CONSTANT_KEY]
[docs] def getObject(self, fst): """ Return the ffio bond block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOBond` :return: The block object for this parameter """ return fst.ffio_block.ffio.addBond()
[docs] def extractData(self, data): """ Extract the data from the database parameters :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. These values in most classes are converted to FFIO cx parameters using the PARAM_ORDER list. """ self.function, pdict = list(data.items())[0] length, force_constant = [pdict[x] for x in self.PARAM_ORDER] if self.ff_type == MARTINI_TYPE: force_constant *= MARTINI_FACTOR self.params = [length, force_constant]
[docs]class FFIOAngle(FFIOParam): """ The class that handles the FFIO angle block """
[docs] def getObject(self, fst): """ Return the ffio angle block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOAngle` :return: The block object for this parameter """ return fst.ffio_block.ffio.addAngle()
[docs] def extractData(self, data): """ Extract the data from the database parameters. Modified from the parent class because the cx parameters must be computed from the database data rather than using the data directly. :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. """ self.function, pdict = list(data.items())[0] if self.function == ANGLES_HARMONIC_KEY: param_order = [ ANGLE_EQ_ANGLE_KEY, ANGLE_HARMONIC_FORCE_CONSTANT_KEY ] angle, force_constant = [pdict[x] for x in param_order] elif self.function == ANGLES_TRIGONOMETRIC_KEY: param_order = [ ANGLE_EQ_ANGLE_KEY, ANGLE_TRIGONOMETRIC_FORCE_CONSTANT_KEY ] angle, force_constant = [pdict[x] for x in param_order] angle = math.cos(math.radians(angle)) if self.ff_type == MARTINI_TYPE: force_constant *= MARTINI_FACTOR self.params = [angle, force_constant]
[docs]class FFIODihedral(FFIOParam): """ The class that handles the FFIO dihedral block """ # FFIO dihedral cx parameters start at 0 for whatever reason CSTART = 0
[docs] def getObject(self, fst): """ Return the ffio dihedral block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIODihedral` :return: The block object for this parameter """ return fst.ffio_block.ffio.addDihedral()
[docs] def extractData(self, data): """ Extract the data from the database parameters. Modified from the parent class because the cx parameters must be computed from the database data rather than using the data directly and to add additional parameters that fill out the required c0-c6 parameter block. :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. """ # see MATSCI-3313 where it was found that values for all indices # must be explicity given otherwise Desmond will exit with an error self.function, pdict = list(data.items())[0] if self.function == DIHEDRALS_OPLS_PROPER_KEY: param_order = [ DIHEDRAL_V1_KEY, DIHEDRAL_V2_KEY, DIHEDRAL_V3_KEY, DIHEDRAL_V4_KEY ] self.params = [pdict[x] for x in param_order] self.params = [0.0] + self.params + [0.0, 0.0] elif self.function == DIHEDRALS_TRIGONOMETRIC_PROPER_KEY: param_order = [ DIHEDRAL_EQ_ANGLE_DUMMY_KEY, DIHEDRAL_FORCE_CONSTANT_DUMMY_KEY ] self.params = [pdict[x] for x in param_order] self.params += 6 * [0.0] multiplicity = pdict[DIHEDRAL_MULT_DUMMY_KEY] self.params[multiplicity + 1] = self.params[1]
[docs]class FFIOImproper(FFIODihedral): """ The class that handles the FFIO improper block, which is just part of the dihedral block """ PARAM_ORDER = [IMPROPER_EQ_ANGLE_KEY, IMPROPER_FORCE_CONSTANT_KEY]
[docs] def extractData(self, data): """ Extract the data from the database parameters. Modified from the parent class to add additional parameters that fill out the required c0-c6 parameter block. :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. These values in most classes are converted to FFIO cx parameters using the PARAM_ORDER list. """ FFIOParam.extractData(self, data) self.params += 5 * [0.0]
[docs]class FFIOExclusion(FFIOParam): """ The class that handles the FFIO exclusion list block """
[docs] def getObject(self, fst): """ Return the ffio exclusion block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOExclusion` :return: The block object for this parameter """ return fst.ffio_block.ffio.addExclusion()
[docs]class FFIOVdwType(FFIOParam): """ The class that handles the FFIO VdW type list block """ NAMED_PROPERTIES = ['name'] FUNCTION = FFIO_VDW_FUNCT
[docs] def getObject(self, fst): """ Return the ffio Vdwtype block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOVdwtype` :return: The block object for this parameter """ return fst.ffio_block.ffio.addVdwtype()
[docs]class FFIOSite(FFIOParam): """ The class that handles the FFIO site block """ # These will be the arguments to the addToFST method NAMED_PROPERTIES = ['site', 'vdwtype', 'mass', 'charge'] CONSTANT_PROPERTIES = {'type': FFIO_SITE_TYPE}
[docs] def getObject(self, fst): """ Return the ffio Site block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOSite` :return: The block object for this parameter """ return fst.ffio_block.ffio.addSite()
[docs] def extractData(self, data): """ Extract the data from the database parameters. Modified from the parent class because unlike other blocks, sites do not have cx parameters but do have named parameters that appear in the database dictionary. :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. """ self.function, pdata = list(data.items())[0] self.mass = pdata[SITE_MASS_KEY] self.charge = pdata[SITE_CHARGE_KEY] self.params = [] # if dielectric ever becomes site dependent then set it here from pdata # and consider redefining charge to always be the scaled charge and # update .getEwaldRealSpaceFFIOConstants and # .getShiftedCoulombFFIOConstants self.dielectric = None
[docs] def setDielectricConstant(self, dielectric): """ Set the dielectric constant for the site. :type dielectric: float :param dielectric: the dielectric constant """ if self.dielectric is None: self.dielectric = dielectric
[docs] def addNamedProperties(self, fobject, named_props): """ Add the parameters passed to the addToFST method to the FFIO block. These parameters generally have custom names. Modified from the parent class to include the mass and charge parameters which have been extracted from the database dictionary. :type fobject: ffiostructure._FFIOSite :param fobject: The object to add the parameters to :type named_props: list :param named_props: Two item list, first item is the name of the site and second is the vdw type of the site """ charge = self.charge / math.sqrt(self.dielectric) named_props = named_props + (self.mass, charge) FFIOParam.addNamedProperties(self, fobject, named_props)
[docs]class FFIOVdwTypesCombined(FFIOParam): """ The class that handles the FFIO combined VdW block """ FUNCTION = FFIO_VDW_FUNCT NAMED_PROPERTIES = ['name1', 'name2'] CONSTANT_PROPERTIES = {'t1': FFIO_VDW_COMBINED_T1}
[docs] def __init__(self, cutoff, dielectric, charge_i, charge_j, *args, **kwargs): """ Create an instance :type cutoff: float :param cutoff: The force field cutoff parameter in Angstroms :type dielectric: float :param dielectric: The force field dielectric constant :type charge_i: float :param charge_i: the charge on the first site :type charge_j: float :param charge_j: the charge on the second site See parent class for additional argument documentation """ self.cutoff = cutoff self.dielectric = dielectric self.charge_i = charge_i self.charge_j = charge_j FFIOParam.__init__(self, *args, **kwargs)
[docs] def getObject(self, fst): """ Return the ffio VdwTypeCombined block :type fst: `CGFFIOStructure` :param fst: The structure to get the block object from :rtype: `schrodinger.application.desmond.ffiostructure._FFIOVdwtypescombined` :return: The block object for this parameter """ return fst.ffio_block.ffio.addVdwtypescombined()
[docs] def extractData(self, data): """ Extract the data from the database parameters. Modified from the parent class because the cx parameters must be computed from the database data rather than using the data directly. :type data: dict :param data: The two-level dictionary of forcefield data for this parameter. The top-level dictionary has one key - the function parameter for this object. The value is another dictionary. The keys in the inner dictionary are force field parameter names and the values are the values for those parameters. """ self.particle_type, pdict = list(data.items())[0] if is_like_lennard_jones(self.particle_type): self.mval = pdict[VDW_M_PARAMETER_KEY] self.nval = pdict[VDW_N_PARAMETER_KEY] self.epsilon = pdict[VDW_EPSILON_KEY] self.sigma = pdict[VDW_SIGMA_KEY] elif self.particle_type == NONBONDS_DISSIPATIVE_PARTICLE_KEY: self.aval = pdict[VDW_A_PARAMETER_KEY] self.params = self.getFFIOConstants()
[docs] def getAlphaParameter(self): """ Return the alpha parameter. :rtype: float :return: the alpha parameter """ mval = float(self.mval) nval = float(self.nval) a1 = (old_div(nval, mval))**(old_div(nval, (mval - nval))) a2 = (old_div(nval, mval))**(old_div(mval, (mval - nval))) return old_div(1, (a1 - a2))
[docs] def getWidthParameter(self, cutoff, ewald_epsilon=EWALD_EPSILON): """ Return the Ewald charge distribution width parameter. :type cutoff: float :param cutoff: The cutoff parameter in Angstroms :type ewald_epsilon: float :param ewald_epsilon: parameter in reciprocal Ang. used to determine the width :rtype: float :return: the width in reciprocal Ang. """ # the Ewald charge distribution rho(r) goes as # (width**3)*exp(-(width*r)**2) # and this width divides up the Coulomb potential into real and # reciprocal space Ewald contributions cxe = cutoff * ewald_epsilon tau = math.sqrt(abs(math.log(cxe))) width = old_div(math.sqrt(abs(math.log(cxe * tau))), cutoff) return width
[docs] def getScaledEwaldConstants(self): """ Scale the given Ewald constant given the cutoff and exponent. :rtype: dict :return: Keys are index and values are scaled Ewald constants to use for the parameter with that index """ # the equation to map polynomical_cij constants from those # in EWALD_FIT for some reference width to those for other widths # is obtained scaling the variable r in the fit to width*r, defining # the constants as functions of the width, and then defining those # functions in terms of the constants in EWALD_FIT and the reference # width width_ratio = (old_div(self.getWidthParameter(self.cutoff), self.getWidthParameter(EWALD_REF_CUTOFF_1))) constants = {} for idx, exponent, constant in EWALD_FIT[EWALD_REF_CUTOFF_1]: constants[idx] = constant * (width_ratio**exponent) return constants
[docs] def getEwaldRealSpaceFFIOConstants(self): """ Return the Ewald real space Coulomb potential FFIO constants. :rtype: dict :return: The real space FFIO constants. Keys are 1-based parameter index as in r_ffio_c<idx>, and values are the constants """ factor = self.charge_i * self.charge_j / self.dielectric real_space_constants = {} for index, constant in self.getScaledEwaldConstants().items(): real_space_constants[index] = constant * factor return real_space_constants
[docs] def getShiftedLennardJonesFFIOConstants(self, epsilon, sigma, key=SHIFTED_LJ_REF_KEY_1): """ Return the shifted Lennard-Jones potential FFIO constants. :type epsilon: float :param epsilon: the epsilon parameter in kcal/mol :type sigma: float :param sigma: the sigma parameter in Ang. :type key: str :param key: the key that indexes the fit constants, should be keys for the sigma dicts used in SHIFTED_LJ_FIT :rtype: dict :return: The shifted Lennard-Jones FFIO constants. Keys are 1-based parameter index as in r_ffio_c<idx>, and values are the constants """ # keying off of the float sigma which comes from a combobox may # look odd but for such potentials (for example Martini) this # value can only take on certain values and so this function # should never be called for any other value but those supported constants = {} for index, exponent, constant in SHIFTED_LJ_FIT[sigma][key]: constants[index] = constant * epsilon return constants
[docs] def getShiftedCoulombFFIOConstants(self, key=SHIFTED_COULOMB_REF_KEY_1): """ Return the shifted Coulomb potential FFIO constants. :type key: str :param key: the key that indexes the fit constants, should be keys used in SHIFTED_COULOMB_FIT :rtype: dict :return: The shifted Coulomb FFIO constants. Keys are 1-based parameter index as in r_ffio_c<idx>, and values are the constants """ factor = self.charge_i * self.charge_j / self.dielectric constants = {} for index, exponent, constant in SHIFTED_COULOMB_FIT[key]: constants[index] = constant * factor return constants
[docs] def getFFIOConstants(self): """ Return the FFIO constants. :rtype: list :return: the FFIO constants as a list in order from 0 to NTERMS. All constants, including those defined as zero, are included in the list. """ # the polynomial_cij = c1/r**0 + c2/r**1 + c3/r**2 + ... + # c13/r**12 + c14*r**1 + c15*r**2 + c16*r**3 # see MATSCI-3266 where it was found that values for all indices # must be explicity given otherwise the Desmond parser will stop # on the first undefined parameter values = [0.0] * FFIO_VDW_NTERMS if self.particle_type == NONBONDS_LENNARD_JONES_KEY: axe = self.getAlphaParameter() * self.epsilon soc = old_div(self.sigma, self.cutoff) c1 = axe * ((soc)**self.nval - (soc)**self.mval) values[0] = c1 values[self.nval] = -axe * self.sigma**self.nval values[self.mval] = axe * self.sigma**self.mval elif self.particle_type == NONBONDS_DISSIPATIVE_PARTICLE_KEY: values[0] = self.aval * self.cutoff**2 values[13] = -2 * self.aval * self.cutoff values[14] = self.aval elif self.particle_type == NONBONDS_SHIFTED_LENNARD_JONES_KEY: constants = self.getShiftedLennardJonesFFIOConstants( self.epsilon, self.sigma, key=SHIFTED_LJ_REF_KEY_1) for one_based_index, constant in constants.items(): values[one_based_index - 1] += constant # handle electrostatics if self.particle_type == NONBONDS_LENNARD_JONES_KEY: ewald_constants = self.getEwaldRealSpaceFFIOConstants() for one_based_index, ewald_constant in ewald_constants.items(): values[one_based_index - 1] += ewald_constant elif self.particle_type == NONBONDS_SHIFTED_LENNARD_JONES_KEY: constants = self.getShiftedCoulombFFIOConstants( key=SHIFTED_COULOMB_REF_KEY_1) for one_based_index, constant in constants.items(): values[one_based_index - 1] += constant return values
def _get_approx_time_step(site_fc_data): """ Return the approximate time step from the given site and force constant data. :type site_fc_data: list :param site_fc_data: contains (first site, second site, force constant) tuples :rtype: float :return: approximate time step in fs """ inv_freqs_squared = [] for site_i, site_j, fc in site_fc_data: mass_i = site_i.mass mass_j = site_j.mass reduced_mass = mass_i * mass_j / (mass_i + mass_j) inv_freqs_squared.append(reduced_mass / fc) return math.sqrt(min(inv_freqs_squared)) def _get_site_fc_data_from_bonds(fst): """ Return site and force constant data from bonds. :type fst: `CGFFIOStructure` :param fst: The FFIO structure containing the required data :rtype: list :return: contains (first site, second site, force constant) tuples """ site_fc_data = {} for bond in fst.ffio_block.ffio.bond: site_i = fst.ffio_block.ffio.site[bond.ai] site_j = fst.ffio_block.ffio.site[bond.aj] names = (site_i.site, site_j.site) if names not in site_fc_data: site_fc_data[names] = (site_i, site_j, bond.c2) return list(site_fc_data.values())
[docs]class MissingParameterError(Exception): """ Raised if the force field is missing required data for a structure """
# Used to store data for a specific VDW site VDWData = namedtuple('VDWData', ['vtype', 'charge']) ShadowData = namedtuple('ShadowData', ['name1', 'name2', 'depth'])
[docs]class ForceField(object): """ Reads in force field data and applies it to a structure """ # Translates database dictionary keys to the class that manages that FFIO # data FP_CLASSES = { SITES_INTERNAL_KEY: FFIOSite, BONDS_INTERNAL_KEY: FFIOBond, ANGLES_INTERNAL_KEY: FFIOAngle, IMPROPERS_INTERNAL_KEY: FFIOImproper }
[docs] def __init__(self, path=None, data=None): """ Create a ForceField instance :type path: str :param path: The path to the force field file. Not used if data is provided :type data: dict :param data: The dictionary that defines the force field. If not provided, the database will be read from path. :raise ValueError: if neither path nor data have been provided """ self.ff_name = None if not data: if not path: msg = ('Either a path or data must be provided.') raise ValueError(msg) data = load_force_field_parameters(path) self.ff_name = get_force_field_name(path) self.data = data
[docs] def getKeyData(self, key, fail_on_missing=True): """ Get the force field data for key :type key: str :param key: A top level key for the force field data dictionary :type fail_on_missing: bool :param fail_on_missing: True if an exception should be raised if the key is not in the force field data. If False, None will be returned for a missing key :return: the value associated with key, or None if key is missing and fail_on_missing is False :raise MissingParameterError: If key is missing and fail_on_missing is True """ try: return self.data[key] except KeyError: if fail_on_missing: msg = 'No force field data found for "%s"' % key raise MissingParameterError(msg) else: return None
[docs] def defineInternals(self, struct): """ Define the internal coordinates (bonds, angles, dihedrals, included and excluded nonbonds, etc) for the structure :type struct: `schrodinger.structure.Structure` :param struct: The structure to use when defining internals. """ exclude = self.getKeyData(EXCLUSIONS_INTERNAL_KEY) self.internals = coarsegrain.Internals(struct, uniqueify=True) self.internals.setSites() self.internals.setBonds() self.internals.setIncludedExcludedNonBonds(exclude_what=exclude) bond_types = [ get_type_name(atype) for atype in self.internals.bonds.keys() ] if not bond_types: return adict = { ANGLES_INTERNAL_KEY: 'setAngles', DIHEDRALS_INTERNAL_KEY: 'setDihedrals', IMPROPERS_INTERNAL_KEY: 'setImpropers' } for attr, setter in adict.items(): internal_types = self.data.get(attr, {}) if not internal_types: continue for internal_type in internal_types.keys(): for bond_type in bond_types: parameterized = bond_type in internal_type if parameterized: break if parameterized: break else: continue getattr(self.internals, setter)()
[docs] def getDatabaseData(self, names, key, fail_on_missing=True): """ Get the database data of type key for the particles with the given names. The names be ordered in the same way they would be ordered by the coarsegrain.Internals class and the proper number of names must be provided for the desired type - two for bonds, three for angles, etc. :type names: tuple :param names: The names of the particles to retrieve data for :type key: str :param key: One of the top level keys in the database dictionary. Should be one of the `*_INTERNAL_KEY` constants :type fail_on_missing: bool :param fail_on_missing: Whether missing database data should cause a MissingParameterError. If False, no error is raised and None is returned :rtype: dict or None :return: A dictionary of database data. The key is the function type and the value is a dictionary of parameter keys and values. None is returned if no data exists and fail_on_missing is False. :raise MissingParameterError: If fail_on_missing is True and there is no data """ tag = get_type_name(names) try: pdata = self.data[key][tag] except KeyError: if fail_on_missing: msg = 'No force field data found for %s %s' % (key, tag) raise MissingParameterError(msg) else: pdata = None return pdata
[docs] def getFFIOParam(self, names, key, fail_on_missing=True): """ Get the FFIOParam object for the given particle names and key :type names: tuple :param names: The names of the particles to retrieve data for :type key: str :param key: One of the top level keys in the database dictionary. Should be one of the `*_INTERNAL_KEY` constants that is a key in the FP_CLASSES dictionary. :type fail_on_missing: bool :param fail_on_missing: Whether missing database data should cause a MissingParameterError. If False, no error is raised and None is returned :rtype: `FFIOParam` or None :return: The FFIOParam object that handles the data for this internal. None is returned if no data exists and fail_on_missing is False. """ pdata = self.getDatabaseData(names, key, fail_on_missing=fail_on_missing) if pdata: ff_type = self.getKeyData(FF_TYPE_INTERNAL_KEY, fail_on_missing=False) if not ff_type: ff_type = GENERAL_TYPE return self.FP_CLASSES[key](pdata, ff_type=ff_type)
def _getTimeStep(self, fst): """ Return the recommended time step for the loaded structure with the loaded FF. :type fst: `CGFFIOStructure` :param fst: The FFIO structure containing the required data :rtype: float :return: time step in fs """ bond_data = _get_site_fc_data_from_bonds(fst) nonbond_type = self.getNonBondType(self.data) if is_like_lennard_jones(nonbond_type): if bond_data: time_step = min(15. * _get_approx_time_step(bond_data), 30.) else: time_step = 30. else: # build deduped site_by_vtype lookup dict sites_by_vtype = {} for site in fst.ffio_block.ffio.site: if site.vdwtype not in sites_by_vtype: sites_by_vtype[site.vdwtype] = site nonbond_data = [] if len(fst.ffio_block.ffio.site) == 1: site_i = list(sites_by_vtype.values())[0] nonbond_data.append( (site_i, site_i, coarsegrain.DEFAULT_A_PARAMETER)) else: for atype in fst.ffio_block.ffio.vdwtypescombined: site_i = sites_by_vtype[atype.name1] site_j = sites_by_vtype[atype.name2] names = (site_i.site, site_j.site) pdata = self.getDatabaseData(names, NONBONDS_INTERNAL_KEY) pfunc, pdict = next(iter(pdata.items())) nonbond_data.append( (site_i, site_j, pdict[VDW_A_PARAMETER_KEY])) time_step = 4. * _get_approx_time_step(bond_data + nonbond_data) return time_step
[docs] def addProperties(self, struct, fst): """ Add some force field properties to the structures. :type struct: `schrodinger.structure.Structure` :param struct: The structure to add properties to :type fst: `CGFFIOStructure` :param fst: The FFIO structure to add properties to """ # Cutoff cutoff = self.getKeyData(CUTOFF_INTERNAL_KEY) struct.property[coarsegrain.FF_CUTOFF_PROPERTY] = cutoff fst.ffio_block.property[coarsegrain.FF_CUTOFF_PROPERTY] = cutoff # Nonbond type (LJ, RH, SLJ) nonbond_type = self.getNonBondType(self.data) coarsegrain.set_nonbond_potential_type(struct, nonbond_type) coarsegrain.set_nonbond_potential_type(fst.ffio_block, nonbond_type) # Angle type adata = self.data.get(ANGLES_INTERNAL_KEY) if adata: angle_type = self.getAngleType(self.data) coarsegrain.set_angle_potential_type(struct, angle_type) coarsegrain.set_angle_potential_type(fst.ffio_block, angle_type) # DPD reduced gamma if nonbond_type == NONBONDS_DISSIPATIVE_PARTICLE_KEY: gamma = self.getKeyData(GAMMA_DISSIPATIVE_PARTICLE_INTERNAL_KEY) struct.property[DPD_REDUCED_GAMMA_PROP_KEY] = gamma fst.ffio_block.property[DPD_REDUCED_GAMMA_PROP_KEY] = gamma # time step time_step = self._getTimeStep(fst) struct.property[coarsegrain.TIME_STEP_PROPERTY] = time_step fst.ffio_block.property[coarsegrain.TIME_STEP_PROPERTY] = time_step
[docs] @staticmethod def getPotentialType(data, key, default): """ Return the most common potential type for this forcefield for the given key. :type data: dict :param data: the FF dictionary :type key: str :param key: the top level key of the block to inspect :type default: str :param default: a default potential type in case one is not present in the given FF dictionary :rtype: str :return: the potential type """ adata = data.get(key) if not adata: atype = default else: adict = {} for aitem in adata.values(): potential = list(aitem.keys())[0] adict.setdefault(potential, []).append(1) atype = max(list(adict.items()), key=lambda x: sum(x[1]))[0] return atype
[docs] @staticmethod def getSiteType(data): """ Determine the site type for this forcefield by inspecting one of the site parameter blocks. We assume all blocks use the same type so an arbitrary block is chosen. :type data: dict :param data: the FF dictionary :rtype: str :return: the site type such as SITES_GENERAL_ATOM_KEY """ return ForceField.getPotentialType(data, SITES_INTERNAL_KEY, SITES_GENERAL_ATOM_KEY)
[docs] @staticmethod def getAngleType(data): """ Determine the angle type for this forcefield by inspecting one of the angle parameter blocks. We assume all blocks use the same type so an arbitrary block is chosen. :type data: dict :param data: the FF dictionary :rtype: str :return: the angle type such as ANGLES_HARMONIC_KEY """ return ForceField.getPotentialType(data, ANGLES_INTERNAL_KEY, ANGLES_HARMONIC_KEY)
[docs] @staticmethod def getDihedralType(data): """ Determine the dihedral type for this forcefield by inspecting one of the dihedral parameter blocks. We assume all blocks use the same type so an arbitrary block is chosen. :type data: dict :param data: the FF dictionary :rtype: str :return: the dihedral type such as DIHEDRALS_OPLS_PROPER_KEY """ return ForceField.getPotentialType(data, DIHEDRALS_INTERNAL_KEY, DIHEDRALS_OPLS_PROPER_KEY)
[docs] @staticmethod def getNonBondType(data): """ Determine the nonbond type for this forcefield by inspecting one of the nonbond parameter blocks. We assume all blocks use the same type so an arbitrary block is chosen. :type data: dict :param data: the FF dictionary :rtype: str :return: the nonbond type such as coarsegrain.LENNARD_JONES """ return ForceField.getPotentialType(data, NONBONDS_INTERNAL_KEY, coarsegrain.LENNARD_JONES)
def _applySites(self, fst): """ Apply the force field site data :type fst: `CGFFIOStructure` :param fst: The FFIO structure to apply the data to :rtype: dict :return: Keys are particle names, values are `VDWData` namedtuples :raise MissingParameterError: If required data is missing from the force field """ ff_type = self.getKeyData(FF_TYPE_INTERNAL_KEY, fail_on_missing=False) if not ff_type: ff_type = GENERAL_TYPE dielectric = self.getKeyData(DIELECTRIC_INTERNAL_KEY) sites = {} vdw_sites = {} for names, (members, unused) in self.internals.sites.items(): # For sites, names is (name, ) name = names[0] # Add a vdw type for this site vtype = get_vdw_type_name(ff_type, name) vt_param = FFIOVdwType(None) vt_param.addToFFST(fst, vtype) fparam = self.getFFIOParam(names, SITES_INTERNAL_KEY) fparam.setDielectricConstant(dielectric) for member in members: sites[member[0]] = (name, fparam ) # name is not an attr of fparam # Save some of this data for later use with VDWTypeCombined blocks vdw_sites[name] = VDWData(vtype=vtype, charge=fparam.charge) # Now add a ffio_sites block, ffio sites and CT particles must have the # same orderings (see MATSCI-4261) for idx in range(1, len(sites) + 1): name, fparam = sites[idx] vtype = get_vdw_type_name(ff_type, name) fparam.addToFFST(fst, name, vtype) return vdw_sites def _applyBAI(self, fst): """ Apply the force field bond, angle and improper data :type fst: `CGFFIOStructure` :param fst: The FFIO structure to apply the data to :raise MissingParameterError: If required bond data is missing from the force field """ # Bonds, angles, and impropers can all be handled in a simple manner gen_types = ['bonds', 'angles', 'impropers'] gen_keys = [ BONDS_INTERNAL_KEY, ANGLES_INTERNAL_KEY, IMPROPERS_INTERNAL_KEY ] data_required = {'bonds'} for gen_type, key in zip(gen_types, gen_keys): needed = gen_type in data_required idata = getattr(self.internals, gen_type) if idata is None: continue for names, (members, unused) in idata.items(): fparam = self.getFFIOParam(names, key, fail_on_missing=needed) if fparam: for member in members: fparam.addToFFST(fst, *member) def _applyDihedrals(self, fst): """ Apply the force field dihedral data :type fst: `CGFFIOStructure` :param fst: The FFIO structure to apply the data to """ # trigonometric is handled differently as here we are potentially # adding multiple instances of the same dihedral quadruple # to the ffio block with different c-parameters (this is interpreted # as a linar combination) idata = self.internals.dihedrals if idata is None: return for names, (members, unused) in idata.items(): pdata = self.getDatabaseData(names, DIHEDRALS_INTERNAL_KEY, fail_on_missing=False) if pdata: function, pdict = list(pdata.items())[0] if function == DIHEDRALS_OPLS_PROPER_KEY: fparam = FFIODihedral(pdata) for member in members: fparam.addToFFST(fst, *member) elif function == DIHEDRALS_TRIGONOMETRIC_PROPER_KEY: for triple in DIHEDRAL_TRIG_TRIPLES: pairs = list(zip(DIHEDRAL_TRIG_DUMMY_TRIPLE, triple)) newpdict = dict([(x, pdict[y]) for x, y in pairs]) if newpdict[DIHEDRAL_FORCE_CONSTANT_DUMMY_KEY] == 0: continue newpdata = {function: newpdict} fparam = FFIODihedral(newpdata) for member in members: fparam.addToFFST(fst, *member) def _applyExclusions(self, fst): """ Apply the exclusions to the fst structure :type fst: `CGFFIOStructure` :param fst: The FFIO structure to apply the data to """ # There is no database info for exclusions excluded_nonbonds = self.internals.excluded_nonbonds for names, (members, unused) in excluded_nonbonds.items(): fparam = FFIOExclusion(None) for member in members: fparam.addToFFST(fst, *member) def _applyVdwCombined(self, fst, vdw_sites): """ Apply the force field vdw combined (nonbond) data :type fst: `CGFFIOStructure` :param fst: The FFIO structure to apply the data to :type vdw_sites: dict :param vdw_sites: Keys are particle names, values are `VDWData` namedtuples :raise MissingParameterError: If required data is missing from the force field """ included_nonbonds = self.internals.included_nonbonds for names, (members, unused) in included_nonbonds.items(): dielectric = self.getKeyData(DIELECTRIC_INTERNAL_KEY) cutoff = self.getKeyData(CUTOFF_INTERNAL_KEY) site1, site2 = [vdw_sites[x] for x in names] type1, charge1 = site1.vtype, site1.charge type2, charge2 = site2.vtype, site2.charge pdata = self.getDatabaseData(names, NONBONDS_INTERNAL_KEY) fparam = FFIOVdwTypesCombined(cutoff, dielectric, charge1, charge2, pdata) fparam.addToFFST(fst, type1, type2)
[docs] @staticmethod def definedShadowBonds(data): """ A generator over all defined shadow bonds in the force field. A shadow bond is a bond that doesn't physically exist in the structure but exists in the forcefield as a bonding attraction between two particles. :type data: dict :param data: The dictionary defining the force field, such as provided by load_force_field_parameters. :rtype: `ShadowData` :return: A ShadowData namedtuple that provides the names of the two particles and the number of bonds that seperate them. """ # Note: Shadow bond info looks like: # ffdata = {SHADOWBOND_TYPE_INTERNAL_KEY: # {'particle_a,particle_b': [3, 5], # 'particle_a,particle_c': [4]}} # Defining 3 shadow bond types - two between a & b separated by either 3 # or 5 bonds, and one between a & c separated by 4 bonds. shadow_info = data.get(SHADOWBOND_TYPE_INTERNAL_KEY, {}) for typename, depths in shadow_info.items(): name1, name2 = split_type_name(typename) for depth in depths: yield ShadowData(name1=name1, name2=name2, depth=depth)
def _findShadowBonds(self): """ Find all the pairs that are not actually bonded but are considered bonded because the user defined them as such in the Force Field panel. """ exclude = self.getKeyData(EXCLUSIONS_INTERNAL_KEY) for sinfo in self.definedShadowBonds(self.data): name1, name2, depth = sinfo.name1, sinfo.name2, sinfo.depth self.internals.findNewBonds(name1, name2, depth, exclude_what=exclude)
[docs] def applyFF(self, struct, internals=None, filename=None, ff_name=None): """ Apply the force field to the structure :type struct: `schrodinger.structure.Structure` :param struct: The structure to apply the force field :type internals: `schrodinger.appliction.matsci.coarsegrain.Internals` :param internals: The set of internal coordinates computed from the structure. If not provided, they will be computed. :type filename: str :param filename: If provided, a Cms structure with the applied force field will be written to this file path :type ff_name: str or None :param ff_name: the name of the forcefield used if self.ff_name is None, if None and self.ff_name is None then the default of FFIO_NAME will be used :rtype: `schrodinger.application.desmond.cms.Cms` :return: The CMS structure with the force field applied :raise MissingParameterError: If required data is missing from the force :raise TypeError: If struct is not a coarse grain structure :raise KeyError: If struct does not have chorus box PBC properties """ if not msutils.is_coarse_grain(struct): raise TypeError('Coarse-grained force fields can only be applied ' 'to coarse-grained structures.') # PBC_POSITION_KEY defines PBC position pbc_position = struct.property.get(xtal.PBC_POSITION_KEY) # Make a copy because we add properties to it struct = struct.copy() if pbc_position: struct.property[xtal.PBC_POSITION_KEY] = pbc_position if self.ff_name is None: if ff_name is None: self.ff_name = FFIO_NAME else: self.ff_name = ff_name # Create the FFIO structure and compute internals fst = CGFFIOStructure(struct) fst.setName(name=self.ff_name) fst.setCombiningRule() description = self.getKeyData(DESCRIPTION_INTERNAL_KEY, fail_on_missing=False) if description: fst.setDescription(description=description) ff_type = self.getKeyData(FF_TYPE_INTERNAL_KEY, fail_on_missing=False) fst.setFFType(ff_type=ff_type) if internals: self.internals = internals else: self.defineInternals(struct) # Detect shadow bonds self._findShadowBonds() # Apply the force field vdw_sites = self._applySites(fst) self._applyBAI(fst) self._applyDihedrals(fst) self._applyExclusions(fst) self._applyVdwCombined(fst, vdw_sites) # Adding properties to the structure before generating the FFIO ensures # that both the full system and solute system structures have those # properties but here some properties need to be created # after the FFIO exists and so it is currently necessary to pass in # both objects as the st in fst necessarily contains an empty ffio # block that should not be written to file self.addProperties(struct, fst) # Write a CMS if requested struct ends up as the full system CT in the # .cms file and fst ends up as the solute CT if filename: fst.writeCMS(struct, filename) cmst = cms.Cms(filename) else: # The tempfilename context manager returns the file name, not the # open file. This was we can use it with structure writers that want # to open the file themselves. On exit, it does a force_remove on # that filename. with fileutils.tempfilename(suffix='.cms') as filename: fst.writeCMS(struct, filename) cmst = cms.Cms(filename) return cmst
[docs]def create_cg_system(struct, ff_path, filename=None, pbc_position=None): """ Apply the force field located at ff_path to struct and return the resulting CMS object :type struct: `schrodinger.structure.Structure` :param struct: The structure to apply the forcefield to :type ff_path: str :param ff_path: The path to the force field file :type filename: str :param filename: If given, the CMS structure will be written to this file :type pbc_position: str or None :param pbc_position: Set the pbc position. If not None, pbc_position must be either xtal.CENTER_PBC_POSITION or xtal.ANCHOR_PBC_POSITION :rtype: `schrodinger.application.desmond.cms.Cms` :return: The Cms object that results from applying the force field to the structure """ if isinstance(struct, cms.Cms): # Since we directly write the passed in structure, we have to avoid # using a Cms because that will write the Cms substructures as well struct = struct.fsys_ct.copy() if pbc_position: struct.property[xtal.PBC_POSITION_KEY] = pbc_position # sync_pbc returns a copy of struct struct = xtal.sync_pbc(struct, create_pbc=True) forcefield = ForceField(path=ff_path) return forcefield.applyFF(struct, filename=filename)
[docs]def validate_options(parser, args, ff_flag, cg_ff_loc_type_flag, input_file_flag=None, input_sts=None): """ Validate options. :type parser: argparse.ArgumentParser :param parser: argument parser to validate :type args: list :param args: the arguments to parse :type ff_flag: str :param ff_flag: the force field flag :type cg_ff_loc_type_flag: str :param cg_ff_loc_type_flag: the coarse-grained force field location type flag :type input_file_flag: str or None :param input_file_flag: the input file flag or None if using input_sts :type input_sts: list or None :param input_sts: the input structures or None if using input_file_flag :rtype: named tuple :return: argparse parsed specified options """ options = parser.parse_args(args) if input_file_flag: input_file = jobutils.get_option(options, input_file_flag) else: input_file = None try: is_coarse_grain = coarsegrain.is_coarse_grain_input( input_file=input_file, input_sts=input_sts) except ValueError as err: parser.error(str(err)) ff = jobutils.get_option(options, ff_flag) if is_coarse_grain: cg_ff_loc_type = jobutils.get_option(options, cg_ff_loc_type_flag) try: get_force_field_file_path(ff, location_type=cg_ff_loc_type, local_type_includes_cwd=True, check_existence=True) except ValueError as err: parser.error(str(err)) else: try: # Validate the FF name and make sure OPLS3e is named correctly jobutils.set_option(options, ff_flag, parserutils.type_forcefield(ff)) except argparse.ArgumentTypeError as err: parser.error(str(err)) return options
[docs]def add_local_type_force_field_to_job_builder_w_check(builder, options, ff_flag, cg_ff_loc_type_flag, input_file_flag=None, input_sts=None): """ Check and set up the Launch API Job Builder to use the given local type force field. :type builder: `schrodinger.job.launchapi.JobSpecificationArgsBuilder` :param builder: the job builder to use for setting the FF as an input file :type options: named tuple :param options: argparse parsed specified options :type ff_flag: str :param ff_flag: the force field flag :type cg_ff_loc_type_flag: str :param cg_ff_loc_type_flag: the coarse-grained force field location type flag :type input_file_flag: str or None :param input_file_flag: the input file flag or None if using input_sts :type input_sts: list or None :param input_sts: the input structures or None if using input_file_flag """ if input_file_flag: input_file = jobutils.get_option(options, input_file_flag) else: input_file = None is_coarse_grain = coarsegrain.is_coarse_grain_input(input_file=input_file, input_sts=input_sts) cg_ff_loc_type = jobutils.get_option(options, cg_ff_loc_type_flag) if is_coarse_grain and \ cg_ff_loc_type == parserutils.LOCAL_CG_FF_LOCATION_TYPE: forcefield = jobutils.get_option(options, ff_flag) add_local_type_force_field_to_job_builder(builder, forcefield)