From e4c08bf130e67006cf28cfd4bb19cb9d100dfce2 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Tue, 17 Dec 2024 22:54:44 +0100 Subject: [PATCH] reverted to original sigma values hydrogen with only h-h h-o interactions updated dihedrals --- .gitignore | 1 + multi-ego-basic.ff/ffbonded.itp | 15 +- src/multiego/ensemble.py | 33 +++-- src/multiego/resources/type_definitions.py | 158 ++++++++++----------- src/multiego/topology.py | 6 +- tools/make_mat/make_mat.py | 26 ++-- 6 files changed, 126 insertions(+), 113 deletions(-) diff --git a/.gitignore b/.gitignore index a55b214f..5863605d 100644 --- a/.gitignore +++ b/.gitignore @@ -137,3 +137,4 @@ outputs/* analysis/* .DS_Store *.ff +!multi-ego-basic.ff diff --git a/multi-ego-basic.ff/ffbonded.itp b/multi-ego-basic.ff/ffbonded.itp index 21d8a28a..ec857a51 100644 --- a/multi-ego-basic.ff/ffbonded.itp +++ b/multi-ego-basic.ff/ffbonded.itp @@ -374,8 +374,8 @@ ; -CHn-NT- 0.9 ; ; OPT November 23, 22 -;#define gd_42 170.0000 2.7000 2 -#define gd_42 -170.0000 2.7000 2 +;#define gd_42 -170.0000 2.7000 2 +#define gd_42 -170.0000 4.0000 2 ; Backbone dihedral angle -N-CA-C-N- ALA/SER/THR ; ;#define gd_43 15.0000 2.500 3 @@ -383,7 +383,7 @@ ; Backbone dihedral angle -C-N-CA-C- ALA/SER/THR ; ;#define gd_44 180.000 4.00 1 -#define gd_44 160.000 4.00 1 +#define gd_44 165.000 4.00 1 ; Backbone dihedral angle -C-N-CA-C- ALA/SER/THR ; ;#define gd_45 -105.00 0.2200 1 @@ -400,7 +400,7 @@ #define gd_44b -110.000 3.90 1 ; Backbone dihedral angle -C-N-CA-C- BULKYAA ; -#define gd_45b 180.000 0.05 1 +#define gd_45b 180.000 0.80 1 ; Backbone dihedral angle -N-CA-C-N- BULKYAA ; ; OPT April, 3 2024 @@ -410,7 +410,7 @@ #define gd_49 180 11 1 ; Backbone dihedral angle -C-N-CA-C GLY ; -#define gd_51 -180 1.7 1 +#define gd_51 -180 2.5 1 ; Backbone dihedral angle -N-CA-C-N GLY ; ; OPT November, 23 2022 @@ -420,10 +420,7 @@ #define gd_53 0.0000 4.000 4 ; Backbone dihedral angle -C-N-CA-C PRO ; -;#define gd_54 -110.4281 0.0000 1 -; Backbone dihedral angle -C-N-CA-C PRO -; -#define gd_55 -60.0000 1.000 1 +#define gd_55 -60.0000 0.400 1 ; Backbone dihedral angle -N-CA-C-N PRO ; #define gd_56 180.000 31.0 2 diff --git a/src/multiego/ensemble.py b/src/multiego/ensemble.py index 3f174b8e..3a2c1c26 100644 --- a/src/multiego/ensemble.py +++ b/src/multiego/ensemble.py @@ -782,7 +782,7 @@ def init_LJ_datasets(meGO_ensemble, matrices, pairs14, exclusion_bonds14, args): pairwise_c12 = np.where( oxygen_mask, - (11.4 / 0.9**12) * np.sqrt(type_ai_mapped.map(type_to_c12) * type_aj_mapped.map(type_to_c12)), + 3e-6, np.sqrt( train_dataset["ai"].map(meGO_ensemble["sbtype_c12_dict"]) * train_dataset["aj"].map(meGO_ensemble["sbtype_c12_dict"]) @@ -792,7 +792,7 @@ def init_LJ_datasets(meGO_ensemble, matrices, pairs14, exclusion_bonds14, args): pairwise_mg_sigma = np.where( oxygen_mask, - ((11.4 / 0.9**12) * np.sqrt(type_ai_mapped.map(type_to_c12) * type_aj_mapped.map(type_to_c12))) ** (1 / 12), + (3e-6) ** (1 / 12), ( train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c12_dict"]) * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c12_dict"]) @@ -807,7 +807,7 @@ def init_LJ_datasets(meGO_ensemble, matrices, pairs14, exclusion_bonds14, args): pairwise_mg_epsilon = np.where( oxygen_mask, - -(11.4 / 0.9**12) * np.sqrt(type_ai_mapped.map(type_to_c12) * type_aj_mapped.map(type_to_c12)), + -3e-6, ( train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c6_dict"]) * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c6_dict"]) @@ -834,16 +834,27 @@ def generate_OO_LJ(meGO_ensemble): sbtype for sbtype, atomtype in meGO_ensemble["sbtype_type_dict"].items() if atomtype == "O" or atomtype == "OM" ] + H_H_sbtype = [sbtype for sbtype, atomtype in meGO_ensemble["sbtype_type_dict"].items() if atomtype == "H"] + + half_matrix = list([(a, b) for a, b in itertools.product(O_OM_sbtype, H_H_sbtype)]) + # Generate all possible combinations combinations = list(itertools.combinations_with_replacement(O_OM_sbtype, 2)) - # Create a DataFrame from the combinations - rc_LJ = pd.DataFrame(combinations, columns=["ai", "aj"]) + OO_LJ = pd.DataFrame(combinations, columns=["ai", "aj"]) + OO_LJ["c12"] = 3e-6 + OO_LJ["c6"] = 0.0 + # Generate all possible combinations + combinations = list(itertools.combinations_with_replacement(H_H_sbtype, 2)) + # Create a DataFrame from the combinations + HH_LJ = pd.DataFrame(combinations, columns=["ai", "aj"]) + HH_LJ["c12"] = 9.14859e-10 + HH_LJ["c6"] = 0.0 + HO_LJ = pd.DataFrame(half_matrix, columns=["ai", "aj"]) + HO_LJ["c12"] = 1.60608e-09 * 0.115 # mg_eps + HO_LJ["c6"] = 8.01519e-05 * 0.115 # mg_eps + rc_LJ = pd.concat([OO_LJ, HO_LJ, HH_LJ], axis=0) rc_LJ["type"] = 1 - rc_LJ["c6"] = 0.0 - rc_LJ["c12"] = (11.4 / 0.9**12) * np.sqrt( - rc_LJ["ai"].map(meGO_ensemble["sbtype_c12_dict"]) * rc_LJ["aj"].map(meGO_ensemble["sbtype_c12_dict"]) - ) rc_LJ["same_chain"] = False rc_LJ["source"] = "mg" rc_LJ["rep"] = rc_LJ["c12"] @@ -1540,9 +1551,7 @@ def make_pairs_exclusion_topology(meGO_ensemble, meGO_LJ_14, args): | (df["aj"].map(meGO_ensemble["sbtype_type_dict"]) == "O") ), "c12", - ] *= ( - 11.4 / 0.9**12 - ) + ] = 3e-6 df["same_chain"] = True df["probability"] = 1.0 df["rc_probability"] = 1.0 diff --git a/src/multiego/resources/type_definitions.py b/src/multiego/resources/type_definitions.py index a8b9f10b..5acf1000 100644 --- a/src/multiego/resources/type_definitions.py +++ b/src/multiego/resources/type_definitions.py @@ -3,11 +3,9 @@ import sys -mg_eps = 0.110 -mg_eps_ch2 = 0.100 -mg_eps_ch1 = 0.090 -sig6 = 0.9**6 -sig12 = 0.9**12 +mg_eps = 0.115 +mg_eps_ch2 = 0.105 +mg_eps_ch1 = 0.095 # Dataframe with GROMOS atom types and associated parameters gromos_atp = pd.DataFrame( { @@ -39,84 +37,82 @@ ], "at.num": [8, 8, 8, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 16, 6, 15, 8, 6, 1, 20], "rc_c12": [ - 2.631580e-07 * sig12, # "O", - 1.724403e-07 * sig12, # "OM", - 5.018430e-07 * sig12, # "OA", - 8.752940e-07 * sig12, # "N", - 2.596154e-06 * sig12, # "NT", - 8.752940e-07 * sig12, # "NL", - 1.506347e-06 * sig12, # "NR", - 8.752940e-07 * sig12, # "NZ", - 8.752940e-07 * sig12, # "NE", - 2.598570e-06 * sig12, # "C", - 2.598570e-06 * sig12, # "CH" - 6.555574e-05 * sig12, # "CH1" - 6.555574e-05 * sig12, # "CAH" - 6.555574e-05 * sig12, # "CH1a" - 1.543890e-05 * sig12, # "CH2" - 8.595562e-06 * sig12, # "CH3" - 1.193966e-05 * sig12, # "CH2r" - 2.724050e-06 * sig12, # "S", - 8.736473e-06 * sig12, # "CH3p" - 3.893600e-06 * sig12, # "P", - 3.558824e-07 * sig12, # "OE", - 6.298560e-06 * sig12, # "CR1", - #1.239247e-09 * sig12, # "H", #1.2392471565152519e-09 * sig12 = 3.5e-10 1KT LJ match - 3.239247e-09 * sig12, # "H", #3.239247e-09 * sig12 = 9.14859e-10 H-H ALA distance match - 2.659360e-07 * sig12, # "C0", + 2.631580e-07, # "O", + 1.724403e-07, # "OM", + 5.018430e-07, # "OA", + 8.752940e-07, # "N", + 2.596154e-06, # "NT", + 8.752940e-07, # "NL", + 1.506347e-06, # "NR", + 8.752940e-07, # "NZ", + 8.752940e-07, # "NE", + 2.598570e-06, # "C", + 2.598570e-06, # "CH" + 6.555574e-05, # "CH1" + 6.555574e-05, # "CAH" + 6.555574e-05, # "CH1a" + 1.543890e-05, # "CH2" + 8.595562e-06, # "CH3" + 1.193966e-05, # "CH2r" + 2.724050e-06, # "S", + 8.736473e-06, # "CH3p" + 3.893600e-06, # "P", + 3.558824e-07, # "OE", + 6.298560e-06, # "CR1", + 9.148590e-10, # "H", + 2.659360e-07, # "C0", ], "mg_c12": [ - 1e-06 / 1.27911 * mg_eps * sig12, # "O", - 7.4149321e-07 / 1.72504 * mg_eps * sig12, # "OM", - 1.505529e-06 / 0.84961 * mg_eps * sig12, # "OA", - 2.319529e-06 / 0.63980 * mg_eps * sig12, # "N", - 5.0625e-06 / 0.29314 * mg_eps * sig12, # "NT", - 2.319529e-06 / 0.63980 * mg_eps * sig12, # "NL", - 3.389281e-06 / 0.43786 * mg_eps * sig12, # "NR", - 2.319529e-06 / 0.63980 * mg_eps * sig12, # "NZ", - 2.319529e-06 / 0.63980 * mg_eps * sig12, # "NE", - 4.937284e-06 / 0.27741 * mg_eps * sig12, # "C", - 4.937284e-06 / 0.27741 * mg_eps * sig12, # "CH" - 9.70225e-05 / 0.09489 * mg_eps_ch1 * sig12, # "CH1" - 9.70225e-05 / 0.09489 * mg_eps_ch1 * sig12, # "CAH" - 9.70225e-05 / 0.09489 * mg_eps_ch1 * sig12, # "CH1a" - 3.3965584e-05 / 0.4105 * mg_eps_ch2 * sig12, # "CH2" - 2.6646244e-05 / 0.8671 * mg_eps * sig12, # "CH3" - 2.8058209e-05 / 0.4792 * mg_eps_ch2 * sig12, # "CH2r" - 1.3075456e-05 / 1.90587 * mg_eps * sig12, # "S", - 2.6646244e-05 / 0.86715 * mg_eps * sig12, # "CH3p" - 2.2193521e-05 / 2.44674 * mg_eps * sig12, # "P", - 1.21e-06 / 1.05711 * mg_eps * sig12, # "OE", - 1.5116544e-05 / 0.50266 * mg_eps * sig12, # "CR1", - #1.239247e-09 * sig12, # "H", #1.2392471565152519e-09 * sig12 = 3.5e-10 1KT LJ match - 3.239247e-09 * sig12, # "H", #3.239247e-09 * sig12 = 9.14859e-10 H-H ALA distance match - 0 * mg_eps * sig12, # "C0", + 1.0000000e-06 / 1.27911 * mg_eps, # "O", + 7.4149321e-07 / 1.72504 * mg_eps, # "OM", + 1.5055290e-06 / 0.84961 * mg_eps, # "OA", + 2.3195290e-06 / 0.63980 * mg_eps, # "N", + 5.0625000e-06 / 0.29314 * mg_eps, # "NT", + 2.3195290e-06 / 0.63980 * mg_eps, # "NL", + 3.3892810e-06 / 0.43786 * mg_eps, # "NR", + 2.3195290e-06 / 0.63980 * mg_eps, # "NZ", + 2.3195290e-06 / 0.63980 * mg_eps, # "NE", + 4.9372840e-06 / 0.27741 * mg_eps, # "C", + 4.9372840e-06 / 0.27741 * mg_eps, # "CH" + 9.7022500e-05 / 0.09489 * mg_eps_ch1, # "CH1" + 9.7022500e-05 / 0.09489 * mg_eps_ch1, # "CAH" + 9.7022500e-05 / 0.09489 * mg_eps_ch1, # "CH1a" + 3.3965584e-05 / 0.41050 * mg_eps_ch2, # "CH2" + 2.6646244e-05 / 0.86710 * mg_eps, # "CH3" + 2.8058209e-05 / 0.47920 * mg_eps_ch2, # "CH2r" + 1.3075456e-05 / 1.90587 * mg_eps, # "S", + 2.6646244e-05 / 0.86715 * mg_eps, # "CH3p" + 2.2193521e-05 / 2.44674 * mg_eps, # "P", + 1.2100000e-06 / 1.05711 * mg_eps, # "OE", + 1.5116544e-05 / 0.50266 * mg_eps, # "CR1", + 0.0000000e-00 / 1.00000 * mg_eps, # "H", + 0.0000000e-00 / 1.00000 * mg_eps, # "C0", ], "mg_c6": [ - 0.0022619536 / 1.27911 * mg_eps * sig6, # "O", - 0.0022619536 / 1.72504 * mg_eps * sig6, # "OM", - 0.0022619536 / 0.84961 * mg_eps * sig6, # "OA", - 0.0024364096 / 0.63980 * mg_eps * sig6, # "N", - 0.0024364096 / 0.29314 * mg_eps * sig6, # "NT", - 0.0024364096 / 0.63980 * mg_eps * sig6, # "NL", - 0.0024364096 / 0.43786 * mg_eps * sig6, # "NR", - 0.0024364096 / 0.63980 * mg_eps * sig6, # "NZ", - 0.0024364096 / 0.63980 * mg_eps * sig6, # "NE", - 0.0023406244 / 0.27741 * mg_eps * sig6, # "C", - 0.0023406244 / 0.27741 * mg_eps * sig6, # "CH" - 0.00606841 / 0.09489 * mg_eps_ch1 * sig6, # "CH1" - 0.00606841 / 0.09489 * mg_eps_ch1 * sig6, # "CAH" - 0.00606841 / 0.09489 * mg_eps_ch1 * sig6, # "CH1a" - 0.0074684164 / 0.41054 * mg_eps_ch2 * sig6, # "CH2" - 0.0096138025 / 0.86715 * mg_eps * sig6, # "CH3" - 0.0073342096 / 0.47928 * mg_eps_ch2 * sig6, # "CH2r" - 0.0099840064 / 1.90587 * mg_eps * sig6, # "S", - 0.0096138025 / 0.86715 * mg_eps * sig6, # "CH3p" - 0.01473796 / 2.44674 * mg_eps * sig6, # "P", - 0.0022619536 / 1.05711 * mg_eps * sig6, # "OE", - 0.0055130625 / 0.50266 * mg_eps * sig6, # "CR1", - 0 * mg_eps * sig6, # "H", # TODO - 0 * mg_eps * sig6, # "C0", + 0.0022619536 / 1.27911 * mg_eps, # "O", + 0.0022619536 / 1.72504 * mg_eps, # "OM", + 0.0022619536 / 0.84961 * mg_eps, # "OA", + 0.0024364096 / 0.63980 * mg_eps, # "N", + 0.0024364096 / 0.29314 * mg_eps, # "NT", + 0.0024364096 / 0.63980 * mg_eps, # "NL", + 0.0024364096 / 0.43786 * mg_eps, # "NR", + 0.0024364096 / 0.63980 * mg_eps, # "NZ", + 0.0024364096 / 0.63980 * mg_eps, # "NE", + 0.0023406244 / 0.27741 * mg_eps, # "C", + 0.0023406244 / 0.27741 * mg_eps, # "CH" + 0.0060684100 / 0.09489 * mg_eps_ch1, # "CH1" + 0.0060684100 / 0.09489 * mg_eps_ch1, # "CAH" + 0.0060684100 / 0.09489 * mg_eps_ch1, # "CH1a" + 0.0074684164 / 0.41054 * mg_eps_ch2, # "CH2" + 0.0096138025 / 0.86715 * mg_eps, # "CH3" + 0.0073342096 / 0.47928 * mg_eps_ch2, # "CH2r" + 0.0099840064 / 1.90587 * mg_eps, # "S", + 0.0096138025 / 0.86715 * mg_eps, # "CH3p" + 0.0147379600 / 2.44674 * mg_eps, # "P", + 0.0022619536 / 1.05711 * mg_eps, # "OE", + 0.0055130625 / 0.50266 * mg_eps, # "CR1", + 0.0000000000 / 1.00000 * mg_eps, # "H", # TODO + 0.0000000000 / 1.00000 * mg_eps, # "C0", ], } ) @@ -164,8 +160,8 @@ def lj14_generator(df): atom_type_combinations = [ # Tuple of atom type combinations for LJ14 pairs ("backbone_carbonyl", "sidechain_cb", 0.275, None, 1), - ("backbone_oxygen", "sidechain_cb", 0.1, None, 0), - ("ct_oxygen", "sidechain_cb", 0.1, None, 0), + ("backbone_oxygen", "sidechain_cb", 0.2, None, 0), + ("ct_oxygen", "sidechain_cb", 0.2, None, 0), ("backbone_nitrogen", "sidechain_cb", 0.65, None, -1), ("first_backbone_nitrogen", "backbone_nitrogen", None, 4.0e-6, 1), ("backbone_nitrogen", "backbone_nitrogen", 0.343, None, 1), diff --git a/src/multiego/topology.py b/src/multiego/topology.py index 03d763e1..54fa04f3 100644 --- a/src/multiego/topology.py +++ b/src/multiego/topology.py @@ -443,7 +443,7 @@ def create_pairs_14_dataframe(atomtype1, atomtype2, c6=0.0, shift=0, prefactor=N if constant is not None: pairs_14_c12.append(constant) if prefactor is not None: - pairs_14_c12.append((prefactor / 0.9**12) * np.sqrt(line_atomtype1["c12"] * line_atomtype2["c12"])) + pairs_14_c12.append(prefactor * np.sqrt(line_atomtype1["c12"] * line_atomtype2["c12"])) pairs_14 = pd.DataFrame( columns=[ @@ -515,7 +515,7 @@ def protein_LJ14(reduced_topology): pairs = pd.concat( [ pairs, - create_pairs_14_dataframe(atomtype1=backbone_oxygen, atomtype2=sidechain_cb, prefactor=0.1), + create_pairs_14_dataframe(atomtype1=backbone_oxygen, atomtype2=sidechain_cb, prefactor=0.2), ], axis=0, sort=False, @@ -525,7 +525,7 @@ def protein_LJ14(reduced_topology): pairs = pd.concat( [ pairs, - create_pairs_14_dataframe(atomtype1=ct_oxygen, atomtype2=sidechain_cb, prefactor=0.1), + create_pairs_14_dataframe(atomtype1=ct_oxygen, atomtype2=sidechain_cb, prefactor=0.2), ], axis=0, sort=False, diff --git a/tools/make_mat/make_mat.py b/tools/make_mat/make_mat.py index a6bd0f7d..91f50827 100644 --- a/tools/make_mat/make_mat.py +++ b/tools/make_mat/make_mat.py @@ -571,16 +571,28 @@ def main_routine(mol_i, mol_j, topology_mego, topology_ref, molecules_name, pref # print(args.bkbn_H) # print(np.array([ protein_ref_i[i].name for i in range(len(protein_ref_i.atoms)) if (protein_ref_i[i].element_name != "H" or protein_ref_i[i].name == args.bkbn_H)])) # print(np.array([ protein_ref_i[i].element_name for i in range(len(protein_ref_i.atoms)) ])) - #exit() - protein_ref_indices_i = np.array([i + 1 for i in range(len(protein_ref_i.atoms)) if (protein_ref_i[i].element_name != "H" or protein_ref_i[i].name == args.bkbn_H)]) - protein_ref_indices_j = np.array([i + 1 for i in range(len(protein_ref_j.atoms)) if (protein_ref_i[i].element_name != "H" or protein_ref_i[i].name == args.bkbn_H)]) + # exit() + protein_ref_indices_i = np.array( + [ + i + 1 + for i in range(len(protein_ref_i.atoms)) + if (protein_ref_i[i].element_name != "H" or protein_ref_i[i].name == args.bkbn_H) + ] + ) + protein_ref_indices_j = np.array( + [ + i + 1 + for i in range(len(protein_ref_j.atoms)) + if (protein_ref_i[i].element_name != "H" or protein_ref_i[i].name == args.bkbn_H) + ] + ) protein_ref_i = [a for a in protein_ref_i.atoms if (a.element_name != "H" or a.name == args.bkbn_H)] - protein_ref_j = [a for a in protein_ref_j.atoms if(a.element_name != "H" or a.name == args.bkbn_H)] + protein_ref_j = [a for a in protein_ref_j.atoms if (a.element_name != "H" or a.name == args.bkbn_H)] print(len(protein_ref_i)) print("\n\n") - #print(len(protein_mego_i)) + # print(len(protein_mego_i)) sorter_i = [str(x.residue.number) + map_if_exists(x.name) for x in protein_ref_i] sorter_mego_i = [str(x.residue.number) + x.name for x in protein_mego_i] @@ -847,9 +859,7 @@ def main_routine(mol_i, mol_j, topology_mego, topology_ref, molecules_name, pref parser.add_argument( "--mode", help="Sets the caculation to be intra/same/cross for histograms processing", default="intra+same+cross" ) - parser.add_argument( - "--bkbn_H", help="Name of backbone hydrogen (default H, charmm HN)", default="H" - ) + parser.add_argument("--bkbn_H", help="Name of backbone hydrogen (default H, charmm HN)", default="H") parser.add_argument("--out", default="./", help="""Sets the output path""") parser.add_argument( "--out_name",