Make-all CG with VS

README.md

@@ -10,7 +10,7 @@ Swarm-CG is designed for automatically optimizing the bonded terms of a coarse-g
 
 ### Publication
 
-> Empereur-mot, C.; Pesce, L.; Bochicchio, D.; Perego, C.; Pavan, G.M. (2020) Swarm-CG: Automatic Parametrization of Bonded Terms in Coarse-Grained Models of Simple to Complex Molecules via Fuzzy Self-Tuning Particle Swarm Optimization. [ChemRxiv. Preprint](https://doi.org/10.26434/chemrxiv.12613427)
+> Empereur-mot, C.; Pesce, L.; Bochicchio, D.; Perego, C.; Pavan, G.M. (2020) Swarm-CG: Automatic Parametrization of Bonded Terms in MARTINI-based Coarse-Grained Models of Simple to Complex Molecules via Fuzzy Self-Tuning Particle Swarm Optimization. [ChemRxiv. Preprint](https://doi.org/10.26434/chemrxiv.12613427)
 
 ### Installation & Usage
 

swarmcg/config.py

@@ -6,7 +6,11 @@
 kB = 0.008314462
 sim_temperature = 300  # Kelvin
 bi_nb_bins = 50  # nb of bins to use for Boltzmann Inversion, will be doubled for dihedrals distributions binning during BI -- this has huge impact on the results of the BI and this value shall STAY AT 50 ! actually I did not try to modify much but this feels like dangerous atm
-bonds_max_range = 20  # nm -- used to define grid for EMD calculations so increasing this only slightly increases computation time, however small bw for bonds has real impact
+bonds_max_range = 15  # nm -- used to define grid for EMD calculations
+# NOTE: increasing bonds_max_range increases computation time, but memory usage increases exponentially
+# TODO: detect when a bond has longer values than bonds_max_range and suggest the user to raise the limit if he really
+#       needs to, but I don't really see what kind of use case would require more than 5 nm (maybe elastic net in
+#       proteins though)
 bw_constraints = 0.002  # nm
 bw_bonds = 0.01  # nm
 bw_angles = 2.5  # degrees
@@ -57,7 +61,7 @@
     'angle': [1, 2],  # tested and verified: 1, 2
     'dihedral': [1, 2, 4],  # tested and verified: 1, 2, 4 -- ongoing: 9 (need to merge the 1+ dihedrals groups on plots)
     'virtual_sites2': [1],  # tested and verified: 1 -- ongoing: 2 (need GMX 2020)
-    'virtual_sites3': [1, 2],  # tested and verified: 1, 2 -- ongoing: 3, 4
+    'virtual_sites3': [1, 2, 3],  # tested and verified: 1, 2, 3 -- ongoing: 4
     'virtual_sites4': [],  # ongoing: 2 -- irrelevant: 1
     'virtual_sitesn': [1, 2, 3]  # tested and verified: 1, 2, 3
 }
@@ -66,7 +70,7 @@
 # plots display parameters
 use_hists = False  # hists are not implemented in a way that they will be displayed with left and right bold borders atm
 line_alpha = 0.6  # line alpha for the density plots
-fill_alpha = 0.35  # fill alpha for the density plots
+fill_alpha = 0.30  # fill alpha for the density plots
 cg_color = '#1f77b4'
 atom_color = '#d62728'
 

swarmcg/evaluate_model.py

@@ -213,9 +213,9 @@ def main():
 								action='store_true', default=False)
 
 	# display help if script was called without arguments
-	if len(sys.argv) == 1:
-		args_parser.print_help()
-		sys.exit()
+	# if len(sys.argv) == 1:
+	# 	args_parser.print_help()
+	# 	sys.exit()
 
 	# arguments handling, display command line if help or no arguments provided
 	ns = args_parser.parse_args()

swarmcg/optimize_model.py

@@ -133,7 +133,10 @@ def run(ns):
         for top_line in top_lines:
             if top_line.startswith('#include'):
                 top_include = top_line.split()[1].replace('"', '').replace("'", '')  # remove potential single and double quotes around filenames
-                top_includes_filenames.append(os.path.dirname(arg_entries[user_provided_filenames[5]]) + '/' + top_include)
+                arg_dirname = os.path.dirname(arg_entries[user_provided_filenames[5]])
+                if arg_dirname == '':
+                    arg_dirname = '.'
+                top_includes_filenames.append(arg_dirname + '/' + top_include)
 
     # check gmx arguments conflicts
     if ns.gmx_args_str != '' and (ns.nb_threads != 0 or ns.gpu_id != ''):
@@ -248,7 +251,7 @@ def run(ns):
     with open(ns.exec_folder+'/'+config.opti_pairwise_distances_file, 'w'):
         pass
 
-    # set these to None to then check the variables have been filled (!= None), so we will do these calculations
+    # set these to None to then check the variables have been filled (is not None), so we will do these calculations
     # one single time in function compare_models that is called at each iteration during optimization
     ns.gyr_aa_mapped, ns.gyr_aa_mapped_std = None, None
     ns.sasa_aa_mapped, ns.sasa_aa_mapped_std = None, None

swarmcg/simulations/vs_functions.py

@@ -74,9 +74,6 @@ def vs3_func_2(ns, traj, vs_def_beads_ids, vs_params, bead_id):
     i, j, k = vs_def_beads_ids
     a, b = vs_params  # weight, nm
     b = b * 10  # retrieve amgstrom for MDA
-    if a < 0 or a > 1:
-        msg = f"Virtual site ID {bead_id + 1} uses an incorrect weight. Expected weight in [0 , 1]."
-        raise exceptions.MissformattedFile(msg)
 
     for ts in ns.aa2cg_universe.trajectory:
         pos_i = ns.aa2cg_universe.atoms[i].position
@@ -94,6 +91,7 @@ def vs3_func_3(ns, traj, vs_def_beads_ids, vs_params):
 
     i, j, k = vs_def_beads_ids
     ang_deg, d = vs_params  # degrees, nm
+    ang_rad = np.deg2rad(ang_deg)  # retrieve radians
     d = d * 10  # retrieve amgstrom for MDA
 
     for ts in ns.aa2cg_universe.trajectory:
@@ -102,10 +100,7 @@ def vs3_func_3(ns, traj, vs_def_beads_ids, vs_params):
         pos_k = ns.aa2cg_universe.atoms[k].position
         r_ij = pos_j - pos_i
         r_jk = pos_k - pos_j
-        ang_rad = np.deg2rad(ang_deg)
-        # comb_ijk = r_jk - (np.dot(r_ij, r_jk) / np.dot(r_ij, r_ij)) * r_ij
-        comb_ijk = r_jk - ((r_ij * r_jk) / (r_ij * r_ij)) * r_ij  # BAD
-        # comb_ijk = r_jk - (np.cross(r_ij, r_jk) / np.cross(r_ij, r_ij)) * r_ij  # BAD
+        comb_ijk = r_jk - (np.dot(r_ij, r_jk) / np.dot(r_ij, r_ij)) * r_ij
         traj[ts.frame] = pos_i + d * np.cos(ang_rad) * (r_ij / mda.lib.mdamath.norm(r_ij)) + d * np.sin(ang_rad) * (comb_ijk / mda.lib.mdamath.norm(comb_ijk))
 
 
@@ -115,22 +110,22 @@ def vs3_func_4(ns, traj, vs_def_beads_ids, vs_params):
 
     i, j, k = vs_def_beads_ids
     a, b, c = vs_params  # weight, weight, nm**(-1)
-    c = c * 10  # retrieve amgstrom**(-1) for MDA
+    c = c / 10  # retrieve amgstrom**(-1) for MDA
 
     for ts in ns.aa2cg_universe.trajectory:
         pos_i = ns.aa2cg_universe.atoms[i].position
         pos_j = ns.aa2cg_universe.atoms[j].position
         pos_k = ns.aa2cg_universe.atoms[k].position
         r_ij = pos_j - pos_i
         r_ik = pos_i - pos_k
-        traj[ts.frame] = pos_i + a * r_ij + b * r_ik + c * (r_ij * r_ik)
+        traj[ts.frame] = pos_i - a * r_ij - b * r_ik + c * (r_ij * r_ik)
 
 
 # Functions for virtual_sites4
 
-# vs_4 func 2 -> Linear combination using 3 reference points -> ?
+# vs_4 func 2 -> Linear combination using 3 reference points
 # NOTE: only function 2 is defined for vs_4 in GROMACS, because it replaces function 1
-#       which still existss for retro compatibility but should be ignored
+#       which still exists for retro compatibility but its usage must be avoided
 def vs4_func_2(ns, traj, vs_def_beads_ids, vs_params):
 
     pass