Skip to content

At dhisto #690

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
brunostega merged 3 commits into from
Jul 29, 2026
Merged

At dhisto #690

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
acdae73
new_atomtypes
brunostega Jul 22, 2026
df443d1
new atomtypes CR,SM and new normalization matrix
brunostega Jul 28, 2026
8d3076e
hb 0.9
brunostega Jul 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions multi-ego-basic.ff/aminoacids.rtp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -499,7 +499,7 @@
H H 0.000 0
CA CAH 0.000 0
CB CH2 0.000 0
CG C 0.000 0
CG CR 0.000 0
ND1 NR 0.000 0
CD2 CH 0.000 0
CE1 CH 0.000 0
Expand Down Expand Up @@ -728,7 +728,7 @@
CA CAH 0.000 0
CB CH2 0.000 0
CG CH2 0.000 0
SD S 0.000 0
SD SM 0.000 0
CE CH3 0.000 0
C C 0.000 0
O O 0.000 0
Expand Down Expand Up @@ -778,7 +778,7 @@
H H 0.000 0
CA CAH 0.000 0
CB CH2 0.000 0
CG CH 0.000 0
CG CR 0.000 0
CD1 CH 0.000 0
CD2 CH 0.000 0
CE1 CH 0.000 0
Expand Down Expand Up @@ -991,11 +991,11 @@
H H 0.000 0
CA CAH 0.000 0
CB CH2 0.000 0
CG CH 0.000 0
CG CR 0.000 0
CD1 CH 0.000 0
CD2 CH 0.000 0
CD2 CR 0.000 0
NE1 NR 0.000 0
CE2 CH 0.000 0
CE2 CR 0.000 0
CE3 CH 0.000 0
CZ2 CH 0.000 0
CZ3 CH 0.000 0
Expand Down Expand Up @@ -1083,7 +1083,7 @@
H H 0.000 0
CA CAH 0.000 0
CB CH2 0.000 0
CG CH 0.000 0
CG CR 0.000 0
CD1 CH 0.000 0
CD2 CH 0.000 0
CE1 CH 0.000 0
Expand Down
2 changes: 2 additions & 0 deletions multi-ego-basic.ff/atomtypes.atp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
C 11.01100 ; bare carbon
CZ 11.01100 ; bare carbon
CH 13.01900 ; aromatic CH
CR 12.01100 ; aromatic C
CAH 12.01100 ; aliphatic or sugar CH-group
CH1 13.01900 ; aliphatic or sugar CH-group
CH1t 13.01900 ; aliphatic or sugar CH-group
Expand All @@ -20,6 +21,7 @@
CH3 15.03500 ; aliphatic CH3-group
CH2r 14.02700 ; CH2-group in a ring
S 32.06000 ; sulfur
SM 32.06000 ; sulfur MET
SH 32.06000 ; sulfur
H 4.03200 ; backbone Hydrogen
; LIPIDS
Expand Down
2 changes: 2 additions & 0 deletions multi-ego-basic.ff/ffnonbonded.itp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
NZ 7 0.000 0.000 A 0.0000000000 8.752940e-07
NE 7 0.000 0.000 A 0.0000000000 8.752940e-07
C 6 0.000 0.000 A 0.0000000000 2.598570e-06
CR 6 0.000 0.000 A 0.0000000000 2.598570e-06
CZ 6 0.000 0.000 A 0.0000000000 2.598570e-06
CH 6 0.000 0.000 A 0.0000000000 2.598570e-06
H 1 0.000 0.000 A 0.0000000000 3.239247e-09
Expand All @@ -21,6 +22,7 @@
CH3 6 0.000 0.000 A 0.0000000000 8.595562e-06
CH2r 6 0.000 0.000 A 0.0000000000 1.193966e-05
S 16 0.000 0.000 A 0.0000000000 2.724050e-06
SM 16 0.000 0.000 A 0.0000000000 2.724050e-06
SH 16 0.000 0.000 A 0.0000000000 2.724050e-06
CH3p 6 0.000 0.000 A 0.0000000000 8.736473e-06
P 15 0.000 0.000 A 0.0000000000 3.893600e-06
Expand Down
47 changes: 47 additions & 0 deletions src/multiego/Minimum_CUTOFF_Manual_modif.csv
View file
Open in desktop

Large diffs are not rendered by default.

Binary file modified src/multiego/atdhisto.pkl
View file
Open in desktop
Binary file not shown.
41 changes: 25 additions & 16 deletions src/multiego/interaction_matrix.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,11 @@ def find_class(self, module, name):

class InteractionMatrix:

def __init__(self, pkl_file, emax=0.16, pth=None, show=False):
def __init__(self, pkl_file, emax=0.16, c12_rep_df=None, pth=None, show=False):
self.atmat = self.read_pickle_file(pkl_file)
self.emax = emax
self.pth = pth
self.c12_rep_df = c12_rep_df
self.define_p_threshold()
self.nonlocal_matrix()
self.plot_energy_matrix(show=show)
Expand All @@ -44,7 +45,7 @@ def regall(self, ps, preference, emax, pth):
emax: maximum energy value
pth: threshold probability to split attractive and repulsive interactions
"""
e = - emax/np.log(1e-4) * np.log(ps/(np.maximum(1e-4, preference*pth)))
e = - emax/np.log(1e-4) * np.log((ps + 1e-10)/(np.maximum(1e-4, preference*pth)))
e /= np.max(e)
e = np.where(e==-np.inf, -2, e)
e = np.where(e > 0.0, 0.07 + e* (emax-0.07), e)
Expand All @@ -59,8 +60,8 @@ def plot_energy_matrix(self, out_path = None, show = False):
# 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", "OA", "CZ", "S", "SH", "CH2r", "CH1t"]
attype_ordering_apolar = ["CH1","CH", "CH2", "CH3"]
attype_ordering_polar = [ "NE", "NR", "NT", "OA", "CZ", "S", "SM", "SH", "CH2r", "CH1t"]
attype_ordering_apolar = ["CH1","CR", "CH", "CH2", "CH3"]



Expand Down Expand Up @@ -157,14 +158,17 @@ def read_pickle_file(self, filename):
tot_reps = []
sum_probs = []
for key in data.keys():
probs.append(data[key].p_repeats)
# probs.append(data[key].p_repeats_n2)
probs.append(data[key].p_density)
# probs.append(data[key].p_mindist_repeats_allpdbs)
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["probability"] = np.array(probs)
# atmat["probability"] = np.where(np.array(cutoffs) > 0 ,np.array(probs)/np.array(cutoffs)**2, np.array(probs))
# atmat["p_water"] =np.zeros(len(probs))
atmat["exp_aver"] = np.array(dists)
atmat["attype1"] = atmat["atom_pair"].str.split("_").str[0]
Expand All @@ -174,6 +178,10 @@ def read_pickle_file(self, filename):
atmat["kde"] = kde
atmat["tot_repeats"] = tot_reps
atmat["sum_probs"] = sum_probs
# atmat["probability"] = atmat["sum_probs"]/atmat["tot_repeats"]**(0.5*1.8) if np.sum(atmat["sum_probs"]) > 0 else 0
# fill nan values in probability with 0
atmat["probability"] = atmat["probability"].fillna(0)

# 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
Expand All @@ -193,8 +201,8 @@ def define_p_threshold(self):
print("Neither atpairref nor pth provided, or both provided. Optimizing pth to best fit KL_SCALE.")
# find the minimum pth such that the following atom-pairs have energy < 0
# OM_OM, O_O, NL_NL, O_OM
pth = 0.001 # Initialize with a default value
dp = 0.001
pth = 0.01 # Initialize with a default value
dp = 0.01
while True:
self.atmat["energy"] = np.array(self.regall(self.atmat["probability"].to_numpy(), 1, self.emax, pth))
if (self.atmat.loc[self.atmat["atom_pair"]=="NL_NL", "energy"].values[0] < 0):
Expand All @@ -209,7 +217,7 @@ def define_p_threshold(self):
pth_max = pth
while True:
self.atmat["energy"] = np.array(self.regall(self.atmat["probability"].to_numpy(), 1, self.emax, pth_max))
if self.atmat.loc[self.atmat["atom_pair"]=="OM_NZ", "energy"].values[0] < 0:
if self.atmat.loc[self.atmat["atom_pair"]=="CH_CH", "energy"].values[0] < 0:
break
pth_max += dp
pth_max -= dp
Expand Down Expand Up @@ -244,18 +252,19 @@ def define_special_nonlocal_dict(self):
"atomtypes": ([atom_pair.split("_")[0]], [atom_pair.split("_")[1]]),
"interaction": "rep",
"sigma": None,
"epsilon": None,
"epsilon": float(self.c12_rep_df.loc[((self.c12_rep_df["atp1"]==atp1) & (self.c12_rep_df["atp2"]==atp2)) | ((self.c12_rep_df["atp1"]==atp2) & (self.c12_rep_df["atp2"]==atp1)), "c12"].values[0]),
})
return special_nonlocal_dict

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.atmat.loc[self.atmat["atom_pair"]=="O_H", "energy"] = 0.9
self.atmat.loc[self.atmat["atom_pair"]=="O_N", "energy"] = 0.9
self.atmat.loc[self.atmat["atom_pair"]=="O_C", "energy"] = 0.9
self.atmat.loc[self.atmat["atom_pair"]=="O_O", "energy"] = 0.08
self.atmat.loc[self.atmat["atom_pair"]=="N_N", "energy"] = 0.08
self.atmat.loc[self.atmat["atom_pair"]=="C_C", "energy"] = 0.08

self.special_nonlocal_dict = self.define_special_nonlocal_dict()


60 changes: 35 additions & 25 deletions src/multiego/type_definitions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@
"NZ",
"NE",
"C",
"CR",
"CZ",
"CH",
"CH1",
Expand All @@ -64,6 +65,7 @@
"CH3",
"CH2r",
"S",
"SM",
"SH",
"CH3p",
"P",
Expand All @@ -72,7 +74,7 @@
"H",
"C0",
],
"at.num": [8, 8, 8, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 16, 16, 6, 15, 8, 6, 1, 20],
"at.num": [8, 8, 8, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 16, 16, 16, 6, 15, 8, 6, 1, 20],
"rc_c12": [
2.5 * 0.262134**12, # "O", 2.631580e-07
2.5 * 0.253061**12, # "OM", 1.724403e-07
Expand All @@ -84,6 +86,7 @@
2.5 * 0.289746**12, # "NZ", 8.752940e-07
2.5 * 0.289746**12, # "NE", 8.752940e-07
2.5 * 0.317248**12, # "C", 2.598570e-06
2.5 * 0.317248**12, # "CR", 2.598570e-06
2.5 * 0.317248**12, # "CZ", 2.598570e-06
2.5 * 0.317248**12, # "CH", 2.598570e-06
2.5 * 0.415167**12, # "CH1", 6.555574e-05
Expand All @@ -94,6 +97,7 @@
2.5 * 0.350505**12, # "CH3", 8.595562e-06
2.5 * 0.360236**12, # "CH2r",1.193966e-05
2.5 * 0.318498**12, # "S", 2.724050e-06
2.5 * 0.318498**12, # "SM", 2.724050e-06
2.5 * 0.318498**12, # "SH", 2.724050e-06
2.5 * 0.350981**12, # "CH3p",8.736473e-06
2.5 * 0.328121**12, # "P", 3.893600e-06
Expand All @@ -114,6 +118,7 @@
4.0 * 0.31365**12 * eps_NZ, # "NZ", sig=0.31365
4.0 * 0.31365**12 * eps_NE, # "NE", sig=0.31365
4.0 * 0.35812**12 * eps_C, # "C", sig=0.35812
4.0 * 0.35812**12 * eps_C, # "CR", sig=0.35812
4.0 * 0.35812**12 * eps_CZ, # "CZ", sig=0.35812
4.0 * 0.35812**12 * eps_CH, # "CH", sig=0.35812
4.0 * 0.44592**12 * eps_CH1, # "CH1", sig=0.50192
Expand All @@ -124,6 +129,7 @@
4.0 * 0.37479**12 * eps_CH3, # "CH3", sig=0.37479
4.0 * 0.39547**12 * eps_CH2r, # "CH2r", sig=0.39547
4.0 * 0.33077**12 * eps_S, # "S", sig=0.33077
4.0 * 0.33077**12 * eps_S, # "S", sig=0.33077
4.0 * 0.33077**12 * eps_SH, # "SH", sig=0.33077
4.0 * 0.37479**12 * eps_CH3p, # "CH3p", sig=0.37479
4.0 * 0.33856**12 * eps_P, # "P", sig=0.33856
Expand All @@ -143,6 +149,7 @@
4.0 * 0.31365**6 * eps_NZ, # "NZ",
4.0 * 0.31365**6 * eps_NE, # "NE",
4.0 * 0.35812**6 * eps_C, # "C",
4.0 * 0.35812**6 * eps_C, # "C",
4.0 * 0.35812**6 * eps_CZ, # "CZ",
4.0 * 0.35812**6 * eps_CH, # "CH"
4.0 * 0.44592**6 * eps_CH1, # "CH1"
Expand All @@ -153,6 +160,7 @@
4.0 * 0.37479**6 * eps_CH3, # "CH3"
4.0 * 0.39547**6 * eps_CH2r, # "CH2r"
4.0 * 0.33077**6 * eps_S, # "S",
4.0 * 0.33077**6 * eps_S, # "S",
4.0 * 0.33077**6 * eps_SH, # "SH",
4.0 * 0.37479**6 * eps_CH3p, # "CH3p"
4.0 * 0.33856**6 * eps_P, # "P",
Expand Down Expand Up @@ -332,12 +340,31 @@
# 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.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
_this_dir = os.path.dirname(os.path.abspath(__file__))
_csv_path = os.path.join(_this_dir, "Minimum_CUTOFF_Manual_modif.csv")
# _csv_path = os.path.join(_this_dir, "Merged_c12_sigma.csv")

# Fallback: older layout had the CSV in src/multiego relative to project root
if not os.path.exists(_csv_path):
_csv_path = os.path.join(_this_dir, "..", "src", "multiego", "Minimum_CUTOFF_Manual_modif.csv")
try:
_c12_df = pd.read_csv(_csv_path, sep="\s+")

except Exception:
# If file not found or parsing fails, expose an empty list so other code
# can import the module without crashing.
raise RuntimeError(
f"Failed to load Minimum_CUTOFF_Manual_modif.csv from {_csv_path}. "
"Please ensure the file exists and is formatted correctly."
)


EMAX = 0.14 #maximum epsilon value for the colorbar in the interaction matrix plot
P_TH = 1.0#0.009#0.008#0.12 # if P_TH is None it will be chosen in the InteractionMatrix class to have NL-NL repulsive
SHOW = False # if SHOW is True the interaction matrix will be plotted and saved in the current directory
# PKL = "atdhisto_density.pkl"
PKL = "atdhisto.pkl"
matrix = InteractionMatrix(pkl_file=PKL, emax = EMAX, pth=P_TH, show=SHOW)
matrix = InteractionMatrix(pkl_file=PKL, emax = EMAX, c12_rep_df = _c12_df, 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])
Expand Down Expand Up @@ -402,26 +429,9 @@
# NTHBOND_C12_OVERRIDES available automatically when the module is imported
# import pandas as pd

_this_dir = os.path.dirname(os.path.abspath(__file__))
_csv_path = os.path.join(_this_dir, "Minimum_CUTOFF_Manual_modif.csv")
# _csv_path = os.path.join(_this_dir, "Merged_c12_sigma.csv")

# Fallback: older layout had the CSV in src/multiego relative to project root
if not os.path.exists(_csv_path):
_csv_path = os.path.join(_this_dir, "..", "src", "multiego", "Minimum_CUTOFF_Manual_modif.csv")

try:
_c12_df = pd.read_csv(_csv_path, sep="\s+")
NTHBOND_C12_OVERRIDES = [
({str(row["atp1"]).strip()}, {str(row["atp2"]).strip()}, float(row["c12"])) for _, row in _c12_df.iterrows()
]
except Exception:
# If file not found or parsing fails, expose an empty list so other code
# can import the module without crashing.
raise RuntimeError(
f"Failed to load Minimum_CUTOFF_Manual_modif.csv from {_csv_path}. "
"Please ensure the file exists and is formatted correctly."
)
NTHBOND_C12_OVERRIDES = [
({str(row["atp1"]).strip()}, {str(row["atp2"]).strip()}, float(row["c12"])) for _, row in _c12_df.iterrows()
]
# NTHBOND_C12_OVERRIDES = [
# ({"O", "OM"}, {"O", "OM"}, mg_OO_c12_rep),
# ({"OM"}, {"OM"}, mg_OMOM_c12_rep),
Expand Down
Loading