Merge pull request multi-ego#502 from multi-ego/mglob_prior

Mglob prior
carlocamilloni · Dec 7, 2024 · 2bf68cb · 2bf68cb
2 parents 99d0a9d + 33ba07e
commit 2bf68cb
Show file tree

Hide file tree

Showing 36 changed files with 194,022 additions and 224,391 deletions.
diff --git a/.gitignore b/.gitignore
@@ -137,5 +137,6 @@ outputs/*
 analysis/*
 .DS_Store
 *.ff
-!multi-ego-basic.ff/
+!multi-rc-prior.ff/
+!multi-mg-prior.ff/
 !multi-ego-ready.ff/
diff --git a/multi-ego-basic.ff/aminoacids.rtp b/multi-ego-basic.ff/aminoacids.rtp
@@ -50,7 +50,7 @@
 [ ALA ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH3       0.000     0
     C     C       0.000     0
     O     O       0.000     0
@@ -84,7 +84,7 @@
 [ ARG ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG   CH2       0.000     0
    CD   CH2       0.000     0
@@ -143,7 +143,7 @@
 [ ASN ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG     C       0.000     0
   OD1     O       0.000     0
@@ -190,7 +190,7 @@
 [ ASP ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG     C       0.000     0
   OD1    OM       0.000     0
@@ -237,7 +237,7 @@
 [ CYS ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    SG     S       0.000     0
     C     C       0.000     0
@@ -275,7 +275,7 @@
 [ CYS2 ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    SG     S       0.000     0
     C     C       0.000     0
@@ -313,7 +313,7 @@
 [ GLN ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG   CH2       0.000     0
    CD     C       0.000     0
@@ -364,7 +364,7 @@
 [ GLU ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG   CH2       0.000     0
    CD     C       0.000     0
@@ -443,7 +443,7 @@
 [ HIS ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG     C       0.000     0
   ND1    NR       0.000     0
@@ -504,7 +504,7 @@
 [ ILE ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH1       0.000     0
   CG1   CH2       0.000     0
   CG2   CH3       0.000     0
@@ -551,7 +551,7 @@
 [ LEU ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG   CH1       0.000     0
   CD1   CH3       0.000     0
@@ -598,7 +598,7 @@
 [ LYS ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG   CH2       0.000     0
    CD   CH2       0.000     0
@@ -648,7 +648,7 @@
 [ MET ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG   CH2       0.000     0
    SD     S       0.000     0
@@ -694,7 +694,7 @@
 [ PHE ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG    CH       0.000     0
   CD1    CH       0.000     0
@@ -759,7 +759,7 @@
 [ PRO ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB  CH2r       0.000     0
    CG  CH2r       0.000     0
    CD  CH2r       0.000     0
@@ -808,7 +808,7 @@
 [ SER ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    OG    OA       0.000     0
     C     C       0.000     0
@@ -846,7 +846,7 @@
 [ THR ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH1       0.000     0
   OG1    OA       0.000     0
   CG2   CH3       0.000     0
@@ -889,7 +889,7 @@
 [ TRP ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG    CH       0.000     0
   CD1    CH       0.000     0
@@ -975,7 +975,7 @@
 [ TYR ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH2       0.000     0
    CG    CH       0.000     0
   CD1    CH       0.000     0
@@ -1045,7 +1045,7 @@
 [ VAL ]
  [ atoms ]
     N     N       0.000     0
-   CA   CH1       0.000     0
+   CA   CH1a      0.000     0
    CB   CH1       0.000     0
   CG1   CH3       0.000     0
   CG2   CH3       0.000     0

diff --git a/multi-ego-basic.ff/atomtypes.atp b/multi-ego-basic.ff/atomtypes.atp
@@ -10,6 +10,7 @@
     C  12.01100 ;     bare carbon
    CH  13.01900 ;     aromatic CH 
   CH1  13.01900 ;     aliphatic or sugar CH-group
+ CH1a  13.01900 ;     aliphatic or sugar CH-group
   CH2  14.02700 ;     aliphatic or sugar CH2-group
   CH3  15.03500 ;     aliphatic CH3-group
  CH2r  14.02700 ;     CH2-group in a ring
@@ -20,4 +21,4 @@
    OE  15.99940 ;     ether or ester oxygen
   CR1  13.01900 ;     aromatic CH-group    
 ; IONS
-   C0  40.08000 ;     calcium (charge 2+) DES-AMBER
+   C0  40.08000 ;     calcium (charge 2+) DES-AMBER
diff --git a/multi-ego-basic.ff/forcefield.itp b/multi-ego-basic.ff/forcefield.itp
@@ -1,5 +1,3 @@
-#define _FF_MAGROS
-
 [ defaults ]
 ; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
   1		1		yes		1.0	1.0

diff --git a/multiego.py b/multiego.py
@@ -1,7 +1,6 @@
 import argparse
 import sys
 import os
-import parmed as pmd
 import pandas as pd
 import time
 import gc
@@ -91,8 +90,17 @@ def meGO_parsing():
     if args.multi_epsilon_intra or args.multi_epsilon_inter_domain or args.multi_epsilon_inter:
         multi_flag = True
 
-    mego_topology = pmd.load_file(f"{args.root_dir}/inputs/{args.system}/topol.top")
-    topol_names = [m for m in mego_topology.molecules]
+    custom_dict = {}
+    if args.custom_dict:
+        custom_dict = parse_json(args.custom_dict)
+        if custom_dict == None:
+            print("ERROR: Custom dictionary was parsed, but the dictionary is empty")
+            sys.exit()
+
+    print(f"Running Multi-eGO: {args.egos}\n")
+    print("- Processing Multi-eGO topology")
+    mego_ensemble = ensemble.init_meGO_ensemble(args, custom_dict)
+    topol_names = [m for m in mego_ensemble["topology"].molecules]
 
     args.names = []
     for name in args.multi_epsilon_intra.keys():
@@ -110,7 +118,7 @@ def meGO_parsing():
     elif not multi_flag:
         args.names = topol_names
 
-    if args.egos != "rc" and not args.reference:
+    if args.egos == "production" and not args.reference:
         args.reference = ["reference"]
 
     if args.epsilon and not args.inter_epsilon:
@@ -127,7 +135,7 @@ def meGO_parsing():
         args.multi_epsilon_inter = {k1: {k2: args.inter_epsilon for k2 in args.names} for k1 in args.names}
 
     # check all epsilons are set and greater than epsilon_min
-    if args.egos != "rc":
+    if args.egos == "production":
         for k, v in args.multi_epsilon_intra.items():
             if v < args.epsilon_min:
                 print("ERROR: epsilon value for " + k + " is less than epsilon_min")
@@ -150,13 +158,6 @@ def meGO_parsing():
     elif args.symmetry:
         args.symmetry = io.parse_symmetry_list(args.symmetry)
 
-    custom_dict = {}
-    if args.custom_dict:
-        custom_dict = parse_json(args.custom_dict)
-        if custom_dict == None:
-            print("ERROR: Custom dictionary was parsed, but the dictionary is empty")
-            sys.exit()
-
     custom_c12_dict = pd.DataFrame()
     if args.custom_c12 is not None:
         custom_c12_dict = io.read_custom_c12_parameters(args.custom_c12)
@@ -169,7 +170,7 @@ def meGO_parsing():
         parser.print_usage()
         sys.exit()
 
-    return args, custom_dict
+    return args, mego_ensemble, custom_dict
 
 
 def main():
@@ -179,20 +180,17 @@ def main():
     """
 
     bt = time.time()
-    print("Multi-eGO\n")
-    args, custom_dict = meGO_parsing()
-
-    if not args.no_header:
-        generate_face.print_welcome()
+    generate_face.print_welcome()
+    args, meGO_ensembles, custom_dict = meGO_parsing()
 
+    st = time.time()
+    elapsed_time = st - bt
+    print("- Done in:", elapsed_time, "seconds")
     print("- Checking for input files and folders")
     io.check_files_existence(args)
-    if args.egos != "rc":
+    if args.egos == "production":
         io.check_matrix_format(args)
 
-    print("- Processing Multi-eGO topology")
-    st = time.time()
-    meGO_ensembles = ensemble.init_meGO_ensemble(args, custom_dict)
     print("\t- Generating bonded interactions")
     meGO_ensembles = ensemble.generate_bonded_interactions(meGO_ensembles)
     print("\t- Generating 1-4 data")
@@ -202,7 +200,7 @@ def main():
     st = et
     print("- Done in:", elapsed_time, "seconds")
 
-    if args.egos != "rc":
+    if args.egos == "production":
         print("- Processing Multi-eGO contact matrices")
         meGO_ensembles, matrices = ensemble.init_meGO_matrices(meGO_ensembles, args, custom_dict)
         et = time.time()
@@ -227,7 +225,7 @@ def main():
         elapsed_time = et - st
         st = et
         print("- Done in:", elapsed_time, "seconds")
-    else:
+    elif args.egos == "mg":
         print("- Generate the LJ dataset")
         meGO_LJ = ensemble.generate_OO_LJ(meGO_ensembles)
         meGO_LJ_14 = pairs14
@@ -237,7 +235,7 @@ def main():
         print("- Done in:", elapsed_time, "seconds")
 
     print("- Finalize pairs and exclusions")
-    meGO_LJ_14 = ensemble.make_pairs_exclusion_topology(meGO_ensembles, meGO_LJ_14)
+    meGO_LJ_14 = ensemble.make_pairs_exclusion_topology(meGO_ensembles, meGO_LJ_14, args)
     et = time.time()
     elapsed_time = et - st
     st = et

diff --git a/src/multiego/arguments.py b/src/multiego/arguments.py
@@ -5,14 +5,10 @@
     },
     "--egos": {
         "type": str,
-        "choices": ["rc", "production"],
-        "help": "rc: creates a force-field for random coil simulations. "
+        "choices": ["mg", "production"],
+        "help": "mg: creates a force-field for molten globule simulations."
         "production: creates a force-field combining random coil simulations and training simulations.",
     },
-    "--epsilon": {
-        "type": float,
-        "help": "Maximum interaction energy per contact. The typical range is 0.2-0.4 kJ/mol",
-    },
     "--reference": {
         "type": lambda x: x.split(","),
         "default": [],
@@ -25,13 +21,17 @@
         "help": "A list of the training simulations to be included in multi-eGO, "
         "corresponding to the subfolders to process and where the contacts are learned.",
     },
+    "--epsilon": {
+        "type": float,
+        "help": "Maximum interaction energy per contact.",
+    },
     "--inter_epsilon": {
         "type": float,
-        "help": "Maximum interaction energy per intermolecular contacts. The typical range is 0.2-0.4 kJ/mol",
+        "help": "Maximum interaction energy per intermolecular contacts.",
     },
     "--inter_domain_epsilon": {
         "type": float,
-        "help": "Maximum interaction energy per interdomain contacts. The typical range is 0.2-0.4 kJ/mol",
+        "help": "Maximum interaction energy per interdomain contacts.",
     },
     "--p_to_learn": {
         "type": float,