docs: def createSystem(self, topology, nonbondedMethod=NoCutoff,
                       nonbondedCutoff=1.0 * u.nanometer, constraints=None,
                       rigidWater=True, removeCMMotion=True, hydrogenMass=None,
                       **args):
          """Construct an OpenMM System representing a Topology with this force field.
  
          Parameters
          ----------
          topology : Topology
              The Topology for which to create a System
          nonbondedMethod : object=NoCutoff
              The method to use for nonbonded interactions.  Allowed values are
              NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
          nonbondedCutoff : distance=1*nanometer
              The cutoff distance to use for nonbonded interactions
          constraints : object=None
              Specifies which bonds and angles should be implemented with constraints.
              Allowed values are None, HBonds, AllBonds, or HAngles.
          rigidWater : boolean=True
              If true, water molecules will be fully rigid regardless of the value
              passed for the constraints argument
          removeCMMotion : boolean=True
              If true, a CMMotionRemover will be added to the System
          hydrogenMass : mass=None
              The mass to use for hydrogen atoms bound to heavy atoms.  Any mass
              added to a hydrogen is subtracted from the heavy atom to keep
              their total mass the same.
          args
               Arbitrary additional keyword arguments may also be specified.
               This allows extra parameters to be specified that are specific to
               particular force fields.
  
          Returns
          -------
          system
              the newly created System
          """
  
          # Overwrite previous _SystemData object
          self._SystemData = app.ForceField._SystemData()
  
          data = self._SystemData
          data.atoms = list(topology.atoms())
          for atom in data.atoms:
              data.excludeAtomWith.append([])
  
          # Make a list of all bonds
          for bond in topology.bonds():
              data.bonds.append(app.ForceField._BondData(bond[0].index, bond[1].index))
  
          # Record which atoms are bonded to each other atom
          bonded_to_atom = []
          for i in range(len(data.atoms)):
              bonded_to_atom.append(set())
              data.atomBonds.append([])
          for i in range(len(data.bonds)):
              bond = data.bonds[i]
              bonded_to_atom[bond.atom1].add(bond.atom2)
              bonded_to_atom[bond.atom2].add(bond.atom1)
              data.atomBonds[bond.atom1].append(i)
              data.atomBonds[bond.atom2].append(i)
  
          # TODO: Better way to lookup nonbonded parameters...?
          nonbonded_params = None
          for generator in self.getGenerators():
              if isinstance(generator, NonbondedGenerator):
                  nonbonded_params = generator.params.paramsForType
                  break
  
          for chain in topology.chains():
              for res in chain.residues():
                  for atom in res.atoms():
                      data.atomType[atom] = atom.id
                      if nonbonded_params:
                          params = nonbonded_params[atom.id]
                          data.atomParameters[atom] = params
  
          # Create the System and add atoms
          sys = mm.System()
          for atom in topology.atoms():
              # Look up the atom type name, returning a helpful error message if it cannot be found.
              if atom not in data.atomType:
                  raise Exception("Could not identify atom type for atom '%s'." % str(atom))
              typename = data.atomType[atom]
  
              # Look up the type name in the list of registered atom types, returning a helpful error message if it cannot be found.
              if typename not in self._atomTypes:
                  msg  = "Could not find typename '%s' for atom '%s' in list of known atom types.\n" % (typename, str(atom))
                  msg += "Known atom types are: %s" % str(self._atomTypes.keys())
                  raise Exception(msg)
  
              # Add the particle to the OpenMM system.
              mass = self._atomTypes[typename].mass
              sys.addParticle(mass)
  
          # Adjust hydrogen masses if requested.
          if hydrogenMass is not None:
              if not u.is_quantity(hydrogenMass):
                  hydrogenMass *= u.dalton
              for atom1, atom2 in topology.bonds():
                  if atom1.element == elem.hydrogen:
                      (atom1, atom2) = (atom2, atom1)
                  if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
                      transfer_mass = hydrogenMass - sys.getParticleMass(atom2.index)
                      sys.setParticleMass(atom2.index, hydrogenMass)
                      mass = sys.getParticleMass(atom1.index) - transfer_mass
                      sys.setParticleMass(atom1.index, mass)
  
          # Set periodic boundary conditions.
          box_vectors = topology.getPeriodicBoxVectors()
          if box_vectors is not None:
              sys.setDefaultPeriodicBoxVectors(box_vectors[0],
                                               box_vectors[1],
                                               box_vectors[2])
          elif nonbondedMethod not in [NoCutoff, CutoffNonPeriodic]:
              raise ValueError('Requested periodic boundary conditions for a '
                               'Topology that does not specify periodic box '
                               'dimensions')
  
          # Make a list of all unique angles
          unique_angles = set()
          for bond in data.bonds:
              for atom in bonded_to_atom[bond.atom1]:
                  if atom != bond.atom2:
                      if atom < bond.atom2:
                          unique_angles.add((atom, bond.atom1, bond.atom2))
                      else:
                          unique_angles.add((bond.atom2, bond.atom1, atom))
              for atom in bonded_to_atom[bond.atom2]:
                  if atom != bond.atom1:
                      if atom > bond.atom1:
                          unique_angles.add((bond.atom1, bond.atom2, atom))
                      else:
                          unique_angles.add((atom, bond.atom2, bond.atom1))
          data.angles = sorted(list(unique_angles))
  
          # Make a list of all unique proper torsions
          unique_propers = set()
          for angle in data.angles:
              for atom in bonded_to_atom[angle[0]]:
                  if atom not in angle:
                      if atom < angle[2]:
                          unique_propers.add((atom, angle[0], angle[1], angle[2]))
                      else:
                          unique_propers.add((angle[2], angle[1], angle[0], atom))
              for atom in bonded_to_atom[angle[2]]:
                  if atom not in angle:
                      if atom > angle[0]:
                          unique_propers.add((angle[0], angle[1], angle[2], atom))
                      else:
                          unique_propers.add((atom, angle[2], angle[1], angle[0]))
          data.propers = sorted(list(unique_propers))
  
          # Make a list of all unique improper torsions
          for atom in range(len(bonded_to_atom)):
              bonded_to = bonded_to_atom[atom]
              if len(bonded_to) > 2:
                  for subset in itertools.combinations(bonded_to, 3):
                      data.impropers.append((atom, subset[0], subset[1], subset[2]))
  
          # Identify bonds that should be implemented with constraints
          if constraints == AllBonds or constraints == HAngles:
              for bond in data.bonds:
                  bond.isConstrained = True
          elif constraints == HBonds:
              for bond in data.bonds:
                  atom1 = data.atoms[bond.atom1]
                  atom2 = data.atoms[bond.atom2]
                  bond.isConstrained = atom1.name.startswith('H') or atom2.name.startswith('H')
          if rigidWater:
              for bond in data.bonds:
                  atom1 = data.atoms[bond.atom1]
                  atom2 = data.atoms[bond.atom2]
                  if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH':
                      bond.isConstrained = True
  
          # Identify angles that should be implemented with constraints
          if constraints == HAngles:
              for angle in data.angles:
                  atom1 = data.atoms[angle[0]]
                  atom2 = data.atoms[angle[1]]
                  atom3 = data.atoms[angle[2]]
                  numH = 0
                  if atom1.name.startswith('H'):
                      numH += 1
                  if atom3.name.startswith('H'):
                      numH += 1
                  data.isAngleConstrained.append(numH == 2 or (numH == 1 and atom2.name.startswith('O')))
          else:
              data.isAngleConstrained = len(data.angles)*[False]
          if rigidWater:
              for i in range(len(data.angles)):
                  angle = data.angles[i]
                  atom1 = data.atoms[angle[0]]
                  atom2 = data.atoms[angle[1]]
                  atom3 = data.atoms[angle[2]]
                  if atom1.residue.name == 'HOH' and atom2.residue.name == 'HOH' and atom3.residue.name == 'HOH':
                      data.isAngleConstrained[i] = True
  
          # Add virtual sites
          for atom in data.virtualSites:
              (site, atoms, excludeWith) = data.virtualSites[atom]
              index = atom.index
              data.excludeAtomWith[excludeWith].append(index)
              if site.type == 'average2':
                  sys.setVirtualSite(index, mm.TwoParticleAverageSite(
                      atoms[0], atoms[1], site.weights[0], site.weights[1]))
              elif site.type == 'average3':
                  sys.setVirtualSite(index, mm.ThreeParticleAverageSite(
                      atoms[0], atoms[1], atoms[2],
                      site.weights[0], site.weights[1], site.weights[2]))
              elif site.type == 'outOfPlane':
                  sys.setVirtualSite(index, mm.OutOfPlaneSite(
                      atoms[0], atoms[1], atoms[2],
                      site.weights[0], site.weights[1], site.weights[2]))
              elif site.type == 'localCoords':
                  local_coord_site = mm.LocalCoordinatesSite(
                      atoms[0], atoms[1], atoms[2],
                      mm.Vec3(site.originWeights[0], site.originWeights[1], site.originWeights[2]),
                      mm.Vec3(site.xWeights[0], site.xWeights[1], site.xWeights[2]),
                      mm.Vec3(site.yWeights[0], site.yWeights[1], site.yWeights[2]),
                      mm.Vec3(site.localPos[0], site.localPos[1], site.localPos[2]))
                  sys.setVirtualSite(index, local_coord_site)
  
          # Add forces to the System
          for force in self._forces:
              force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
          if removeCMMotion:
              sys.addForce(mm.CMMotionRemover())
  
          # Let force generators do postprocessing
          for force in self._forces:
              if 'postprocessSystem' in dir(force):
                  force.postprocessSystem(sys, data, args)
  
          # Execute scripts found in the XML files.
          for script in self._scripts:
              exec(script, locals())
  
          return sys