diff --git a/G1_DATA/cg_topol.tpr b/G1_DATA/cg_topol.tpr index 0ce7463..e69de29 100644 Binary files a/G1_DATA/cg_topol.tpr and b/G1_DATA/cg_topol.tpr differ diff --git a/G1_DATA/cg_traj.xtc b/G1_DATA/cg_traj.xtc index d0f2ad3..e69de29 100644 Binary files a/G1_DATA/cg_traj.xtc and b/G1_DATA/cg_traj.xtc differ diff --git a/README.md b/README.md index 5aafd7e..f572f95 100644 --- a/README.md +++ b/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 diff --git a/swarmcg/config.py b/swarmcg/config.py index f27da0d..462ea18 100644 --- a/swarmcg/config.py +++ b/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' diff --git a/swarmcg/evaluate_model.py b/swarmcg/evaluate_model.py index 0f43a34..c3835cc 100644 --- a/swarmcg/evaluate_model.py +++ b/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() diff --git a/swarmcg/optimize_model.py b/swarmcg/optimize_model.py index 13c9cda..b50fe58 100644 --- a/swarmcg/optimize_model.py +++ b/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 diff --git a/swarmcg/simulations/vs_functions.py b/swarmcg/simulations/vs_functions.py index 61fe9ff..8f76163 100644 --- a/swarmcg/simulations/vs_functions.py +++ b/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,7 +110,7 @@ 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 @@ -123,14 +118,14 @@ def vs3_func_4(ns, traj, vs_def_beads_ids, vs_params): 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 diff --git a/swarmcg/swarmCG.py b/swarmcg/swarmCG.py index 9e99c8b..5fff0a9 100644 --- a/swarmcg/swarmCG.py +++ b/swarmcg/swarmCG.py @@ -135,12 +135,12 @@ def section_switch(section_read, section_active): for section_current in section_read: section_read[section_current] = False - if section_active != None: + if section_active is not None: section_read[section_active] = True # vs_type in [2, 3, 4, n], then they each have specific functions to define their positions -def vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, line_nb): +def vs_error_control(ns, bead_id, vs_type, func, line_nb, vs_def_beads_ids=None): if bead_id >= len(ns.cg_itp['atoms']): msg = ( @@ -149,13 +149,14 @@ def vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, line_nb): ) raise exceptions.MissformattedFile(msg) - for bid in vs_def_beads_ids: - if bid >= len(ns.cg_itp['atoms']): - msg = ( - f"The definition of virtual site ID {bead_id + 1} makes use of ID {bid + 1}, while this ID exceeds" - f" the number of atoms defined in the CG ITP file." - ) - raise exceptions.MissformattedFile(msg) + if vs_def_beads_ids is not None: + for bid in vs_def_beads_ids: + if bid >= len(ns.cg_itp['atoms']): + msg = ( + f"The definition of virtual site ID {bead_id + 1} makes use of ID {bid + 1}, while this ID exceeds" + f" the number of atoms defined in the CG ITP file." + ) + raise exceptions.MissformattedFile(msg) if not ns.cg_itp['atoms'][bead_id]['bead_type'].startswith('v'): msg = ( @@ -400,7 +401,7 @@ def read_cg_itp_file(ns): bead_id = int(sp_itp_line[0])-1 vs_def_beads_ids = [int(bid)-1 for bid in sp_itp_line[1:3]] func = sp_itp_line[3] # will be casted to int in the verification below (for factorizing checks) - func = vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, i) # i is the line number + func = vs_error_control(ns, bead_id, vs_type, func, i, vs_def_beads_ids) # i is the line number vs_params = float(sp_itp_line[4]) ns.cg_itp['atoms'][bead_id]['vs_type'] = vs_type ns.cg_itp['virtual_sites2'][bead_id] = {'bead_id': bead_id, 'func': func, 'vs_def_beads_ids': vs_def_beads_ids, 'vs_params': vs_params} @@ -411,7 +412,7 @@ def read_cg_itp_file(ns): bead_id = int(sp_itp_line[0])-1 vs_def_beads_ids = [int(bid)-1 for bid in sp_itp_line[1:4]] func = sp_itp_line[4] # will be casted to int in the verification below (for factorizing checks) - func = vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, i) # i is the line number + func = vs_error_control(ns, bead_id, vs_type, func, i, vs_def_beads_ids) # i is the line number if func in [1, 2, 3]: vs_params = [float(param) for param in sp_itp_line[5:7]] elif func == 4: @@ -425,7 +426,7 @@ def read_cg_itp_file(ns): bead_id = int(sp_itp_line[0]) - 1 vs_def_beads_ids = [int(bid) - 1 for bid in sp_itp_line[1:5]] func = sp_itp_line[5] # will be casted to int in the verification below (for factorizing checks) - func = vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, i) # i is the line number + func = vs_error_control(ns, bead_id, vs_type, func, i, vs_def_beads_ids) # i is the line number vs_params = [float(param) for param in sp_itp_line[6:9]] ns.cg_itp['atoms'][bead_id]['vs_type'] = vs_type ns.cg_itp['virtual_sites4'][bead_id] = {'bead_id': bead_id, 'func': func, 'vs_def_beads_ids': vs_def_beads_ids, 'vs_params': vs_params} @@ -436,14 +437,14 @@ def read_cg_itp_file(ns): bead_id = int(sp_itp_line[0])-1 func = sp_itp_line[1] # will be casted to int in verification below (for factorizing checks) # here we do the check in 2 steps, because the reading of beads_ids depends on the function - func = vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, i) # i is the line number + func = vs_error_control(ns, bead_id, vs_type, func, i, vs_def_beads_ids=None) # i is the line number if func == 3: vs_def_beads_ids = [int(sp_itp_line[2:][i])-1 for i in range(0, len(sp_itp_line[2:]), 2)] vs_params = [float(sp_itp_line[2:][i]) for i in range(1, len(sp_itp_line[2:]), 2)] else: vs_def_beads_ids = [int(bid) - 1 for bid in sp_itp_line[2:]] vs_params = None - func = vs_error_control(ns, bead_id, vs_type, func, vs_def_beads_ids, i) # i is the line number + func = vs_error_control(ns, bead_id, vs_type, func, i, vs_def_beads_ids) # i is the line number ns.cg_itp['atoms'][bead_id]['vs_type'] = vs_type ns.cg_itp['virtual_sitesn'][bead_id] = {'bead_id': bead_id, 'func': func, 'vs_def_beads_ids': vs_def_beads_ids, 'vs_params': vs_params} @@ -817,7 +818,7 @@ def write_cg_itp_file(itp_obj, out_path_itp, print_sections=['constraint', 'bond for i in range(len(itp_obj['atoms'])): # if the ITP did NOT contain masses, they are set at 0 in this field during ITP reading - if itp_obj['atoms'][i]['mass'] != None: + if itp_obj['atoms'][i]['mass'] is not None: fp.write('{0:<4} {1:>4} {6:>2} {2:>6} {3:>6} {4:<4} {5:9.5f} {7:<5.2f}\n'.format( itp_obj['atoms'][i]['bead_id']+1, itp_obj['atoms'][i]['bead_type'], itp_obj['atoms'][i]['residue'], itp_obj['atoms'][i]['atom'], i+1, itp_obj['atoms'][i]['charge'], @@ -1374,7 +1375,7 @@ def get_AA_bonds_distrib(ns, beads_ids, grp_type, grp_nb): print(' Ref. AA-mapped distrib. rescaled to avg', bond_avg_final, 'nm for', grp_type, grp_nb+1, '(initially', bond_avg_init, 'nm)') # or if specific lengths were provided for constraints and/or bonds - elif ns.bonds_scaling_specific != None: + elif ns.bonds_scaling_specific is not None: if grp_type.startswith('constraint'): geom_id_full = 'C'+str(grp_nb+1) @@ -1728,8 +1729,17 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ns.aa2cg_universe._topology.masses.values = np.array(masses) # create fake bonds in the CG MDA universe, that will be used only for making the molecule whole + # we make bonds between each VS and their beads definition, so we retrieve the connectivity + # iteratively towards the real CG beads, that are all connected if len(ns.vs_beads_ids) > 0: - fake_bonds = [[vs_bid, ns.real_beads_ids[0]] for vs_bid in ns.vs_beads_ids] + fake_bonds = [] + for vs_type in ['2', '3', '4', 'n']: + try: + for bead_id in ns.cg_itp['virtual_sites'+vs_type]: + for vs_def_bead_id in ns.cg_itp['virtual_sites'+vs_type][bead_id]['vs_def_beads_ids']: + fake_bonds.append([bead_id, vs_def_bead_id]) + except (IndexError, ValueError): + pass ns.cg_universe.add_bonds(fake_bonds, guessed=False) # select the whole molecule as an MDA atomgroup and make its coordinates whole, inplace, across the complete trajectory @@ -1780,7 +1790,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False constraint_avg, constraint_hist, _ = get_AA_bonds_distrib(ns, beads_ids=ns.cg_itp['constraint'][grp_constraint]['beads'], grp_type='constraints group', grp_nb=grp_constraint) constraints[grp_constraint]['AA']['avg'] = constraint_avg constraints[grp_constraint]['AA']['hist'] = constraint_hist - else: # use atomistic reference that was loaded by the optimization routines + else: # use atomistic reference that was loaded by the optimization routines constraints[grp_constraint]['AA']['avg'] = ns.cg_itp['constraint'][grp_constraint]['avg'] constraints[grp_constraint]['AA']['hist'] = ns.cg_itp['constraint'][grp_constraint]['hist'] @@ -1839,7 +1849,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False bond_avg, bond_hist, _ = get_AA_bonds_distrib(ns, beads_ids=ns.cg_itp['bond'][grp_bond]['beads'], grp_type='bonds group', grp_nb=grp_bond) bonds[grp_bond]['AA']['avg'] = bond_avg bonds[grp_bond]['AA']['hist'] = bond_hist - else: # use atomistic reference that was loaded by the optimization routines + else: # use atomistic reference that was loaded by the optimization routines bonds[grp_bond]['AA']['avg'] = ns.cg_itp['bond'][grp_bond]['avg'] bonds[grp_bond]['AA']['hist'] = ns.cg_itp['bond'][grp_bond]['hist'] @@ -1898,7 +1908,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False angle_avg, angle_hist, _, _ = get_AA_angles_distrib(ns, beads_ids=ns.cg_itp['angle'][grp_angle]['beads']) angles[grp_angle]['AA']['avg'] = angle_avg angles[grp_angle]['AA']['hist'] = angle_hist - else: # use atomistic reference that was loaded by the optimization routines + else: # use atomistic reference that was loaded by the optimization routines angles[grp_angle]['AA']['avg'] = ns.cg_itp['angle'][grp_angle]['avg'] angles[grp_angle]['AA']['hist'] = ns.cg_itp['angle'][grp_angle]['hist'] @@ -1949,7 +1959,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False dihedral_avg, dihedral_hist, _, _ = get_AA_dihedrals_distrib(ns, beads_ids=ns.cg_itp['dihedral'][grp_dihedral]['beads']) dihedrals[grp_dihedral]['AA']['avg'] = dihedral_avg dihedrals[grp_dihedral]['AA']['hist'] = dihedral_hist - else: # use atomistic reference that was loaded by the optimization routines + else: # use atomistic reference that was loaded by the optimization routines dihedrals[grp_dihedral]['AA']['avg'] = ns.cg_itp['dihedral'][grp_dihedral]['avg'] dihedrals[grp_dihedral]['AA']['hist'] = ns.cg_itp['dihedral'][grp_dihedral]['hist'] @@ -1998,12 +2008,12 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False if larger_group > ncols: hidden_cols = larger_group - ncols if ns.atom_only: - print('Displaying max '+str(ncols)+' distributions per row using the CG ITP file ordering of distributions groups ('+str(hidden_cols)+' more are hidden)', flush=True) + print(f'Displaying max {ncols} distributions per row using the CG ITP file ordering of distributions groups ({hidden_cols} more are hidden)') else: if not ns.mismatch_order: - print(styling.header_warning + 'Displaying max ' + str(ncols) + ' distributions groups per row and this can be MISLEADING because ordering by pairwise AA-mapped vs. CG distributions mismatch is DISABLED (' + str(hidden_cols) + ' more are hidden)', flush=True) + print(f'{styling.header_warning}Displaying max {ncols} distributions groups per row and this can be MISLEADING because ordering by pairwise AA-mapped vs. CG distributions mismatch is DISABLED ({hidden_cols} more are hidden)') else: - print('Displaying max '+str(ncols)+' distributions groups per row ordered by pairwise AA-mapped vs. CG distributions difference ('+str(hidden_cols)+' more are hidden)', flush=True) + print(f'Displaying max {ncols} distributions groups per row ordered by pairwise AA-mapped vs. CG distributions difference ({hidden_cols} more are hidden)') else: print() if not ns.mismatch_order: @@ -2012,7 +2022,8 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False print('Distributions groups will be displayed using ranked mismatch score between pairwise AA-mapped and CG distributions') nrows -= sum([ns.nb_constraints == 0, ns.nb_bonds == 0, ns.nb_angles == 0, ns.nb_dihedrals == 0]) - # fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(ncols*3, nrows*3), squeeze=False) # this fucking line was responsible of the big memory leak (figures were not closing) + # fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(ncols*3, nrows*3), squeeze=False) + # this fucking line was responsible of the big memory leak (figures were not closing) so I let this here for memory fig = plt.figure(figsize=(ncols*3, nrows*3)) ax = fig.subplots(nrows=nrows, ncols=ncols, squeeze=False) @@ -2037,7 +2048,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(constraints[grp_constraint]['AA']['avg'], 0, color=config.atom_color, marker='D') if not ns.atom_only: - ax[nrow][i].set_title('Constraint grp '+str(grp_constraint+1)+' - EMD Δ '+str(round(avg_diff_grp_constraints[grp_constraint], 3))) + ax[nrow][i].set_title(f'Constraint grp {grp_constraint + 1} - EMD Δ {round(avg_diff_grp_constraints[grp_constraint], 3)}') if config.use_hists: ax[nrow][i].step(constraints[grp_constraint]['CG']['x'], constraints[grp_constraint]['CG']['y'], label='CG', color=config.cg_color, where='mid', alpha=config.line_alpha) ax[nrow][i].fill_between(constraints[grp_constraint]['CG']['x'], constraints[grp_constraint]['CG']['y'], color=config.cg_color, step='mid', alpha=config.fill_alpha) @@ -2045,13 +2056,11 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(constraints[grp_constraint]['CG']['x'], constraints[grp_constraint]['CG']['y'], label='CG', color=config.cg_color, alpha=config.line_alpha) ax[nrow][i].fill_between(constraints[grp_constraint]['CG']['x'], constraints[grp_constraint]['CG']['y'], color=config.cg_color, alpha=config.fill_alpha) ax[nrow][i].plot(constraints[grp_constraint]['CG']['avg'], 0, color=config.cg_color, marker='D') - # if ns.verbose: - print('Constraint '+str(grp_constraint+1)+' -- AA Avg: '+str(round(constraints[grp_constraint]['AA']['avg'], 3))+' nm -- CG Avg: '+str(round(constraints[grp_constraint]['CG']['avg'], 3))+' nm') + print(f"Constraint {grp_constraint + 1} -- AA Avg: {round(constraints[grp_constraint]['AA']['avg'], 3)} nm -- CG Avg: {round(constraints[grp_constraint]['CG']['avg'], 3)}") else: ax[nrow][i].set_title('Constraint grp '+str(grp_constraint+1)+' - Avg '+str(round(avg_diff_grp_constraints[grp_constraint], 3))+' nm') - print('Constraint '+str(grp_constraint+1)+' -- AA Avg: '+str(round(constraints[grp_constraint]['AA']['avg'], 3))) + print(f"Constraint {grp_constraint + 1} -- AA Avg: {round(constraints[grp_constraint]['AA']['avg'], 3)}") ax[nrow][i].grid(zorder=0.5) - # ax[nrow][i].set_ylim(bottom=0) if ns.row_x_scaling: ax[nrow][i].set_xlim(np.mean(row_wise_ranges['constraints'][grp_constraint])-row_wise_ranges['max_range_constraints']/2*1.1, np.mean(row_wise_ranges['constraints'][grp_constraint])+row_wise_ranges['max_range_constraints']/2*1.1) if i % 2 == 0: @@ -2081,7 +2090,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(bonds[grp_bond]['AA']['avg'], 0, color=config.atom_color, marker='D') if not ns.atom_only: - ax[nrow][i].set_title('Bond grp '+str(grp_bond+1)+' - EMD Δ '+str(round(avg_diff_grp_bonds[grp_bond], 3))) + ax[nrow][i].set_title(f'Bond grp {grp_bond+1} - EMD Δ {round(avg_diff_grp_bonds[grp_bond], 3)}') if config.use_hists: ax[nrow][i].step(bonds[grp_bond]['CG']['x'], bonds[grp_bond]['CG']['y'], label='CG', color=config.cg_color, where='mid', alpha=config.line_alpha) ax[nrow][i].fill_between(bonds[grp_bond]['CG']['x'], bonds[grp_bond]['CG']['y'], color=config.cg_color, step='mid', alpha=config.fill_alpha) @@ -2089,13 +2098,11 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(bonds[grp_bond]['CG']['x'], bonds[grp_bond]['CG']['y'], label='CG', color=config.cg_color, alpha=config.line_alpha) ax[nrow][i].fill_between(bonds[grp_bond]['CG']['x'], bonds[grp_bond]['CG']['y'], color=config.cg_color, alpha=config.fill_alpha) ax[nrow][i].plot(bonds[grp_bond]['CG']['avg'], 0, color=config.cg_color, marker='D') - # if ns.verbose: - print('Bond '+str(grp_bond+1)+' -- AA Avg: '+str(round(bonds[grp_bond]['AA']['avg'], 3))+' nm -- CG Avg: '+str(round(bonds[grp_bond]['CG']['avg'], 3))+' nm') + print(f"Bond {grp_bond + 1} -- AA Avg: {round(bonds[grp_bond]['AA']['avg'], 3)} nm -- CG Avg: {round(bonds[grp_bond]['CG']['avg'], 3)} nm") else: - ax[nrow][i].set_title('Bond grp '+str(grp_bond+1)+' - Avg '+str(round(avg_diff_grp_bonds[grp_bond], 3))+' nm') - print('Bond '+str(grp_bond+1)+' -- AA Avg: '+str(round(bonds[grp_bond]['AA']['avg'], 3))) + ax[nrow][i].set_title(f"Bond grp {grp_bond+1} - Avg {round(avg_diff_grp_bonds[grp_bond], 3)} nm") + print(f"Bond {grp_bond+1} -- AA Avg: {round(bonds[grp_bond]['AA']['avg'], 3)}") ax[nrow][i].grid(zorder=0.5) - # ax[nrow][i].set_ylim(bottom=0) if ns.row_x_scaling: ax[nrow][i].set_xlim(np.mean(row_wise_ranges['bonds'][grp_bond])-row_wise_ranges['max_range_bonds']/2*1.1, np.mean(row_wise_ranges['bonds'][grp_bond])+row_wise_ranges['max_range_bonds']/2*1.1) if i % 2 == 0: @@ -2125,7 +2132,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(angles[grp_angle]['AA']['avg'], 0, color=config.atom_color, marker='D') if not ns.atom_only: - ax[nrow][i].set_title('Angle grp '+str(grp_angle+1)+' - EMD Δ '+str(round(avg_diff_grp_angles[grp_angle], 3))) + ax[nrow][i].set_title(f'Angle grp {grp_angle+1} - EMD Δ {round(avg_diff_grp_angles[grp_angle], 3)}') if config.use_hists: ax[nrow][i].step(angles[grp_angle]['CG']['x'], angles[grp_angle]['CG']['y'], label='CG', color=config.cg_color, where='mid', alpha=config.line_alpha) ax[nrow][i].fill_between(angles[grp_angle]['CG']['x'], angles[grp_angle]['CG']['y'], color=config.cg_color, step='mid', alpha=config.fill_alpha) @@ -2133,13 +2140,11 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(angles[grp_angle]['CG']['x'], angles[grp_angle]['CG']['y'], label='CG', color=config.cg_color, alpha=config.line_alpha) ax[nrow][i].fill_between(angles[grp_angle]['CG']['x'], angles[grp_angle]['CG']['y'], color=config.cg_color, alpha=config.fill_alpha) ax[nrow][i].plot(angles[grp_angle]['CG']['avg'], 0, color=config.cg_color, marker='D') - # if ns.verbose: - print('Angle '+str(grp_angle+1)+' -- AA Avg: '+str(round(angles[grp_angle]['AA']['avg'], 1))+'° -- CG Avg: '+str(round(angles[grp_angle]['CG']['avg'], 1))+'°') + print(f"Angle {grp_angle+1} -- AA Avg: {round(angles[grp_angle]['AA']['avg'], 1)}° -- CG Avg: {round(angles[grp_angle]['CG']['avg'], 1)}°") else: - ax[nrow][i].set_title('Angle grp '+str(grp_angle+1)+' - Avg '+str(round(avg_diff_grp_angles[grp_angle], 1))+'°') - print('Angle '+str(grp_angle+1)+' -- AA Avg: '+str(round(angles[grp_angle]['AA']['avg'], 1))) + ax[nrow][i].set_title(f"Angle grp {grp_angle+1} - Avg {round(avg_diff_grp_angles[grp_angle], 1)}°") + print(f"Angle {grp_angle+1} -- AA Avg: {round(angles[grp_angle]['AA']['avg'], 1)}") ax[nrow][i].grid(zorder=0.5) - # ax[nrow][i].set_ylim(bottom=0) if ns.row_x_scaling: ax[nrow][i].set_xlim(np.mean(row_wise_ranges['angles'][grp_angle])-row_wise_ranges['max_range_angles']/2*1.1, np.mean(row_wise_ranges['angles'][grp_angle])+row_wise_ranges['max_range_angles']/2*1.1) if i % 2 == 0: @@ -2169,7 +2174,7 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(dihedrals[grp_dihedral]['AA']['avg'], 0, color=config.atom_color, marker='D') if not ns.atom_only: - ax[nrow][i].set_title('Dihedral grp '+str(grp_dihedral+1)+' - EMD Δ '+str(round(avg_diff_grp_dihedrals[grp_dihedral], 3))) + ax[nrow][i].set_title(f'Dihedral grp {grp_dihedral+1} - EMD Δ {round(avg_diff_grp_dihedrals[grp_dihedral], 3)}') if config.use_hists: ax[nrow][i].step(dihedrals[grp_dihedral]['CG']['x'], dihedrals[grp_dihedral]['CG']['y'], label='CG', color=config.cg_color, where='mid', alpha=config.line_alpha) ax[nrow][i].fill_between(dihedrals[grp_dihedral]['CG']['x'], dihedrals[grp_dihedral]['CG']['y'], color=config.cg_color, step='mid', alpha=config.fill_alpha) @@ -2177,13 +2182,11 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False ax[nrow][i].plot(dihedrals[grp_dihedral]['CG']['x'], dihedrals[grp_dihedral]['CG']['y'], label='CG', color=config.cg_color, alpha=config.line_alpha) ax[nrow][i].fill_between(dihedrals[grp_dihedral]['CG']['x'], dihedrals[grp_dihedral]['CG']['y'], color=config.cg_color, alpha=config.fill_alpha) ax[nrow][i].plot(dihedrals[grp_dihedral]['CG']['avg'], 0, color=config.cg_color, marker='D') - # if ns.verbose: - print('Dihedral '+str(grp_dihedral+1)+' -- AA Avg: '+str(round(dihedrals[grp_dihedral]['AA']['avg'], 1))+'° -- CG Avg: '+str(round(dihedrals[grp_dihedral]['CG']['avg'], 1))+'°') + print(f"Dihedral {grp_dihedral+1} -- AA Avg: {round(dihedrals[grp_dihedral]['AA']['avg'], 1)}° -- CG Avg: {round(dihedrals[grp_dihedral]['CG']['avg'], 1)}°") else: - ax[nrow][i].set_title('Dihedral grp '+str(grp_dihedral+1)+' - Avg '+str(round(avg_diff_grp_dihedrals[grp_dihedral], 1))+'°') - print('Dihedral '+str(grp_dihedral+1)+' -- AA Avg: '+str(round(dihedrals[grp_dihedral]['AA']['avg'], 1))) + ax[nrow][i].set_title(f'Dihedral grp {grp_dihedral+1} - Avg {round(avg_diff_grp_dihedrals[grp_dihedral], 1)}°') + print(f"Dihedral {grp_dihedral+1} -- AA Avg: {round(dihedrals[grp_dihedral]['AA']['avg'], 1)}") ax[nrow][i].grid(zorder=0.5) - # ax[nrow][i].set_ylim(bottom=0) if ns.row_x_scaling: ax[nrow][i].set_xlim(np.mean(row_wise_ranges['dihedrals'][grp_dihedral])-row_wise_ranges['max_range_dihedrals']/2*1.1, np.mean(row_wise_ranges['dihedrals'][grp_dihedral])+row_wise_ranges['max_range_dihedrals']/2*1.1) if i % 2 == 0: @@ -2224,7 +2227,6 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False fit_score_total, fit_score_constraints_bonds, fit_score_angles, fit_score_dihedrals = 0, 0, 0, 0 for i in range(ns.nb_constraints): - # dist_pairwise = np.sqrt(avg_diff_grp_constraints[diff_ordered_grp_constraints[i]]) dist_pairwise = avg_diff_grp_constraints[diff_ordered_grp_constraints[i]] all_dist_pairwise += str(dist_pairwise)+' ' all_emd_dist_geoms['constraints'].append(dist_pairwise) @@ -2239,7 +2241,6 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False fit_score_constraints_bonds += dist_pairwise for i in range(ns.nb_bonds): - # dist_pairwise = np.sqrt(avg_diff_grp_bonds[diff_ordered_grp_bonds[i]]) dist_pairwise = avg_diff_grp_bonds[diff_ordered_grp_bonds[i]] all_dist_pairwise += str(dist_pairwise)+' ' all_emd_dist_geoms['bonds'].append(dist_pairwise) @@ -2254,7 +2255,6 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False fit_score_constraints_bonds += dist_pairwise for i in range(ns.nb_angles): - # dist_pairwise = np.sqrt(avg_diff_grp_angles[diff_ordered_grp_angles[i]]) dist_pairwise = avg_diff_grp_angles[diff_ordered_grp_angles[i]] all_dist_pairwise += str(dist_pairwise)+' ' all_emd_dist_geoms['angles'].append(dist_pairwise) @@ -2270,7 +2270,6 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False # dihedrals_dist_pairwise = 0 for i in range(ns.nb_dihedrals): - # dist_pairwise = np.sqrt(avg_diff_grp_dihedrals[diff_ordered_grp_dihedrals[i]]) dist_pairwise = avg_diff_grp_dihedrals[diff_ordered_grp_dihedrals[i]] all_dist_pairwise += str(dist_pairwise)+' ' all_emd_dist_geoms['dihedrals'].append(dist_pairwise) @@ -2283,7 +2282,6 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False dist_pairwise = dist_pairwise ** 2 fit_score_dihedrals += dist_pairwise - # dihedrals_dist_pairwise += dist_pairwise fit_score_constraints_bonds = np.sqrt(fit_score_constraints_bonds) fit_score_angles = np.sqrt(fit_score_angles) @@ -2302,11 +2300,10 @@ def compare_models(ns, manual_mode=True, ignore_dihedrals=False, calc_sasa=False print(' Dihedrals constribution to fitness score:', fit_score_dihedrals, flush=True) plt.tight_layout(rect=[0, 0, 1, 0.9]) - # plt.suptitle('FITNESS SCORE\nTotal: '+str(fit_score_total)+' -- Constraints/Bonds: '+str(fit_score_constraints_bonds)+' -- Angles: '+str(fit_score_angles)+' -- Dihedrals: '+str(fit_score_dihedrals)) eval_score = fit_score_total if ignore_dihedrals and ns.nb_dihedrals > 0: eval_score -= fit_score_dihedrals - sup_title = 'FITNESS SCORE\nTotal: '+str(round(eval_score, 3))+' -- Constraints/Bonds: '+str(fit_score_constraints_bonds)+' -- Angles: '+str(fit_score_angles)+' -- Dihedrals: '+str(fit_score_dihedrals) + sup_title = f'FITNESS SCORE\nTotal: {round(eval_score, 3)} -- Constraints/Bonds: {fit_score_constraints_bonds} -- Angles: {fit_score_angles} -- Dihedrals: {fit_score_dihedrals}' if ignore_dihedrals and ns.nb_dihedrals > 0: sup_title += ' (ignored)' plt.suptitle(sup_title) @@ -2393,7 +2390,7 @@ def modify_mdp(mdp_filename, sim_time=None, nb_frames=1500, log_write_freq=5000, nstxout_compressed_line = i # adjust simulation time according to timestep - if sim_time != None: + if sim_time is not None: if dt_line != -1 and nsteps_line != -1: nsteps = int(sim_time*1000 / dt) mdp_lines_in[nsteps_line] = sp_nsteps_line[0]+'= '+str(nsteps)+' ; automatically modified by Swarm-CG' @@ -2411,22 +2408,22 @@ def modify_mdp(mdp_filename, sim_time=None, nb_frames=1500, log_write_freq=5000, raise exceptions.MissformattedFile(msg) # force NOT writting coordinates data, as this can only slow the simulation and we don't need it + nstxout = nsteps if nstxout_line != -1: - nstxout = nsteps mdp_lines_in[nstxout_line] = sp_nstxout_line[0]+'= '+str(nstxout)+' ; automatically modified by Swarm-CG' else: mdp_lines_in.append(f'nstxout = {nstxout} ; automatically added by Swarm-CG') # force NOT writting velocities data, as this can only slow the simulation and we don't need it + nstvout = nsteps if nstvout_line != -1: - nstvout = nsteps mdp_lines_in[nstvout_line] = sp_nstvout_line[0]+'= '+str(nstvout)+' ; automatically modified by Swarm-CG' else: mdp_lines_in.append(f'nstvout = {nstvout} ; automatically added by Swarm-CG') # force NOT writting forces data, as this can only slow the simulation and we don't need it + nstfout = nsteps if nstfout_line != -1: - nstfout = nsteps mdp_lines_in[nstfout_line] = sp_nstfout_line[0]+'= '+str(nstfout)+' ; automatically modified by Swarm-CG' else: mdp_lines_in.append(f'nstfout = {nstfout} ; automatically added by Swarm-CG') @@ -2664,7 +2661,7 @@ def eval_function(parameters_set, ns): ns.total_model_eval_time += datetime.now().timestamp() - start_model_eval_ts # if gmx sasa failed to compute, it's most likely because there were inconsistent shifts across PBC in the trajectory = failed run - if ns.sasa_cg != None: + if ns.sasa_cg is not None: # store the distributions for each evaluation step shutil.move('distributions.png', '../'+config.distrib_plots_all_evals_dirname+'/distributions_eval_step_'+str(ns.nb_eval)+'.png') diff --git a/tests/data/cg_model.itp b/tests/data/cg_model.itp index e5bfcd4..e33ba59 100644 --- a/tests/data/cg_model.itp +++ b/tests/data/cg_model.itp @@ -52,8 +52,6 @@ G1 1 29 vT 1 G1 V29 29 0.00000 30 vT 1 G1 V30 30 0.00000 31 vT 1 G1 V31 31 0.00000 -32 vT 1 G1 V32 32 0.00000 -33 vT 1 G1 V33 33 0.00000 [ bonds ] @@ -93,25 +91,10 @@ G1 1 25 26 1 0 0 ; B4 ; bond group 5 - 1 28 1 0 0 ; B5 + 2 28 1 0 0 ; B5 -; bond group 6 - 2 28 1 0 0 ; B6 - -; bond group 7 - 16 32 1 0 0 ; B7 - -; bond group 8 - 3 29 1 0 0 ; B8 - -; bond group 9 - 7 27 1 0 0 ; B9 - -; bond group 10 - 1 33 1 0 0 ; B10 - -; bond group 11 - 2 33 1 0 0 ; B11 +; bond group 5 + 2 29 1 0 0 ; B6 [ angles ] @@ -157,19 +140,6 @@ G1 1 22 23 24 2 180 0 ; A5 22 25 26 2 180 0 ; A5 -; angle group 6 - 1 3 29 2 0 0 ; A6 - - -[ dihedrals ] -; i j k l funct dihedral force.c. mult. - -; dihedral group 1 - 3 1 2 21 1 0 0 2 ; D1 - 15 2 1 9 1 0 0 2 ; D1 - 21 2 1 3 1 0 0 2 ; D1 - 9 1 2 15 1 0 0 2 ; D1 - [ virtual_sites2 ] ; vs i j func param @@ -178,14 +148,16 @@ G1 1 [ virtual_sites3 ] ; vs i j k func params - 28 1 3 9 1 -0.3 -0.7 - 33 1 3 9 2 0.2 0.35 - 29 3 1 9 3 120 0.25 + 28 1 3 9 1 0.5 0.5 + 29 1 3 9 1 -0.6 -0.4 + +; 29 1 3 9 2 0.6 0.4 +; 30 1 3 9 1 -0.5 -0.5 +; 31 1 3 9 3 120 0.5 [ virtual_sitesn ] ; vs func def - 32 1 16 20 - 30 3 23 0.1 24 0.1 25 0.3 26 0.3 + 30 1 23 24 25 26 31 2 17 18 19 20 diff --git a/tests/data/system.top b/tests/data/system.top index 74c8bc0..9bd4bf1 100644 --- a/tests/data/system.top +++ b/tests/data/system.top @@ -10,5 +10,5 @@ G1 (PAMAM) in water [ molecules ] ; Compound #mols G1 1 -W 1241 +W 1243 CL- 8 \ No newline at end of file