virtual-sites-matches.py

from openff.toolkit.tests.create_molecules import create_dinitrogen
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.units import unit

xml_snippet = """<?xml version="1.0" encoding="utf-8"?>
<SMIRNOFF version="0.3" aromaticity_model="OEAroModel_MDL">
    <VirtualSites version="0.3">
        <VirtualSite
            type="BondCharge"
            name="EP1"
            smirks="[*:1]~[*:2]"
            distance="0.111*angstrom"
            charge_increment1="0.1*elementary_charge"
            charge_increment2="0.1*elementary_charge"
            sigma="0.1*angstrom"
            epsilon="0.1*kilocalories_per_mole"
            match="once" >
        </VirtualSite>
        <VirtualSite
            type="BondCharge"
            name="EP2"
            smirks="[#7:1]~[#7:2]"
            distance="0.222*angstrom"
            charge_increment1="0.2*elementary_charge"
            charge_increment2="0.2*elementary_charge"
            sigma="0.2*angstrom"
            epsilon="0.2*kilocalories_per_mole"
            match="all_permutations" >
        </VirtualSite>
        <VirtualSite
            type="BondCharge"
            name="EP3"
            smirks="[#7:1]~[#7:2]"
            distance="0.333*angstrom"
            charge_increment1="0.2*elementary_charge"
            charge_increment2="0.2*elementary_charge"
            sigma="0.2*angstrom"
            epsilon="0.2*kilocalories_per_mole"
            match="once" >
        </VirtualSite>
    </VirtualSites>
</SMIRNOFF>
"""

box_vectors = unit.Quantity([[4, 0, 0], [0, 4, 0], [0, 0, 4]], units=unit.nanometer)

molecule = create_dinitrogen()

topology = molecule.to_topology()
topology.box_vectors = box_vectors

force_field = ForceField("openff-2.0.0.offxml", xml_snippet)

print("names and distances of bond charge parameters:")
for parameter in force_field["VirtualSites"].parameters:
    print(parameter.name, parameter.distance)

original_matches = force_field["VirtualSites"].find_matches(topology)

print("atom tuples and parameter distances after find_matches")
for key, val in original_matches.items():
    print(key, val[0].parameter_type.name, val[0].parameter_type.distance)
# Some later post-processing (including applying `match="once"` is done here,
# but we're already in a state in which (0, 1) and (1, 0) have different
# *parameters*, so just picking the first pair of atoms could result in assigning
# different parameters besed on the order of the `matches` dict

reduced_matches = force_field["VirtualSites"]._reduce_virtual_particles_to_sites(
    original_matches
)