Support for neutral point mutation

transformato/mutate.py

@@ -358,13 +358,13 @@ def __init__(
             self.dummy_region_cc2: DummyRegion
 
             self.asfe: bool = False
-            self._check_cgenff_versions()
-
-        except:
+            # self._check_cgenff_versions()
 
+        except AttributeError:
             logger.info(
                 "Only information about one structure, assume an ASFE simulation is requested"
             )
+
             mol1_name: str = "m1"
             self.system: dict = {"system1": s1}
             self.mols: dict = {mol1_name: s1.mol}
@@ -602,10 +602,9 @@ def _check_for_lp(
         ]
         lp_dict_dummy_region = defaultdict(list)
         lp_dict_common_core = defaultdict(list)
-        print(f"die version {pm.__version__}")
+
         for atom in psf.view[f":{tlc}"].atoms:
             if atom.name.find("LP") == False:
-                print(f"die Atome {atom}")
                 if atom.frame_type.atom1.idx in flat_ordered_connected_dummy_regions:
                     lp_dict_dummy_region[atom.frame_type.atom1.idx].append(atom.idx)
 
@@ -1222,8 +1221,8 @@ def _transform_common_core(self) -> list:
             logger.warning(f"Charge parameters mutation: {charge_mutation}.")
 
             t = CommonCoreTransformation(
-                self.get_common_core_idx_mol1() + self.dummy_region_cc1.lj_default,
-                self.get_common_core_idx_mol2() + self.dummy_region_cc2.lj_default,
+                self.get_common_core_idx_mol1(),
+                self.get_common_core_idx_mol2(),
                 self.psfs["m1"],
                 self.psfs["m2"],
                 self.s1_tlc,
@@ -1458,8 +1457,55 @@ def _get_atom_mapping(self) -> dict:
             cc_names_struc1.append(ligand1_atom.name)
             cc_names_struc2.append(ligand2_atom.name)
 
-        print(f"CC Struc1: {cc_names_struc1}")
-        print(f"CC Struc2: {cc_names_struc2}")
+        logger.info(f"CC Struc1: {cc_names_struc1}")
+        logger.info(f"CC Struc2: {cc_names_struc2}")
+
+        ### DANGER ONLY FOR ONE MUTATION NEEEDS TO BE FIXED #####
+        # match_atom_names_cc1_to_cc2 = {
+        #     "O5'": "O5'",
+        #     "H5T": "H5T",
+        #     "C5'": "C5'",
+        #     "H5''": "H5''",
+        #     "H5'": "H5'",
+        #     "C4'": "C4'",
+        #     "H4'": "H4'",
+        #     "O4'": "O4'",
+        #     "C1'": "C1'",
+        #     "C2'": "C2'",
+        #     "C3'": "C3'",
+        #     "H3'": "H3'",
+        #     "O3'": "O3'",
+        #     "P": "P",
+        #     "O2P": "O2P",
+        #     "O1P": "O1P",
+        #     "H1'": "H1'",
+        #     "H2''": "H2''",
+        #     "O2'": "O2'",
+        #     "H2'": "H2'",
+        #     "N1": "N1",
+        #     "C6": "C6",
+        #     "H6": "H6",
+        #     "C5": "C5",
+        #     "H5": "H5",
+        #     "C4": "C4",
+        #     "N3": "N3",
+        #     "C2": "C2",
+        #     "O2": "O2",
+        #     "N4": "N4",
+        #     "H42": "H42",
+        #     "H41": "H41",
+        #     "N9": "N9",
+        #     "N7": "N7",
+        #     "C8": "C8",
+        #     "H8": "H8",
+        #     "N2": "N2",
+        #     "H21": "H21",
+        #     "H22": "H22",
+        #     "H1": "H1",
+        #     "O6": "O6",
+        #     "H3T": "H3T",
+        # }
+
         return match_atom_names_cc1_to_cc2
 
     def _mutate_charges(self, psf: pm.charmm.CharmmPsfFile, scale: float):
@@ -1472,7 +1518,15 @@ def _mutate_charges(self, psf: pm.charmm.CharmmPsfFile, scale: float):
 
             # compare to charge compenstated psf 2
             for ligand2_atom in self.charge_compensated_ligand2_psf:
-                if self.atom_names_mapping[ligand1_atom.name] == ligand2_atom.name:
+                # Assure that only atoms from the same resdiue are compared, and only atoms belonging to the same chain!
+                if (
+                    self.atom_names_mapping[ligand1_atom.name] == ligand2_atom.name
+                    and ligand1_atom.residue.name == ligand2_atom.residue.name
+                    and ligand1_atom.type == ligand2_atom.type
+                    and len(ligand1_atom.residue.atoms)
+                    == len(ligand2_atom.residue.atoms)
+                    and ligand1_atom.residue.chain == ligand2_atom.residue.chain
+                ):
                     found = True
                     # are the atoms different?
                     logger.debug(f"Modifying atom: {ligand1_atom}")
@@ -1488,7 +1542,34 @@ def _mutate_charges(self, psf: pm.charmm.CharmmPsfFile, scale: float):
                     ligand1_atom.charge = modified_charge
 
             if not found:
-                raise RuntimeError("No corresponding atom in cc2 found")
+                try:
+                    for ligand2_atom in self.charge_compensated_ligand2_psf:
+                        # make sure resdiue PSU which is like CYT are nevertheless found
+                        if (
+                            self.atom_names_mapping[ligand1_atom.name]
+                            == ligand2_atom.name
+                            and len(ligand1_atom.residue.atoms)
+                            == len(ligand2_atom.residue.atoms)
+                            and ligand1_atom.residue.chain == ligand2_atom.residue.chain
+                        ):
+                            found = True
+                            # are the atoms different?
+                            logger.debug(f"Modifying atom: {ligand1_atom}")
+                            logger.debug(f"Template atom: {ligand2_atom}")
+
+                            # scale epsilon
+                            modified_charge = (
+                                scale * ligand1_atom.charge
+                                + (1 - scale) * ligand2_atom.charge
+                            )
+                            logger.debug(
+                                f"Current charge: {ligand1_atom.charge}; target charge: {ligand2_atom.charge}; modified charge: {modified_charge}"
+                            )
+                            ligand1_atom.charge = modified_charge
+                except:
+                    raise RuntimeError(
+                        f"No corresponding atom for {ligand1_atom} in cc2 found"
+                    )
 
     def _mutate_atoms(self, psf: pm.charmm.CharmmPsfFile, lambda_value: float):
         """
@@ -1516,8 +1597,10 @@ def _mutate_atoms(self, psf: pm.charmm.CharmmPsfFile, lambda_value: float):
                 # is there a match up?
                 if self.atom_names_mapping[ligand1_atom.name] == ligand2_atom.name:
                     found = True
-                    # are the atoms different?
-                    if ligand1_atom.type != ligand2_atom.type:
+                    # are the atoms different? and assure that only atomtypes in the same residue are compared
+                    if ligand1_atom.type != ligand2_atom.type and len(
+                        ligand1_atom.residue.atoms
+                    ) == len(ligand2_atom.residue.atoms):
                         if "DDX" in ligand1_atom.type:
                             logger.warning(
                                 "This is the terminal LJ atom. If everything went correct, this does not have to change atom types."
@@ -1845,6 +1928,9 @@ def _modify_type_in_cc(atom: pm.Atom, psf: pm.charmm.CharmmPsfFile):
         if hasattr(atom, "initial_type"):
             # only change parameters
             pass
+        elif atom.residue.chain == "RNAB":
+            # only change parameters
+            pass
         else:
             logger.info(f"Setting RRR atomtype for atom: {atom}.")
             atom.initial_type = atom.type
@@ -2045,6 +2131,9 @@ def _modify_type(atom, psf, atom_type_suffix: str):
         if hasattr(atom, "initial_type"):
             # only change parameters
             pass
+        elif atom.residue.chain == "RNAB":
+            # only change parameters
+            pass
         else:
             atom.initial_type = atom.type
             if atom_type_suffix == "DDD":

transformato/state.py

@@ -132,7 +132,7 @@ def write_state(
                     )
 
                 else:
-                    # mutation_type.print_details()
+
                     logger.debug(f"Lambda electrostatics: {lambda_value_electrostatic}")
                     logger.debug(f"Lambda vdw: {lambda_value_vdw}")
 
@@ -628,12 +628,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 os.path.exists(f"{basedir}/waterbox/{self.system.tlc.lower()}/"):
+            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
         self._copy_crd_file((intermediate_state_file_path))
@@ -1031,35 +1032,30 @@ def _init_base_dir(self):
 
     def _write_toppar_str(self, output_file_base):
 
-        toppar_format = f"""
-../../toppar/top_all36_prot.rtf
-../../toppar/par_all36m_prot.prm
-../../toppar/top_all36_na.rtf
-../../toppar/par_all36_na.prm
-../../toppar/top_all36_carb.rtf
-../../toppar/par_all36_carb.prm
-../../toppar/top_all36_lipid.rtf
-../../toppar/par_all36_lipid.prm
-../../toppar/top_all36_cgenff.rtf
-../../toppar/par_all36_cgenff.prm
-../../toppar/toppar_water_ions.str
-../../toppar/toppar_dum_noble_gases.str
-../../toppar/toppar_all36_prot_c36m_d_aminoacids.str
-../../toppar/toppar_all36_prot_fluoro_alkanes.str
-../../toppar/toppar_all36_prot_heme.str
-../../toppar/toppar_all36_prot_na_combined.str
-../../toppar/toppar_all36_prot_retinol.str
-../../toppar/toppar_all36_na_nad_ppi.str
-../../toppar/toppar_all36_lipid_bacterial.str
-../../toppar/toppar_all36_lipid_cardiolipin.str
-../../toppar/toppar_all36_lipid_cholesterol.str
-../../toppar/toppar_all36_lipid_inositol.str
-../../toppar/toppar_all36_lipid_lps.str
-../../toppar/toppar_all36_lipid_miscellaneous.str
-../../toppar/toppar_all36_lipid_model.str
-../../toppar/toppar_all36_lipid_prot.str
-../../toppar/toppar_all36_lipid_sphingo.str
-"""
+        from transformato.system import parameter_files
+
+        base = "../../toppar"
+        toppar_format = ""
+        # Only necessary because the variable self.configuration["system"]["structure2"]["tlc"] gives a keyError for asfe
+        # since its not defined there
+        try:
+            if self.configuration["system"]["structure2"]["tlc"].lower():
+                struct2 = self.configuration["system"]["structure2"]["tlc"].lower()
+            else:
+                struct2 = "emptyString"
+
+        except KeyError:
+            struct2 = "emptyString"
+
+        for i in parameter_files:
+            # the parameterfiles for customized ligands are in another directory and are considered below
+            if self.configuration["system"]["structure1"]["tlc"].lower() in i:
+                pass
+            if struct2 in i:
+                pass
+            else:
+                toppar_format += f"{base}/{i.split('/')[-1]} \n"
+
         if os.path.isfile(
             f"{self.system.charmm_gui_base}/waterbox/{self.system.tlc.lower()}/{self.system.tlc.lower()}_g.rtf"
         ):
@@ -1068,10 +1064,16 @@ def _write_toppar_str(self, output_file_base):
 dummy_atom_definitions.rtf
 dummy_parameters.prm
 """
-        else:
+        elif os.path.isfile(
+            f"{self.system.charmm_gui_base}/waterbox/{self.system.tlc.lower()}/{self.system.tlc.lower()}.str"
+        ):
             toppar_format += f"""{self.system.tlc.lower()}.str
 dummy_atom_definitions.rtf
 dummy_parameters.prm
+"""
+        else:
+            toppar_format += f"""dummy_atom_definitions.rtf
+dummy_parameters.prm
 """
 
         f = open(f"{output_file_base}/toppar.str", "w+")
@@ -1100,7 +1102,7 @@ def _write_psf(psf, output_file_base: str, env: str):
         psf.write_psf(string_object)
         # read in psf and correct some aspects of the file not suitable for CHARMM
         corrected_psf = psf_correction(string_object)
-        with open(f"{output_file_base}/lig_in_{env}_corr.psf", "w+") as f:
+        with open(f"{output_file_base}/lig_in_{env}.psf", "w+") as f:
             f.write(corrected_psf)
         # write pdb
         psf.write_pdb(f"{output_file_base}/lig_in_{env}.pdb")