Merge branch 'master' into fix-order-smirnoff-cache

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 = []

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":

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

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

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")

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,15 +365,22 @@ 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'):
  self.offxml = [self.FF.offxml]
  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 = []
@@ -509,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)
@@ -524,7 +534,7 @@ def _update_positions(self, X1, disable_vsite):
  else:
  X1 = (X1 + [Vec3(0.0, 0.0, 0.0)] * n_v_sites) * angstrom
 
- self.simulation.context.setPositions(X1)
+ self.simulation.context.setPositions(X1 * angstrom)
  self.simulation.context.computeVirtualSites()
 
  def interaction_energy(self, fraga, fragb):