Net charge abfe -fix alchemical ion selection heuristic

[aqemia-benedict-tan]

Fixed some problems with the alchemical ion selection in the net_charge_abfe feature branch for an experimental implementation of ABFE calculations for charged ligands.

Issue 1: the previous implementation searches for ions by element name:

openfe/openfe/protocols/openmm_afe/equil_binding_afe_method.py

Line 1070 in d55d90c

if at.element in IONS[total_charge]

For larger ligands, this can end up selecting Cl atoms that have a distance >10 Angstrom to another atom in the ligand

Fix 1: select alchemical ions based on residue name ('NA', 'CL', 'K', et.c) - this excludes the ligand atoms as these have the residue name 'UNK'

Issue 2: the previous search volume for candidate ions uses a search shell with 1.0-1.1 nm distance radius. For a large protein, the closest ion can often be >1.1 nm away from the ligand - this returns no candidate ions.

Fix 2: we can keep the 1.0 nm minimum distance but expand the maximum limit to a PBC safe distance. The safe limit is obtained by reading in the unit cell vectors and computing the cell inradius. The inradius minus a small buffer should ensure that we only select ions within the unit cell. The implementation should work for any triclinic cell.

Testing:
The new _get_alchemical_ion method was tested using
test_ion_selection.py . Ion selection was performed for the ejm_52charg and jmc_32charg ligands from the OpenFE Benchmark (equilibration performed using the default AbsoluteBindingProtocol. The test verified the following:

• A suitable counterion is chosen (NA for a negative ligand, CL for a positive ligand)
• The selected alchemical ion is not a ligand atom
• The selected ion is within the distance cutoff, i.e. within the unit cell inradius

### ejm_52charg
  ligand total_charge = +1  ->  expecting counter-ion 'CL'
  unit cell [A,B,C,α,β,γ] : [78.31, 78.31, 78.31, 60.0, 60.0, 90.0]  (Å, deg)
  box vectors (Å):
      a = [  78.310    0.000    0.000]
      b = [   0.000   78.310    0.000]
      c = [  39.155   39.155   55.373]
  centre->edge (inradius) : 27.69 Å = 2.769 nm  (min-image radius; ion-search max = this - 0.1 nm)
  selected atom index : 34084
  resname / resid     : CL 10119
  distance to ligand  : 1.81 nm
  CHECK is expected counter-ion (CL): True
  CHECK not a ligand atom (UNK)               : True

### jmc_32charg
  ligand total_charge = -1  ->  expecting counter-ion 'NA'
  unit cell [A,B,C,α,β,γ] : [78.12, 78.12, 78.12, 60.0, 60.0, 90.0]  (Å, deg)
  box vectors (Å):
      a = [  78.124    0.000    0.000]
      b = [   0.000   78.124    0.000]
      c = [  39.062   39.062   55.242]
  centre->edge (inradius) : 27.62 Å = 2.762 nm  (min-image radius; ion-search max = this - 0.1 nm)
  selected atom index : 34049
  resname / resid     : NA 10089
  distance to ligand  : 0.87 nm
  CHECK is expected counter-ion (NA): True
  CHECK not a ligand atom (UNK)               : True

Note: the current 1.0 nm minimum distance applies to the pairwise ligand-counterion distances over all ligand atoms scanned using FindHostAtoms. This means that ions with a distance of <1.0 nm can be selected, as is the case for jmc_32charg (0.87 nm). This may be confusing for the user. We could easily fix this by filtering the atoms pairs explicitly with a minimum distance cutoff.

Developers certificate of origin

• I certify that this contribution is covered by the MIT License here and the Developer Certificate of Origin at https://developercertificate.org/.