From 605884ec67169407dc484f7fdcc2990d0cc37548 Mon Sep 17 00:00:00 2001 From: SimonBoothroyd Date: Thu, 1 Apr 2021 11:53:05 +0100 Subject: [PATCH 1/7] Load SMIRNOFF plugins by default --- src/evaluator_io.py | 4 ++-- src/forcefield.py | 2 +- src/recharge_io.py | 4 +++- src/smirnoffio.py | 2 +- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/evaluator_io.py b/src/evaluator_io.py index 8d5ef27c1..d244a9a0f 100644 --- a/src/evaluator_io.py +++ b/src/evaluator_io.py @@ -438,13 +438,13 @@ def submit_jobs(self, mvals, AGrad=True, AHess=True): self.FF.make(mvals) force_field = smirnoff.ForceField( - self.FF.offxml, allow_cosmetic_attributes=True + self.FF.offxml, allow_cosmetic_attributes=True, load_plugins=True ) # strip out cosmetic attributes with tempfile.NamedTemporaryFile(mode="w", suffix=".offxml") as file: force_field.to_file(file.name, discard_cosmetic_attributes=True) - force_field = smirnoff.ForceField(file.name) + force_field = smirnoff.ForceField(file.name, load_plugins=True) # Determine which gradients (if any) we should be estimating. parameter_gradient_keys = [] diff --git a/src/forcefield.py b/src/forcefield.py index 0625a5ccc..47f9a6339 100644 --- a/src/forcefield.py +++ b/src/forcefield.py @@ -458,7 +458,7 @@ def addff(self,ffname,xmlScript=False): from openff.toolkit.typing.engines.smirnoff import ForceField as OpenFF_ForceField self.offxml = ffname self.openff_forcefield = OpenFF_ForceField(os.path.join(self.root, self.ffdir, self.offxml), - allow_cosmetic_attributes=True) + allow_cosmetic_attributes=True, load_plugins=True) self.amber_mol2 = [] if fftype == "mol2": diff --git a/src/recharge_io.py b/src/recharge_io.py index 65ec397c4..7efb8dabe 100644 --- a/src/recharge_io.py +++ b/src/recharge_io.py @@ -83,7 +83,9 @@ def _initialize(self): # Determine which BCC parameters are being optimized. force_field = smirnoff.ForceField( - os.path.join(self.FF.ffdir, self.FF.offxml), allow_cosmetic_attributes=True + os.path.join(self.FF.ffdir, self.FF.offxml), + allow_cosmetic_attributes=True, + load_plugins=True ) bcc_handler = force_field.get_parameter_handler("ChargeIncrementModel") diff --git a/src/smirnoffio.py b/src/smirnoffio.py index ef7ef4176..92bdc232a 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -324,7 +324,7 @@ def prepare(self, pbc=False, mmopts={}, **kwargs): self.forcefield = self.FF.openff_forcefield else: self.offxml = listfiles(kwargs.get('offxml'), 'offxml', err=True) - self.forcefield = OpenFF_ForceField(*self.offxml) + self.forcefield = OpenFF_ForceField(*self.offxml, load_plugins=True) ## Load mol2 files for smirnoff topology openff_mols = [] From 89e5a08b2631d7fc796d0b176c44f4bf7dd8967a Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Tue, 4 May 2021 16:06:17 +0100 Subject: [PATCH 2/7] add localcoordinate site support and a test --- src/openmmio.py | 28 +++ src/tests/files/forcefield/vs_mol.xml | 276 ++++++++++++++++++++++++++ src/tests/files/vs_mol.pdb | 39 ++++ src/tests/test_openmmio.py | 27 +++ 4 files changed, 370 insertions(+) create mode 100644 src/tests/files/forcefield/vs_mol.xml create mode 100644 src/tests/files/vs_mol.pdb diff --git a/src/openmmio.py b/src/openmmio.py index 63edd339f..c6235ac40 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -151,6 +151,34 @@ def vsfunc(pos, idx_, wt_): v2 = pos[idx_[2]] - pos[idx_[0]] cross = np.array([v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0]]) return pos[idx_[0]] + wt_[0]*v1 + wt_[1]*v2 + wt_[2]*cross + elif isinstance(vs, LocalCoordinatesSite): + vsidx = [_openmm.VirtualSite_getParticle(vs, i) for i in range(_openmm.VirtualSite_getNumParticles(vs))] + vswt = [np.array(_openmm.LocalCoordinatesSite_getOriginWeights(vs)), np.array(_openmm.LocalCoordinatesSite_getXWeights(vs)), np.array(_openmm.LocalCoordinatesSite_getYWeights(vs)), np.array(_openmm.LocalCoordinatesSite_getLocalPosition(vs))] + def vsfunc(pos, idx_, wt_): + """Calculate the vsite position within a orthonormal coordinate system described here + http://docs.openmm.org/latest/api-c++/generated/OpenMM.LocalCoordinatesSite.html#localcoordinatessite""" + # origin weights + ows = wt_[0] + # xdir weights + xws = wt_[1] + # ydir weights + yws = wt_[2] + # vs position in local coordinates + vpos = wt_[3] + # dependent atom positions + dpos = np.array([pos[j] for j in idx_]) + origin = np.array(dpos * ows[:, None]).sum(axis=0) + xdir = np.array(dpos * xws[:, None]).sum(axis=0) + ydir = np.array(dpos * yws[:, None]).sum(axis=0) + zdir = np.cross(xdir, ydir) + ydir = np.cross(zdir, xdir) + xdir /= np.linalg.norm(xdir) + ydir /= np.linalg.norm(ydir) + zdir /= np.linalg.norm(zdir) + return origin + np.array(np.array([xdir, ydir, zdir]) * vpos[:, None]).sum(axis=0) + else: + raise NotImplementedError(f"The virtual site type {vs.__class__.__name__} is not currently supported.") + else: isvsites.append(0) vsfunc = None diff --git a/src/tests/files/forcefield/vs_mol.xml b/src/tests/files/forcefield/vs_mol.xml new file mode 100644 index 000000000..207cfb289 --- /dev/null +++ b/src/tests/files/forcefield/vs_mol.xml @@ -0,0 +1,276 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/tests/files/vs_mol.pdb b/src/tests/files/vs_mol.pdb new file mode 100644 index 000000000..feb67973a --- /dev/null +++ b/src/tests/files/vs_mol.pdb @@ -0,0 +1,39 @@ +COMPND mol21 +HETATM 1 C1 UNL 1 4.129 0.259 -0.327 1.00 0.00 C +HETATM 2 C2 UNL 1 2.835 1.066 -0.137 1.00 0.00 C +HETATM 3 C3 UNL 1 1.722 0.328 0.630 1.00 0.00 C +HETATM 4 C4 UNL 1 1.138 -0.887 -0.120 1.00 0.00 C +HETATM 5 N1 UNL 1 0.055 -1.622 0.550 1.00 0.00 N +HETATM 6 C5 UNL 1 -1.205 -0.894 0.771 1.00 0.00 C +HETATM 7 C6 UNL 1 -1.998 -0.691 -0.528 1.00 0.00 C +HETATM 8 C7 UNL 1 -3.422 -0.145 -0.317 1.00 0.00 C +HETATM 9 C8 UNL 1 -3.490 1.286 0.239 1.00 0.00 C +HETATM 10 H1 UNL 1 4.898 0.864 -0.820 1.00 0.00 H +HETATM 11 H2 UNL 1 3.973 -0.634 -0.941 1.00 0.00 H +HETATM 12 H3 UNL 1 4.534 -0.068 0.638 1.00 0.00 H +HETATM 13 H4 UNL 1 2.452 1.377 -1.119 1.00 0.00 H +HETATM 14 H5 UNL 1 3.074 1.993 0.400 1.00 0.00 H +HETATM 15 H6 UNL 1 2.112 -0.006 1.604 1.00 0.00 H +HETATM 16 H7 UNL 1 0.919 1.042 0.849 1.00 0.00 H +HETATM 17 H8 UNL 1 1.933 -1.614 -0.318 1.00 0.00 H +HETATM 18 H9 UNL 1 0.775 -0.565 -1.104 1.00 0.00 H +HETATM 19 H10 UNL 1 0.398 -1.962 1.447 1.00 0.00 H +HETATM 20 H11 UNL 1 -1.062 0.074 1.279 1.00 0.00 H +HETATM 21 H12 UNL 1 -1.809 -1.512 1.449 1.00 0.00 H +HETATM 22 H13 UNL 1 -1.455 -0.012 -1.199 1.00 0.00 H +HETATM 23 H14 UNL 1 -2.050 -1.659 -1.040 1.00 0.00 H +HETATM 24 H15 UNL 1 -3.974 -0.822 0.350 1.00 0.00 H +HETATM 25 H16 UNL 1 -3.948 -0.175 -1.279 1.00 0.00 H +HETATM 26 H17 UNL 1 -2.942 1.985 -0.405 1.00 0.00 H +HETATM 27 H18 UNL 1 -4.527 1.635 0.298 1.00 0.00 H +HETATM 28 H19 UNL 1 -3.062 1.357 1.244 1.00 0.00 H +CONECT 1 2 10 11 12 +CONECT 2 3 13 14 +CONECT 3 4 15 16 +CONECT 4 5 17 18 +CONECT 5 6 19 +CONECT 6 7 20 21 +CONECT 7 8 22 23 +CONECT 8 9 24 25 +CONECT 9 26 27 28 +END diff --git a/src/tests/test_openmmio.py b/src/tests/test_openmmio.py index 6f4392b53..97eccc100 100644 --- a/src/tests/test_openmmio.py +++ b/src/tests/test_openmmio.py @@ -1,5 +1,11 @@ from __future__ import absolute_import import forcebalance +from forcebalance.openmmio import PrepareVirtualSites, ResetVirtualSites_fast +import simtk.openmm as mm +from simtk.openmm import app +import numpy as np +from simtk import unit +import os import shutil from .test_target import TargetTests # general targets tests defined in test_target.py import pytest @@ -58,3 +64,24 @@ def teardown_method(self): shutil.rmtree('temp') super(TestInteraction_OpenMM, self).teardown_method() + +def test_local_coord_sites(): + """Make sure that the internal prep of vs positions matches that given by OpenMM.""" + # make a system + mol = app.PDBFile(os.path.join("files", "vs_mol.pdb")) + modeller = app.Modeller(topology=mol.topology, positions=mol.positions) + forcefield = app.ForceField(os.path.join("files", "forcefield", "vs_mol.xml")) + modeller.addExtraParticles(forcefield) + system = forcefield.createSystem(modeller.topology, constraints=None) + integrator = mm.VerletIntegrator(1.0 * unit.femtoseconds) + platform = mm.Platform.getPlatformByName("Reference") + simulation = app.Simulation(modeller.topology, system, integrator, platform) + simulation.context.setPositions(modeller.positions) + # update the site positions + simulation.context.computeVirtualSites() + # one vs and it should be last + vs_pos = simulation.context.getState(getPositions=True).getPositions(asNumpy=True)[-1] + # now compute and compare + vsinfo = PrepareVirtualSites(system=system) + new_pos = ResetVirtualSites_fast(positions=modeller.positions, vsinfo=vsinfo)[-1] + assert np.allclose(vs_pos._value, np.array([new_pos.x, new_pos.y, new_pos.z])) From 593791866e622ab4eae23ce29a0bed27499a118d Mon Sep 17 00:00:00 2001 From: Mateusz Bieniek Date: Fri, 28 May 2021 14:44:19 +0100 Subject: [PATCH 3/7] the files are compressed This was throwing the UnicodeDecodeError: 'utf-8' codec can't decode byte 0x80 in position 0: invalid start byte, when loading the file. --- src/optimizer.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/optimizer.py b/src/optimizer.py index 25176be8c..2582fe738 100644 --- a/src/optimizer.py +++ b/src/optimizer.py @@ -21,7 +21,7 @@ from copy import deepcopy import forcebalance from forcebalance.parser import parse_inputs -from forcebalance.nifty import col, flat, row, printcool, printcool_dictionary, pvec1d, pmat2d, warn_press_key, invert_svd, wopen, bak, est124 +from forcebalance.nifty import col, flat, row, printcool, printcool_dictionary, pvec1d, pmat2d, warn_press_key, invert_svd, wopen, bak, est124, lp_load from forcebalance.finite_difference import f1d7p, f1d5p, fdwrap from collections import OrderedDict import random @@ -1568,7 +1568,7 @@ def readchk(self): if self.rchk_fnm is not None: absfnm = os.path.join(self.root,self.rchk_fnm) if os.path.exists(absfnm): - self.chk = pickle.load(open(absfnm)) + self.chk = lp_load(absfnm) else: logger.info("\x1b[40m\x1b[1;92mWARNING:\x1b[0m read_chk is set to True, but checkpoint file not loaded (wrong filename or doesn't exist?)\n") return self.chk From 2949bc25827921ec4745f7a68425b9ddc974d4fd Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Fri, 4 Jun 2021 15:41:01 +0100 Subject: [PATCH 4/7] unify openmm and smirnoff prepare and readsrc functions --- src/smirnoffio.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/smirnoffio.py b/src/smirnoffio.py index 90b4c0295..410a0a67c 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -281,13 +281,14 @@ def readsrc(self, **kwargs): This could be the same as the pdb argument from above. """ - pdbfnm = kwargs.get('pdb') + pdbfnm = None # Determine the PDB file name. - if not pdbfnm: - raise RuntimeError('Name of PDB file not provided.') - elif not os.path.exists(pdbfnm): - logger.error("%s specified but doesn't exist\n" % pdbfnm) - raise RuntimeError + if 'pdb' in kwargs and os.path.exists(kwargs['pdb']): + # Case 1. The PDB file name is provided explicitly + pdbfnm = kwargs['pdb'] + if not os.path.exists(pdbfnm): + logger.error("%s specified but doesn't exist\n" % pdbfnm) + raise RuntimeError if 'mol' in kwargs: self.mol = kwargs['mol'] @@ -314,10 +315,11 @@ def readsrc(self, **kwargs): else: logger.error("Must provide a list of .mol2 files.\n") - self.abspdb = os.path.abspath(pdbfnm) - mpdb = Molecule(pdbfnm) - for i in ["chain", "atomname", "resid", "resname", "elem"]: - self.mol.Data[i] = mpdb.Data[i] + if pdbfnm is not None: + self.abspdb = os.path.abspath(pdbfnm) + mpdb = Molecule(pdbfnm) + for i in ["chain", "atomname", "resid", "resname", "elem"]: + self.mol.Data[i] = mpdb.Data[i] # Store a separate copy of the molecule for reference restraint positions. self.ref_mol = deepcopy(self.mol) @@ -363,7 +365,14 @@ def prepare(self, pbc=False, mmopts={}, **kwargs): but we are calling ForceField() from the OpenFF toolkit and ignoring AMOEBA stuff. """ - self.pdb = PDBFile(self.abspdb) + + if hasattr(self, 'abspdb'): + self.pdb = PDBFile(self.abspdb) + else: + pdb1 = "%s-1.pdb" % os.path.splitext(os.path.basename(self.mol.fnm))[0] + self.mol[0].write(pdb1) + self.pdb = PDBFile(pdb1) + os.unlink(pdb1) # Create the OpenFF ForceField object. if hasattr(self, 'FF'): @@ -371,7 +380,7 @@ def prepare(self, pbc=False, mmopts={}, **kwargs): self.forcefield = self.FF.openff_forcefield else: self.offxml = listfiles(kwargs.get('offxml'), 'offxml', err=True) - self.forcefield = OpenFF_ForceField(*self.offxml) + self.forcefield = OpenFF_ForceField(*self.offxml, load_plugins=True) ## Load mol2 files for smirnoff topology openff_mols = [] From d0e5298c31f1190d1c78c47435e202b751419e1a Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Fri, 4 Jun 2021 17:31:38 +0100 Subject: [PATCH 5/7] fix error when padding positions when there are no virtual sites. --- src/smirnoffio.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/smirnoffio.py b/src/smirnoffio.py index 410a0a67c..fbaa73889 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -526,9 +526,9 @@ def _update_positions(self, X1, disable_vsite): n_v_sites = ( self.mod.getTopology().getNumAtoms() - self.pdb.topology.getNumAtoms() ) - - # Add placeholder positions for an v-sites. - X1 = (X1 + [Vec3(0.0, 0.0, 0.0)] * n_v_sites) * angstrom + if n_v_sites > 0: + # Add placeholder positions for any v-sites. + X1 = (X1 + [Vec3(0.0, 0.0, 0.0)] * n_v_sites) * angstrom self.simulation.context.setPositions(X1) self.simulation.context.computeVirtualSites() From 2e902c7dcc67474626dfc5b202a2c0d94b9e08c7 Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 7 Jun 2021 10:49:15 +0100 Subject: [PATCH 6/7] make sure units are correct on updated positions --- src/smirnoffio.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/smirnoffio.py b/src/smirnoffio.py index fbaa73889..0240869f3 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -528,9 +528,9 @@ def _update_positions(self, X1, disable_vsite): ) if n_v_sites > 0: # Add placeholder positions for any v-sites. - X1 = (X1 + [Vec3(0.0, 0.0, 0.0)] * n_v_sites) * angstrom + X1 = (X1 + [Vec3(0.0, 0.0, 0.0)] * n_v_sites) - self.simulation.context.setPositions(X1) + self.simulation.context.setPositions(X1 * angstrom) self.simulation.context.computeVirtualSites() def interaction_energy(self, fraga, fragb): From 308b333c5ec6ea8556c689a2ddcfe7cdceefd80c Mon Sep 17 00:00:00 2001 From: Josh Horton Date: Mon, 7 Jun 2021 18:13:20 +0100 Subject: [PATCH 7/7] fix the virtual site padding to work with numpy ndarray --- src/smirnoffio.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/smirnoffio.py b/src/smirnoffio.py index 0240869f3..6c0da03b7 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -518,6 +518,7 @@ def update_simulation(self, **kwargs): self.create_simulation(**self.simkwargs) def _update_positions(self, X1, disable_vsite): + # X1 is a numpy ndarray not vec3 if disable_vsite: super(SMIRNOFF, self)._update_positions(X1, disable_vsite) @@ -526,9 +527,10 @@ def _update_positions(self, X1, disable_vsite): n_v_sites = ( self.mod.getTopology().getNumAtoms() - self.pdb.topology.getNumAtoms() ) - if n_v_sites > 0: - # Add placeholder positions for any v-sites. - X1 = (X1 + [Vec3(0.0, 0.0, 0.0)] * n_v_sites) + + # Add placeholder positions for any v-sites. + sites = np.zeros((n_v_sites, 3)) + X1 = np.vstack((X1, sites)) self.simulation.context.setPositions(X1 * angstrom) self.simulation.context.computeVirtualSites()