Merge pull request #124 from cbc-univie/amber_rbfe

reading in the two structures

transformato/bin/openmm_run.py

@@ -41,7 +41,7 @@
     top = gen_box(top, crd)
 
 # Build system
-if env == "waterbox":
+if env == "waterbox" or env == "complex":
     nboptions = dict(
         nonbondedMethod=inputs.coulomb,
         nonbondedCutoff=inputs.r_off * unit.nanometers,
@@ -53,6 +53,7 @@
         nonbondedMethod=NoCutoff,
         constraints=inputs.cons,
     )
+print(f"Applying the following nonbonded options {nboptions}")
 
 if inputs.vdw == "Switch" and env != "vacuum":
     print(f"Setting the vdw switching function to the defalut Openmm Switch")

transformato/mutate.py

@@ -871,6 +871,13 @@ def _prepare_cc_for_charge_transfer(self):
             [self.get_common_core_idx_mol1(), self.get_common_core_idx_mol2()],
             [self.dummy_region_cc1, self.dummy_region_cc2],
         ):
+            ## We need this for point mutations, because if we give a resid, the mol here
+            ## consists only of on residue which resid is always 1
+            try:
+                int(tlc)
+                tlc = "1"
+            except ValueError:
+                tlc = tlc
             # set `initial_charge` parameter for Mutation
             for atom in psf.view[f":{tlc}"].atoms:
                 # charge, epsilon and rmin are directly modiefied
@@ -1474,8 +1481,13 @@ def _mutate_to_common_core(
         """
 
         mutations = defaultdict(list)
-        tlc = self.s1_tlc
-
+        ## We need this for point mutations, because if we give a resid, the mol here
+        ## consists only of on residue which resid is always 1
+        try:
+            int(self.s1_tlc)
+            tlc = "1"
+        except ValueError:
+            tlc = self.s1_tlc
         if self.asfe:
             psf = self.psf1["waterbox"]
             cc_idx = []  # no CC in ASFE
@@ -1496,7 +1508,13 @@ def _mutate_to_common_core(
             logger.info(f"Terminal dummy atoms: {list_termin_dummy_atoms}")
 
             if mol_name == "m2":
-                tlc = self.s2_tlc
+                ## We need this for point mutations, because if we give a resid, the mol here
+                ## consists only of on residue which resid is always 1
+                try:
+                    int(self.s2_tlc)
+                    tlc = "1"
+                except ValueError:
+                    tlc = self.s2_tlc
 
         # iterate through atoms and select atoms that need to be mutated
         atoms_to_be_mutated = []
@@ -1997,7 +2015,16 @@ def _mutate_torsions(self, psf: pm.charmm.CharmmPsfFile, lambda_value: float):
                     f = max((1 - ((1 - lambda_value) * 2)), 0.0)
 
                     if f > 0.0 or lambda_value == 0.5:
-                        for torsion_t in original_torsion.type:
+                        ## Necessary, because in Amber topologies this is not a list
+                        if (
+                            type(original_torsion.type)
+                            == pm.topologyobjects.DihedralType
+                        ):
+                            orig_torsion_as_list = [original_torsion.type]
+                        else:
+                            orig_torsion_as_list = original_torsion.type
+
+                        for torsion_t in orig_torsion_as_list:
                             modified_phi_k = torsion_t.phi_k * f
                             mod_types.append(
                                 mod_type(
@@ -2012,7 +2039,14 @@ def _mutate_torsions(self, psf: pm.charmm.CharmmPsfFile, lambda_value: float):
                     # torsion present at cc2 needs to be fully turned on at lambda_value == 0.0
                     f = 1 - min((lambda_value) * 2, 1.0)
                     if f > 0.0:
-                        for torsion_t in new_torsion.type:
+
+                        ## Necessary, because in Amber topologies this is not a list
+                        if type(new_torsion.type) == pm.topologyobjects.DihedralType:
+                            new_torsion_as_list = [new_torsion.type]
+                        else:
+                            new_torsion_as_list = new_torsion.type
+
+                        for torsion_t in new_torsion_as_list:
                             modified_phi_k = torsion_t.phi_k * f
                             if modified_phi_k >= 0.0:
                                 mod_types.append(
@@ -2199,7 +2233,12 @@ def mutate(
         logger.debug(f"LJ scaling factor: {lambda_value_electrostatic}")
         logger.debug(f"VDW scaling factor: {lambda_value_vdw}")
 
-        offset = min([a.idx for a in psf.view[f":{self.tlc.upper()}"].atoms])
+        try:
+            offset = min([a.idx for a in psf.view[f":{self.tlc.upper()}"].atoms])
+        ### This give a ValueErrror for point mutation, where a resid is specified
+        ### but here we have only one ligand or the residue, which should be mutated left
+        except ValueError:
+            offset = min([a.idx for a in psf.view[f":1"].atoms])
 
         if lambda_value_electrostatic < 1.0:
             self._mutate_charge(psf, lambda_value_electrostatic, offset)
@@ -2284,7 +2323,7 @@ def _compensate_charge(
             [atom.charge for atom in psf.view[f":{self.tlc.upper()}"].atoms]
         )
 
-        if not (np.isclose(new_charge, total_charge, rtol=1e-4)):
+        if not (np.isclose(round(new_charge, 3), round(total_charge, 3), rtol=1e-4)):
             raise RuntimeError(
                 f"Charge compensation failed. Introducing non integer total charge: {new_charge}. Target total charge: {total_charge}."
             )

transformato/state.py

@@ -101,6 +101,9 @@ def _write_amber_files(
         if env == "waterbox":
             psf.write_parm(f"{output_file_base}/lig_in_{env}.parm7")
             psf.write_rst7(f"{output_file_base}/lig_in_{env}.rst7")
+        if env == "complex":
+            psf.write_parm(f"{output_file_base}/lig_in_{env}.parm7")
+            psf.write_rst7(f"{output_file_base}/lig_in_{env}.rst7")
         elif env == "vacuum":
             psf[f":{tlc}"].write_parm(f"{output_file_base}/lig_in_{env}.parm7")
             psf[f":{tlc}"].write_rst7(f"{output_file_base}/lig_in_{env}.rst7")
@@ -381,7 +384,7 @@ def _copy_charmm_files(self, intermediate_state_file_path: str):
 
         # copy diverse set of helper files for CHARMM
         for env in self.system.envs:
-            if env != "vacuum":
+            if env != "vacuum" and self.system.ff.lower() != "amber":
                 FILES = [
                     "crystal_image.str",
                     "step3_pbcsetup.str",
@@ -666,12 +669,13 @@ def _copy_files(self, intermediate_state_file_path: str):
 
         basedir = self.system.charmm_gui_base
 
-        try:
-            self._copy_ligand_specific_top_and_par(
-                basedir, intermediate_state_file_path
-            )
-        except:
-            self._copy_ligand_specific_str(basedir, intermediate_state_file_path)
+        if self.system.ff.lower() == "charmm":
+            try:
+                self._copy_ligand_specific_top_and_par(
+                    basedir, intermediate_state_file_path
+                )
+            except:
+                self._copy_ligand_specific_str(basedir, intermediate_state_file_path)
 
         # copy crd file
         try:

transformato/system.py

@@ -47,18 +47,33 @@ def __init__(self, configuration: dict, structure: str):
         # running a binding-free energy calculation?
         if configuration["simulation"]["free-energy-type"] == "rbfe":
             self.envs = set(["complex", "waterbox"])
-            for env in self.envs:
-                parameter = self._read_parameters(env)
-                # set up system
-                self.psfs[env] = self._initialize_system(configuration, env)
-                # load parameters
-                self.psfs[env].load_parameters(parameter)
-                # get offset
-                self.offset[env] = (
-                    self._determine_offset_and_set_possible_dummy_properties(
-                        self.psfs[env]
+
+            if self.ff == "charmm":
+                for env in self.envs:
+                    parameter = self._read_parameters(env)
+                    # set up system
+                    self.psfs[env] = self._initialize_system(configuration, env)
+                    # load parameters
+                    self.psfs[env].load_parameters(parameter)
+                    # get offset
+                    self.offset[env] = (
+                        self._determine_offset_and_set_possible_dummy_properties(
+                            self.psfs[env]
+                        )
+                    )
+            elif self.ff == "amber":
+                for env in self.envs:
+                    self.psfs[env] = pm.load_file(
+                        f"{self.charmm_gui_base}/{env}/openmm/step3_input.parm7"
+                    )
+                    self.psfs[env].load_rst7(
+                        f"{self.charmm_gui_base}/{env}/openmm/step3_input.rst7"
+                    )
+                    self.offset[env] = (
+                        self._determine_offset_and_set_possible_dummy_properties(
+                            self.psfs[env]
+                        )
                     )
-                )
 
             # generate rdkit mol object of small molecule
             self.mol: Chem.Mol = self._generate_rdkit_mol(
@@ -378,6 +393,40 @@ def _return_small_molecule(self, env: str) -> Chem.rdchem.Mol:
                 "Could not load small molecule sdf file in {charmm_gui_env}. Aborting."
             )
 
+    def _create_sdf_from_pdb(self, psf: pm.topologyobjects):
+        """
+        This function creates a sdf file from the file provided in openmm/step3_input.pdb
+        and reads in this sdf file. This should now have the same order as the psf topology object
+        and is from now on used as the rdkit mol object.
+
+        Parameters
+        ----------
+        psf: pm.topologyobjects
+            a topology provided by parmed (either CHARMM type or AMBER)
+        Returns
+        ----------
+        mol : rdkit.Chem.rdchem.Mol
+            mol object
+        """
+
+        file_path = f"{self.charmm_gui_base}/waterbox/openmm/"
+        psf.save(f"{file_path}/step3_input_reordered.pdb", overwrite=True)
+
+        from openbabel import openbabel
+
+        obConversion = openbabel.OBConversion()
+        obConversion.SetInAndOutFormats("pdb", "sdf")
+        mol = openbabel.OBMol()
+        obConversion.ReadFile(mol, f"{file_path}/step3_input_reordered.pdb")
+        obConversion.WriteFile(mol, f"{file_path}/step3_input_reduced.sdf")
+
+        suppl = Chem.SDMolSupplier(
+            f"{file_path}/step3_input_reduced.sdf", removeHs=False
+        )
+        mol = next(suppl)
+
+        return mol
+
     def _generate_rdkit_mol(
         self, env: str, psf: pm.charmm.CharmmPsfFile
     ) -> Chem.rdchem.Mol:
@@ -393,53 +442,72 @@ def _generate_rdkit_mol(
         mol: rdkit.Chem.rdchem.Mol
         """
 
-        mol = self._return_small_molecule(env)
-        (
-            atom_idx_to_atom_name,
-            _,
-            atom_name_to_atom_type,
-            atom_idx_to_atom_partial_charge,
-        ) = self.generate_atom_tables_from_psf(psf)
+        try:
+            mol = self._return_small_molecule(env)
 
-        for atom in mol.GetAtoms():
-            atom.SetProp("atom_name", atom_idx_to_atom_name[atom.GetIdx()])
-            atom.SetProp(
-                "atom_type",
-                atom_name_to_atom_type[atom_idx_to_atom_name[atom.GetIdx()]],
-            )
-            atom.SetProp("atom_index", str(atom.GetIdx()))
-            atom.SetProp(
-                "atom_charge", str(atom_idx_to_atom_partial_charge[atom.GetIdx()])
-            )
+        except FileNotFoundError:
+
+            mol = self._create_sdf_from_pdb(psf)
+
+        mol = self.generate_atom_tables_from_psf(psf, mol)
 
         # check if psf and sdf have same indeces
+        try:
+            for a in mol.GetAtoms():
+                if str(psf[a.GetIdx()].element_name) == str(a.GetSymbol()):
+                    pass
+                else:
+                    raise RuntimeError(
+                        "First PSF to mol conversion did not work, trying it with the pdb file"
+                    )
+        except RuntimeError:
+            # It's possible, that there is a sdf file in the TCL folder but the order does not match
+            # in thuis case, we try to create a matching sdf file by using the pdb file in the openmm folder
+            logger.info(
+                f"It seems like the sdf file provided in the {self.tlc} does not have the correct order, trying it again by using the pdb file"
+            )
+            mol = self._create_sdf_from_pdb(psf)
+            mol = self.generate_atom_tables_from_psf(psf, mol)
+
+        # final check if psf and sdf have same indeces
         for a in mol.GetAtoms():
             if str(psf[a.GetIdx()].element_name) == str(a.GetSymbol()):
                 pass
             else:
-                raise RuntimeError("PSF to mol conversion did not work! Aborting.")
+                raise RuntimeError(
+                    "PSF to mol conversion did not work! Remove the TLC.sdf file in the tlc folder and try again!"
+                )
 
         return mol
 
     def generate_atom_tables_from_psf(
-        self, psf: pm.charmm.CharmmPsfFile
-    ) -> Tuple[dict, dict, dict, dict]:
+        self, psf: pm.charmm.CharmmPsfFile, mol: Chem.rdchem.Mol
+    ) -> Chem.rdchem.Mol:
         """
         Generate mapping dictionaries for a molecule in a psf.
         Parameters
         ----------
-        psf: pm.charmm.CharmmPsfFile
+        psf: pm.charmm.CharmmPsfFile or parmed.amber._amberparm.AmberParm
         Returns
         ----------
-        dict's
+        mol: Chem.rdchem.Mol
+            mol object with atom name and type properties
         """
 
         atom_idx_to_atom_name = dict()
         atom_name_to_atom_idx = dict()
         atom_name_to_atom_type = dict()
         atom_idx_to_atom_partial_charge = dict()
 
-        for atom in psf.view[f":{self.tlc}"].atoms:
+        ## We need this for point mutations, because if we give a resid, the mol here
+        ## consists only of on residue which resid is always 1
+        try:
+            int(self.tlc)
+            tlc = "1"
+        except ValueError:
+            tlc = self.tlc
+
+        for atom in psf.view[f":{tlc}"].atoms:
             atom_name = atom.name
             atom_index = atom.idx
             atom_type = atom.type
@@ -450,9 +518,15 @@ def generate_atom_tables_from_psf(
             atom_name_to_atom_type[atom_name] = atom_type
             atom_idx_to_atom_partial_charge[atom_index] = atom_charge
 
-        return (
-            atom_idx_to_atom_name,
-            atom_name_to_atom_idx,
-            atom_name_to_atom_type,
-            atom_idx_to_atom_partial_charge,
-        )
+        for atom in mol.GetAtoms():
+            atom.SetProp("atom_name", atom_idx_to_atom_name[atom.GetIdx()])
+            atom.SetProp(
+                "atom_type",
+                atom_name_to_atom_type[atom_idx_to_atom_name[atom.GetIdx()]],
+            )
+            atom.SetProp("atom_index", str(atom.GetIdx()))
+            atom.SetProp(
+                "atom_charge", str(atom_idx_to_atom_partial_charge[atom.GetIdx()])
+            )
+
+        return mol