Skip to content

Commit

Permalink
reverted to original sigma values
Browse files Browse the repository at this point in the history 
hydrogen with only h-h h-o interactions
updated dihedrals
  • Loading branch information
@carlocamilloni
carlocamilloni committed Dec 17, 2024
1 parent f742489 commit e4c08bf
Show file tree
Hide file tree
Showing 6 changed files with 126 additions and 113 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,3 +137,4 @@ outputs/*
analysis/*
.DS_Store
*.ff
!multi-ego-basic.ff
15 changes: 6 additions & 9 deletions multi-ego-basic.ff/ffbonded.itp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -374,16 +374,16 @@
; -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
#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 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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
33 changes: 21 additions & 12 deletions src/multiego/ensemble.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"])
Expand All @@ -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"])
Expand All @@ -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"])
Expand All @@ -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"]
Expand Down Expand Up @@ -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
Expand Down
158 changes: 77 additions & 81 deletions src/multiego/resources/type_definitions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
{
Expand Down Expand Up @@ -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",
],
}
)
Expand Down Expand Up @@ -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),
Expand Down
6 changes: 3 additions & 3 deletions src/multiego/topology.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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=[
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down
26 changes: 18 additions & 8 deletions tools/make_mat/make_mat.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand Down Expand Up @@ -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",
Expand Down

0 comments on commit e4c08bf

Please sign in to comment.