questions on RBFE examples

[blakemertz]

1. How is the displacement vector determined for each respective system? I noticed that for eralpha it is displacement=("22.0" "22.0" "22.0") whereas with CDK2 it is displacement=("-20.0" "0.0" "-20.0"). Your 2023 JCIM paper with Gianni states at least three layers of water molecules in between, which I guess is the general rule thumb. I am assuming the starting position of the protein-ligand complex should be centered at the origin to make the translation of the ligand in bulk solvent easier as well.
2. Is it sufficient to modify the setup-atm.sh script to utilize sdf files for GAFF-based force field generation or is it necessary to modify the provided python scripts for preparing ligands (make_atm_system_from_Amber.py, make_atm_system_from_pdb.py, and make_atm_system_from_rcpt_lig.py)? In the eralpha example .mol2 files are used for the ligands. I can convert an sdf to mol2, but most of my current workflows utilize sdf.
3. What tool(s) are used to construct the FEP map? That is one of the biggest advantages of schrodingers FEP+ -- the mapping tool does a pretty good job of a first approximation of what series of ligands should be paired together to get reasonable transformations between the two ligands. I know the Mobley group has developed HiMAP, but I didnt know if your group had used something else.

[egallicc]

1. Please see this about general considerations regarding the choice of displacement vector. Some corporate users developed specialized automated algorithms to choose the displacement to minimize the system's size.
2. The examples are just that, examples. AToM can work with any system set up for MD. The ERalpha example uses an Amber-based setup starting with mol2 files. If your workflow is based on sdf files, and you want to go with GAFF, I recommend using OpenFF and openmmforcefields as in the CDK2 example. Simply replace the openff-2.0.0 default ligand force field argument (--ligandForceField) of make_atm_system_from_rcpt_lig.py with GAFF.
3. Right, sorry, I do not have useful suggestions, as we do not routinely use automated FEP map generators. We are aware of corporate users who deployed their automated tools in conjunction with the AToM backend.

[blakemertz]

@egallicc thanks for the quick response -- apologies for not checking your group website for the documentation on displacement vectors. And of course for my third concern that is the beauty of using open-source tools -- gotta be willing to get your hands dirty on the coding side. Truly appreciate the AToM backend that is already in place, it definitely has a viable place in the landscape of free energy tools.

[blakemertz]

Reopening this issue as I had not been working on a complete solution until now. I followed the suggestion to substitute GAFF for SMIRNOFF in the --ligandForceField argument of make_atm_system_from_rcpt_lig.py. However, this also required me to make a call to the GAFFTemplateGenerator at the beginning of the python script (from openmmforcefields.generators import GAFFTemplateGenerator). However, when I went to execute setup-atm.sh, the call to make_atm_system_from_rcpt_lig.py crashes:

Call ForceField for protein and water
Receptor in PDB format
Number of atoms in receptor: 4717
Call Modeller: include receptor
Calculating receptor bounding box:
Areas of faces [Quantity(value=23.313064640000004, unit=nanometer**2), Quantity(value=20.937380080000004, unit=nanometer**2), Quantity(value=26.60136272, unit=nanometer**2)]
Direction of smallest area dimention: 1
Read ligand 1:
Number of atoms in ligand 1: 59
Call Modeller: include ligand 1
Read ligand 2:
Number of atoms in ligand 1: 59
Call Modeller: include ligand 2
Indexes of ligand 2 (starting from 0): [4776, 4777, 4778, 4779, 4780, 4781, 4782, 4783, 4784, 4785, 4786, 4787, 4788, 4789, 4790, 4791, 4792, 4793, 4794, 4795, 4796, 4797, 4798, 4799, 4800, 4801, 4802, 4803, 4804, 4805, 4806, 4807, 4808, 4809, 4810, 4811, 4812, 4813, 4814, 4815, 4816, 4817, 4818, 4819, 4820, 4821, 4822, 4823, 4824, 4825, 4826, 4827, 4828, 4829, 4830, 4831, 4832, 4833, 4834]
Calculating system bounding box:
boxVectors: (Quantity(value=Vec3(x=7.36895, y=0.0, z=0.0), unit=nanometer), Quantity(value=Vec3(x=0.0, y=7.4424, z=0.0), unit=nanometer), Quantity(value=Vec3(x=0.0, y=0.0, z=6.283600000000001), unit=nanometer))

Set up the combined protein + ligand + water system for simulation
Call the GAFFTemplateGenerator function for ligand
Traceback (most recent call last):
  File "/media/Data/binaries/github/AToM-OpenMM/make_atm_system_from_rcpt_lig.py", line 380, in <module>
    gaff = GAFFTemplateGenerator(molecules=ligandmolecules, forcefield=ligandforcefield, cache=ffcachefile )
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/Data/binaries/miniconda3/envs/atm8.1.1/lib/python3.12/site-packages/openmmforcefields/generators/template_generators.py", line 510, in __init__
    raise GAFFNotSupportedError(
openmmforcefields.generators.template_generators.GAFFNotSupportedError: This release (0.13.x) of openmmforcefields temporarily drops GAFF support and thereby the GAFFTemplateGenerator class. Support will be re-introduced in future releases (0.14.x). To use this class, install version 0.12.0 or older.

version 0.14.1 was recently released. Has the code base for AToM been updated to reflect this change in openmmforcefields, or will this require downgrading my environment to v. 0.12? Thanks.

[egallicc]

The examples are not tied to a particular version of openmmforcefields. Unless the API has changed, either downgrading or upgrading should work.

Please note that the examples and workflows we distribute with AToM-OpenMM are only examples of ways users can prepare inputs for AToM-OpenMM. They are not meant to be exhaustive or exclusive.