diff --git a/multiego.py b/multiego.py index 7a24a825..2766e134 100644 --- a/multiego.py +++ b/multiego.py @@ -103,10 +103,8 @@ def meGO_parsing(): print(f"ERROR: --epsilon_min ({ref['epsilon_min']}) must be greater than 0.") sys.exit() - if ref["epsilon"] < ref["epsilon_min"]: - raise ValueError( - f"ERROR: in {ref}. \nEpsilon value epsi = {ref['epsilon']} is below the minimum meaningful value." - ) + if ref["epsilon"] <= 0.0: + raise ValueError(f"ERROR: in {ref}. \nEpsilon value epsi = {ref['epsilon']} must be greater than 0.") custom_dict = {} if args.custom_dict: diff --git a/src/multiego/arguments.py b/src/multiego/arguments.py index 2c7adca6..cbfdc58a 100644 --- a/src/multiego/arguments.py +++ b/src/multiego/arguments.py @@ -83,7 +83,7 @@ "help": "partition function normalization", }, "--relative_c12d": { - "default": 0.001, + "default": 0.01, "type": float, "help": "Relative deviation from default to set new replulsive c12", }, diff --git a/src/multiego/ensemble.py b/src/multiego/ensemble.py index 7bc7ee02..594bf81f 100644 --- a/src/multiego/ensemble.py +++ b/src/multiego/ensemble.py @@ -171,9 +171,9 @@ def initialize_molecular_contacts(contact_matrix, prior_matrix, args, reference) 1 - (-(reference["epsilon_min"] - np.maximum(0, prior_matrix["epsilon_prior"])) / contact_matrix["epsilon_0"]) ) - #modify limit_rc_att in the cases where epsilon_prior is negative and limit_rc_att is below 1 == epsilon_0 < epsilon_min) - contact_matrix.loc[(contact_matrix["limit_rc_att"] < 1 ) & (prior_matrix["epsilon_prior"] < 0), "limit_rc_att"] = 1 - + # modify limit_rc_att in the cases where epsilon_prior is negative and limit_rc_att is below 1 == epsilon_0 < epsilon_min) + contact_matrix.loc[(contact_matrix["limit_rc_att"] < 1) & (prior_matrix["epsilon_prior"] < 0), "limit_rc_att"] = 1 + contact_matrix["limit_rc_rep"] = contact_matrix["rc_threshold"] ** ( np.maximum(0, prior_matrix["epsilon_prior"]) / contact_matrix["epsilon_0"] ) * contact_matrix["zf"] ** (1 + (np.maximum(0, prior_matrix["epsilon_prior"]) / contact_matrix["epsilon_0"])) @@ -964,7 +964,7 @@ def set_sig_epsilon(meGO_LJ, parameters): meGO_LJ["epsilon"] = np.where( meGO_LJ["epsilon_prior"] < 0, -meGO_LJ["rep"] * (meGO_LJ["distance"] / meGO_LJ["rc_distance"]) ** 12, # If epsilon_prior == 0 - meGO_LJ["epsilon_prior"], # Otherwise, set to epsilon_prior + 0, # meGO_LJ["epsilon_prior"], # Otherwise, set to epsilon_prior ) # Attractive interactions @@ -985,19 +985,15 @@ def set_sig_epsilon(meGO_LJ, parameters): # negative epsilon are used to identify non-attractive interactions meGO_LJ.loc[ ( - ( - np.maximum(meGO_LJ["probability"], meGO_LJ["rc_threshold"]) - < meGO_LJ["limit_rc_rep"] * np.maximum(meGO_LJ["rc_probability"], meGO_LJ["rc_threshold"]) - ) + (meGO_LJ["probability"] < meGO_LJ["limit_rc_att"] * np.maximum(meGO_LJ["rc_probability"], meGO_LJ["rc_threshold"])) & (meGO_LJ["probability"] > meGO_LJ["md_threshold"]) ), "epsilon", - ] = -(meGO_LJ["epsilon_0"] / (np.log(meGO_LJ["zf"] * meGO_LJ["rc_threshold"]))) * meGO_LJ["distance"] ** 12 * ( - np.log( - np.maximum(meGO_LJ["probability"], meGO_LJ["rc_threshold"]) - / (meGO_LJ["zf"] * np.maximum(meGO_LJ["rc_probability"], meGO_LJ["rc_threshold"])) - ) - - np.maximum(0.0, meGO_LJ["epsilon_prior"]) + ] = -np.maximum( + 0, + (meGO_LJ["distance"] ** 12) + * (meGO_LJ["epsilon_0"] / (np.log(meGO_LJ["zf"] * meGO_LJ["rc_threshold"]))) + * np.log(meGO_LJ["probability"] / (meGO_LJ["zf"] * np.maximum(meGO_LJ["rc_probability"], meGO_LJ["rc_threshold"]))), ) - ( meGO_LJ["rep"] * (meGO_LJ["distance"] / meGO_LJ["rc_distance"]) ** 12 ) @@ -1011,32 +1007,32 @@ def set_sig_epsilon(meGO_LJ, parameters): meGO_LJ = meGO_LJ[meGO_LJ.epsilon != 0] # lower value for repulsion - meGO_LJ.loc[ - (meGO_LJ["1-4"] != "1_4") & (meGO_LJ["epsilon"] < 0.0) & (-meGO_LJ["epsilon"] < 0.1 * meGO_LJ["rep"]), - "epsilon", - ] = ( - -0.1 * meGO_LJ["rep"] - ) - # higher value for repulsion - meGO_LJ.loc[ - (meGO_LJ["1-4"] != "1_4") & (meGO_LJ["epsilon"] < 0.0) & (-meGO_LJ["epsilon"] > 20.0 * meGO_LJ["rep"]), - "epsilon", - ] = ( - -20.0 * meGO_LJ["rep"] - ) + # meGO_LJ.loc[ + # (meGO_LJ["1-4"] != "1_4") & (meGO_LJ["epsilon"] < 0.0) & (-meGO_LJ["epsilon"] < 0.1 * meGO_LJ["rep"]), + # "epsilon", + # ] = ( + # -0.1 * meGO_LJ["rep"] + # ) + ## higher value for repulsion + # meGO_LJ.loc[ + # (meGO_LJ["1-4"] != "1_4") & (meGO_LJ["epsilon"] < 0.0) & (-meGO_LJ["epsilon"] > 20.0 * meGO_LJ["rep"]), + # "epsilon", + # ] = ( + # -20.0 * meGO_LJ["rep"] + # ) # but within a lower - meGO_LJ.loc[ - (meGO_LJ["1-4"] == "1_4") & (-meGO_LJ["epsilon"] < 0.1 * meGO_LJ["rep"]), - "epsilon", - ] = ( - -0.1 * meGO_LJ["rep"] - ) - # and an upper value - meGO_LJ.loc[ - (meGO_LJ["1-4"] == "1_4") & (-meGO_LJ["epsilon"] > meGO_LJ["rep"]), - "epsilon", - ] = -meGO_LJ["rep"] + # meGO_LJ.loc[ + # (meGO_LJ["1-4"] == "1_4") & (-meGO_LJ["epsilon"] < 0.1 * meGO_LJ["rep"]), + # "epsilon", + # ] = ( + # -0.1 * meGO_LJ["rep"] + # ) + ## and an upper value + # meGO_LJ.loc[ + # (meGO_LJ["1-4"] == "1_4") & (-meGO_LJ["epsilon"] > meGO_LJ["rep"]), + # "epsilon", + # ] = -meGO_LJ["rep"] # Sigma is set from the estimated interaction length meGO_LJ = meGO_LJ.assign(sigma=meGO_LJ["distance"] / 2 ** (1.0 / 6.0))