`PDBFixer.addSolvent` often incorrectly neutralises systems with charged ligands as though they were neutral

[Yoshanuikabundi]

The addSolvent method uses a force field to assign charges, which it then uses to add ions that neutralise the net charge of the box. Unfortunately, the force field seems to assign formal charges of zero to many charged groups in small molecules; for example the citrate ion, which has a -3 formal charge:

from openff.toolkit import Molecule

top = Molecule.from_smiles("C(C(=O)[O-])C(CC(=O)[O-])(C(=O)[O-])O").to_topology()
top.molecule(0).generate_conformers(n_conformers=1)
top.to_file("citrate.pdb")

from pdbfixer import PDBFixer
import openmm

fixer = PDBFixer("citrate.pdb")

fixer.addSolvent(
    padding=2.0 * openmm.unit.nanometer,
    ionicStrength=0.1 * openmm.unit.molar,
    boxShape="dodecahedron",
)

counts = {"NA": 0, "UNK": 0, "CL": 0, "HOH": 0}
for residue in fixer.topology.residues():
    counts[residue.name] += 1
counts # 
{'NA': 3, 'UNK': 1, 'CL': 3, 'HOH': 1457}

The citrate.pdb file generated above is available here: citrate.pdb.zip

With the current main branch (commit 92044f3), you can also see this behavior with a PDB file like 1PO0, which has 2 citrate ions (each with -3 formal charge).

from pdbfixer import PDBFixer
import openmm

fixer = PDBFixer("1PO0.pdb")

fixer.findMissingResidues()

# This section adds caps; leave it commented to rebuild terminal loops
chains_to_cap = {chain for chain, resi in fixer.missingResidues}
for chainidx in chains_to_cap:
    chain = [*fixer.topology.chains()][chainidx]
    last_resi = len([*chain.residues()])
    fixer.missingResidues[chainidx, 0] = ['ACE']
    fixer.missingResidues[chainidx, last_resi] = ['NME']

fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()
# fixer.removeHeterogens(keepWater=True)
fixer.findMissingAtoms()
fixer.addMissingAtoms()

fixer.addMissingHydrogens(pH=7.4)

fixer.addSolvent(
    padding=2.0 * openmm.unit.nanometer,
    ionicStrength=0.1 * openmm.unit.molar,
    boxShape="dodecahedron",
)

charges = {
 'ASP': -1,
 'ARG': +1,
 'GLU': -1,
 'LYS': +1,
 'FLC': -3,
 'NA': 1,
 'CL': -1
} # I think this is right!
charge = 0
for residue in fixer.topology.residues():
    if residue.name in charges:
        charge += charges[residue.name]
charge # gives -6 (-3 for each citrate residue)

It would be great if the correct formal charges could be stored when the template is downloaded! Then they could be used when neutralising the system. This might also give an opportunity for users to manually specify the formal charges of UNK residues.

Thanks!

(also the ACE and NME residues in the second example above don't seem to be getting hydrogens - let me know if you'd like me to open a second issue for that!)