-
Notifications
You must be signed in to change notification settings - Fork 19
Description
I attempted to generate the .top and .gro files for a relative FEP calculation in GROMACS, where the perturbation is between two ligands that differ only in their stereochemistry (one is the R-enantiomer and the other is the S-enantiomer). I used the following Python code:
"""
import BioSimSpace as BSS
import os
define all the folder locations
main_folder = os.getcwd()
prot_p = BSS.IO.readMolecules(["RORa.gas.receptor.prmtop", "RORa.gas.receptor.rst7"])
print(prot_p)
ligand3S = BSS.IO.readMolecules(["RORa_ligand3S_bcc.acpype/RORa_ligand3S_bcc_AC.prmtop", "RORa_ligand3S_bcc.acpype/RORa_ligand3S_bcc_AC.inpcrd"])
ligand3S_prep = ligand3S.getMolecules()[0]
print(ligand3S_prep)
ligand3R = BSS.IO.readMolecules(["RORa_ligand3R_bcc.acpype/RORa_ligand3R_bcc_AC.prmtop", "RORa_ligand3R_bcc.acpype/RORa_ligand3R_bcc_AC.inpcrd"])
ligand3R_prep = ligand3R.getMolecules()[0]
print(ligand3R_prep)
mapping = BSS.Align.matchAtoms(ligand3R_prep, ligand3S_prep)
print(mapping)
inv_mapping = {v: k for k, v in mapping.items()}
ligand3S_aligned = BSS.Align.rmsdAlign(ligand3S_prep, ligand3R_prep, inv_mapping)
merged = BSS.Align.merge(ligand3R_prep, ligand3S_aligned, mapping)
merged_system = prot_p + merged
BSS.IO.saveMolecules("new_inputs_ligand3R_ligand3S/holo_ligand3R_ligand3S", merged_system, ["gro87", "grotop"])
"""
However, the topology file ends up with no dummy atoms, so no atoms are actually being perturbed. Additionally, the .gro file contains only the atoms from one state. Which function should I use in this situation? Is there a way to perform distance-based alignment between atoms, or to indicate—perhaps by using different atom names—which atoms I want to be perturbed?
Thanks in advance for your time.