diff --git a/multi-ego-basic.ff/aminoacids.rtp b/multi-ego-basic.ff/aminoacids.rtp index d92da36a..220febd8 100644 --- a/multi-ego-basic.ff/aminoacids.rtp +++ b/multi-ego-basic.ff/aminoacids.rtp @@ -202,6 +202,7 @@ ; ai aj ak al gromos type -CA -C N CA gd_14 -C N CA C gd_44N + -C N CA C gd_44Nb -C N CA C gd_43N N CA CB CG gd_34 N CA C +N gd_45N @@ -552,6 +553,7 @@ ; ai aj ak al gromos type -CA -C N CA gd_14 -C N CA C gd_44H + -C N CA C gd_44Hb -C N CA C gd_43H N CA CB CG gd_34 N CA C +N gd_45H @@ -710,6 +712,7 @@ ; ai aj ak al gromos type -CA -C N CA gd_14 -C N CA C gd_44K + -C N CA C gd_44Kb -C N CA C gd_43K N CA CB CG gd_34 N CA C +N gd_45K @@ -927,6 +930,7 @@ ; ai aj ak al gromos type -CA -C N CA gd_14 -C N CA C gd_44S + -C N CA C gd_44Sb -C N CA C gd_43S N CA CB OG gd_34 N CA C +N gd_45S diff --git a/multi-ego-basic.ff/ffbonded.itp b/multi-ego-basic.ff/ffbonded.itp index eac40834..b4d4c66a 100644 --- a/multi-ego-basic.ff/ffbonded.itp +++ b/multi-ego-basic.ff/ffbonded.itp @@ -391,16 +391,19 @@ ;;;;;;;;;;;;;;;;;;;; ; ASN ;;;;;;;;;;;;;;;;;;;; -#define gd_42N 130.0000 1.5000 2 +#define gd_42N 100.0000 1.5000 2 ; Backbone dihedral angle -N-CA-C-N- ASN ; -#define gd_43N 15.0000 0.600 3 +#define gd_43N -35.0000 1.700 3 ; Backbone dihedral angle -C-N-CA-C- ASN ; -#define gd_44N -110.000 3.3 1 +#define gd_44N -85.000 2.7 1 ; Backbone dihedral angle -C-N-CA-C- ASN ; -#define gd_45N 180.00 2.6 1 +#define gd_44Nb 30.000 2.2 4 +; Backbone dihedral angle -C-N-CA-C- ASN +; +#define gd_45N 0.00 1.6 1 ; Backbone dihedral angle -N-CA-C-N- ASN ; ;;;;;;;;;;;;;;;;;;;; @@ -421,31 +424,32 @@ ;;;;;;;;;;;;;;;;;;;; ; CYS ;;;;;;;;;;;;;;;;;;;; -#define gd_42C 98.0000 2.0000 2 +#define gd_42C 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- CYS ; -#define gd_43C 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43C 70.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- CYS ; -#define gd_43Cb -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_44C 68.000 3 1 ; Backbone dihedral angle -C-N-CA-C- CYS ; -#define gd_44C 68.000 3.0 1 +#define gd_43Cb -10.0000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- CYS ; -#define gd_45C -125.00 1.0 1 +; +#define gd_45C 0.00 5.50 1 ; Backbone dihedral angle -N-CA-C-N- CYS ; ;;;;;;;;;;;;;;;;;;;; ;GLU ;;;;;;;;;;;;;;;;;;;; -#define gd_42E 120.0000 2.5000 2 +#define gd_42E 90.0000 2.5000 2 ; Backbone dihedral angle -N-CA-C-N- GLU ; #define gd_43E 15.0000 2.60 3 ;heigher k to shift down J coupling and have a cleaner distinction between betas ; Backbone dihedral angle -C-N-CA-C- GLU ; -#define gd_44E -120.000 1.50 1 +#define gd_44E -160.000 1.50 1 ; Backbone dihedral angle -C-N-CA-C- GLU ; #define gd_45E 180.00 2.900 1 @@ -454,28 +458,30 @@ ;;;;;;;;;;;;;;;;;;;; ; HIS ;;;;;;;;;;;;;;;;;;;; -#define gd_42H 150.0000 2.5000 2 +#define gd_42H 80.0000 2.5000 2 ; Backbone dihedral angle -N-CA-C-N- HIS ; -#define gd_43H 15.0000 0.20 3 +#define gd_43H -50.0000 2.70 3 ; Backbone dihedral angle -C-N-CA-C- HIS ; -#define gd_44H -110.000 3.10 1 +#define gd_44H -25.000 1.60 1 +; Backbone dihedral angle -C-N-CA-C- HIS +#define gd_44Hb 40.000 2.40 4 ; Backbone dihedral angle -C-N-CA-C- HIS ; -#define gd_45H 160.00 0.5 1 +#define gd_45H 0.00 3.0 1 ; Backbone dihedral angle -N-CA-C-N- HIS ; ;;;;;;;;;;;;;;;;;;;; ; LEU ;;;;;;;;;;;;;;;;;;;; -#define gd_42L 120.0000 2.5000 2 +#define gd_42L 90.0000 2.5000 2 ; Backbone dihedral angle -N-CA-C-N- LEU ; #define gd_43L 15.0000 0.600 3 ; Backbone dihedral angle -C-N-CA-C- LEU ; -#define gd_44L -140.000 2.5 1 +#define gd_44L -160.000 2.5 1 ; Backbone dihedral angle -C-N-CA-C- LEU ; #define gd_45L 180.00 2.0 1 @@ -484,7 +490,7 @@ ;;;;;;;;;;;;;;;;;;;; ; MET ;;;;;;;;;;;;;;;;;;;; -#define gd_42M 120.0000 2.5000 2 +#define gd_42M 90.0000 2.5000 2 ; Backbone dihedral angle -N-CA-C-N- MET ; #define gd_43M 15.0000 0.600 3 @@ -499,13 +505,13 @@ ;;;;;;;;;;;;;;;;;;;; ; ARG ;;;;;;;;;;;;;;;;;;;; -#define gd_42R 120.0000 2.5000 2 +#define gd_42R 95.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- ARG ; #define gd_43R 15.0000 1.600 3 ; Backbone dihedral angle -C-N-CA-C- ARG ; -#define gd_44R -110.000 2.50 1 +#define gd_44R -150.000 2.50 1 ; Backbone dihedral angle -C-N-CA-C- ARG ; #define gd_45R 180.000 2.50 1 @@ -514,31 +520,31 @@ ;;;;;;;;;;;;;;;;;;;; ; TRP ;;;;;;;;;;;;;;;;;;;; -#define gd_42W 98.0000 2.0000 2 +#define gd_42W 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- TRP ; -#define gd_43W 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43W 80.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- TRP ; -#define gd_43Wb -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Wb -10.000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- TRP ; #define gd_44W 68.000 3.0 1 ; Backbone dihedral angle -C-N-CA-C- TRP ; -#define gd_45W -125.00 2.0 1 +#define gd_45W 0.00 6.0 1 ; Backbone dihedral angle -N-CA-C-N- TRP ; ;;;;;;;;;;;;;;;;;;;; ; GLN ;;;;;;;;;;;;;;;;;;;; -#define gd_42Q 120.0000 1.5000 2 +#define gd_42Q 90.0000 1.5000 2 ; Backbone dihedral angle -N-CA-C-N- GLN ; #define gd_43Q 15.0000 1.000 3 ; Backbone dihedral angle -C-N-CA-C- GLN ; -#define gd_44Q -140.000 2.60 1 +#define gd_44Q 180.000 2.60 1 ; Backbone dihedral angle -C-N-CA-C- GLN ; #define gd_45Q 180.000 2.50 1 @@ -547,122 +553,128 @@ ;;;;;;;;;;;;;;;;;;;; ; LYS ;;;;;;;;;;;;;;;;;;;; -#define gd_42K 150.0000 1.5000 2 +#define gd_42K 100.0000 2.3000 2 ; Backbone dihedral angle -N-CA-C-N- LYS ; -#define gd_43K 15.0000 2.000 3 +#define gd_43K 30.0000 2.000 3 ; Backbone dihedral angle -C-N-CA-C- LYS ; -#define gd_44K -120.000 1.00 1 +#define gd_44K 0.000 1.50 1 ; Backbone dihedral angle -C-N-CA-C- LYS ; -#define gd_45K 180.000 1.40 1 +#define gd_44Kb 65.000 1.80 4 +; Backbone dihedral angle -C-N-CA-C- LYS +; +#define gd_45K 0.000 3.0 1 ; Backbone dihedral angle -N-CA-C-N- LYS ; ;;;;;;;;;;;;;;;;;;;; ; ILE ;;;;;;;;;;;;;;;;;;;; -#define gd_42I 98.0000 2.0000 2 +#define gd_42I 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- ILE ; -#define gd_43I 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43I 80.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- ILE ; -#define gd_43Ib -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Ib -10.0000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- ILE ; #define gd_44I 68.000 3.0 1 ; Backbone dihedral angle -C-N-CA-C- ILE ; -#define gd_45I -125.00 1.0 1 +#define gd_45I 0.00 6.0 1 ; Backbone dihedral angle -N-CA-C-N- ILE ; ;;;;;;;;;;;;;;;;;;;; ; VAL ;;;;;;;;;;;;;;;;;;;; -#define gd_42V 98.0000 2.0000 2 +#define gd_42V 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- VAL ; -#define gd_43V 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43V 95.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- VAL ; -#define gd_43Vb -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Vb -10.0000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- VAL ; -#define gd_44V 68.000 3.0 1 +#define gd_44V 50.000 3.0 1 ; Backbone dihedral angle -C-N-CA-C- VAL ; -#define gd_45V -125.00 0.0 1 +#define gd_45V 0.00 6.0 1 ; Backbone dihedral angle -N-CA-C-N- VAL ;;;;;;;;;;;;;;;;;;;; ; PHE ;;;;;;;;;;;;;;;;;;;; -#define gd_42F 98.0000 2.0000 2 +#define gd_42F 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- PHE ; -#define gd_43F 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43F 80.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- PHE ; -#define gd_43Fb -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Fb -10.0000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- PHE ; #define gd_44F 68.000 3.0 1 ; Backbone dihedral angle -C-N-CA-C- PHE ; -#define gd_45F -125.00 0.7 1 +#define gd_45F 0.00 5.7 1 ; Backbone dihedral angle -N-CA-C-N- PHE ;;;;;;;;;;;;;;;;;;;; ;SER ;;;;;;;;;;;;;;;;;;;; -#define gd_42S 120.0000 2.5000 2 +#define gd_42S 110.0000 2.9000 2 ; Backbone dihedral angle -N-CA-C-N- SER ; -#define gd_43S 15.0000 2.000 3 +#define gd_43S -15.0000 2.000 3 +; Backbone dihedral angle -C-N-CA-C- SER +; +#define gd_44S -10.000 1.50 1 ; Backbone dihedral angle -C-N-CA-C- SER ; -#define gd_44S -140.000 1.50 1 +#define gd_44Sb 40.000 1.70 4 ; Backbone dihedral angle -C-N-CA-C- SER ; -#define gd_45S 180.00 0.60 1 +#define gd_45S 0.00 3.00 1 ; Backbone dihedral angle -N-CA-C-N- SER ; ;;;;;;;;;;;;;;;;;;;; ;THR ;;;;;;;;;;;;;;;;;;;; -#define gd_42T 98.0000 2.0000 2 +#define gd_42T 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- THR ; -#define gd_43T 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43T 80.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- THR ; -#define gd_43Tb -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Tb -10.0000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- THR ; #define gd_44T 68.000 3.0 1 ; Backbone dihedral angle -C-N-CA-C- THR ; -#define gd_45T -125.00 1.5 1 +#define gd_45T 0.00 6.0 1 ; Backbone dihedral angle -N-CA-C-N- THR ;;;;;;;;;;;;;;;;;;;; ; TYR ;;;;;;;;;;;;;;;;;;;; -#define gd_42Y 98.0000 2.0000 2 +#define gd_42Y 30.0000 2.0000 2 ; Backbone dihedral angle -N-CA-C-N- TYR ; -#define gd_43Y 70.0000 4.300 3 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Y 80.0000 2.200 3 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- TYR ; -#define gd_43Yb -24.0000 5.500 2 ; lower k for higher J coupling and more population in intermediate beta region +#define gd_43Yb -10.0000 2.700 2 ; lower k for higher J coupling and more population in intermediate beta region ; Backbone dihedral angle -C-N-CA-C- TYR ; #define gd_44Y 68.000 3.0 1 ; Backbone dihedral angle -C-N-CA-C- TYR ; -#define gd_45Y -125.00 1.5 1 +#define gd_45Y 0.00 5.5 1 ; Backbone dihedral angle -N-CA-C-N- TYR diff --git a/src/multiego/atdhisto.pkl b/src/multiego/atdhisto.pkl index 70c9b55d..ee4c8957 100644 Binary files a/src/multiego/atdhisto.pkl and b/src/multiego/atdhisto.pkl differ diff --git a/src/multiego/interaction_matrix.py b/src/multiego/interaction_matrix.py index ce295117..b6fa2153 100644 --- a/src/multiego/interaction_matrix.py +++ b/src/multiego/interaction_matrix.py @@ -21,13 +21,13 @@ def find_class(self, module, name): class InteractionMatrix: - def __init__(self, emax=0.16, pth=None): - self.atmat = self.read_pickle_file("atdhisto.pkl") + def __init__(self, pkl_file, emax=0.16, pth=None, show=False): + self.atmat = self.read_pickle_file(pkl_file) self.emax = emax self.pth = pth self.define_p_threshold() self.nonlocal_matrix() - self.plot_energy_matrix() + self.plot_energy_matrix(show=show) def roundPartial(self, value, resolution=0.005): """ @@ -51,15 +51,17 @@ def regall(self, ps, preference, emax, pth): e = np.array([self.roundPartial(ei) for ei in e]) return e - def plot_energy_matrix(self, out_path = None): + def plot_energy_matrix(self, out_path = None, show = False): """ Plot the energy matrix as a heatmap. """ # TODO read them by gromos to avoid hardcoding them here + bkbnd_atoms = ["C", "O", "N", "CAH", "CAH2", "H"] attype_ordering_charged = [ "OM", "NL", "NZ"] - attype_ordering_polar = [ "NE", "NR", "NT", "O", "OA", "H","N", "CZ", "S", "SH", "CAH", "CH2", "CAH2", "CH1t"] - attype_ordering_apolar = ["C", "CH", "CH1", "CH3"] + attype_ordering_polar = [ "NE", "NR", "NT", "OA", "CZ", "S", "SH", "CH2r", "CH1t"] + attype_ordering_apolar = ["CH1","CH", "CH2", "CH3"] + list_atom1 = np.sort(np.array(self.atmat["attype1"].unique())) @@ -72,8 +74,8 @@ def plot_energy_matrix(self, out_path = None): list_atom2 = [at for at in list_atom2 if at not in noplot] # create ordered list of atom types based on attype_ordering - order_list_atom1 = [at for at in attype_ordering_charged if at in list_atom1] + [at for at in attype_ordering_polar if at in list_atom1] + [at for at in attype_ordering_apolar if at in list_atom1] - order_list_atom2 = [at for at in attype_ordering_charged if at in list_atom2] + [at for at in attype_ordering_polar if at in list_atom2] + [at for at in attype_ordering_apolar if at in list_atom2] + order_list_atom1 = [at for at in bkbnd_atoms if at in list_atom1] + [at for at in attype_ordering_charged if at in list_atom1] + [at for at in attype_ordering_polar if at in list_atom1] + [at for at in attype_ordering_apolar if at in list_atom1] + order_list_atom2 = [at for at in bkbnd_atoms if at in list_atom2] + [at for at in attype_ordering_charged if at in list_atom2] + [at for at in attype_ordering_polar if at in list_atom2] + [at for at in attype_ordering_apolar if at in list_atom2] energies = [] for at1 in order_list_atom1: @@ -121,7 +123,8 @@ def plot_energy_matrix(self, out_path = None): ax.grid() if out_path is not None: plt.savefig(out_path) - plt.show() + if show: + plt.show() def read_pickle_file(self, filename): @@ -149,16 +152,36 @@ def read_pickle_file(self, filename): probs = [] dists = [] cutoffs = [] + bins = [] + kde = [] + tot_reps = [] + sum_probs = [] for key in data.keys(): probs.append(data[key].p_repeats) dists.append(data[key].exp_aver) cutoffs.append(data[key].cutoff) + bins.append(data[key].xbins) + kde.append(data[key].kde) + tot_reps.append(data[key].tot_repeats) + sum_probs.append(np.sum(data[key].sum_probs)) atmat["probability"] = np.array(probs)#/np.array(cutoffs)**2 # atmat["p_water"] =np.zeros(len(probs)) atmat["exp_aver"] = np.array(dists) atmat["attype1"] = atmat["atom_pair"].str.split("_").str[0] atmat["attype2"] = atmat["atom_pair"].str.split("_").str[1] atmat["cutoff"] = np.array(cutoffs) + atmat["bins"] = bins + atmat["kde"] = kde + atmat["tot_repeats"] = tot_reps + atmat["sum_probs"] = sum_probs + # where O-H set distance to 0.195 + # atmat.loc[atmat["atom_pair"]=="O_H", "exp_aver"] = 0.195 + # atmat.loc[atmat["atom_pair"]=="O_N", "exp_aver"] = 0.29 + # atmat.loc[atmat["atom_pair"]=="C_N", "exp_aver"] = 0.41 + ps = [] + for i in range(len(atmat)): + ps.append((np.sum(kde[i]*np.diff(bins[i])[0]/bins[i]**2/ tot_reps[i] )))#* cutoffs[i]**3 )) ) + # atmat["probability"] = np.array(ps) return atmat @@ -177,7 +200,7 @@ def define_p_threshold(self): if (self.atmat.loc[self.atmat["atom_pair"]=="NL_NL", "energy"].values[0] < 0): # self.atmat.loc[self.atmat["atom_pair"]=="OM_OM", "energy"].values[0] < 0 and # self.atmat.loc[self.atmat["atom_pair"]=="O_O", "energy"].values[0] < 0 and - # self.atmat.loc[self.atmat["atom_pair"]=="O_OM", "energy"].values[0] < 0) + # self.atmat.loc[self.atmat["atom_pair"]=="OA_OA", "energy"].values[0] < 0): break pth += dp # find also the maximum pth such that the following atom-pairs have energy > 0 @@ -227,6 +250,12 @@ def define_special_nonlocal_dict(self): def nonlocal_matrix(self): self.atmat["energy"] = np.array(self.regall(self.atmat["probability"].to_numpy(), 1, self.emax, self.pth)) + self.atmat.loc[self.atmat["atom_pair"]=="O_H", "energy"] = 0.45 + self.atmat.loc[self.atmat["atom_pair"]=="O_N", "energy"] = 0.45 + self.atmat.loc[self.atmat["atom_pair"]=="O_C", "energy"] = 0.45 + self.atmat.loc[self.atmat["atom_pair"]=="O_O", "energy"] = 0.09 + self.atmat.loc[self.atmat["atom_pair"]=="N_N", "energy"] = 0.09 + self.special_nonlocal_dict = self.define_special_nonlocal_dict() diff --git a/src/multiego/type_definitions.py b/src/multiego/type_definitions.py index b66af1f5..ca5873a6 100644 --- a/src/multiego/type_definitions.py +++ b/src/multiego/type_definitions.py @@ -332,9 +332,12 @@ # from .atdhist_matrix import special_nonlocal as special_non_local from .interaction_matrix import InteractionMatrix -EMAX = 0.15 # maximum epsilon value for the colorbar in the interaction matrix plot -P_TH = 0.15 # if P_TH is None it will be chosen in the InteractionMatrix class to have NL-NL repulsive -matrix = InteractionMatrix(emax = EMAX, pth=P_TH) +EMAX = 0.15 #maximum epsilon value for the colorbar in the interaction matrix plot +P_TH = 0.006#None#1.0#0.3 # 0.005 # if P_TH is None it will be chosen in the InteractionMatrix class to have NL-NL repulsive +SHOW = False +# PKL = "atdhisto_density.pkl" +PKL = "atdhisto.pkl" +matrix = InteractionMatrix(pkl_file=PKL, emax = EMAX, pth=P_TH, show=SHOW) special_non_local = matrix.special_nonlocal_dict # print(len(special_non_local_A), len(special_non_local)) # print(special_non_local_A[0]) diff --git a/tools/make_mat/make_mat.py b/tools/make_mat/make_mat.py index a2dfab9a..c71c3146 100644 --- a/tools/make_mat/make_mat.py +++ b/tools/make_mat/make_mat.py @@ -495,17 +495,33 @@ def generate_c12_values(df, types, combinations, molecule_type): # ].values[0] # print(f"c12 for {a}-{b}: {aaa}") - # TODO vectorize it as in lj.py? - c12_matrix = np.array( + print("Mapping c12 values for all atom-type combinations...") + # c12_matrix = np.array( + # [ + # type_definitions._c12_df.loc[ + # ((type_definitions._c12_df["atp1"] == a) & (type_definitions._c12_df["atp2"] == b)) + # | ((type_definitions._c12_df["atp1"] == b) & (type_definitions._c12_df["atp2"] == a)), + # "c12", + # ].values[0] + # for a, b in all_combinations + # ] + # ).reshape(len(mego_types_i), len(mego_types_j)) + # vectorized version + c12_lookup = pd.concat( [ - type_definitions._c12_df.loc[ - ((type_definitions._c12_df["atp1"] == a) & (type_definitions._c12_df["atp2"] == b)) - | ((type_definitions._c12_df["atp1"] == b) & (type_definitions._c12_df["atp2"] == a)), - "c12", - ].values[0] - for a, b in all_combinations + type_definitions._c12_df.set_index( + type_definitions._c12_df["atp1"] + "_" + type_definitions._c12_df["atp2"] + )["c12"], + type_definitions._c12_df.set_index( + type_definitions._c12_df["atp2"] + "_" + type_definitions._c12_df["atp1"] + )["c12"], ] - ).reshape(len(mego_types_i), len(mego_types_j)) + ) + # if a pair appears in both directions with different values, keep one consistently + c12_lookup = c12_lookup[~c12_lookup.index.duplicated(keep="first")] + # select only the relevant pairs which are present in the current molecule: all_combinations + all_combinations_str = [f"{a}_{b}" for a, b in all_combinations] + c12_matrix = c12_lookup.reindex(all_combinations_str).to_numpy().reshape(len(mego_types_i), len(mego_types_j)) all_c12 = np.sqrt(df["c12"].to_numpy() * df["c12"].to_numpy()[:, np.newaxis]) all_c12 = c12_matrix.copy() @@ -705,7 +721,6 @@ def main_routine(mol_i, mol_j, topology_mego, topology_ref, molecules_name, pref c12_lookup.update(zip(custom_c12_dict.name, custom_c12_dict.rc_c12)) topology_df_i["c12"] = topology_df_i["mego_type"].map(c12_lookup) - print(topology_df_i["mego_type"]) # preparing topology of molecule j topology_df_j["ref_ai"] = protein_ref_indices_j topology_df_j["ref_type"] = [a.name for a in protein_ref_j]