diff --git a/docs/releasehistory.md b/docs/releasehistory.md
index 3711106a0..3b89bcfef 100644
--- a/docs/releasehistory.md
+++ b/docs/releasehistory.md
@@ -15,6 +15,14 @@ Releases follow the `major.minor.micro` scheme recommended by [PEP440](https://w
 - [PR #1284](https://github.com/openforcefield/openforcefield/pull/1284): Fixes
   [Issue #1283](https://github.com/openforcefield/openff-toolkit/issues/1283) - force fields containing `BondCharge` virtual sites cannot be loaded due to an issue with how `outOfPlaneAngle` and `inPlaneAngle` keywords are validated
 
+### Minor bugfixes
+- [PR #1290](https://github.com/openforcefield/openforcefield/pull/1290): Fixes
+  [Issue #1216](https://github.com/openforcefield/openff-toolkit/issues/1216) by adding internal logic to handle
+  the possibility that multiple vsites share the same parent atom, and makes the return value of 
+  `VirtualSiteHandler.find_matches` be closer to the base class.
+
+
+
 ## 0.10.5 Bugfix release
 
 - [PR #1252](https://github.com/openforcefield/openforcefield/pull/1252): Refactors virtual 
diff --git a/openff/toolkit/tests/test_forcefield.py b/openff/toolkit/tests/test_forcefield.py
index e7cf3d7f7..fd4fc324d 100644
--- a/openff/toolkit/tests/test_forcefield.py
+++ b/openff/toolkit/tests/test_forcefield.py
@@ -2181,6 +2181,49 @@ def test_hash_strip_ids(self):
 
         assert hash(ff_with_id) == hash(ff_without_id)
 
+    def test_issue_1216(self):
+        """Test that the conflict between vsitehandler and librarychargehandler
+        identified in issue #1216 remains resolved.
+
+        https://github.com/openforcefield/openff-toolkit/issues/1216
+        """
+        force_field = ForceField()
+        force_field.get_parameter_handler("Electrostatics")
+
+        vsite_handler: VirtualSiteHandler = force_field.get_parameter_handler(
+            "VirtualSites"
+        )
+        vsite_handler.add_parameter(
+            {
+                "smirks": "[#6:1][#9:2]",
+                "name": "EP",
+                "type": "BondCharge",
+                "distance": 1.0 * unit.angstrom,
+                "match": "all_permutations",
+                "charge_increment1": 0.2 * unit.elementary_charge,
+                "charge_increment2": 0.1 * unit.elementary_charge,
+                "sigma": 1.0 * unit.angstrom,
+                "epsilon": 0.0 * unit.kilocalorie_per_mole,
+            }
+        )
+
+        library_handler: LibraryChargeHandler = force_field.get_parameter_handler(
+            "LibraryCharges"
+        )
+        library_handler.add_parameter(
+            {
+                "smirks": "[F:2][C:1]([H:3])([H:4])([H:5])",
+                "charge": [
+                    0.3 * unit.elementary_charge,
+                    -0.15 * unit.elementary_charge,
+                    -0.05 * unit.elementary_charge,
+                    -0.05 * unit.elementary_charge,
+                    -0.05 * unit.elementary_charge,
+                ],
+            }
+        )
+        force_field.label_molecules(Molecule.from_smiles("CF").to_topology())
+
 
 def generate_monatomic_ions():
     return (
diff --git a/openff/toolkit/tests/test_parameters.py b/openff/toolkit/tests/test_parameters.py
index 6ce07f4df..be53cdd24 100644
--- a/openff/toolkit/tests/test_parameters.py
+++ b/openff/toolkit/tests/test_parameters.py
@@ -2372,11 +2372,11 @@ def test_find_matches(
 
         matched_smirks = defaultdict(set)
 
-        for match in matches:
-
-            matched_smirks[match.environment_match.topology_atom_indices].add(
-                (match.parameter_type.smirks, match.parameter_type.name)
-            )
+        for match_list in matches.values():
+            for match in match_list:
+                matched_smirks[match.environment_match.topology_atom_indices].add(
+                    (match.parameter_type.smirks, match.parameter_type.name)
+                )
 
         assert {**matched_smirks} == expected_matches
 
diff --git a/openff/toolkit/typing/engines/smirnoff/forcefield.py b/openff/toolkit/typing/engines/smirnoff/forcefield.py
index e508de916..4c3b48e31 100644
--- a/openff/toolkit/typing/engines/smirnoff/forcefield.py
+++ b/openff/toolkit/typing/engines/smirnoff/forcefield.py
@@ -1486,14 +1486,13 @@ def label_molecules(self, topology):
         for molecule_idx, molecule in enumerate(topology.reference_molecules):
             top_mol = Topology.from_molecules([molecule])
             current_molecule_labels = dict()
-            param_is_list = False
             for tag, parameter_handler in self._parameter_handlers.items():
+                param_is_list = False
 
                 if type(parameter_handler) == VirtualSiteHandler:
                     param_is_list = True
 
                 matches = parameter_handler.find_matches(top_mol)
-                print(matches)
                 # Remove the chemical environment matches from the
                 # matched results.
 
diff --git a/openff/toolkit/typing/engines/smirnoff/parameters.py b/openff/toolkit/typing/engines/smirnoff/parameters.py
index 231a34314..46a7f217e 100644
--- a/openff/toolkit/typing/engines/smirnoff/parameters.py
+++ b/openff/toolkit/typing/engines/smirnoff/parameters.py
@@ -54,7 +54,7 @@
 import re
 from collections import OrderedDict, defaultdict
 from enum import Enum
-from typing import Any, List, Optional, Union, cast
+from typing import Any, List, Optional, Tuple, Union, cast
 
 try:
     import openmm
@@ -5151,13 +5151,12 @@ def _find_matches(
         entity: Topology,
         transformed_dict_cls=dict,
         unique=False,
-    ) -> List[ParameterHandler._Match]:
+    ) -> Dict[Tuple[int], List[ParameterHandler._Match]]:
 
         assigned_matches_by_parent = self._find_matches_by_parent(entity)
-        assigned_matches = []
-
+        return_dict = {}
         for parent_index, assigned_parameters in assigned_matches_by_parent.items():
-
+            assigned_matches = []
             for assigned_parameter, match_orientations in assigned_parameters:
 
                 for match in match_orientations:
@@ -5165,8 +5164,9 @@ def _find_matches(
                     assigned_matches.append(
                         ParameterHandler._Match(assigned_parameter, match)
                     )
+            return_dict[(parent_index,)] = assigned_matches
 
-        return assigned_matches
+        return return_dict
 
     def create_force(self, system: openmm.System, topology: Topology, **kwargs):