Merge pull request #226 from jthorton/polar_fix

Fix #225

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'):
@@ -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)
@@ -518,10 +528,11 @@ def _update_positions(self, X1, disable_vsite):
  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
+ # Add placeholder positions for any v-sites.
+ sites = np.zeros((n_v_sites, 3))
+ X1 = np.vstack((X1, sites))
 
- self.simulation.context.setPositions(X1)
+ self.simulation.context.setPositions(X1 * angstrom)
  self.simulation.context.computeVirtualSites()
 
  def interaction_energy(self, fraga, fragb):