Repository: qubekit/QUBEKit (language: python)

Issue:
XML combiner 
Now we know the xml combiner [here](https://github.com/CRingrose94/QUBEKit2Testing/blob/master/scripts/xml_combiner.py) works it would be good to incorporate this into the qubekit CLI. I think it might make sense to expose the available fitting parameters as options on the CLI so users who want to fit the Rfree for elements H, C, N, O for a selection of molecules could use a command like 
`qubekit combine --output master.xml -p c -p h -p n -p o`
where p is parameterise?

All Rfree would be options along with alpha and beta and they would start the optimisation from the current values, ie the ones used in the initial protocol.

## Questions
- [ ] should we have an option to reset the rfree and alpha and beta to some other values or just alpha and beta as each optimisation in the paper started from `alpha=1, beta=0.5`.


Discussion / hints:
We could also make the parameterise option in the command optional meaning by default it would add parameterise tags to all Rfree present in the set of molecules and alpha and beta by default?


Fix patch:
diff --git a/devtools/conda-envs/basic.yaml b/devtools/conda-envs/basic.yaml
index 2a8ede41..e4e772f0 100644
--- a/devtools/conda-envs/basic.yaml
+++ b/devtools/conda-envs/basic.yaml
@@ -10,6 +10,7 @@ dependencies:
     # Testing
   - pytest
   - pytest-cov
+  - deepdiff
 
   # Core deps
   - openff-toolkit-base >=0.9.2
diff --git a/devtools/conda-envs/everything.yaml b/devtools/conda-envs/everything.yaml
index 9e91f89b..cb8ecabc 100644
--- a/devtools/conda-envs/everything.yaml
+++ b/devtools/conda-envs/everything.yaml
@@ -10,6 +10,7 @@ dependencies:
   # Testing
   - pytest
   - pytest-cov
+  - deepdiff
 
     # Core deps
   - openff-toolkit-base >=0.9.2
diff --git a/qubekit/cli/cli.py b/qubekit/cli/cli.py
index fcc0cdf1..46dcf7b6 100644
--- a/qubekit/cli/cli.py
+++ b/qubekit/cli/cli.py
@@ -2,6 +2,7 @@
 
 import qubekit
 from qubekit.cli.bulk import bulk
+from qubekit.cli.combine import combine
 from qubekit.cli.config import config
 from qubekit.cli.progress import progress
 from qubekit.cli.run import restart, run
@@ -18,6 +19,7 @@ def cli():
 cli.add_command(restart)
 cli.add_command(bulk)
 cli.add_command(progress)
+cli.add_command(combine)
 
 if __name__ == "__main__":
     cli()
diff --git a/qubekit/cli/combine.py b/qubekit/cli/combine.py
new file mode 100644
index 00000000..268add52
--- /dev/null
+++ b/qubekit/cli/combine.py
@@ -0,0 +1,379 @@
+import os
+import xml.etree.ElementTree as ET
+from datetime import datetime
+from typing import TYPE_CHECKING, Dict, List, Optional, Tuple, Union
+from xml.dom.minidom import parseString
+
+import click
+
+import qubekit
+from qubekit.utils import constants
+from qubekit.workflow import WorkFlowResult
+
+if TYPE_CHECKING:
+    from qubekit.molecules import Atom, Ligand
+
+
+elements = ["H", "C", "N", "O", "X", "Cl", "S", "F", "Br", "P", "I", "B", "Si", "AB"]
+parameters_to_fit = click.Choice(
+    elements,
+    case_sensitive=True,
+)
+
+
+@click.command()
+@click.argument("filename", type=click.STRING)
+@click.option(
+    "-p",
+    "--parameters",
+    type=parameters_to_fit,
+    help="The elements whose Rfree is to be optimised, if not provided all will be fit.",
+    multiple=True,
+)
+def combine(filename: str, parameters: Optional[List[str]] = None):
+    """
+    Combine a list of molecules together and create a single master XML force field file.
+    """
+
+    molecules, rfrees = _find_molecules_and_rfrees()
+    if not parameters:
+        # fit everything
+        parameters = elements
+
+    xml_data = _combine_molecules(
+        molecules=molecules, rfree_data=rfrees, parameters=parameters
+    ).getroot()
+    messy = ET.tostring(xml_data, "utf-8")
+
+    pretty_xml = parseString(messy).toprettyxml(indent="")
+
+    with open(filename, "w") as xml_doc:
+        xml_doc.write(pretty_xml)
+
+
+def _combine_molecules(
+    molecules: List["Ligand"],
+    parameters: List[str],
+    rfree_data: Dict[str, Dict[str, Union[str, float]]],
+) -> ET.ElementTree:
+    """
+    Main worker function used to combine  the molecule force fields.
+    """
+    fit_ab = False
+    # if alpha and beta should be fit
+    if "AB" in parameters:
+        fit_ab = True
+
+    rfree_codes = set()  # keep track of all rfree codes used by these molecules
+
+    root = ET.Element("ForceField")
+    ET.SubElement(
+        root,
+        "QUBEKit",
+        attrib={
+            "Version": qubekit.__version__,
+            "Date": datetime.now().strftime("%Y_%m_%d"),
+        },
+    )
+    AtomTypes = ET.SubElement(root, "AtomTypes")
+    Residues = ET.SubElement(root, "Residues")
+
+    force_by_type = {}
+
+    increment = 0
+
+    # loop over molecules and create forces as they are found
+    for molecule in molecules:
+
+        resiude = ET.SubElement(Residues, "Residue", name=molecule.name)
+
+        for atom in molecule.atoms:
+            atom_type = f"QUBE_{atom.atom_index + increment}"
+            ET.SubElement(
+                AtomTypes,
+                "Type",
+                attrib={
+                    "name": atom_type,
+                    "class": str(atom.atom_index + increment),
+                    "element": atom.atomic_symbol,
+                    "mass": str(atom.atomic_mass),
+                },
+            )
+
+            ET.SubElement(
+                resiude,
+                "Atom",
+                attrib={"name": str(atom.atom_index + increment), "type": atom_type},
+            )
+
+        for (
+            i,
+            site,
+        ) in enumerate(molecule.extra_sites, start=1 + increment):
+            site_name = f"v-site{i}"
+            site_class = f"X{i}"
+            ET.SubElement(
+                AtomTypes,
+                "Type",
+                attrib={"name": site_name, "class": site_class, "mass": "0"},
+            )
+            ET.SubElement(
+                resiude, "Atom", attrib={"name": site_class, "type": site_name}
+            )
+
+        if molecule.BondForce.openmm_group() not in force_by_type:
+            new_bond_force = ET.SubElement(
+                root,
+                molecule.BondForce.openmm_group(),
+                attrib=molecule.BondForce.xml_data(),
+            )
+            # add to dict
+            force_by_type[molecule.BondForce.openmm_group()] = new_bond_force
+
+        BondForce = force_by_type[molecule.BondForce.openmm_group()]
+        for parameter in molecule.BondForce:
+            ET.SubElement(
+                BondForce,
+                parameter.openmm_type(),
+                attrib=_update_increment(
+                    force_data=parameter.xml_data(), increment=increment
+                ),
+            )
+            ET.SubElement(
+                resiude,
+                "Bond",
+                attrib={"from": str(parameter.atoms[0]), "to": str(parameter.atoms[1])},
+            )
+
+        if molecule.AngleForce.openmm_group() not in force_by_type:
+            new_angle_force = ET.SubElement(
+                root,
+                molecule.AngleForce.openmm_group(),
+                attrib=molecule.AngleForce.xml_data(),
+            )
+            force_by_type[molecule.AngleForce.openmm_group()] = new_angle_force
+
+        AngleForce = force_by_type[molecule.AngleForce.openmm_group()]
+        for parameter in molecule.AngleForce:
+            ET.SubElement(
+                AngleForce,
+                parameter.openmm_type(),
+                attrib=_update_increment(
+                    force_data=parameter.xml_data(), increment=increment
+                ),
+            )
+
+        if molecule.TorsionForce.openmm_group() not in force_by_type:
+            new_torsion_force = ET.SubElement(
+                root,
+                molecule.TorsionForce.openmm_group(),
+                attrib=molecule.TorsionForce.xml_data(),
+            )
+            force_by_type[molecule.TorsionForce.openmm_group()] = new_torsion_force
+
+        TorsionForce = force_by_type[molecule.TorsionForce.openmm_group()]
+
+        for parameter in molecule.TorsionForce:
+            ET.SubElement(
+                TorsionForce,
+                parameter.openmm_type(),
+                attrib=_update_increment(
+                    force_data=parameter.xml_data(), increment=increment
+                ),
+            )
+        for parameter in molecule.ImproperTorsionForce:
+            ET.SubElement(
+                TorsionForce,
+                parameter.openmm_type(),
+                attrib=_update_increment(
+                    force_data=parameter.xml_data(), increment=increment
+                ),
+            )
+
+        # not common so check if we need the section
+        if (
+            molecule.RBTorsionForce.n_parameters > 0
+            or molecule.ImproperRBTorsionForce.n_parameters > 0
+        ):
+            if molecule.RBTorsionForce.openmm_group() not in force_by_type:
+                new_rb_torsion_force = ET.SubElement(
+                    root,
+                    molecule.RBTorsionForce.openmm_group(),
+                    attrib=molecule.RBTorsionForce.xml_data(),
+                )
+                force_by_type[
+                    molecule.RBTorsionForce.openmm_group()
+                ] = new_rb_torsion_force
+
+            RBTorsion = force_by_type[molecule.RBTorsionForce.openmm_group()]
+            for parameter in molecule.RBTorsionForce:
+                ET.SubElement(
+                    RBTorsion,
+                    parameter.openmm_type(),
+                    attrib=_update_increment(
+                        force_data=parameter.xml_data(), increment=increment
+                    ),
+                )
+            for parameter in molecule.ImproperRBTorsionForce:
+                ET.SubElement(
+                    RBTorsion,
+                    parameter.openmm_type(),
+                    attrib=_update_increment(
+                        force_data=parameter.xml_data(), increment=increment
+                    ),
+                )
+
+        for i, site in enumerate(molecule.extra_sites):
+            site_data = site.xml_data()
+            site_data["index"] = str(i + molecule.n_atoms)
+            ET.SubElement(resiude, site.openmm_type(), attrib=site_data)
+
+        if molecule.NonbondedForce.openmm_group() not in force_by_type:
+            new_nb_force = ET.SubElement(
+                root,
+                molecule.NonbondedForce.openmm_group(),
+                attrib=molecule.NonbondedForce.xml_data(),
+            )
+            force_by_type[molecule.NonbondedForce.openmm_group()] = new_nb_force
+
+        NonbondedForce = force_by_type[molecule.NonbondedForce.openmm_group()]
+        for parameter in molecule.NonbondedForce:
+            # work out if the atom is being optimised
+            rfree_code = _get_parameter_code(
+                molecule=molecule, atom_index=parameter.atoms[0]
+            )
+            atom_data = {
+                "charge": str(parameter.charge),
+                "sigma": str(parameter.sigma),
+                "epsilon": str(parameter.epsilon),
+                "type": f"QUBE_{str(parameter.atoms[0] + increment)}",
+            }
+            if rfree_code in parameters or fit_ab:
+                # keep track of present codes to optimise
+                rfree_codes.add(rfree_code)
+                # this is to be refit
+                atom = molecule.atoms[parameter.atoms[0]]
+                eval_string = _get_eval_string(
+                    atom=atom,
+                    rfree_data=rfree_data[rfree_code],
+                    a_and_b=fit_ab,
+                    rfree_code=rfree_code if rfree_code in parameters else None,
+                )
+                atom_data["parameter_eval"] = eval_string
+                atom_data["volume"] = str(atom.aim.volume)
+                atom_data["bfree"] = str(rfree_data[rfree_code]["b_free"])
+                atom_data["vfree"] = str(rfree_data[rfree_code]["v_free"])
+
+            ET.SubElement(NonbondedForce, parameter.openmm_type(), attrib=atom_data)
+
+        for i, site in enumerate(molecule.extra_sites, start=1 + increment):
+            site_name = f"v-site{i}"
+            ET.SubElement(
+                NonbondedForce,
+                "Atom",
+                attrib={
+                    "charge": str(site.charge),
+                    "epsilon": "0",
+                    "sigma": "1",
+                    "type": site_name,
+                },
+            )
+        increment += molecule.n_atoms
+
+    # add the ForceBalance tags
+    ForceBalance = ET.SubElement(root, "ForceBalance")
+    for parameter_to_fit in parameters:
+        if parameter_to_fit != "AB" and parameter_to_fit in rfree_codes:
+            fit_data = {
+                f"{parameter_to_fit.lower()}free": str(
+                    rfree_data[parameter_to_fit]["r_free"]
+                ),
+                "parameterize": f"{parameter_to_fit.lower()}free",
+            }
+            ET.SubElement(ForceBalance, f"{parameter_to_fit}Element", attrib=fit_data)
+    if fit_ab:
+        ET.SubElement(
+            ForceBalance, "xalpha", alpha=str(rfree_data["alpha"]), parameterize="alpha"
+        )
+        ET.SubElement(
+            ForceBalance, "xbeta", beta=str(rfree_data["beta"]), parameterize="beta"
+        )
+
+    return ET.ElementTree(root)
+
+
+def _get_eval_string(
+    atom: "Atom",
+    rfree_data: Dict[str, str],
+    a_and_b: bool,
+    rfree_code: Optional[str] = None,
+) -> str:
+    """
+    Create the parameter eval string used by ForceBalance to calculate the sigma and epsilon parameters.
+    """
+
+    if a_and_b:
+        alpha = "PARM['xalpha/alpha']*"
+        beta = f"*({atom.aim.volume}/{rfree_data['v_free']})**PARM['xbeta/beta']"
+    else:
+        alpha, beta = "", ""
+    if rfree_code is not None:
+        rfree = f"PARM['{rfree_code}Element/{rfree_code.lower()}free']"
+    else:
+        rfree = f"{rfree_data['r_free']}"
+
+    eval_string = (
+        f"epsilon=({alpha}{rfree_data['b_free']}{beta})/(128*{rfree}**6)*{constants.EPSILON_CONVERSION}, "
+        f"sigma=2**(5/6)*({atom.aim.volume}/{rfree_data['v_free']})**(1/3)*{rfree}*{constants.SIGMA_CONVERSION}"
+    )
+
+    return eval_string
+
+
+def _find_molecules_and_rfrees() -> Tuple[
+    List["Ligand"], Dict[str, Dict[str, Union[str, float]]]
+]:
+    """
+    Loop over the local directories looking for qubekit WorkFlowResults and extract the ligands and a list of all
+    unique free params used to parameterise the molecules.
+    """
+    molecules = []
+    element_params = {}
+
+    for folder in os.listdir("."):
+        if os.path.isdir(folder) and folder.startswith("QUBEKit_"):
+            workflow_result = WorkFlowResult.parse_file(
+                os.path.join(folder, "workflow_result.json")
+            )
+            molecules.append(workflow_result.current_molecule)
+            lj_stage = workflow_result.results["non_bonded"].stage_settings
+            element_params.update(lj_stage["free_parameters"])
+            element_params["alpha"] = lj_stage["alpha"]
+            element_params["beta"] = lj_stage["beta"]
+
+    print(f"{len(molecules)} molecules found, combining...")
+    return molecules, element_params
+
+
+def _update_increment(force_data: Dict[str, str], increment: int) -> Dict[str, str]:
+    """
+    Update the increment on some force data.
+    """
+    for key, value in force_data.items():
+        if key.startswith("class"):
+            force_data[key] = str(int(value) + increment)
+
+    return force_data
+
+
+def _get_parameter_code(molecule: "Ligand", atom_index: int) -> str:
+    """
+    Work out the rfree code for an atom in a molecule, mainly used for special cases.
+    """
+    atom = molecule.atoms[atom_index]
+    rfree_code = atom.atomic_symbol
+    if rfree_code == "H":
+        # check if polar
+        if molecule.atoms[atom.bonds[0]].atomic_symbol in ["O", "N", "S"]:
+            rfree_code = "X"
+    return rfree_code
diff --git a/qubekit/data/pyridine.xml b/qubekit/data/pyridine.xml
index 7f7a2974..ffb6f77c 100644
--- a/qubekit/data/pyridine.xml
+++ b/qubekit/data/pyridine.xml
@@ -99,7 +99,7 @@
 <Improper class1="H809" class2="C804" class3="N803" class4="C805" k1="0.000000" k2="10.460000" k3="0.000000" k4="0.000000" periodicity1="1" periodicity2="2" periodicity3="3" periodicity4="4" phase1="0.00" phase2="3.141592653589793" phase3="0.00" phase4="3.141592653589793"/>
 <Improper class1="C800" class2="C805" class3="H810" class4="C804" k1="0.000000" k2="10.460000" k3="0.000000" k4="0.000000" periodicity1="1" periodicity2="2" periodicity3="3" periodicity4="4" phase1="0.00" phase2="3.141592653589793" phase3="0.00" phase4="3.141592653589793"/>
 </PeriodicTorsionForce>
-<NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
+<NonbondedForce coulomb14scale="0.8333333333" lj14scale="0.5">
 <Atom type="opls_800" charge="0.025900" sigma="0.343870" epsilon="0.258990" />
 <Atom type="opls_801" charge="-0.268800" sigma="0.370010" epsilon="0.258990" />
 <Atom type="opls_802" charge="0.237100" sigma="0.328250" epsilon="0.258990" />


Test patch:
diff --git a/qubekit/tests/cli/test_combine.py b/qubekit/tests/cli/test_combine.py
new file mode 100644
index 00000000..714a7046
--- /dev/null
+++ b/qubekit/tests/cli/test_combine.py
@@ -0,0 +1,273 @@
+import os
+import xml.etree.ElementTree as ET
+from xml.dom.minidom import parseString
+
+import pytest
+import xmltodict
+from deepdiff import DeepDiff
+from openmm import XmlSerializer, app
+
+import qubekit
+from qubekit.cli.combine import (
+    _combine_molecules,
+    _find_molecules_and_rfrees,
+    _get_eval_string,
+    _get_parameter_code,
+    _update_increment,
+    combine,
+    elements,
+)
+from qubekit.molecules import Ligand
+from qubekit.nonbonded import get_protocol
+from qubekit.utils import constants
+from qubekit.utils.file_handling import get_data
+
+
+@pytest.mark.parametrize(
+    "molecule, atom_index, expected",
+    [
+        pytest.param("water", 0, "O", id="O"),
+        pytest.param("water", 1, "X", id="Polar H"),
+        pytest.param("acetone", 4, "H", id="Normal H"),
+    ],
+)
+def test_get_param_code(molecule, atom_index, expected, request):
+    """
+    Test getting the rfree param code including polar h.
+    """
+    molecule = request.getfixturevalue(molecule)
+    rfree_code = _get_parameter_code(molecule=molecule, atom_index=atom_index)
+    assert rfree_code == expected
+
+
+def test_update_increment():
+    """
+    Make sure force data has the class incremented correctly.
+    """
+    force_data = {"class1": "1", "class2": "2", "k": "300"}
+    new_data = _update_increment(force_data=force_data, increment=100)
+    assert new_data["class1"] == "101"
+    assert new_data["class2"] == "102"
+
+
+def test_find_molecules_and_rfree(acetone, tmpdir, rdkit_workflow):
+    """
+    Test looping over folders to find rfree
+    """
+
+    with tmpdir.as_cwd():
+        # make some example folders
+        os.mkdir("QUBEKit_mol1")
+        result = rdkit_workflow._build_initial_results(molecule=acetone)
+        result.to_file(os.path.join("QUBEKit_mol1", "workflow_result.json"))
+
+        molecules, rfree = _find_molecules_and_rfrees()
+        assert len(molecules) == 1
+        assert len(rfree) == 7
+
+
+@pytest.mark.parametrize(
+    "atom, a_and_b, rfree_code, expected",
+    [
+        pytest.param(
+            1,
+            True,
+            "C",
+            f"epsilon=(PARM['xalpha/alpha']*46.6*(21.129491/34.4)**PARM['xbeta/beta'])/(128*PARM['CElement/cfree']**6)*{constants.EPSILON_CONVERSION}, sigma=2**(5/6)*(21.129491/34.4)**(1/3)*PARM['CElement/cfree']*{constants.SIGMA_CONVERSION}",
+            id="AB Rfree",
+        ),
+        pytest.param(
+            1,
+            False,
+            "C",
+            f"epsilon=(46.6)/(128*PARM['CElement/cfree']**6)*{constants.EPSILON_CONVERSION}, sigma=2**(5/6)*(21.129491/34.4)**(1/3)*PARM['CElement/cfree']*{constants.SIGMA_CONVERSION}",
+            id="No AB Rfree",
+        ),
+        pytest.param(
+            1,
+            True,
+            None,
+            f"epsilon=(PARM['xalpha/alpha']*46.6*(21.129491/34.4)**PARM['xbeta/beta'])/(128*2**6)*{constants.EPSILON_CONVERSION}, sigma=2**(5/6)*(21.129491/34.4)**(1/3)*2*{constants.SIGMA_CONVERSION}",
+            id="AB No Rfree",
+        ),
+        pytest.param(
+            1,
+            False,
+            None,
+            f"epsilon=(46.6)/(128*2**6)*{constants.EPSILON_CONVERSION}, sigma=2**(5/6)*(21.129491/34.4)**(1/3)*2*{constants.SIGMA_CONVERSION}",
+            id="No AB No Rfree",
+        ),
+    ],
+)
+def test_get_eval_string(atom, a_and_b, rfree_code, expected, coumarin):
+    """
+    Test generating the eval string for ForceBalance with different settings.
+    """
+    rfree_data = {"v_free": 34.4, "b_free": 46.6, "r_free": 2}  # carbon
+
+    eval_string = _get_eval_string(
+        atom=coumarin.atoms[atom],
+        rfree_data=rfree_data,
+        a_and_b=a_and_b,
+        rfree_code=rfree_code,
+    )
+    assert eval_string == expected
+
+
+def test_combine_molecules_deepdiff(acetone, openff, coumarin, tmpdir, rfree_data):
+
+    with tmpdir.as_cwd():
+        openff.run(acetone)
+        acetone_ref_system = xmltodict.parse(open("serialised.xml").read())
+        openff.run(coumarin)
+        coumarin_ref_system = xmltodict.parse(open("serialised.xml").read())
+
+        combined_xml = _combine_molecules(
+            molecules=[acetone, coumarin], parameters=elements, rfree_data=rfree_data
+        ).getroot()
+        messy = ET.tostring(combined_xml, "utf-8")
+
+        pretty_xml = parseString(messy).toprettyxml(indent="")
+
+        with open("combined.xml", "w") as xml_doc:
+            xml_doc.write(pretty_xml)
+        root = combined_xml.find("QUBEKit")
+        assert qubekit.__version__ == root.get("Version")
+
+        # load up new systems and compare
+        combinded_ff = app.ForceField("combined.xml")
+
+        acetone_combine_system = xmltodict.parse(
+            XmlSerializer.serialize(
+                combinded_ff.createSystem(
+                    acetone.to_openmm_topology(),
+                    nonbondedCutoff=0.9,
+                    removeCMMotion=False,
+                )
+            )
+        )
+        # slight differences in masses we need to ignore
+        acetone_diff = DeepDiff(
+            acetone_ref_system,
+            acetone_combine_system,
+            ignore_order=True,
+            significant_digits=6,
+            exclude_regex_paths="mass",
+        )
+        assert len(acetone_diff) == 0
+
+        coumarin_combine_system = xmltodict.parse(
+            XmlSerializer.serialize(
+                combinded_ff.createSystem(
+                    coumarin.to_openmm_topology(),
+                    nonbondedCutoff=0.9,
+                    removeCMMotion=False,
+                )
+            )
+        )
+        coumarin_diff = DeepDiff(
+            coumarin_ref_system,
+            coumarin_combine_system,
+            ignore_order=True,
+            significant_digits=6,
+            exclude_regex_paths="mass",
+        )
+        assert len(coumarin_diff) == 0
+
+
+def test_combine_molecules_sites_deepdiff(openff, xml, acetone, rfree_data, tmpdir):
+    """
+    Test combining molecules with virtual sites and ensure they are correctly applied and the energy break down matches.
+    """
+    with tmpdir.as_cwd():
+        openff.run(acetone)
+        acetone_ref_system = xmltodict.parse(open("serialised.xml").read())
+        pyridine = Ligand.from_file(file_name=get_data("pyridine.sdf"))
+        xml.run(molecule=pyridine, input_files=[get_data("pyridine.xml")])
+        pyridine_ref_system = xmltodict.parse(open("serialised.xml").read())
+
+        combined_xml = _combine_molecules(
+            molecules=[acetone, pyridine], parameters=elements, rfree_data=rfree_data
+        ).getroot()
+        messy = ET.tostring(combined_xml, "utf-8")
+
+        pretty_xml = parseString(messy).toprettyxml(indent="")
+        print(pretty_xml)
+        with open("combined.xml", "w") as xml_doc:
+            xml_doc.write(pretty_xml)
+        root = combined_xml.find("QUBEKit")
+        assert qubekit.__version__ == root.get("Version")
+
+        # load up new systems and compare
+        combinded_ff = app.ForceField("combined.xml")
+
+        acetone_combine_system = xmltodict.parse(
+            XmlSerializer.serialize(
+                combinded_ff.createSystem(
+                    acetone.to_openmm_topology(),
+                    nonbondedCutoff=0.9,
+                    removeCMMotion=False,
+                )
+            )
+        )
+        # slight differences in masses we need to ignore
+        acetone_diff = DeepDiff(
+            acetone_ref_system,
+            acetone_combine_system,
+            ignore_order=True,
+            significant_digits=6,
+            exclude_regex_paths="mass",
+        )
+        assert len(acetone_diff) == 0
+        # add v-site
+        pyridine_mod = app.Modeller(
+            pyridine.to_openmm_topology(), pyridine.openmm_coordinates()
+        )
+        pyridine_mod.addExtraParticles(combinded_ff)
+        pyridine_combine_system = xmltodict.parse(
+            XmlSerializer.serialize(combinded_ff.createSystem(pyridine_mod.topology))
+        )
+        pyridine_diff = DeepDiff(
+            pyridine_ref_system,
+            pyridine_combine_system,
+            ignore_order=True,
+            significant_digits=6,
+            exclude_regex_paths="mass",
+        )
+        assert len(acetone_diff) == 0
+
+
+@pytest.mark.parametrize(
+    "parameters, expected",
+    [
+        pytest.param(["-p", "C"], 1, id="Only C"),
+        pytest.param(
+            ["-p", "C", "-p", "AB", "-p", "N", "-p", "O"], 5, id="CNO alpha beta"
+        ),
+        pytest.param([], 7, id="All present"),
+    ],
+)
+def test_combine_cli_all(
+    run_cli, tmpdir, acetone, coumarin, openff, rdkit_workflow, parameters, expected
+):
+    """
+    Test combining via the cli with different parameters.
+    """
+    workflow = rdkit_workflow.copy(deep=True)
+    workflow.non_bonded = get_protocol(protocol_name="5b")
+    with tmpdir.as_cwd():
+        openff.run(coumarin)
+        openff.run(acetone)
+        # make some dummy dirs
+        os.mkdir("QUBEKit_acetone")
+        os.mkdir("QUBEKit_coumarin")
+        result = workflow._build_initial_results(molecule=acetone)
+        result.to_file(os.path.join("QUBEKit_acetone", "workflow_result.json"))
+        result = workflow._build_initial_results(molecule=coumarin)
+        result.to_file(os.path.join("QUBEKit_coumarin", "workflow_result.json"))
+        output = run_cli.invoke(combine, args=["combined.xml", *parameters])
+        assert output.exit_code == 0
+        assert "2 molecules found, combining..." in output.output
+
+        data = xmltodict.parse(open("combined.xml").read())
+        assert len(data["ForceField"]["ForceBalance"]) == expected
diff --git a/qubekit/tests/conftest.py b/qubekit/tests/conftest.py
index 5e04b552..a644d529 100644
--- a/qubekit/tests/conftest.py
+++ b/qubekit/tests/conftest.py
@@ -50,3 +50,20 @@ def rdkit_workflow():
     rdkit_spec = QCOptions(program="rdkit", method="uff", basis=None)
     workflow = WorkFlow(qc_options=rdkit_spec)
     return workflow
+
+
+@pytest.fixture()
+def rfree_data():
+    return {
+        "H": {"v_free": 7.6, "b_free": 6.5, "r_free": 1.738},
+        "X": {"v_free": 7.6, "b_free": 6.5, "r_free": 1.083},
+        "C": {"v_free": 34.4, "b_free": 46.6, "r_free": 2.008},
+        "N": {"v_free": 25.9, "b_free": 24.2, "r_free": 1.765},
+        "O": {"v_free": 22.1, "b_free": 15.6, "r_free": 1.499},
+        "F": {"v_free": 18.2, "b_free": 9.5, "r_free": 1.58},
+        "Cl": {"v_free": 65.1, "b_free": 94.6, "r_free": 1.88},
+        "Br": {"v_free": 95.7, "b_free": 162.0, "r_free": 1.96},
+        "S": {"v_free": 75.2, "b_free": 134.0, "r_free": 2.0},
+        "alpha": 1,
+        "beta": 0.5,
+    }