From d353b77593096906d0fb9e47ed7c9f7a89244122 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Thu, 16 Jan 2025 01:41:02 +0100 Subject: [PATCH] updated mg model and fixes --- multi-ego-basic.ff/ffbonded.itp | 22 ++---- src/multiego/ensemble.py | 87 +++++++++++++++------- src/multiego/resources/type_definitions.py | 1 + 3 files changed, 68 insertions(+), 42 deletions(-) diff --git a/multi-ego-basic.ff/ffbonded.itp b/multi-ego-basic.ff/ffbonded.itp index ec857a51..a22d80b4 100644 --- a/multi-ego-basic.ff/ffbonded.itp +++ b/multi-ego-basic.ff/ffbonded.itp @@ -64,7 +64,7 @@ #define gb_20 0.143500 2.500000e+05 ; ; CHn - OA (sugar) 600 ; -#define gb_21 0.147000 1.000000e+05 ; +#define gb_21 0.147000 9.000000e+04 ; ; CHn - N, NT, NL, NZ, NE 900 ; #define gb_22 0.148000 2.500000e+05 ; @@ -373,24 +373,18 @@ #define gd_41 0.000 3.77 6 ; -CHn-NT- 0.9 ; -; OPT November 23, 22 -;#define gd_42 -170.0000 2.7000 2 -#define gd_42 -170.0000 4.0000 2 +#define gd_42 170.0000 2.5000 2 ; Backbone dihedral angle -N-CA-C-N- ALA/SER/THR ; -;#define gd_43 15.0000 2.500 3 #define gd_43 15.0000 2.500 3 ; Backbone dihedral angle -C-N-CA-C- ALA/SER/THR ; -;#define gd_44 180.000 4.00 1 -#define gd_44 165.000 4.00 1 +#define gd_44 180.000 5.50 1 ; Backbone dihedral angle -C-N-CA-C- ALA/SER/THR ; -;#define gd_45 -105.00 0.2200 1 -#define gd_45 -180.00 0.40 1 +#define gd_45 180.00 0.20 1 ; Backbone dihedral angle -N-CA-C-N- ALA/SER/THR ; -; OPT November 24, 22 #define gd_42b 130.0000 2.5000 2 ; Backbone dihedral angle -N-CA-C-N- BULKYAA ; @@ -400,27 +394,25 @@ #define gd_44b -110.000 3.90 1 ; Backbone dihedral angle -C-N-CA-C- BULKYAA ; -#define gd_45b 180.000 0.80 1 +#define gd_45b 180.000 0.20 1 ; Backbone dihedral angle -N-CA-C-N- BULKYAA ; -; OPT April, 3 2024 #define gd_48 180.0 5 2 ; Backbone dihedral angle -N-CA-C-N GLY ; #define gd_49 180 11 1 ; Backbone dihedral angle -C-N-CA-C GLY ; -#define gd_51 -180 2.5 1 +#define gd_51 -180 1.8 1 ; Backbone dihedral angle -N-CA-C-N GLY ; -; OPT November, 23 2022 #define gd_52 150.0000 4.0000 2 ; Backbone dihedral angle -N-CA-C-N PRO ; #define gd_53 0.0000 4.000 4 ; Backbone dihedral angle -C-N-CA-C PRO ; -#define gd_55 -60.0000 0.400 1 +#define gd_55 -60.0000 1.200 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 03307794..ab0c3654 100644 --- a/src/multiego/ensemble.py +++ b/src/multiego/ensemble.py @@ -784,19 +784,33 @@ def init_LJ_datasets(meGO_ensemble, matrices, pairs14, exclusion_bonds14, args): symmetrize=True, ) - pairwise_c12 = np.where( - oxygen_mask, - 3e-6, + ho_mask = masking.create_linearized_mask( + type_ai_mapped.to_numpy(), + type_aj_mapped.to_numpy(), + [("H", "O"), ("H", "OM")], + symmetrize=True, + ) + + # Define condition where only ai or aj (but not both) starts with "H" + h_condition = ( + train_dataset["ai"].str.startswith("H") ^ train_dataset["aj"].str.startswith("H") + ) + + hh_condition = ( + train_dataset["ai"].str.startswith("H") & train_dataset["aj"].str.startswith("H") + ) + + pairwise_c12 = ( np.sqrt( train_dataset["ai"].map(meGO_ensemble["sbtype_c12_dict"]) * train_dataset["aj"].map(meGO_ensemble["sbtype_c12_dict"]) - ), + ) ) train_dataset["rep"] = train_dataset["rep"].fillna(pd.Series(pairwise_c12)) + train_dataset.loc[oxygen_mask, "rep"] = 3e-6 + print(train_dataset.to_string()) - pairwise_mg_sigma = np.where( - oxygen_mask, - (3e-6) ** (1 / 12), + pairwise_mg_sigma = ( ( train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c12_dict"]) * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c12_dict"]) @@ -805,26 +819,37 @@ def init_LJ_datasets(meGO_ensemble, matrices, pairs14, exclusion_bonds14, args): * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c6_dict"]) ) ) - ** (1 / 12), + ** (1 / 12) ) train_dataset["mg_sigma"] = pd.Series(pairwise_mg_sigma) - - pairwise_mg_epsilon = np.where( - oxygen_mask, - -3e-6, - ( - train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c6_dict"]) - * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c6_dict"]) + train_dataset.loc[oxygen_mask, "mg_sigma"] = (3e-6) ** (1 / 12) + # Apply the specific value for this condition + #train_dataset.loc[h_condition, "mg_sigma"] = 0. + train_dataset.loc[hh_condition, "mg_sigma"] = train_dataset["rep"] ** (1 /12) + train_dataset.loc[ho_mask, "mg_sigma"] = 0.169500 + + # Generate the default pairwise_mg_epsilon + pairwise_mg_epsilon = ( + train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c6_dict"]) + * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c6_dict"]) + ) / ( + 4 + * np.sqrt( + train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c12_dict"]) + * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c12_dict"]) ) - / ( - 4 - * np.sqrt( - train_dataset["ai"].map(meGO_ensemble["sbtype_mg_c12_dict"]) - * train_dataset["aj"].map(meGO_ensemble["sbtype_mg_c12_dict"]) - ) - ), ) + + # Initialize pairwise_mg_epsilon with default values train_dataset["mg_epsilon"] = pd.Series(pairwise_mg_epsilon) + train_dataset.loc[oxygen_mask, "mg_epsilon"] = -3e-6 + + # Apply the specific value for this condition + train_dataset.loc[h_condition, "mg_epsilon"] = 0.0 + train_dataset.loc[hh_condition, "mg_epsilon"] = -train_dataset["rep"] + train_dataset.loc[ho_mask, "mg_epsilon"] = 0.11 + + train_dataset.dropna(subset=["mg_sigma"], inplace=True) return train_dataset @@ -839,7 +864,7 @@ def generate_OO_LJ(meGO_ensemble): ] H_H_sbtype = [sbtype for sbtype, atomtype in meGO_ensemble["sbtype_type_dict"].items() if atomtype == "H"] - full_matrix = list(itertools.product(H_H_sbtype, O_OM_sbtype)) + list(itertools.product(O_OM_sbtype, H_H_sbtype)) + full_matrix_OH = list(itertools.product(H_H_sbtype, O_OM_sbtype)) + list(itertools.product(O_OM_sbtype, H_H_sbtype)) # Generate all possible combinations combinations = list(itertools.product(O_OM_sbtype, repeat=2)) @@ -861,9 +886,7 @@ def generate_OO_LJ(meGO_ensemble): HH_LJ["sigma"] = HH_LJ["c12"] ** (1.0 / 12.0) HH_LJ["mg_sigma"] = HH_LJ["c12"] ** (1 / 12) HH_LJ["mg_epsilon"] = -HH_LJ["c12"] - HO_LJ = pd.DataFrame(full_matrix, columns=["ai", "aj"]) - # HO_LJ["c12"] = 1.153070e-08 * type_definitions.mg_eps - # HO_LJ["c6"] = 2.147622e-04 * type_definitions.mg_eps + HO_LJ = pd.DataFrame(full_matrix_OH, columns=["ai", "aj"]) HO_LJ["c12"] = 2.249554e-09 * type_definitions.mg_eps HO_LJ["c6"] = 9.485893e-05 * type_definitions.mg_eps HO_LJ["epsilon"] = type_definitions.mg_eps @@ -1548,8 +1571,18 @@ def make_pairs_exclusion_topology(meGO_ensemble, meGO_LJ_14, args): filtered_combinations.append((sbtype_with_residue[i][0], sbtype_with_residue[j][0])) j += 1 + # Filter out invalid combinations + valid_combinations = [ + (ai, aj) + for ai, aj in filtered_combinations + if not ( + (meGO_ensemble["sbtype_type_dict"][ai] == "H" and meGO_ensemble["sbtype_type_dict"][aj] not in {"H", "O", "OM"}) + or (meGO_ensemble["sbtype_type_dict"][aj] == "H" and meGO_ensemble["sbtype_type_dict"][ai] not in {"H", "O", "OM"}) + ) + ] + # Create a DataFrame from the filtered combinations - df = pd.DataFrame(filtered_combinations, columns=["ai", "aj"]) + df = pd.DataFrame(valid_combinations, columns=["ai", "aj"]) df["c6"] = 0.0 df["c12"] = np.sqrt( df["ai"].map(meGO_ensemble["sbtype_c12_dict"]) * df["aj"].map(meGO_ensemble["sbtype_c12_dict"]) diff --git a/src/multiego/resources/type_definitions.py b/src/multiego/resources/type_definitions.py index 79cb3362..82d89e24 100644 --- a/src/multiego/resources/type_definitions.py +++ b/src/multiego/resources/type_definitions.py @@ -120,6 +120,7 @@ # Dictionary mapping atom types from a force field to multiego representation from_ff_to_multiego = { + "HN": "H", "OC1": "O1", "OC2": "O2", "OT1": "O1",