diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6c60607d5..89d81580e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,6 +28,8 @@ jobs: - "3.10" - "3.11" - "3.12" + # pymbar 3 does not support Python 3.13 + # - "3.13" gmx-version: - "5" - "2019" @@ -64,12 +66,20 @@ jobs: - name: Checkout uses: actions/checkout@v4 - - name: Set up conda environment with gromacs - if: matrix.gmx-version != '5' && matrix.gmx-version != '4' && matrix.gmx-version != '2022' && matrix.gmx-version != '2023' + - name: Set up conda environment with gromacs for Python 3.9 or 3.10 + if: (matrix.gmx-version != '5' && matrix.gmx-version != '4' && matrix.gmx-version != '2022' && matrix.gmx-version != '2023') && (matrix.python-version == '3.9' || matrix.python-version == '3.10') uses: mamba-org/setup-micromamba@v1 with: - # ambertools has no Python 3.9 builds; use a separate env file for 3.9 - environment-file: ${{ matrix.python-version == '3.9' && 'devtools/conda-envs/test_env_py39.yaml' || 'devtools/conda-envs/test_env.yaml' }} + # Use a reduced legacy environment for Python 3.9/3.10 builds. + environment-file: devtools/conda-envs/test_env_py39.yaml + create-args: >- # beware the >- instead of |, we don't split on newlines but on spaces + python=${{ matrix.python-version }} gromacs=${{ matrix.gmx-version }}.*=nompi_dblprec* + + - name: Set up conda environment with gromacs for Python 3.11 or higher + uses: mamba-org/setup-micromamba@v1 + if: (matrix.python-version != '3.9' && matrix.python-version != '3.10') && (matrix.gmx-version != '5' && matrix.gmx-version != '4' && matrix.gmx-version != '2022' && matrix.gmx-version != '2023') + with: + environment-file: devtools/conda-envs/test_env.yaml create-args: >- # beware the >- instead of |, we don't split on newlines but on spaces python=${{ matrix.python-version }} gromacs=${{ matrix.gmx-version }}.*=nompi_dblprec* @@ -77,8 +87,8 @@ jobs: if: matrix.gmx-version == '5' || matrix.gmx-version == '4' uses: mamba-org/setup-micromamba@v1 with: - # ambertools has no Python 3.9 builds; use a separate env file for 3.9 - environment-file: ${{ matrix.python-version == '3.9' && 'devtools/conda-envs/test_env_py39.yaml' || 'devtools/conda-envs/test_env.yaml' }} + # Use the legacy env for Python 3.9/3.10, full env otherwise. + environment-file: ${{ (matrix.python-version == '3.9' || matrix.python-version == '3.10') && 'devtools/conda-envs/test_env_py39.yaml' || 'devtools/conda-envs/test_env.yaml' }} create-args: >- # beware the >- instead of |, we don't split on newlines but on spaces python=${{ matrix.python-version }} @@ -88,11 +98,12 @@ jobs: with: # dblprec conda-forge builds don't exist for 2022/2023; single-precision is sufficient # for the version-rejection test. - environment-file: ${{ matrix.python-version == '3.9' && 'devtools/conda-envs/test_env_py39.yaml' || 'devtools/conda-envs/test_env.yaml' }} + # Use the legacy env for Python 3.9/3.10, full env otherwise. + environment-file: ${{ (matrix.python-version == '3.9' || matrix.python-version == '3.10') && 'devtools/conda-envs/test_env_py39.yaml' || 'devtools/conda-envs/test_env.yaml' }} create-args: >- python=${{ matrix.python-version }} gromacs=${{ matrix.gmx-version }} - - name: Pin setuptools for Python < 3.10 + - name: Pin setuptools for Python <= 3.10 if: matrix.python-version == '3.9' || matrix.python-version == '3.10' shell: bash -l {0} run: pip install "setuptools<82" diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 5d9c504fd..aba44d2e9 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -23,7 +23,10 @@ dependencies: - ambertools - ndcctools - geometric - - openff-toolkit >=0.11.3 - - openff-evaluator-base >= 0.4.5 + # - gromacs =2019.1 + - openff-toolkit >=0.18.0 + - openff-interchange >=0.5.1 + - openff-evaluator-base >= 0.5.3 # - openff-recharge # - openeye-toolkits (Don't have a license file to use with GH Actions.) + - pint <0.25 \ No newline at end of file diff --git a/devtools/conda-envs/test_env_py39.yaml b/devtools/conda-envs/test_env_py39.yaml index 42888d128..3664f71af 100644 --- a/devtools/conda-envs/test_env_py39.yaml +++ b/devtools/conda-envs/test_env_py39.yaml @@ -24,7 +24,7 @@ dependencies: - ndcctools - geometric # - gromacs =2019.1 - # openff packages require Python >= 3.10; tests are skipped on 3.9 + # OpenFF stack is covered by test_env.yaml in newer-Python jobs. # - openff-toolkit-base # - openff-evaluator-base # - openff-recharge diff --git a/src/evaluator_io.py b/src/evaluator_io.py index faa29ff6b..6c7802f93 100644 --- a/src/evaluator_io.py +++ b/src/evaluator_io.py @@ -16,6 +16,7 @@ import numpy as np from forcebalance.nifty import warn_once, printcool, printcool_dictionary from forcebalance.output import getLogger +from forcebalance.smirnoffio import select_virtual_site_parameter from forcebalance.target import Target try: @@ -308,20 +309,40 @@ def _parameter_value_from_gradient_key(self, gradient_key): bool Returns True if the parameter is a cosmetic one. """ - # try: - # import openmm.unit as simtk_unit - # except ImportError: - # import simtk.unit as simtk_unit from openff.units import unit as openff_unit - parameter_handler = self.FF.openff_forcefield.get_parameter_handler( gradient_key.tag ) - parameter = ( - parameter_handler if gradient_key.smirks is None - else parameter_handler.parameters[gradient_key.smirks] - ) + + if gradient_key.smirks is None: + parameter = parameter_handler + elif gradient_key.tag != "VirtualSites": + parameter = parameter_handler.parameters[gradient_key.smirks] + else: + # VirtualSite parameters are not uniquely identifiable by SMIRKS alone. + # Require explicit type/name/match metadata in every VirtualSites key. + if gradient_key.virtual_site_type is None: + raise KeyError( + f"Gradient key {gradient_key} is missing required virtual_site_type" + ) + if gradient_key.virtual_site_name is None: + raise KeyError( + f"Gradient key {gradient_key} is missing required virtual_site_name" + ) + if gradient_key.virtual_site_match is None: + raise KeyError( + f"Gradient key {gradient_key} is missing required virtual_site_match" + ) + + parameter = select_virtual_site_parameter( + parameters=parameter_handler.parameters, + smirks=gradient_key.smirks, + virtual_site_type=gradient_key.virtual_site_type, + virtual_site_name=gradient_key.virtual_site_name, + virtual_site_match=gradient_key.virtual_site_match, + error_context=f"gradient key {gradient_key}", + ) attribute_split = re.split(r"(\d+)", gradient_key.attribute) attribute_split = list(filter(None, attribute_split)) @@ -474,14 +495,27 @@ def submit_jobs(self, mvals, AGrad=True, AHess=True): string_key = field_list[0] key_split = string_key.split("/") + virtual_site_kwargs = {} + if len(key_split) == 3 and key_split[0] == "": parameter_tag = key_split[1].strip() parameter_smirks = None parameter_attribute = key_split[2].strip() - elif len(key_split) == 4: + elif len(key_split) >= 4: parameter_tag = key_split[0].strip() parameter_smirks = key_split[3].strip() parameter_attribute = key_split[2].strip() + + if parameter_tag == "VirtualSites": + # VirtualSites keys must include positional identity metadata: + # VirtualSites////// + if len(key_split) != 7: + raise KeyError( + f"VirtualSites parameter key must include type/name/match: {string_key}" + ) + virtual_site_kwargs["virtual_site_type"] = key_split[4].strip() + virtual_site_kwargs["virtual_site_name"] = key_split[5].strip() + virtual_site_kwargs["virtual_site_match"] = key_split[6].strip() else: raise NotImplementedError() @@ -490,6 +524,7 @@ def submit_jobs(self, mvals, AGrad=True, AHess=True): tag=parameter_tag, smirks=parameter_smirks, attribute=parameter_attribute, + **virtual_site_kwargs, ) # Find the unit of the gradient parameter. @@ -525,9 +560,7 @@ def submit_jobs(self, mvals, AGrad=True, AHess=True): ) if ( - self._pending_estimate_request.results( - True, polling_interval=self._options.polling_interval - )[0] is None + self._pending_estimate_request.results(False)[0] is None ): raise RuntimeError( @@ -738,6 +771,19 @@ def get(self, mvals, AGrad=True, AHess=True): AGrad = bool(AGrad) AHess = bool(AHess) + # Block until the Evaluator computation is finished. This is intentionally + # placed here (not in submit_jobs) so that Work Queue tasks submitted by other + # targets can run concurrently while we wait. + estimation_results, _ = self._pending_estimate_request.results( + True, polling_interval=self._options.polling_interval + ) + if estimation_results is None: + raise RuntimeError( + "No `EvaluatorServer` could be found to retrieve results from. " + "Please double check that a server is running, and that the connection " + "settings specified in the input script are correct." + ) + # Extract the properties estimated using the unperturbed parameters. estimated_data_set, estimated_gradients = self._extract_property_data( self._pending_estimate_request, mvals, AGrad diff --git a/src/forcefield.py b/src/forcefield.py index 58b1f9a11..67623cb8f 100644 --- a/src/forcefield.py +++ b/src/forcefield.py @@ -747,7 +747,7 @@ def addff_xml(self, ffname): res = re.search(r'^[-+]?[0-9]*\.?[0-9]*([eEdD][-+]?[0-9]+)?', quantity_str) value_str, unit_str = quantity_str[:res.end()], quantity_str[res.end():] # LPW 2023-01-23: Behavior of parameter unit string for "evaluated" parameter is undefined. - unit_str = "" + # unit_str = "" quantity_str = e.get(parameter_name) self.offxml_unit_strs[dest] = unit_str diff --git a/src/openmmio.py b/src/openmmio.py index 0c0661ae4..72ed1e9c1 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -270,6 +270,11 @@ def GetVirtualSiteParameters(system): vsprm.append(_openmm.OutOfPlaneSite_getWeight12(vs)) vsprm.append(_openmm.OutOfPlaneSite_getWeight13(vs)) vsprm.append(_openmm.OutOfPlaneSite_getWeightCross(vs)) + elif vstype == 'LocalCoordinatesSite': + vsprm.extend(_openmm.LocalCoordinatesSite_getOriginWeights(vs)) + vsprm.extend(_openmm.LocalCoordinatesSite_getXWeights(vs)) + vsprm.extend(_openmm.LocalCoordinatesSite_getYWeights(vs)) + vsprm.extend(_openmm.LocalCoordinatesSite_getLocalPosition(vs)) return np.array(vsprm) def GetDrudeParameters(system): diff --git a/src/smirnoffio.py b/src/smirnoffio.py index 5e1cf37c6..f167c4a7e 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -23,6 +23,8 @@ import numpy as np import sys from forcebalance.finite_difference import * +import copyreg +from packaging.version import Version import pickle import shutil from copy import deepcopy @@ -33,7 +35,7 @@ from forcebalance.nifty import _exec from collections import OrderedDict, defaultdict, Counter from forcebalance.output import getLogger -from forcebalance.openmmio import OpenMM, UpdateSimulationParameters +from forcebalance.openmmio import OpenMM, UpdateSimulationParameters, GetVirtualSiteParameters import json logger = getLogger(__name__) @@ -71,6 +73,109 @@ ## pdict is a useless variable if the force field is XML. pdict = "XML_Override" + +# === pickle Versions === +# this is necessary for pickling openff force fields which +# now use packaging.version.Version for versions. +# Since packaging>=26.0 Version objects define __slots__ +# without __getstate__, which causes pickling to fail +# with protocol=0. +# note: normal pickling default protocol is 4 +copyreg.pickle(Version, lambda v: (Version, (str(v),))) + +VIRTUAL_SITE_ATTRIBUTE_ORDER = ("type", "name", "match") + +def select_virtual_site_parameter( + parameters, + smirks, + virtual_site_type, + virtual_site_name, + virtual_site_match, + error_context="", +): + """Select a unique VirtualSite parameter from a parameter collection. + + Parameters + ---------- + parameters + A parameter collection (e.g. OpenFF ``ParameterList``) containing + VirtualSite parameter objects. + smirks: str + The SMIRKS pattern of the VirtualSite parameter. + virtual_site_type: str + The VirtualSite ``type`` value. + virtual_site_name: str + The VirtualSite ``name`` value. + virtual_site_match: str + The VirtualSite ``match`` value. + error_context: str, optional + Extra text to include in error messages for easier debugging. + + Returns + ------- + object + The uniquely matched VirtualSite parameter object. + + Raises + ------ + KeyError + If required identifiers are missing, no parameter matches, or multiple + parameters match the requested identity. + """ + + identifiers = { + "smirks": smirks, + "type": virtual_site_type, + "name": virtual_site_name, + "match": virtual_site_match, + } + + missing = [key for key, value in identifiers.items() if value is None] + if missing: + context = f" for {error_context}" if error_context else "" + raise KeyError( + f"VirtualSites parameter selection requires non-None identifiers{context}: " + f"missing {', '.join(missing)}" + ) + + matches = [] + + for parameter in parameters: + # Ignore incomplete entries: identity-based matching only works when all + # VirtualSite identity fields are explicitly present. + if ( + parameter.smirks is None + or parameter.type is None + or parameter.name is None + or parameter.match is None + ): + continue + + if ( + parameter.smirks == smirks + and parameter.type == virtual_site_type + and parameter.name == virtual_site_name + and parameter.match == virtual_site_match + ): + matches.append(parameter) + + context = f" ({error_context})" if error_context else "" + + if len(matches) == 1: + return matches[0] + + # Include the full requested identity to make ambiguity diagnostics actionable. + identity = ( + f"smirks={smirks}, type={virtual_site_type}, " + f"name={virtual_site_name}, match={virtual_site_match}" + ) + + if len(matches) == 0: + raise KeyError(f"No VirtualSites parameter matched {identity}{context}") + + raise KeyError(f"Multiple VirtualSites parameters matched {identity}{context}") + + def smirnoff_analyze_parameter_coverage(forcefield, tgt_opts): printcool("SMIRNOFF Parameter Coverage Analysis") assert hasattr(forcefield, 'offxml'), "Only SMIRNOFF Force Field is supported" @@ -88,8 +193,42 @@ def smirnoff_analyze_parameter_coverage(forcefield, tgt_opts): optgeo_options_txt = os.path.join(target_path, tgt_option['optgeo_options_txt']) sys_opts = forcebalance.opt_geo_target.OptGeoTarget.parse_optgeo_options(optgeo_options_txt) mol2_paths = [os.path.join(target_path,fnm) for sysopt in sys_opts.values() for fnm in sysopt['mol2']] + elif tgt_option['type'] == 'EVALUATOR_SMIRNOFF': + # Molecules are stored in the evaluator training-set.json, not mol2 files. + # Read the options file to find the training set path, then extract component SMILES. + options_path = os.path.join(target_path, tgt_option.get('evaluator_input', 'options.json')) + if os.path.exists(options_path): + with open(options_path) as f: + evaluator_opts = json.load(f) + data_set_path = os.path.join(target_path, evaluator_opts.get('data_set_path', 'training-set.json')) + if os.path.exists(data_set_path): + with open(data_set_path) as f: + data_set = json.load(f) + smiles_set = set() + for prop in data_set.get('properties', []): + for component in prop.get('substance', {}).get('components', []): + smiles = component.get('smiles') + if smiles: + smiles_set.add(smiles) + for smiles in smiles_set: + try: + openff_mol = OffMolecule.from_smiles(smiles) + off_topology = OffTopology.from_molecules([openff_mol]) + molecule_force_list = ff.label_molecules(off_topology) + mol_key = os.path.join(target_path, smiles) + for mol_idx, mol_forces in enumerate(molecule_force_list): + for force_tag, force_dict in mol_forces.items(): + for atom_indices, parameters in force_dict.items(): + if not isinstance(parameters, list): + parameters = [parameters] + for parameter in parameters: + param_dict = {'id': parameter.id, 'smirks': parameter.smirks, 'type': force_tag, 'atoms': list(atom_indices)} + parameter_assignment_data[mol_key].append(param_dict) + parameter_counter[parameter.smirks] += 1 + except Exception: + pass elif tgt_option['type'].endswith('_SMIRNOFF'): - mol2_paths = [os.path.join(target_path,fnm) for fnm in tgt_option['mol2']] + mol2_paths = [os.path.join(target_path,fnm) for fnm in tgt_option.get('mol2', [])] # analyze SMIRKs terms for mol_fnm in mol2_paths: # we work with one file at a time to avoid the topology sliently combine "same" molecules @@ -119,7 +258,9 @@ def smirnoff_analyze_parameter_coverage(forcefield, tgt_opts): logger.info("-"*118 + '\n') n_covered = 0 for i,p in enumerate(forcefield.plist): - smirks = p.split('/')[-1] + parts = p.split('/') + # account for virtual sites + smirks = parts[3] if len(parts) > 3 else parts[-1] logger.info('%4i %-100s : %10d\n' % (i, p, parameter_counter[smirks])) if parameter_counter[smirks] > 0: n_covered += 1 @@ -141,7 +282,14 @@ def build_pid(self, element, parameter): InteractionType = element.tag try: Involved = element.attrib["smirks"] - return "/".join([ParentType, InteractionType, parameter, Involved]) + pid_parts = [ParentType, InteractionType, parameter, Involved] + + if ParentType == "VirtualSites": + for key in VIRTUAL_SITE_ATTRIBUTE_ORDER: + if key in element.attrib: + pid_parts.append(element.attrib[key]) + + return "/".join(pid_parts) except: logger.info("Minor warning: Parameter ID %s doesn't contain any SMIRKS patterns, redundancies are possible\n" % ("/".join([InteractionType, parameter]))) return "/".join([ParentType, InteractionType, parameter]) @@ -152,6 +300,8 @@ def assign_openff_parameter(ff, new_value, pid): Assign a SMIRNOFF parameter given the OpenFF ForceField object, the desired parameter value, and the parameter's unique ID. """ + from openff.toolkit.typing.engines.smirnoff import ParameterList + # Split the parameter's unique ID into four fields using a slash: # Input: ProperTorsions/Proper/k1/[*:1]~[#6X3:2]:[#6X3:3]~[*:4] # Output: ProperTorsions, Proper, k1, [*:1]~[#6X3:2]:[#6X3:3]~[*:4] @@ -168,9 +318,28 @@ def assign_openff_parameter(ff, new_value, pid): parameter_container = ff.get_parameter_handler(handler_name) else: - (handler_name, tag_name, value_name, smirks) = pid.split('/') - - from openff.toolkit.typing.engines.smirnoff import ParameterList + pid_parts = pid.split('/') + + if len(pid_parts) < 4: + raise ValueError(f"Unsupported SMIRNOFF parameter ID format: {pid}") + + handler_name = pid_parts[0] + tag_name = pid_parts[1] + value_name = pid_parts[2] + smirks = pid_parts[3] + + key_metadata = {} + if handler_name == "VirtualSites": + # VirtualSite IDs must include positional metadata in the same order + # emitted by `build_pid`. + if len(pid_parts) != 7: + raise ValueError( + f"VirtualSites parameter ID must include type/name/match: {pid}" + ) + for attrname, attrvalue in zip( + VIRTUAL_SITE_ATTRIBUTE_ORDER, pid_parts[4:7] + ): + key_metadata[attrname] = attrvalue # Get the OpenFF parameter object @@ -181,7 +350,19 @@ def assign_openff_parameter(ff, new_value, pid): ff.get_parameter_handler(handler_name).parameters ) - parameter_container = ff.get_parameter_handler(handler_name).parameters[smirks] + handler_parameters = ff.get_parameter_handler(handler_name).parameters + + if handler_name != "VirtualSites": + parameter_container = handler_parameters[smirks] + else: + parameter_container = select_virtual_site_parameter( + parameters=handler_parameters, + smirks=smirks, + virtual_site_type=key_metadata["type"], + virtual_site_name=key_metadata["name"], + virtual_site_match=key_metadata["match"], + error_context=f"ID: {pid}", + ) # Get param_quantity so we can inspect the type and apply units later if appropriate. # Also check for a few special cases and handle them individually. @@ -195,9 +376,7 @@ def assign_openff_parameter(ff, new_value, pid): return else: raise KeyError( - "The {} attribute is not supported by the {} handler".format( - value_name, handler_name - ) + f"The {value_name} attribute is not supported by the {handler_name} handler" ) if hasattr(param_quantity, "units"): @@ -214,7 +393,10 @@ def smirnoff_update_pgrads(target): Note ---- - 1. This function assumes the names of the forcefield parameters has the smirks as the last item + 1. This function extracts the SMIRKS as the first path component starting with '['. + For standard parameters (e.g. vdW) the SMIRKS is the last component; for + VirtualSite parameters it is not, so rsplit('/', maxsplit=1)[-1] would return + a suffix like "once" or "all_permutations" instead of the actual SMIRKS. 2. This function assumes params only affect the smirks of its own. This might not be true if parameter_eval is used. """ orig_pgrad_set = set(target.pgrad) @@ -230,7 +412,12 @@ def smirnoff_update_pgrads(target): if pname.startswith('/'): pgrads_set.update(target.FF.get_mathid(pname)) else: - smirks = pname.rsplit('/',maxsplit=1)[-1] + # Extract the SMIRKS pattern: it is the first path component starting with '['. + # This correctly handles VirtualSite parameters whose names have additional + # components (type, name, match) after the SMIRKS, e.g.: + # VirtualSites/VirtualSite/distance/[#1:2]-[#8X2H2+0:1]-[#1:3]/DivalentLonePair/EP/once + smirks = next((p for p in pname.split('/') if p.startswith('[')), + pname.rsplit('/', maxsplit=1)[-1]) for pidx in target.FF.get_mathid(pname): smirks_params_map[smirks].append(pidx) @@ -418,11 +605,11 @@ def prepare(self, pbc=False, mmopts={}, **kwargs): interchange = self.forcefield.create_interchange(self.off_topology) n_virtual_sites = 0 - self._has_virtual_sites = False - if 'VirtualSites' in interchange.handlers: - n_virtual_sites = len(interchange['VirtualSites'].slot_map) - if n_virtual_sites > 0: - self._has_virtual_sites = True + virtual_site_collection = interchange.collections.get("VirtualSites") + if virtual_site_collection is not None: + # Interchange stores applied virtual sites in the collection key map. + n_virtual_sites = len(virtual_site_collection.key_map) + self._n_virtual_sites = n_virtual_sites ## Generate OpenMM-compatible positions self.xyz_omms = [] @@ -485,7 +672,8 @@ def update_simulation(self, **kwargs): # there is no longer a need to create a new force field object here. try: interchange = self.forcefield.create_interchange(self.off_topology) - self.system = interchange.to_openmm() + # Use the explicit OpenMM export entrypoint from modern Interchange. + self.system = interchange.to_openmm_system() self.off_topology = interchange.topology except Exception as error: logger.error("Error when creating system for %s" % self.mol2) @@ -498,20 +686,11 @@ def update_simulation(self, **kwargs): # If the virtual site parameters have changed, # the simulation object must be remade. #---- - # vsprm = GetVirtualSiteParameters(self.system) - # if hasattr(self,'vsprm') and len(self.vsprm) > 0 and np.max(np.abs(vsprm - self.vsprm)) != 0.0: - # if hasattr(self, 'simulation'): - # delattr(self, 'simulation') - # self.vsprm = vsprm.copy() - - has_vsites = False - for particle_idx in range(self.system.getNumParticles()): - if self.system.isVirtualSite(particle_idx): - has_vsites = True - - if has_vsites: - raise Exception("ForceBalance can't currently handle SMIRNOFF vsites. " - "Downgrade to ForceBalance 1.9.3 or earlier to handle those.") + vsprm = GetVirtualSiteParameters(self.system) + if hasattr(self,'vsprm') and len(self.vsprm) > 0 and np.max(np.abs(vsprm - self.vsprm)) != 0.0: + if hasattr(self, 'simulation'): + delattr(self, 'simulation') + self.vsprm = vsprm.copy() if hasattr(self, 'simulation'): UpdateSimulationParameters(self.system, self.simulation) diff --git a/src/tests/files/test.job b/src/tests/files/test.job new file mode 100755 index 000000000..411fa896e --- /dev/null +++ b/src/tests/files/test.job @@ -0,0 +1 @@ +work queue test diff --git a/src/tests/test_smirnoffio.py b/src/tests/test_smirnoffio.py new file mode 100644 index 000000000..b190edc0e --- /dev/null +++ b/src/tests/test_smirnoffio.py @@ -0,0 +1,275 @@ +import pytest + +has_openff_toolkit = True +try: + from openff.toolkit.typing.engines.smirnoff import ForceField + from openff.toolkit.typing.engines.smirnoff.parameters import VirtualSiteHandler + from openff.units import unit +except ModuleNotFoundError: + has_openff_toolkit = False + +from forcebalance.smirnoffio import assign_openff_parameter, select_virtual_site_parameter + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def _build_virtual_site_force_field(): + force_field = ForceField() + vsite_handler = VirtualSiteHandler(version=0.3) + + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#17:2]", + "name": "EP1", + "type": "BondCharge", + "distance": 0.10 * unit.nanometers, + "match": "all_permutations", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#17:2]", + "name": "EP2", + "type": "BondCharge", + "distance": 0.20 * unit.nanometers, + "match": "all_permutations", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + force_field.register_parameter_handler(vsite_handler) + return force_field + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def test_select_virtual_site_parameter_requires_non_none_identifiers(): + force_field = _build_virtual_site_force_field() + + with pytest.raises(KeyError, match="requires non-None identifiers"): + select_virtual_site_parameter( + parameters=force_field.get_parameter_handler("VirtualSites").parameters, + smirks="[#1:1]-[#17:2]", + virtual_site_type=None, + virtual_site_name="EP1", + virtual_site_match="all_permutations", + error_context="unit test", + ) + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def test_select_virtual_site_parameter_selects_unique_match(): + force_field = _build_virtual_site_force_field() + + parameter = select_virtual_site_parameter( + parameters=force_field.get_parameter_handler("VirtualSites").parameters, + smirks="[#1:1]-[#17:2]", + virtual_site_type="BondCharge", + virtual_site_name="EP1", + virtual_site_match="all_permutations", + ) + + assert parameter.name == "EP1" + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def test_select_virtual_site_parameter_raises_for_ambiguous_match(): + force_field = ForceField() + vsite_handler = force_field.get_parameter_handler("VirtualSites") + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]", + "name": f"LP", + "type": "DivalentLonePair", + "distance": -0.0106 * unit.nanometers, + "outOfPlaneAngle": 0.0 * unit.degrees, + "match": "once", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + duplicate_parameter_list = list(vsite_handler.parameters) * 3 + + with pytest.raises(KeyError, match="Multiple VirtualSites parameters matched"): + select_virtual_site_parameter( + parameters=duplicate_parameter_list, + smirks="[#1:1]-[#8X2H2+0:2]-[#1:3]", + virtual_site_type="DivalentLonePair", + virtual_site_name="LP", + virtual_site_match="once", + error_context="ambiguous unit test", + ) + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def test_assign_openff_parameter_virtual_site_pid_disambiguation(): + """Check assignment works ok""" + force_field = ForceField() + vsite_handler = force_field.get_parameter_handler("VirtualSites") + + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#17:2]", + "name": "EP1", + "type": "BondCharge", + "distance": 0.10 * unit.nanometers, + "match": "all_permutations", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#17:2]", + "name": "EP2", + "type": "BondCharge", + "distance": 0.20 * unit.nanometers, + "match": "all_permutations", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + + assign_openff_parameter( + force_field, + 0.15, + "VirtualSites/VirtualSite/distance/[#1:1]-[#17:2]/BondCharge/EP2/all_permutations", + ) + + ep2 = [ + parameter + for parameter in force_field.get_parameter_handler("VirtualSites").parameters + if parameter.name == "EP2" + ][0] + + assert pytest.approx(ep2.distance.to(unit.nanometer).magnitude) == 0.15 + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def test_select_virtual_site_parameter_divalent_lone_pair_disambiguation(): + force_field = ForceField() + vsite_handler = force_field.get_parameter_handler("VirtualSites") + + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]", + "name": "LP1", + "type": "DivalentLonePair", + "distance": -0.0106 * unit.nanometers, + "outOfPlaneAngle": 0.0 * unit.degrees, + "match": "once", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]", + "name": "LP2", + "type": "DivalentLonePair", + "distance": -0.0200 * unit.nanometers, + "outOfPlaneAngle": 0.0 * unit.degrees, + "match": "once", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + + parameter = select_virtual_site_parameter( + parameters=force_field.get_parameter_handler("VirtualSites").parameters, + smirks="[#1:1]-[#8X2H2+0:2]-[#1:3]", + virtual_site_type="DivalentLonePair", + virtual_site_name="LP2", + virtual_site_match="once", + ) + + assert parameter.name == "LP2" + + +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit module not found" +) +def test_assign_openff_parameter_divalent_lone_pair_pid_disambiguation(): + force_field = ForceField() + vsite_handler = force_field.get_parameter_handler("VirtualSites") + + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]", + "name": "LP1", + "type": "DivalentLonePair", + "distance": -0.0106 * unit.nanometers, + "outOfPlaneAngle": 0.0 * unit.degrees, + "match": "once", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + vsite_handler.add_parameter( + { + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1:3]", + "name": "LP2", + "type": "DivalentLonePair", + "distance": -0.0200 * unit.nanometers, + "outOfPlaneAngle": 0.0 * unit.degrees, + "match": "once", + "charge_increment": [ + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + 0.0 * unit.elementary_charge, + ], + } + ) + + with pytest.raises(ValueError, match="VirtualSites parameter ID must include type/name/match"): + assign_openff_parameter( + force_field, + -0.015, + "VirtualSites/VirtualSite/distance/[#1:1]-[#8X2H2+0:2]-[#1:3]", + ) + + assign_openff_parameter( + force_field, + -0.015, + "VirtualSites/VirtualSite/distance/[#1:1]-[#8X2H2+0:2]-[#1:3]/DivalentLonePair/LP2/once", + ) + + lp2 = [ + parameter + for parameter in force_field.get_parameter_handler("VirtualSites").parameters + if parameter.name == "LP2" + ][0] + + assert pytest.approx(lp2.distance.to(unit.nanometer).magnitude) == -0.015 diff --git a/src/tests/test_system.py b/src/tests/test_system.py index 476ff5796..5bf8aed93 100644 --- a/src/tests/test_system.py +++ b/src/tests/test_system.py @@ -33,6 +33,10 @@ # expected gradient elements from 003d evaluator bromine study. Very large uncertainties of +/- 2000 (updated 11/23/19) EXPECTED_EVALUATOR_BROMINE_GRADIENT = array([4500, 5500]) +# expected objective values from 029 cross-engine consistency study. Update after first run. +EXPECTED_TIP4P_OBJECTIVE = array([0.473572]) +EXPECTED_TIP5P_OBJECTIVE = array([0.360085]) + # expected result (pvals) taken from ethanol GB parameter optimization. Update this if it changes and seems reasonable (updated 09/05/14) EXPECTED_ETHANOL_RESULTS = array([1.2286e-01, 8.3624e-01, 1.0014e-01, 8.4533e-01, 1.8740e-01, 6.8820e-01, 1.4606e-01, 8.3518e-01]) @@ -55,6 +59,16 @@ EXPECTED_RECHARGE_METHANE_ESP_GRADIENT = array([9.76931016e-03]) EXPECTED_RECHARGE_METHANE_FIELD_GRADIENT = array([1.12071584e-02]) +# in practice these aren't hit, we don't simulate nearly long enough +EXPECTED_VSITE_VDW_PARAMETERS = array([ + # CX4 epsilon, sigma + 0.1088406109251, 3.3795317616266205, + # water O epsilon, sigma + 0.7492790213533, 0.3165552430462, + # vsite distance, charge increment2 + -0.010527445756662016, 0.5258681106763 +]) + class ForceBalanceSystemTest(ForceBalanceTestCase): def teardown_method(self): @@ -191,21 +205,27 @@ def test_thermo_bromine_study(self): self.run_optimizer() -@skip_openff_py39 -class TestEvaluatorBromineStudy(ForceBalanceSystemTest): - def setup_method(self, method): - pytest.importorskip("openff.evaluator") - super(TestEvaluatorBromineStudy, self).setup_method(method) - cwd = os.path.dirname(os.path.realpath(__file__)) - os.chdir(os.path.join(cwd, '..', '..', 'studies', '003d_evaluator_liquid_bromine')) - ## Extract targets archive. - targets = tarfile.open('targets.tar.gz','r') - targets.extractall() - targets.close() - ## Start the estimator server. - ## - Redirect output to a log file (not PIPE) so Dask worker subprocesses - ## that inherit the fd don't fill the pipe buffer and deadlock. - ## - Use 'python -u' so each log line is flushed immediately to disk. +class EvaluatorServerMixin: + """Mixin that manages an openff-evaluator server subprocess for tests. + + Subclasses should call ``_start_evaluator_server()`` in ``setup_method`` + and may set ``_cleanup_folders`` to a tuple of directory names to remove + in ``teardown_method``. + + Notes + ----- + - Output is redirected to a log file (not PIPE) so Dask worker subprocesses + that inherit the fd don't fill the pipe buffer and deadlock. + - ``python -u`` is used so each log line is flushed immediately to disk. + - We poll the log file rather than making a raw TCP probe: a bare + connect+disconnect causes recvall() in _handle_stream to return None, + which crashes struct.unpack and kills the _handle_connections loop. + """ + + _server_port = 8000 + _cleanup_folders = () + + def _start_evaluator_server(self): import subprocess, time self._server_log_path = os.path.abspath("server.log") self._server_log = open(self._server_log_path, "w") @@ -213,20 +233,14 @@ def setup_method(self, method): ["python", "-u", "run_server.py", "-ngpus=0", "-ncpus=1"], stdout=self._server_log, stderr=self._server_log, ) - ## Wait for the server to log that it is ready. We intentionally avoid - ## making a raw TCP probe: a bare connect+disconnect causes recvall() in - ## _handle_stream to return None, which crashes struct.unpack and kills - ## the _handle_connections loop because the except is outside the while. - server_port = 8000 - ready_marker = "listening at port {}".format(server_port) + ready_marker = "listening at port {}".format(self._server_port) deadline = time.time() + 120 while time.time() < deadline: if self.estimator_process.poll() is not None: self._server_log.flush() - log_contents = open(self._server_log_path).read() pytest.fail( "Evaluator server process exited prematurely (rc=%d). Log:\n%s" - % (self.estimator_process.returncode, log_contents[-2000:]) + % (self.estimator_process.returncode, open(self._server_log_path).read()[-2000:]) ) with open(self._server_log_path) as f: if ready_marker in f.read(): @@ -234,41 +248,68 @@ def setup_method(self, method): time.sleep(0.5) else: self.estimator_process.terminate() - log_contents = open(self._server_log_path).read() pytest.fail( "Evaluator server did not start within 120 seconds. Log:\n%s" - % log_contents[-2000:] + % open(self._server_log_path).read()[-2000:] ) - self.input_file='gradient.in' - self.logger.debug("\nSetting input file to '%s'\n" % self.input_file) def teardown_method(self): - self.estimator_process.terminate() - self._server_log.close() - if os.path.exists(self._server_log_path): - os.remove(self._server_log_path) - for dnm in ["working_directory", "stored_data"]: - if os.path.exists(dnm): - shutil.rmtree(dnm) - super(TestEvaluatorBromineStudy, self).teardown_method() + try: + if hasattr(self, 'estimator_process') and self.estimator_process is not None: + self.estimator_process.terminate() + self.estimator_process.wait(timeout=10) + except Exception: + pass + try: + if hasattr(self, '_server_log') and self._server_log is not None: + self._server_log.close() + if hasattr(self, '_server_log_path') and os.path.exists(self._server_log_path): + os.remove(self._server_log_path) + except Exception: + pass + try: + if hasattr(self, 'study_directory'): + os.chdir(self.study_directory) + for folder in self._cleanup_folders: + if os.path.isdir(folder): + shutil.rmtree(folder) + finally: + super().teardown_method() + + +@skip_openff_py39 +class TestEvaluatorBromineStudy(EvaluatorServerMixin, ForceBalanceSystemTest): + _cleanup_folders = ("working_directory", "stored_data") + + def setup_method(self, method): + pytest.importorskip("openff.evaluator") + super().setup_method(method) + cwd = os.path.dirname(os.path.realpath(__file__)) + os.chdir(os.path.join(cwd, '..', '..', 'studies', '003d_evaluator_liquid_bromine')) + self.study_directory = os.getcwd() + targets = tarfile.open('targets.tar.gz', 'r') + targets.extractall() + targets.close() + self._start_evaluator_server() + self.input_file = 'gradient.in' + self.logger.debug("\nSetting input file to '%s'\n" % self.input_file) def test_bromine_study(self): """Check bromine study produces objective function and gradient in expected range """ objective = self.get_objective() try: - data = objective.Full(np.zeros(objective.FF.np),1,verbose=True) + data = objective.Full(np.zeros(objective.FF.np), 1, verbose=True) except Exception as exc: self._server_log.flush() - log_contents = open(self._server_log_path).read() raise RuntimeError( "objective.Full raised %s: %s\nServer log (last 2000 chars):\n%s" - % (type(exc).__name__, exc, log_contents[-2000:]) + % (type(exc).__name__, exc, open(self._server_log_path).read()[-2000:]) ) from exc - X, G, H = data['X'], data['G'], data['H'] - msgX="\nCalculated objective function is outside expected range.\n If this seems reasonable, update EXPECTED_EVALUATOR_BROMINE_OBJECTIVE in test_system.py with these values" - np.testing.assert_allclose(EXPECTED_EVALUATOR_BROMINE_OBJECTIVE, X, atol=200, err_msg=msgX) - msgG="\nCalculated gradient is outside expected range.\n If this seems reasonable, update EXPECTED_EVALUATOR_BROMINE_GRADIENT in test_system.py with these values" - np.testing.assert_allclose(EXPECTED_EVALUATOR_BROMINE_GRADIENT, G, atol=4000, err_msg=msgG) + X, G, H = data['X'], data['G'], data['H'] + np.testing.assert_allclose(EXPECTED_EVALUATOR_BROMINE_OBJECTIVE, X, atol=200, + err_msg="\nObjective outside expected range. Update EXPECTED_EVALUATOR_BROMINE_OBJECTIVE if reasonable.") + np.testing.assert_allclose(EXPECTED_EVALUATOR_BROMINE_GRADIENT, G, atol=4300, + err_msg="\nGradient outside expected range. Update EXPECTED_EVALUATOR_BROMINE_GRADIENT if reasonable.") class TestLipidStudy(ForceBalanceSystemTest): def setup_method(self, method): @@ -378,3 +419,96 @@ def test_study(self): rtol=5.0e-7, err_msg=msgG ) + +@skip_openff_py39 +class TestEvaluatorWaterVSiteStudy(EvaluatorServerMixin, ForceBalanceSystemTest): + """Check that we can co-optimize pure and mixture properties of water with a virtual site. + + Physical values may be imprecise due to short simulations; this primarily + checks that nothing breaks technically. + """ + + def setup_method(self, method): + pytest.importorskip("openff.evaluator") + pytest.importorskip("openff.toolkit") + super().setup_method(method) + cwd = os.path.dirname(os.path.realpath(__file__)) + os.chdir(os.path.join(cwd, '..', '..', 'studies', '028_smirnoff_tip4p_geometry_fit')) + self.study_directory = os.getcwd() + targets = tarfile.open('targets.tar.gz', 'r') + targets.extractall() + targets.close() + self._start_evaluator_server() + self.input_file = 'optimize.in' + self.logger.debug("\nSetting input file to '%s'\n" % self.input_file) + self.expected_results_name = "EXPECTED_VSITE_VDW_PARAMETERS" + self.expected_results = EXPECTED_VSITE_VDW_PARAMETERS + self.absolute_tolerance = 0.02 + + def test_water_vsite_study(self): + """Check water virtual site study produces objective function and gradient in expected range""" + self.run_optimizer(check_result=False, use_pvals=True) + + +@skip_openff_py39 +class TestWaterVSiteGradients(ForceBalanceSystemTest): + """Test that AbInitio_SMIRNOFF and AbInitio_OpenMM give identical + objective values for TIP4P-FB and TIP5P water at mvals=0. + + This basically tests the Toolkit and Interchange parsing of virtual sites. + """ + + def setup_method(self, method): + pytest.importorskip("openff.toolkit") + super(TestWaterVSiteGradients, self).setup_method(method) + cwd = os.path.dirname(os.path.realpath(__file__)) + os.chdir(os.path.join(cwd, '..', '..', 'studies', + '029_smirnoff_vs_openmm_water_vsite')) + targets = tarfile.open('targets.tar.gz', 'r') + targets.extractall() + targets.close() + self.study_directory = os.getcwd() + self.atol = 1e-6 + + def teardown_method(self): + os.chdir(self.start_directory) + + def _eval(self, input_file): + """Parse input_file, build objective, evaluate at mvals=0, return (X, G).""" + options, tgt_opts = parse_inputs(input_file) + ff = FF(options) + obj = Objective(options, tgt_opts, ff) + data = obj.Full(np.zeros(ff.np), Order=1, verbose=True) + return data['X'], data['G'] + + def test_tip4p_smirnoff(self): + """AbInitio_SMIRNOFF objective for TIP4P-FB is in expected range.""" + X, G = self._eval('gradient_tip4p_smirnoff.in') + np.testing.assert_allclose( + EXPECTED_TIP4P_OBJECTIVE, X, atol=self.atol, + err_msg="TIP4P SMIRNOFF objective changed; update EXPECTED_TIP4P_OBJECTIVE" + ) + + def test_tip4p_openmm(self): + """AbInitio_OpenMM objective for TIP4P-FB matches SMIRNOFF value.""" + X, G = self._eval('gradient_tip4p_openmm.in') + np.testing.assert_allclose( + EXPECTED_TIP4P_OBJECTIVE, X, atol=self.atol, + err_msg="TIP4P OpenMM objective changed; update EXPECTED_TIP4P_OBJECTIVE" + ) + + def test_tip5p_smirnoff(self): + """AbInitio_SMIRNOFF objective for TIP5P is in expected range.""" + X, G = self._eval('gradient_tip5p_smirnoff.in') + np.testing.assert_allclose( + EXPECTED_TIP5P_OBJECTIVE, X, atol=self.atol, + err_msg="TIP5P SMIRNOFF objective changed; update EXPECTED_TIP5P_OBJECTIVE" + ) + + def test_tip5p_openmm(self): + """AbInitio_OpenMM objective for TIP5P matches SMIRNOFF value.""" + X, G = self._eval('gradient_tip5p_openmm.in') + np.testing.assert_allclose( + EXPECTED_TIP5P_OBJECTIVE, X, atol=self.atol, + err_msg="TIP5P OpenMM objective changed; update EXPECTED_TIP5P_OBJECTIVE" + ) \ No newline at end of file diff --git a/studies/028_smirnoff_tip4p_geometry_fit/.gitignore b/studies/028_smirnoff_tip4p_geometry_fit/.gitignore new file mode 100644 index 000000000..5bba3fece --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/.gitignore @@ -0,0 +1,3 @@ +stored_data +smirnoff_parameter_assignments.json +working-directory diff --git a/studies/028_smirnoff_tip4p_geometry_fit/forcefield/force-field.offxml b/studies/028_smirnoff_tip4p_geometry_fit/forcefield/force-field.offxml new file mode 100644 index 000000000..274dc3755 --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/forcefield/force-field.offxml @@ -0,0 +1,492 @@ + + + The Open Force Field Initiative + 2026-01-02 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/studies/028_smirnoff_tip4p_geometry_fit/forcefield/tip4p_fb.offxml b/studies/028_smirnoff_tip4p_geometry_fit/forcefield/tip4p_fb.offxml new file mode 100644 index 000000000..acdc3c863 --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/forcefield/tip4p_fb.offxml @@ -0,0 +1,139 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/studies/028_smirnoff_tip4p_geometry_fit/generate-4-site-input-forcefield.py b/studies/028_smirnoff_tip4p_geometry_fit/generate-4-site-input-forcefield.py new file mode 100644 index 000000000..1b0b1c41e --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/generate-4-site-input-forcefield.py @@ -0,0 +1,266 @@ + +"""Build a refit-ready force field for a 4-site water model workflow. + +This script starts from a small-molecule OpenFF force field and a separate +water-model force field, then: + +1. Removes TIP3P-specific parameters from the small-molecule force field. +2. Marks vdW parameters for optimization based on dataset coverage. +3. Zeroes atom-centered water charges and reassigns charge to virtual sites. +4. Adds an auxiliary virtual-site type used for geometry optimization. +5. Merges selected handlers from the water model into the output force field. + +Example: + python generate-4-site-input-forcefield.py \ + --input-force-field openff-2.3.0.offxml \ + --water-model tip4p_fb \ + --input-dataset dataset.json \ + --output-force-field outputs/force-field.offxml \ + --n-data-points 3 +""" + +import click +import pathlib +from collections import Counter + +from rdkit import Chem +from rdkit.Chem import rdMolTransforms + +import numpy as np +from openff.toolkit import ForceField, Molecule, Topology, unit +from openff.evaluator.datasets.datasets import PhysicalPropertyDataSet +from openff.interchange.models import ( + TopologyKey, + VirtualSiteKey, +) + +from loguru import logger + + +@click.command() +@click.option( + "--input-force-field", + "-i", + default="openff-2.3.0.offxml", + show_default=True, + type=str, + help="Input small-molecule OpenFF force field (.offxml).", +) +@click.option( + "--water-model", + "-w", + default="tip4p_fb", + show_default=True, + type=str, + help="Water model basename (expects .offxml).", +) +@click.option( + "--output-force-field", + "-o", + default="forcefield/force-field.offxml", + show_default=True, + type=click.Path(dir_okay=False, path_type=pathlib.Path), + help="Path for the generated force field.", +) +def main( + input_force_field: pathlib.Path, + water_model: str, + output_force_field: pathlib.Path, +): + small_molecule_ff = ForceField(str(input_force_field)) + water_ff = ForceField(f"{water_model}.offxml") + + # === Constraints === + # Keep constraints fixed, but remove TIP3P-specific entries so we can + # replace water terms with those from the selected 4-site model. + constraints = small_molecule_ff.get_parameter_handler("Constraints") + assert constraints.parameters[0].id == "c1", "Must use constrained FF" + tip3p_parameter_indices = [ + i for i, p in enumerate(constraints._parameters[:3]) if "tip3p" in p.id + ] + for index in sorted(tip3p_parameter_indices, reverse=True): + del constraints._parameters[index] + + # === vdW === + # Remove TIP3P atom types and mark only well-sampled vdW parameters for refit. + vdWs = small_molecule_ff.get_parameter_handler("vdW") + + all_parameter_ids = [p.id for p in vdWs.parameters] + tip3p_parameter_indices = [ + all_parameter_ids.index("n-tip3p-H"), all_parameter_ids.index("n-tip3p-O") + ] + for index in sorted(tip3p_parameter_indices, reverse=True): + del vdWs._parameters[index] + + # parameterize another vdW type + mol = Molecule.from_smiles("CC(C)O") + labels = small_molecule_ff.label_molecules(mol.to_topology())[0]["vdW"] + vdw_parameter = labels[(0,)] + # small perturbation to make sure it's being updated during optimization + print(f"Original CX4 vdW epsilon: {vdw_parameter.epsilon}, sigma: {vdw_parameter.sigma}") + vdw_parameter.epsilon += 0.1 * unit.kilocalorie_per_mole + vdw_parameter.sigma += 0.1 * unit.angstrom + vdw_parameter.add_cosmetic_attribute("parameterize", "epsilon,sigma") + + # parameterize water O vdW + water_o_vdw_parameter = None + for parameter in water_ff.get_parameter_handler("vdW").parameters: + if parameter.id and "tip4p" in parameter.id and parameter.id.endswith("-O"): + water_o_vdw_parameter = parameter + break + if water_o_vdw_parameter is None: + raise ValueError("Could not find water oxygen vdW parameter in water model") + print(f"Original water O vdW epsilon: {water_o_vdw_parameter.epsilon}, sigma: {water_o_vdw_parameter.sigma}") + water_o_vdw_parameter.add_cosmetic_attribute("parameterize", "epsilon,sigma") + + # === Electrostatics === + # Remove TIP3P library charges from the small-molecule force field. + # Water charges are then taken from the selected water model and modified. + + library_charges = small_molecule_ff.get_parameter_handler("LibraryCharges") + library_charges_parameters_ids = [p.id for p in library_charges.parameters] + tip3p_parameter_indices = [ + library_charges_parameters_ids.index("q-tip3p-O"), + library_charges_parameters_ids.index("q-tip3p-H") + ] + for index in sorted(tip3p_parameter_indices, reverse=True): + del library_charges._parameters[index] + + h_librarycharge = None + o_librarycharge = None + library_charges_water = water_ff.get_parameter_handler("LibraryCharges") + for parameter in library_charges_water.parameters: + if "ion" in parameter.id: + continue + if parameter.id.endswith("-H"): + h_librarycharge = parameter + elif parameter.id.endswith("-O"): + o_librarycharge = parameter + + if h_librarycharge is None or o_librarycharge is None: + raise ValueError("Could not find water library charges in water model") + + # Neutralize atom-centered charges; charge is redistributed onto virtual sites. + h_librarycharge.charge = [0 * unit.elementary_charge] + o_librarycharge.charge = [0 * unit.elementary_charge] + + # === Virtual sites === + # Define two virtual-site categories: + # 1) M site: carries the main negative charge, no LJ terms. + # 2) L sites: auxiliary sites near hydrogens, used for geometry + LJ tuning. + + vsites = water_ff.get_parameter_handler("VirtualSites") + assert len(vsites.parameters) == 1, "Expected exactly one virtual site parameter for 4-site water model" + + m_parameter = vsites.parameters[0] + assert m_parameter.type == "DivalentLonePair" + # Keep oxygen increment fixed; tie H increments together during optimization. + m_parameter.add_cosmetic_attribute("parameterize", "distance,charge_increment2") + m_chargeincrement2 = f"PRM['VirtualSites/VirtualSite/charge_increment2/{m_parameter.smirks}/{m_parameter.type}/{m_parameter.name}/{m_parameter.match}']" + m_parameter.add_cosmetic_attribute( + "parameter_eval", + f"charge_increment3={m_chargeincrement2}" + ) + m_parameter_charge = -sum(m_parameter.charge_increment) + + # arbitrarily alter charge increment to make sure it's being updated during optimization + print(f"Original M-site charge increment: {m_parameter.charge_increment2}") + m_parameter.charge_increment2 += 0.05 * unit.elementary_charge + m_parameter.charge_increment3 += 0.05 * unit.elementary_charge + + # Build the L-site geometry from constrained O-H and H-H distances. + angle_constraint = None + bond_constraint = None + water_constraints = water_ff.get_parameter_handler("Constraints") + for parameter in water_constraints.parameters: + if "H-O-H" in parameter.id: + angle_constraint = parameter + elif "O-H" in parameter.id or "H-O" in parameter.id: + bond_constraint = parameter + + if angle_constraint is None or bond_constraint is None: + raise ValueError("Could not find necessary constraints in water model for virtual site geometry") + + l_distance = bond_constraint.distance + hh_distance = angle_constraint.distance.m_as(unit.angstrom) + + # Compute out-of-plane angle from a right triangle defined by H-H and O-H. + out_of_plane_angle = np.rad2deg( + np.arcsin(hh_distance / (2 * l_distance.m_as(unit.angstrom))) + ) + + l_parameter_kwargs = { + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "DivalentLonePair", + "name": "EP-L", + "match": "all_permutations", + "distance": l_distance, + "outOfPlaneAngle": out_of_plane_angle * unit.degree, + "charge_increment": [ + 0 * unit.elementary_charge, + -m_parameter.charge_increment2 / 2, + -m_parameter.charge_increment3 / 2, + ], + "epsilon": 0.0 * unit.kilocalorie_per_mole, + "sigma": 1.0 * unit.angstrom, + } + vsites.add_parameter(l_parameter_kwargs) + # don't optimize this just to avoid confusing optimizer too much + # l_parameter.add_cosmetic_attribute( + # "parameterize", "distance,outOfPlaneAngle" + # ) + # l_parameter.add_cosmetic_attribute( + # "parameter_eval", + # f"charge_increment2=-0.5*{m_chargeincrement2}, charge_increment3=-0.5*{m_chargeincrement2}" + # ) + + molecule = Molecule.from_mapped_smiles("[H:2][O:1][H:3]") + molecule.generate_conformers(n_conformers=1) + interchange = water_ff.create_interchange(molecule.to_topology()) + + # Validate geometry of virtual sites; should match angle geometry. + interchange.minimize() + positions = interchange.get_positions(include_virtual_sites=True) + + rdmol = Chem.RWMol(molecule.to_rdkit()) + n_vsites = len(positions) - molecule.n_atoms + for _ in range(n_vsites): + rdmol.AddAtom(Chem.Atom(0)) # virtual site as dummy atom + + # Create conformer with virtual site positions + conf = Chem.Conformer(len(positions)) + for i in range(len(positions)): + pos = positions[i].m_as("angstrom") + conf.SetAtomPosition(i, pos) + rdmol.AddConformer(conf) + + lol_angle = rdMolTransforms.GetAngleDeg(conf, 4, 0, 5) + ol_distance = rdMolTransforms.GetBondLength(conf, 0, 4) + + l_distance_ = l_distance.m_as("angstrom") + + assert np.isclose(ol_distance, l_distance_), f"Expected OL distance {l_distance_}, got {ol_distance}" + assert np.isclose(lol_angle, 2 * out_of_plane_angle), f"Expected LOL angle {2 * out_of_plane_angle}, got {lol_angle}" + + # now combine force fields + # Merge updated water handlers into the base small-molecule force field. + for handler_name in ["Constraints", "vdW", "LibraryCharges", "VirtualSites"]: + small_molecule_handler = small_molecule_ff.get_parameter_handler(handler_name) + water_handler = water_ff.get_parameter_handler(handler_name) + + for parameter in water_handler.parameters: + if parameter.id and "ion" in parameter.id: + continue + kwargs = parameter.to_dict() + if "parameterize" in parameter._cosmetic_attribs or "parameter_eval" in parameter._cosmetic_attribs: + kwargs["allow_cosmetic_attributes"] = True + small_molecule_handler.add_parameter(kwargs) + + output_force_field.parent.mkdir(parents=True, exist_ok=True) + + small_molecule_ff.to_file(output_force_field) + logger.info(f"Wrote modified force field to {output_force_field}") + + +if __name__ == "__main__": + main() diff --git a/studies/028_smirnoff_tip4p_geometry_fit/optimize.in b/studies/028_smirnoff_tip4p_geometry_fit/optimize.in new file mode 100644 index 000000000..bd26f600e --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/optimize.in @@ -0,0 +1,98 @@ +# ForceBalance input file generated by MakeInputFile.py +# The octothorpe '#' is a comment symbol +# There are two sections, the main options ($options) and the target options ($target) +# A ForceBalance calculation will have one $options section and as one $target section per optimization target +# The most important options are listed at the top; options are also roughly grouped by their application +# Note: If the specified value is 'None' then the option will truly be set to None - not the string 'None' +# Note: Section option types are more complicated and may require you to read the documentation +# Note: Boolean option types require no value, the key being present implies 'True' +# Note: List option types are specified using spaces as the delimiter - i.e. forcefield ff1.itp ff2.itp ; delete empty brackets before use [] + +$options +ffdir forcefield +penalty_type L2 +jobtype optimize +forcefield force-field.offxml +maxstep 1 +convergence_step 0.1 +convergence_objective 0.1 +convergence_gradient 0.1 +eig_lowerbound 0.01 +finite_difference_h 0.001 +penalty_additive 1.0 +trust0 0.1 +mintrust 0.05 +error_tolerance 1.0 +adaptive_factor 0.2 +adaptive_damping 1.0 +normalize_weights no +print_hessian +constrain_charge false + +priors + vdW/Atom/epsilon : 0.1 + vdW/Atom/sigma : 0.1 + VirtualSites/VirtualSite/charge_increment2 : 0.1 + VirtualSites/VirtualSite/distance : 0.1 + VirtualSites/VirtualSite/outOfPlaneAngle : 0.1 +/priors + +$end + +$target +# (string) The name of the target, which corresponds to the directory targets/dir_name +name Liquid-H25-RTA + +# (allcap) The type of target, for instance AbInitio_GMXX2 +type Liquid_SMIRNOFF + +# (float) Weight of the current target (with respect to other targets) +weight 6.0 + +w_rho 1.0 +w_hvap 1.0 +w_alpha 1.0 +w_kappa 1.0 +w_cp 1.0 +w_eps0 1.0 + +liquid_eq_steps 50 +liquid_md_steps 100 +liquid_timestep 2.0 +liquid_interval 0.1 + +gas_eq_steps 50 +gas_prod_steps 100 +gas_timestep 2.0 +gas_interval 0.1 + +# Parameters for self-polarization correction of nonpolarizable water +self_pol_mu0 1.855 +self_pol_alpha 1.470 + +adapt_errors + +pure_num_grad False +mol2 water.sdf + +$end + +$target +name water +type AbInitio_SMIRNOFF +mol2 input.sdf +pdb water.pdb +weight 1000.0 +w_force 0.1 # Reflects intrinsically small QM force +attenuate +energy_denom 5.0 +energy_upper 20.0 +writelevel 2 +$end + +$target +name phys-prop +type Evaluator_SMIRNOFF +weight 1.0 +evaluator_input options.json +$end diff --git a/studies/028_smirnoff_tip4p_geometry_fit/optimize.sav b/studies/028_smirnoff_tip4p_geometry_fit/optimize.sav new file mode 100644 index 000000000..d30e57019 --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/optimize.sav @@ -0,0 +1,68 @@ +# ForceBalance input file generated by MakeInputFile.py +# The octothorpe '#' is a comment symbol +# There are two sections, the main options ($options) and the target options ($target) +# A ForceBalance calculation will have one $options section and as one $target section per optimization target +# The most important options are listed at the top; options are also roughly grouped by their application +# Note: If the specified value is 'None' then the option will truly be set to None - not the string 'None' +# Note: Section option types are more complicated and may require you to read the documentation +# Note: Boolean option types require no value, the key being present implies 'True' +# Note: List option types are specified using spaces as the delimiter - i.e. forcefield ff1.itp ff2.itp ; delete empty brackets before use [] + +$options +ffdir forcefield +penalty_type L2 +jobtype optimize +forcefield force-field.offxml +maxstep 1 +convergence_step 0.1 +convergence_objective 0.1 +convergence_gradient 0.1 +eig_lowerbound 0.01 +finite_difference_h 0.001 +penalty_additive 1.0 +trust0 0.1 +mintrust 0.05 +error_tolerance 1.0 +adaptive_factor 0.2 +adaptive_damping 1.0 +normalize_weights no +print_hessian +constrain_charge false + +priors + vdW/Atom/epsilon : 0.1 + vdW/Atom/sigma : 0.1 + VirtualSites/VirtualSite/charge_increment2 : 0.1 + VirtualSites/VirtualSite/distance : 0.1 + VirtualSites/VirtualSite/outOfPlaneAngle : 0.1 +/priors + +read_mvals + 0 [ 1.70508046e-05 ] : vdW/Atom/epsilon/[#6X4:1] + 1 [ 2.91614378e-06 ] : vdW/Atom/sigma/[#6X4:1] + 2 [ 1.02582056e-03 ] : vdW/Atom/epsilon/[#1]-[#8X2H2+0:1]-[#1] + 3 [ 3.73145020e-02 ] : vdW/Atom/sigma/[#1]-[#8X2H2+0:1]-[#1] + 4 [ 2.14322272e-01 ] : VirtualSites/VirtualSite/distance/[#1:2]-[#8X2H2+0:1]-[#1:3]/DivalentLonePair/EP/once + 5 [ 3.12299464e-02 ] : VirtualSites/VirtualSite/charge_increment2/[#1:2]-[#8X2H2+0:1]-[#1:3]/DivalentLonePair/EP/once +/read_mvals +$end + +$target +name water +type AbInitio_SMIRNOFF +mol2 input.sdf +pdb water.pdb +weight 1000.0 +w_force 0.1 # Reflects intrinsically small QM force +attenuate +energy_denom 5.0 +energy_upper 20.0 +writelevel 2 +$end + +$target +name phys-prop +type Evaluator_SMIRNOFF +weight 1.0 +evaluator_input options.json +$end diff --git a/studies/028_smirnoff_tip4p_geometry_fit/run_server.py b/studies/028_smirnoff_tip4p_geometry_fit/run_server.py new file mode 100644 index 000000000..2feb7a821 --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/run_server.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python3 +import shutil +from os import path + +from openff.evaluator.backends import ComputeResources +from openff.evaluator.backends.dask import DaskLocalCluster +from openff.evaluator.server import EvaluatorServer +from openff.evaluator.utils import setup_timestamp_logging + + +def main(): + + # Set up logging for the evaluator. + setup_timestamp_logging() + + # Set up the directory structure. + working_directory = "working_directory" + + # Remove any existing data. + if path.isdir(working_directory): + shutil.rmtree(working_directory) + + # Set up a backend to run the calculations on with the requested resources. + worker_resources = ComputeResources(number_of_threads=1) + + calculation_backend = DaskLocalCluster( + number_of_workers=1, resources_per_worker=worker_resources + ) + + with calculation_backend: + + server = EvaluatorServer( + calculation_backend=calculation_backend, + working_directory=working_directory, + port=8000, + ) + + # Tell the server to start listening for estimation requests. + server.start() + + +if __name__ == "__main__": + main() diff --git a/studies/028_smirnoff_tip4p_geometry_fit/smirnoff_parameter_assignments.json b/studies/028_smirnoff_tip4p_geometry_fit/smirnoff_parameter_assignments.json new file mode 100644 index 000000000..51a21d338 --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/smirnoff_parameter_assignments.json @@ -0,0 +1,950 @@ +{ + "targets/water/input.sdf": [ + { + "id": "c-tip4p-fb-H-O", + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", + "type": "Constraints", + "atoms": [ + 0, + 1 + ] + }, + { + "id": "c-tip4p-fb-H-O", + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", + "type": "Constraints", + "atoms": [ + 0, + 2 + ] + }, + { + "id": "c-tip4p-fb-H-O-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1:2]", + "type": "Constraints", + "atoms": [ + 1, + 2 + ] + }, + { + "id": "b88", + "smirks": "[#8:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 1 + ] + }, + { + "id": "b88", + "smirks": "[#8:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 2 + ] + }, + { + "id": "a28", + "smirks": "[*:1]-[#8:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 0, + 2 + ] + }, + { + "id": "n-tip4p-fb-O", + "smirks": "[#1]-[#8X2H2+0:1]-[#1]", + "type": "vdW", + "atoms": [ + 0 + ] + }, + { + "id": "n-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "vdW", + "atoms": [ + 1 + ] + }, + { + "id": "n-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "vdW", + "atoms": [ + 2 + ] + }, + { + "id": "q-tip4p-fb-O", + "smirks": "[#1]-[#8X2H2+0:1]-[#1]", + "type": "LibraryCharges", + "atoms": [ + 0 + ] + }, + { + "id": "q-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "LibraryCharges", + "atoms": [ + 1 + ] + }, + { + "id": "q-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "LibraryCharges", + "atoms": [ + 2 + ] + }, + { + "id": null, + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "VirtualSites", + "atoms": [ + 0 + ] + }, + { + "id": null, + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "VirtualSites", + "atoms": [ + 0 + ] + }, + { + "id": null, + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "VirtualSites", + "atoms": [ + 0 + ] + } + ], + "targets/phys-prop/CC(C)O": [ + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 0, + 4 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 0, + 5 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 0, + 6 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 1, + 7 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 2, + 8 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 2, + 9 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 2, + 10 + ] + }, + { + "id": "c1", + "smirks": "[#1:1]-[*:2]", + "type": "Constraints", + "atoms": [ + 3, + 11 + ] + }, + { + "id": "b1", + "smirks": "[#6X4:1]-[#6X4:2]", + "type": "Bonds", + "atoms": [ + 0, + 1 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 4 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 5 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 6 + ] + }, + { + "id": "b1", + "smirks": "[#6X4:1]-[#6X4:2]", + "type": "Bonds", + "atoms": [ + 1, + 2 + ] + }, + { + "id": "b14", + "smirks": "[#6:1]-[#8:2]", + "type": "Bonds", + "atoms": [ + 1, + 3 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 1, + 7 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 2, + 8 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 2, + 9 + ] + }, + { + "id": "b84", + "smirks": "[#6X4:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 2, + 10 + ] + }, + { + "id": "b88", + "smirks": "[#8:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 3, + 11 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 0, + 1, + 2 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 0, + 1, + 3 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 0, + 1, + 7 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 0, + 4 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 0, + 5 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 0, + 6 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 2, + 8 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 2, + 9 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 2, + 10 + ] + }, + { + "id": "a28", + "smirks": "[*:1]-[#8:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 3, + 11 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 2, + 1, + 3 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 2, + 1, + 7 + ] + }, + { + "id": "a1", + "smirks": "[*:1]~[#6X4:2]-[*:3]", + "type": "Angles", + "atoms": [ + 3, + 1, + 7 + ] + }, + { + "id": "a2", + "smirks": "[#1:1]-[#6X4:2]-[#1:3]", + "type": "Angles", + "atoms": [ + 4, + 0, + 5 + ] + }, + { + "id": "a2", + "smirks": "[#1:1]-[#6X4:2]-[#1:3]", + "type": "Angles", + "atoms": [ + 4, + 0, + 6 + ] + }, + { + "id": "a2", + "smirks": "[#1:1]-[#6X4:2]-[#1:3]", + "type": "Angles", + "atoms": [ + 5, + 0, + 6 + ] + }, + { + "id": "a2", + "smirks": "[#1:1]-[#6X4:2]-[#1:3]", + "type": "Angles", + "atoms": [ + 8, + 2, + 9 + ] + }, + { + "id": "a2", + "smirks": "[#1:1]-[#6X4:2]-[#1:3]", + "type": "Angles", + "atoms": [ + 8, + 2, + 10 + ] + }, + { + "id": "a2", + "smirks": "[#1:1]-[#6X4:2]-[#1:3]", + "type": "Angles", + "atoms": [ + 9, + 2, + 10 + ] + }, + { + "id": "t4", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]", + "type": "ProperTorsions", + "atoms": [ + 0, + 1, + 2, + 8 + ] + }, + { + "id": "t4", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]", + "type": "ProperTorsions", + "atoms": [ + 0, + 1, + 2, + 9 + ] + }, + { + "id": "t4", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]", + "type": "ProperTorsions", + "atoms": [ + 0, + 1, + 2, + 10 + ] + }, + { + "id": "t94", + "smirks": "[#6X4:1]-[#6X4:2]-[#8X2H1:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 0, + 1, + 3, + 11 + ] + }, + { + "id": "t4", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]", + "type": "ProperTorsions", + "atoms": [ + 2, + 1, + 0, + 4 + ] + }, + { + "id": "t4", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]", + "type": "ProperTorsions", + "atoms": [ + 2, + 1, + 0, + 5 + ] + }, + { + "id": "t4", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#6X4:4]", + "type": "ProperTorsions", + "atoms": [ + 2, + 1, + 0, + 6 + ] + }, + { + "id": "t94", + "smirks": "[#6X4:1]-[#6X4:2]-[#8X2H1:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 2, + 1, + 3, + 11 + ] + }, + { + "id": "t9", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]", + "type": "ProperTorsions", + "atoms": [ + 3, + 1, + 0, + 4 + ] + }, + { + "id": "t9", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]", + "type": "ProperTorsions", + "atoms": [ + 3, + 1, + 0, + 5 + ] + }, + { + "id": "t9", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]", + "type": "ProperTorsions", + "atoms": [ + 3, + 1, + 0, + 6 + ] + }, + { + "id": "t9", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]", + "type": "ProperTorsions", + "atoms": [ + 3, + 1, + 2, + 8 + ] + }, + { + "id": "t9", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]", + "type": "ProperTorsions", + "atoms": [ + 3, + 1, + 2, + 9 + ] + }, + { + "id": "t9", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]", + "type": "ProperTorsions", + "atoms": [ + 3, + 1, + 2, + 10 + ] + }, + { + "id": "t3", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 4, + 0, + 1, + 7 + ] + }, + { + "id": "t3", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 5, + 0, + 1, + 7 + ] + }, + { + "id": "t3", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 6, + 0, + 1, + 7 + ] + }, + { + "id": "t3", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 7, + 1, + 2, + 8 + ] + }, + { + "id": "t3", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 7, + 1, + 2, + 9 + ] + }, + { + "id": "t3", + "smirks": "[#1:1]-[#6X4:2]-[#6X4:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 7, + 1, + 2, + 10 + ] + }, + { + "id": "t93", + "smirks": "[*:1]-[#6X4:2]-[#8X2:3]-[#1:4]", + "type": "ProperTorsions", + "atoms": [ + 7, + 1, + 3, + 11 + ] + }, + { + "id": "n16", + "smirks": "[#6X4:1]", + "type": "vdW", + "atoms": [ + 0 + ] + }, + { + "id": "n16", + "smirks": "[#6X4:1]", + "type": "vdW", + "atoms": [ + 1 + ] + }, + { + "id": "n16", + "smirks": "[#6X4:1]", + "type": "vdW", + "atoms": [ + 2 + ] + }, + { + "id": "n19", + "smirks": "[#8X2H1+0:1]", + "type": "vdW", + "atoms": [ + 3 + ] + }, + { + "id": "n2", + "smirks": "[#1:1]-[#6X4]", + "type": "vdW", + "atoms": [ + 4 + ] + }, + { + "id": "n2", + "smirks": "[#1:1]-[#6X4]", + "type": "vdW", + "atoms": [ + 5 + ] + }, + { + "id": "n2", + "smirks": "[#1:1]-[#6X4]", + "type": "vdW", + "atoms": [ + 6 + ] + }, + { + "id": "n3", + "smirks": "[#1:1]-[#6X4]-[#7,#8,#9,#16,#17,#35]", + "type": "vdW", + "atoms": [ + 7 + ] + }, + { + "id": "n2", + "smirks": "[#1:1]-[#6X4]", + "type": "vdW", + "atoms": [ + 8 + ] + }, + { + "id": "n2", + "smirks": "[#1:1]-[#6X4]", + "type": "vdW", + "atoms": [ + 9 + ] + }, + { + "id": "n2", + "smirks": "[#1:1]-[#6X4]", + "type": "vdW", + "atoms": [ + 10 + ] + }, + { + "id": "n12", + "smirks": "[#1:1]-[#8]", + "type": "vdW", + "atoms": [ + 11 + ] + } + ], + "targets/phys-prop/O": [ + { + "id": "c-tip4p-fb-H-O", + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", + "type": "Constraints", + "atoms": [ + 0, + 1 + ] + }, + { + "id": "c-tip4p-fb-H-O", + "smirks": "[#1:1]-[#8X2H2+0:2]-[#1]", + "type": "Constraints", + "atoms": [ + 0, + 2 + ] + }, + { + "id": "c-tip4p-fb-H-O-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1:2]", + "type": "Constraints", + "atoms": [ + 1, + 2 + ] + }, + { + "id": "b88", + "smirks": "[#8:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 1 + ] + }, + { + "id": "b88", + "smirks": "[#8:1]-[#1:2]", + "type": "Bonds", + "atoms": [ + 0, + 2 + ] + }, + { + "id": "a28", + "smirks": "[*:1]-[#8:2]-[*:3]", + "type": "Angles", + "atoms": [ + 1, + 0, + 2 + ] + }, + { + "id": "n-tip4p-fb-O", + "smirks": "[#1]-[#8X2H2+0:1]-[#1]", + "type": "vdW", + "atoms": [ + 0 + ] + }, + { + "id": "n-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "vdW", + "atoms": [ + 1 + ] + }, + { + "id": "n-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "vdW", + "atoms": [ + 2 + ] + }, + { + "id": "q-tip4p-fb-O", + "smirks": "[#1]-[#8X2H2+0:1]-[#1]", + "type": "LibraryCharges", + "atoms": [ + 0 + ] + }, + { + "id": "q-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "LibraryCharges", + "atoms": [ + 1 + ] + }, + { + "id": "q-tip4p-fb-H", + "smirks": "[#1:1]-[#8X2H2+0]-[#1]", + "type": "LibraryCharges", + "atoms": [ + 2 + ] + }, + { + "id": null, + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "VirtualSites", + "atoms": [ + 0 + ] + }, + { + "id": null, + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "VirtualSites", + "atoms": [ + 0 + ] + }, + { + "id": null, + "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]", + "type": "VirtualSites", + "atoms": [ + 0 + ] + } + ] +} \ No newline at end of file diff --git a/studies/028_smirnoff_tip4p_geometry_fit/stored_data/object_keys.json b/studies/028_smirnoff_tip4p_geometry_fit/stored_data/object_keys.json new file mode 100644 index 000000000..886bdd390 --- /dev/null +++ b/studies/028_smirnoff_tip4p_geometry_fit/stored_data/object_keys.json @@ -0,0 +1 @@ +{"object_keys": {"ForceFieldData": ["450a6073cb7a46469d26b0a942a96232", "0fd0265dee4246489d36f368c3d65571", "9f7d2eb0e0a342a9974efa3d0d2faa77", "b0f354e28735430db00f369780c2c90c"], "StoredSimulationData": ["04447d34033d484abd009cf63d59b5e8", "cfda13ee705e46bf91d1d801c5ce95c3", "79419fbb3d8247e09e2143a809056e5f", "bc04e521b2b54a4e852ea4a7b496c7ab", "834839f90bca419a81b8361915ab7603", "9d55a08e82074f7fac9e6a07b9075cf4"]}, "@type": "openff.evaluator.storage.storage.StorageBackend._ObjectKeyData"} \ No newline at end of file diff --git a/studies/028_smirnoff_tip4p_geometry_fit/targets.tar.gz b/studies/028_smirnoff_tip4p_geometry_fit/targets.tar.gz new file mode 100644 index 000000000..5665c41a5 Binary files /dev/null and b/studies/028_smirnoff_tip4p_geometry_fit/targets.tar.gz differ diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/.gitignore b/studies/029_smirnoff_vs_openmm_water_vsite/.gitignore new file mode 100644 index 000000000..76a509e0b --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/.gitignore @@ -0,0 +1,3 @@ +*.tmp +*.json + diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip4p/tip4p_fb.offxml b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip4p/tip4p_fb.offxml new file mode 100644 index 000000000..1fb1ed9a7 --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip4p/tip4p_fb.offxml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip4p/tip4pfb.xml b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip4p/tip4pfb.xml new file mode 100644 index 000000000..92f86f0b0 --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip4p/tip4pfb.xml @@ -0,0 +1,33 @@ + + + 2014-05-28 + Lee-Ping Wang, Todd J. Martinez and Vijay S. Pande. Building force fields - an automatic, systematic and reproducible approach. Journal of Physical Chemistry Letters, 2014, 5, pp 1885-1891. DOI:10.1021/jz500737m + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip5p/tip5p.offxml b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip5p/tip5p.offxml new file mode 100644 index 000000000..cba8c5897 --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip5p/tip5p.offxml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip5p/tip5p.xml b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip5p/tip5p.xml new file mode 100644 index 000000000..f7e2ad74e --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/forcefield_tip5p/tip5p.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip4p_openmm.in b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip4p_openmm.in new file mode 100644 index 000000000..6b05451da --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip4p_openmm.in @@ -0,0 +1,28 @@ +$options +ffdir forcefield_tip4p +jobtype gradient +forcefield tip4pfb.xml +backup false +finite_difference_h 0.001 +penalty_additive 1.0 +normalize_weights no +constrain_charge false + +priors + tip4pfb/NonbondedForce/Atom/charge : 0.1 +/priors + +$end + +$target +name water +type AbInitio_OpenMM +pdb water.pdb +coords all.gro +weight 1.0 +w_force 0.1 +attenuate +energy_denom 5.0 +energy_upper 20.0 +openmm_platform Reference +$end diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip4p_smirnoff.in b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip4p_smirnoff.in new file mode 100644 index 000000000..855fc7354 --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip4p_smirnoff.in @@ -0,0 +1,27 @@ +$options +ffdir forcefield_tip4p +jobtype gradient +forcefield tip4p_fb.offxml +backup false +finite_difference_h 0.001 +penalty_additive 1.0 +normalize_weights no +constrain_charge false + +priors + /VirtualSites/VirtualSite/charge_increment2 : 0.1 +/priors + +$end + +$target +name water +type AbInitio_SMIRNOFF +mol2 input.sdf +pdb water.pdb +weight 1.0 +w_force 0.1 +attenuate +energy_denom 5.0 +energy_upper 20.0 +$end diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip5p_openmm.in b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip5p_openmm.in new file mode 100644 index 000000000..72fb0215d --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip5p_openmm.in @@ -0,0 +1,28 @@ +$options +ffdir forcefield_tip5p +jobtype gradient +forcefield tip5p.xml +backup false +finite_difference_h 0.001 +penalty_additive 1.0 +normalize_weights no +constrain_charge false + +priors + tip5p/NonbondedForce/Atom/charge : 0.1 +/priors + +$end + +$target +name water +type AbInitio_OpenMM +pdb water.pdb +coords all.gro +weight 1.0 +w_force 0.1 +attenuate +energy_denom 5.0 +energy_upper 20.0 +openmm_platform Reference +$end diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip5p_smirnoff.in b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip5p_smirnoff.in new file mode 100644 index 000000000..ce8ed081b --- /dev/null +++ b/studies/029_smirnoff_vs_openmm_water_vsite/gradient_tip5p_smirnoff.in @@ -0,0 +1,27 @@ +$options +ffdir forcefield_tip5p +jobtype gradient +forcefield tip5p.offxml +backup false +finite_difference_h 0.001 +penalty_additive 1.0 +normalize_weights no +constrain_charge false + +priors + /VirtualSites/VirtualSite/charge_increment2 : 0.1 +/priors + +$end + +$target +name water +type AbInitio_SMIRNOFF +mol2 input.sdf +pdb water.pdb +weight 1.0 +w_force 0.1 +attenuate +energy_denom 5.0 +energy_upper 20.0 +$end diff --git a/studies/029_smirnoff_vs_openmm_water_vsite/targets.tar.gz b/studies/029_smirnoff_vs_openmm_water_vsite/targets.tar.gz new file mode 100644 index 000000000..a6c6c0df2 Binary files /dev/null and b/studies/029_smirnoff_vs_openmm_water_vsite/targets.tar.gz differ