Add ChargeIncrementModelHandler

The-SMIRNOFF-force-field-format.md

@@ -231,21 +231,20 @@ In keeping with the AMBER force field philosophy, especially as implemented in s
 
 Here is an example:
 ```XML
-<ChargeIncrementModel version="0.3" number_of_conformers="10" quantum_chemical_method="AM1" partial_charge_method="CM2">
+<ChargeIncrementModel version="0.3" number_of_conformers="1" partial_charge_method="AM1-Mulliken">
   <!-- A fractional charge can be moved along a single bond -->
   <ChargeIncrement smirks="[#6X4:1]-[#6X3a:2]" charge_increment1="-0.0073*elementary_charge" charge_increment2="0.0073*elementary_charge"/>
-  <ChargeIncrement smirks="[#6X4:1]-[#6X3a:2]-[#7]" charge_increment1=+0.0943*elementary_charge" charge_increment2="-0.0943*elementary_charge"/>
+  <ChargeIncrement smirks="[#6X4:1]-[#6X3a:2]-[#7]" charge_increment1="0.0943*elementary_charge" charge_increment2="-0.0943*elementary_charge"/>
   <ChargeIncrement smirks="[#6X4:1]-[#8:2]" charge_increment1="-0.0718*elementary_charge" charge_increment2="0.0718*elementary_charge"/>
   <!--- Alternatively, fractional charges can be redistributed among any number of bonded atoms -->
-  <ChargeIncrement smirks="[N:1](H:2)(H:3)" charge_increment1="0.02*elementary_charge" charge_increment2="-0.01*elementary_charge" charge_increment3="-0.01*elementary_charge"/>
+  <ChargeIncrement smirks="[N:1]([H:2])([H:3])" charge_increment1="0.02*elementary_charge" charge_increment2="-0.01*elementary_charge" charge_increment3="-0.01*elementary_charge"/>
 </ChargeIncrementModel>
 ```
 The sum of formal charges for the molecule or fragment will be used to determine the total charge the molecule or fragment will possess.
 
 `<ChargeIncrementModel>` provides several optional attributes to control its behavior:
-* The `number_of_conformers` attribute (default: `"10"`) is used to specify how many conformers will be generated for the molecule (or capped fragment) prior to charging.
-* The `quantum_chemical_method` attribute (default: `"AM1"`) is used to specify the quantum chemical method applied to the molecule or capped fragment.
-* The `partial_charge_method` attribute (default: `"Mulliken"`) is used to specify how uncorrected partial charges are to be generated from the quantum chemical wavefunction. Later additions will add restrained electrostatic potential fitting (RESP) capabilities.   
+* The `number_of_conformers` attribute (default: `"1"`) is used to specify how many conformers will be generated for the molecule (or capped fragment) prior to charging.
+* The `partial_charge_method` attribute (default: `"AM1-Mulliken"`) is used to specify how uncorrected partial charges are to be generated. Later additions will add restrained electrostatic potential fitting (RESP) capabilities.   
 
 The `<ChargeIncrement>` tags specify how the quantum chemical derived charges are to be corrected to produce the final charges.
 The `charge_increment#` attributes specify how much the charge on the associated tagged atom index (replacing `#`) should be modified.

openforcefield/tests/test_forcefield.py

@@ -227,7 +227,7 @@
 
 xml_spec_docs_charge_increment_model_xml = '''
 <SMIRNOFF version="0.3" aromaticity_model="OEAroModel_MDL">
-  <ChargeIncrementModel version="0.3" number_of_conformers="10" quantum_chemical_method="AM1" partial_charge_method="CM2">
+  <ChargeIncrementModel version="0.3" number_of_conformers="1" partial_charge_method="AM1-Mulliken">
     <!-- A fractional charge can be moved along a single bond -->
     <ChargeIncrement smirks="[#6X4:1]-[#6X3a:2]" charge_increment1="-0.0073*elementary_charge" charge_increment2="0.0073*elementary_charge"/>
     <ChargeIncrement smirks="[#6X4:1]-[#6X3a:2]-[#7]" charge_increment1="0.0943*elementary_charge" charge_increment2="-0.0943*elementary_charge"/>
@@ -284,6 +284,36 @@ def create_ethanol():
     ethanol.partial_charges = charges
     return ethanol
 
+
+def create_reversed_ethanol():
+    """
+    Creates an openforcefield.topology.Molecule representation of
+    ethanol without the use of a cheminformatics toolkit. This function
+    reverses the atom indexing of create_ethanol
+    """
+    # Create an ethanol molecule without using a toolkit
+    ethanol = Molecule()
+    ethanol.add_atom(1, 0, False)  # H0
+    ethanol.add_atom(1, 0, False)  # H1
+    ethanol.add_atom(1, 0, False)  # H2
+    ethanol.add_atom(1, 0, False)  # H3
+    ethanol.add_atom(1, 0, False)  # H4
+    ethanol.add_atom(1, 0, False)  # H5
+    ethanol.add_atom(8, 0, False)  # O6
+    ethanol.add_atom(6, 0, False)  # C7
+    ethanol.add_atom(6, 0, False)  # C8
+    ethanol.add_bond(8, 7, 1, False)  # C8 - C7
+    ethanol.add_bond(7, 6, 1, False)  # C7 - O6
+    ethanol.add_bond(8, 5, 1, False)  # C8 - H5
+    ethanol.add_bond(8, 4, 1, False)  # C8 - H4
+    ethanol.add_bond(8, 3, 1, False)  # C8 - H3
+    ethanol.add_bond(7, 2, 1, False)  # C7 - H2
+    ethanol.add_bond(7, 1, 1, False)  # C7 - H1
+    ethanol.add_bond(6, 0, 1, False)  # O6 - H0
+    charges = unit.Quantity(np.array([0.4, 0.3, 0.2, 0.1, 0.00001, -0.1, -0.2, -0.3, -0.4]), unit.elementary_charge)
+    ethanol.partial_charges = charges
+    return ethanol
+
 def create_cyclohexane():
     """
     Creates an openforcefield.topology.Molecule representation of
@@ -1155,6 +1185,30 @@ def test_parse_charge_increment_model_from_spec_docs(self):
         #       We should implement something like doctests for the XML snippets on the SMIRNOFF spec page.
         ff = ForceField(xml_spec_docs_charge_increment_model_xml)
 
+    def test_charge_increment_model_zero_net_charge(self):
+        """Test application of ChargeIncrements"""
+        test_charge_increment_model_ff = '''
+<SMIRNOFF version="0.3" aromaticity_model="OEAroModel_MDL">
+  <Electrostatics version="0.3" method="PME" scale12="0.0" scale13="0.0" scale14="0.833333" cutoff="9.0 * angstrom"/>
+  <ChargeIncrementModel version="0.3" number_of_conformers="1" partial_charge_method="formal_charge">
+    <ChargeIncrement smirks="[#6X4:1]-[#8:2]" charge_increment1="-0.05*elementary_charge" charge_increment2="0.05*elementary_charge"/>
+    <ChargeIncrement smirks="[C:1][C:2][O:3]" charge_increment1="0.2*elementary_charge" charge_increment2="-0.1*elementary_charge" charge_increment3="-0.1*elementary_charge"/>
+  </ChargeIncrementModel>
+</SMIRNOFF>
+        '''
+        file_path = get_data_file_path('test_forcefields/smirnoff99Frosst.offxml')
+
+        ff = ForceField(file_path, test_charge_increment_model_ff)
+        del ff._parameter_handlers['ToolkitAM1BCC']
+        top = Topology.from_molecules([create_ethanol(), create_reversed_ethanol()])
+        sys = ff.create_openmm_system(top)
+        nonbonded_force = [force for force in sys.getForces() if isinstance(force, openmm.NonbondedForce)][0]
+        expected_charges = [0.2, -0.15, -0.05, 0., 0., 0., 0., 0., 0.,
+                            0., 0., 0., 0., 0., 0., -0.05, -0.15, 0.2] * unit.elementary_charge
+        for idx, expected_charge in enumerate(expected_charges):
+            charge, _, _ = nonbonded_force.getParticleParameters(idx)
+            assert abs(charge - expected_charge) < 1.e-6 * unit.elementary_charge
+
     @pytest.mark.parametrize("inputs", partial_charge_method_resolution_matrix)
     def test_partial_charge_resolution(self, inputs):
         """Check that the proper partial charge methods are available, and that unavailable partial charge methods
@@ -1223,8 +1277,8 @@ def test_library_charges_to_two_waters(self):
             q, sigma, epsilon = nonbondedForce.getParticleParameters(particle_index)
             assert q == expected_charge
 
-    def test_library_charges_to_two_ethanols_different_atom_ordering(self):
-        """Test assigning charges to two ethanols with different atom orderings"""
+    def test_library_charges_to_three_ethanols_different_atom_ordering(self):
+        """Test assigning charges to three ethanols with different atom orderings"""
         from simtk.openmm import NonbondedForce
 
         # Define a library charge parameter for ethanol (C1-C2-O3) where C1 has charge -0.2, and its Hs have -0.02,
@@ -1245,18 +1299,29 @@ def test_library_charges_to_two_ethanols_different_atom_ordering(self):
         # H6 - C1 - C3 - O2 - H4
         #      |    |
         #      H7   H9
+        #
+        # create_reversed_ethanol()
+        #      H5   H2
+        #      |    |
+        # H4 - C8 - C7 - O6 - H0
+        #      |    |
+        #      H3   H1
+
 
         molecules = [Molecule.from_file(get_data_file_path('molecules/ethanol.sdf')),
-                     Molecule.from_file(get_data_file_path('molecules/ethanol_reordered.sdf'))]
+                     Molecule.from_file(get_data_file_path('molecules/ethanol_reordered.sdf')),
+                     create_reversed_ethanol()]
         top = Topology.from_molecules(molecules)
         omm_system = ff.create_openmm_system(top)
         nonbondedForce = [f for f in omm_system.getForces() if type(f) == NonbondedForce][0]
-        expected_charges = [-0.2, -0.1, 0.3, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01, -0.2,
-                            0.3, -0.1, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01] * unit.elementary_charge
+        expected_charges = [-0.2, -0.1, 0.3, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01,
+                            -0.2, 0.3, -0.1, 0.08, -0.02, -0.02, -0.02, -0.01, -0.01,
+                            0.08, -0.01, -0.01, -0.02, -0.02, -0.02, 0.3, -0.1, -0.2] * unit.elementary_charge
         for particle_index, expected_charge in enumerate(expected_charges):
             q, sigma, epsilon = nonbondedForce.getParticleParameters(particle_index)
             assert q == expected_charge
 
+
     def test_charge_method_hierarchy(self):
         """Ensure that molecules are parameterized by charge_from_molecules first, then library charges
         if not applicable, then AM1BCC otherwise"""

openforcefield/typing/engines/smirnoff/parameters.py

@@ -2823,8 +2823,6 @@ def find_matches(self, entity):
             matching the tuple of particle indices in ``entity``.
         """
 
-        # TODO: Right now, this method is only ever called with an entity that is a Topoogy.
-        #  Should we reduce its scope and have a check here to make sure entity is a Topology?
         return self._find_matches(entity, transformed_dict_cls=dict)
 
     def create_force(self, system, topology, **kwargs):
@@ -3021,12 +3019,12 @@ class ChargeIncrementType(ParameterType):
 
     _TAGNAME = 'ChargeIncrementModel'  # SMIRNOFF tag name to process
     _INFOTYPE = ChargeIncrementType  # info type to store
-    # TODO: The structure of this is still undecided
+    _DEPENDENCIES = [vdWHandler, ElectrostaticsHandler, LibraryChargeHandler]
 
     number_of_conformers = ParameterAttribute(default=1, converter=int)
 
     partial_charge_method = ParameterAttribute(
-        default='AM1-Mulliken'
+        default='AM1-Mulliken', converter=str
     )
 
     def check_handler_compatibility(self,
@@ -3047,11 +3045,32 @@ def check_handler_compatibility(self,
         """
 
         int_attrs_to_compare = ['number_of_conformers']
-        string_attrs_to_compare = ['quantum_chemical_method', 'partial_charge_method']
+        string_attrs_to_compare = ['partial_charge_method']
 
         self._check_attributes_are_equal(other_handler,
                                          identical_attrs=string_attrs_to_compare+int_attrs_to_compare)
 
+    def find_matches(self, entity):
+        """Find the elements of the topology/molecule matched by a parameter type.
+
+        Parameters
+        ----------
+        entity : openforcefield.topology.Topology
+            Topology to search.
+
+        Returns
+        ---------
+        matches : ValenceDict[Tuple[int], ParameterHandler._Match]
+            ``matches[particle_indices]`` is the ``ParameterType`` object
+            matching the tuple of particle indices in ``entity``.
+        """
+        # It's necessary to overwrite the default behavior of this method to prevent it
+        # from using transformed_dict_class=ValenceDict, which could change the order
+        # in which chargeincrements match topology particles
+
+        return self._find_matches(entity, transformed_dict_cls=dict)
+
+
 
     def create_force(self, system, topology, **kwargs):
 
@@ -3066,6 +3085,8 @@ def create_force(self, system, topology, **kwargs):
         else:
             force = existing[0]
 
+
+
         for ref_mol in topology.reference_molecules:
 
             # If charges were already assigned, skip this molecule
@@ -3079,14 +3100,7 @@ def create_force(self, system, topology, **kwargs):
             temp_mol.generate_conformers(n_conformers=self.number_of_conformers)
             temp_mol.compute_partial_charges(partial_charge_method=self.partial_charge_method)
 
-            bond_matches = self.find_matches(temp_mol)
-
-            for (atoms, charge_increment_match) in bond_matches.items():
-                charge_increment = charge_increment_match.parameter_type
-
-                for ref_mol_atom_idx, charge_increment in zip(atoms, charge_increment.charge_increment):
-                    temp_mol.partial_charges[ref_mol_atom_idx] += charge_increment
-
+            charges_to_assign = {}
 
             # Assign charges to relevant atoms
             for topology_molecule in topology._reference_molecule_to_topology_molecules[ref_mol]:
@@ -3098,14 +3112,30 @@ def create_force(self, system, topology, **kwargs):
                         ref_mol_particle_index = topology_particle.virtual_site.molecule_particle_index
                     particle_charge = temp_mol._partial_charges[ref_mol_particle_index]
 
-                    # Retrieve nonbonded parameters for reference atom (charge not set yet)
-                    _, sigma, epsilon = force.getParticleParameters(topology_particle_index)
-                    # Set the nonbonded force with the partial charge
-                    force.setParticleParameters(topology_particle_index,
-                                                particle_charge, sigma,
-                                                epsilon)
+                    charges_to_assign[topology_particle_index] = particle_charge
+                    # # Retrieve nonbonded parameters for reference atom (charge not set yet)
+                    # _, sigma, epsilon = force.getParticleParameters(topology_particle_index)
+                    # # Set the nonbonded force with the partial charge
+                    # force.setParticleParameters(topology_particle_index,
+                    #                     #                             particle_charge, sigma,
+                    #                     #                             epsilon)
 
 
+            charge_increment_matches = self.find_matches(topology)
+
+            for (atoms, charge_increment_match) in charge_increment_matches.items():
+                charge_increment = charge_increment_match.parameter_type
+
+                for top_particle_idx, charge_increment in zip(atoms, charge_increment.charge_increment):
+                    #print(top_particle_idx, charge_increment)
+                    if top_particle_idx in charges_to_assign:
+                        charges_to_assign[top_particle_idx] += charge_increment
+            #1/0
+            for topology_particle_index, charge_to_assign in charges_to_assign.items():
+                _, sigma, epsilon = force.getParticleParameters(topology_particle_index)
+                force.setParticleParameters(topology_particle_index,
+                                            charge_to_assign, sigma,
+                                            epsilon)
 
             # Finally, mark that charges were assigned for this reference molecule
             self.mark_charges_assigned(ref_mol, topology)