Ligand parameterization for GROMACS constant pH simulations

This is a notebook to prepare input files for running GROMACS constant pH (1, 2) simulations for ionizable ligands. The notebook should be run in GoogleColab, but can be adapted for jupyter, etc.

Implemented pipeline

1. Defines the protonation states of the ligand (provided as SMILES string) and allows the user to select states to be used in the simulation. Protonation states are defined using Dimorphite tool, which creates SMILES for each possible protonation state. Dimorphite is running a substructure search over the small molecule using RdKit. The idea is based on searching each of the ionizable sites from a predefined library in the molecule and assigning all possible protonation states to each site found. As a result of this step, the notebook outputs all possible protonation states for the provided SMILES string and asks the user to select the states to be used in the simulation.
2. For the selected states all possible symmetrical states are created. Symmetrical states can be of two types:
   • Sites with multiple hydrogens (e.g. -NH$_3$). Such groups might have several hydrogens detached in the deprotonated state. While it is not usually a problem for SMILES representation of a molecule, in 3D MD simulations (especially in a protein environment) it might be quite important (e.g. one of the hydrogens might be involved in the hydrogen bonding, and making it disappear can significantly affect the results.) Therefore, those states can be treated independently.
   • Actually symmetrical sites: e.g. first and third oxygen atoms in SMILES [O:1]-[C:2]-[O:3]. At the SMILES level, deprotonating the first or the third oxygen makes no difference, leading to the same molecule. But as in the case of multiple hydrogens bonded to the site, in the more complex environment, one site can accept/donate proton much easier than the second one. Again, those states can be treated independently. After detecting all the symmetrical states, the user is asked to select whether or not symmetry has to be considered for each of the states with symmetrical variants. Symmetry search is implemented using RdKit.
3. After determining all the states to be used in simulations transitions between states are defined and the microscopic p$K_{\text{a}}$ values for each state have to be assigned. We recommend to use experimental p$K_{\text{a}}$ values. If those are not available, QM calculations can be useful in determining p$K_{\text{a}}$ values, as reported in SAMPL challenge.
4. When p$K_{\text{a}}$ values are assigned, the $\lambda$-group parameters should be assigned. This is done automatically for you. All the computations are done using SyMpy package.
5. Finally, the parameters have to be assigned to prepare lambdagrouptypes.dat file for an ionizable ligand of interest. lambdagrouptypes.dat is an essential input file for a phbuilder tool, which automates the preparation of input files for GROMACS constant pH simulations. At the moment, the parameters are assigned using Open Force Field. TODO Support for user-provided parameters to be added soon.