From 4c45a16cee102b5f4e22b33b8c9fb9c9c5c385eb Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Tue, 14 Dec 2021 08:26:46 -0500 Subject: [PATCH] Combined fixes and updates to gen_multimap.py Allow different deformation amplitude between leaflets, small fixes Fix reading input from --deformation-map-mask option Issue identified by @1004parky, mask map was initialized to 1 in all cases. Silence NumPy warnings Small fixes Fix handling of PDB files in generated VMD input Make single-map CVs optional Allow customizing selection for COM and orientation restraints Make generated NAMD and VMD scripts more uniform Unify writing of Colvars input between NAMD and VMD Add support for mapFile keyword --- colvartools/gen_multimap.py | 382 ++++++++++++++++++++++-------------- 1 file changed, 234 insertions(+), 148 deletions(-) diff --git a/colvartools/gen_multimap.py b/colvartools/gen_multimap.py index ba5050f70..7aa820fad 100644 --- a/colvartools/gen_multimap.py +++ b/colvartools/gen_multimap.py @@ -5,7 +5,7 @@ # Used in this tutorial: https://colvars.github.io/multi-map/multi-map.html # Download link: https://github.com/Colvars/colvars/blob/master/colvartools/gen_multimap.py?raw=true -# Last updated: 2020-10-28 + # Contact: giacomo.fiorin@gmail.com @@ -16,14 +16,13 @@ import scipy.ndimage as spnd from scipy.interpolate import InterpolatedUnivariateSpline - try: import gridData except ImportError: - print("Module gridData not found: " - "it is usually bundled with MDAnalysis.") - if (__name__ == '__main__'): - sys.exit(1) + raise ImportError("Error: module gridData not found, please install MDAnalysis.") + + +__version__ = '2024-10-24' def format_array(a, fmt="%.3f"): @@ -257,22 +256,22 @@ def get_cell_dimensions(args): return cell -def create_grid_object(cell, args): +def create_grid_object(cell, npoints): """ Create a new 3D map using gridData + + Parameters + ---------- + cell : array_like + Array with 3 dimensions + + npoints : array_like + Number of points along each of the three dimensions """ - # Set up 3D map cell_origin = -0.5*cell - cell_dx = np.zeros(shape=(3)) - cell_n = np.zeros(shape=(3), dtype=int) - if not args.grid_npoints is None: - cell_dx = np.divide(cell, np.array(args.grid_npoints)) - cell_n = np.array(args.grid_npoints, dtype=int) - else: - raise Exception("Undefined number of grid points.") - + cell_dx = np.divide(cell, np.array(npoints)) + cell_n = np.array(npoints, dtype=int) cell_origin += 0.5*cell_dx - map_obj = gridData.Grid(np.ones(shape=cell_n, dtype=float), origin=cell_origin, delta=cell_dx) @@ -290,7 +289,7 @@ def init_z_profiles(origin, length, dz, args): for input_file in args.input_profiles: # Load 1D profiles from file z_profile = np.loadtxt(input_file, - dtype=[('z', np.float), ('p', np.float)]) + dtype=[('z', np.float64), ('p', np.float64)]) z_profiles += [z_profile] elif (len(args.input_centers)): @@ -483,7 +482,8 @@ def get_vmd_load_map_cmds(map_file): def multimap_colvar_def_tcl(name, map_labels, coefficients=[], map_norms=[], indices=None, pdb_files=None, use_pdb_weights=False, - scripted_function=None, width=None): + map_files=None, + scripted_function=None, width=None, pdb_file_cache=[]): """ Write a Tcl script to define the collective variable in VMD or NAMD. @@ -513,6 +513,9 @@ def multimap_colvar_def_tcl(name, map_labels, coefficients=[], map_norms=[], use_pdb_weights : bool If True, set the atomWeights keyword using PDB columns + map_files : list of strings + If given, Colvars will load each maps internally using mapFile + scripted_function : string If defined, use it as the value of the scriptedFunction keyword @@ -520,7 +523,12 @@ def multimap_colvar_def_tcl(name, map_labels, coefficients=[], map_norms=[], If defined, set the width parameter of the colvar (the string may represent either a number or a Tcl variable expansion) + pdb_file_cache : list + If define, keep track of what PDB files are being loaded for atom + selections and reuse each file + """ + if len(coefficients) == 0 and scripted_function is None: raise Exception('Need either a set of coefficients or a scripting ' 'function.') @@ -566,22 +574,38 @@ def multimap_colvar_def_tcl(name, map_labels, coefficients=[], map_norms=[], componentCoeff %13.6e """ % (label, coefficients[i] * scale) - if indices is None: - conf += """\ + if map_files is None: + if indices is None: + conf += """\ mapName %s """ % label - else: - conf += """\ + else: + conf += """\ mapID %d """ % indices[i] - - if not pdb_files is None: + else: conf += """\ + mapFile %s +""" % map_files[i] + + + if (indices is not None) or args.load_maps_internally: + if not pdb_files[i] in pdb_file_cache: + conf += """\ atoms { + name %s atomsFile %s atomsCol O } +""" % (pdb_files[i], pdb_files[i]) + pdb_file_cache += [pdb_files[i]] + else: + conf += """\ + atoms { + atomsOfGroup %s + } """ % pdb_files[i] + if use_pdb_weights: conf += """\ atomWeights $weights(%s) @@ -599,10 +623,11 @@ def multimap_colvar_def_tcl(name, map_labels, coefficients=[], map_norms=[], def singlemap_colvars_def_tcl(map_labels, coefficients=[], map_norms=[], - indices=None, pdb_files=None, - use_pdb_weights=False): + indices=None, pdb_files=None, use_pdb_weights=False, + map_files=None, pdb_file_cache=[]): + conf = """ -# Define single-map variables for diagnostics (no extra computational cost) +# Define single-map variables for diagnostics """ for i, label in enumerate(map_labels): @@ -618,29 +643,43 @@ def singlemap_colvars_def_tcl(map_labels, coefficients=[], map_norms=[], mapTotal { name %s """ % (label, label) - if indices is None: - # Map labels for NAMD - conf += """\ - mapName %s + + if map_files is None: + if indices is None: + conf += """\ + mapName %s """ % label + else: + conf += """\ + mapID %d +""" % indices[i] else: - # Numeric IDs for VMD conf += """\ - mapID %d -""" % indices[i] + mapFile %s +""" % map_files[i] if scale != 1.0: conf += """\ componentCoeff %13.6e """ % scale - if not pdb_files is None: - conf += """\ + if (pdb_files is not None) or args.load_maps_internally: + if not pdb_files[i] in pdb_file_cache: + conf += """\ atoms { + name %s atomsFile %s atomsCol O } +""" % (pdb_files[i], pdb_files[i]) + pdb_file_cache += [pdb_files[i]] + else: + conf += """\ + atoms { + atomsOfGroup %s + } """ % pdb_files[i] + if use_pdb_weights: conf += """\ atomWeights $weights(%s) @@ -656,6 +695,7 @@ def singlemap_colvars_def_tcl(map_labels, coefficients=[], map_norms=[], def namd_com_restraint_def(pdb_file, name='com_dist', force_constant=5.0): + conf = """ # Define a center-of-mass restraint on all requested atoms @@ -684,11 +724,12 @@ def namd_com_restraint_def(pdb_file, name='com_dist', force_constant=5.0): def namd_ori_restraint_def(pdb_file, name='ori', force_constant=5.0): + conf = """ # Define an orientation restraint on requested atoms # WARNING: unlike the center of mass, this computation is not parallelized, -# and reducing the size of the selection (e.g. only C-alpha atoms for +# and reducing the size of the selection (e.g. only C-alpha atoms for # proteins) is highly recommended to maintain reasonable parallel performance. cv config \" colvar { @@ -714,8 +755,8 @@ def namd_ori_restraint_def(pdb_file, name='ori', force_constant=5.0): return conf -def namd_com_z_restraint_def(pdb_files, name='com_dist', force_constant=5.0): - unique_files = list(set(pdb_files)) +def namd_com_z_restraint_def(pdb_file, name='com_dist', force_constant=5.0): + conf = """ # Define a center-of-mass restraint on all requested atoms @@ -723,22 +764,17 @@ def namd_com_z_restraint_def(pdb_files, name='com_dist', force_constant=5.0): colvar { name %s width 0.1 -""" % name - - for pdb_file in unique_files: - conf += """\ distanceZ { - componentCoeff %f main { atomsFile %s atomsCol O } ref { dummyAtom (0.0, 0.0, 0.0) } } -""" % (1.0/len(unique_files), pdb_file) +} +""" % (name, pdb_file) conf += """\ -} harmonic { name r_%s @@ -902,7 +938,7 @@ def parse_cmdline_args(namespace=None): 'in the same format as --deformation-map. ' 'May be used together with other masking functions.') group.add_argument('--mask-by-amplitude', - type=np.float, + type=float, help='Mask out regions that move less than this ' 'fraction of the maximum amplitude. ' 'May be used together with other masking functions.', @@ -979,17 +1015,31 @@ def parse_cmdline_args(namespace=None): 'PDB columns to define the weights of each atom within ' 'each map; when false, all weights are equal to 1. ' 'Defaults to True if --input-pdb-files is used.', - default=None) + default=False) + group.add_argument('--load-maps-internally', + action='store_true', + help="When true, Colvars will load the maps internally using " + "\"mapFile\" and the loop over atoms is done inside Colvars " + "instead of NAMD; " + "does not have an effect in VMD", + default=False) + group.add_argument('--define-single-maps', + action='store_true', + help='Define single-map variables in addition to Multi-Map', + default=False) group.add_argument('--com-restraint', action='store_true', help='Add the definition of a center-of-mass restraint ' 'to the NAMD input', default=True) - + group.add_argument('--com-restraint-pdb-file', + type=str, + help="PDB file-based selection for the COM restraint") group.add_argument('--ori-restraint', action='store_true', help='Add the definition of an orientational restraint ' - 'to the NAMD input; it is highly recommended to reduce ' + 'to the NAMD input; uses the same file as --com-restraint to select atoms; ' + 'it is highly recommended to reduce ' 'the size of the selection (e.g. only C-alpha atoms ' 'for proteins) to maintain reasonable parallel ' 'performance. Not valid for bilayers.', @@ -1000,11 +1050,11 @@ def parse_cmdline_args(namespace=None): type=str, required=True, help='Prefix for output files') group.add_argument('--namd-script', - type=str, default=None, + type=str, help='Write NAMD commands to define the Multi-Map ' 'collective variable into this file.') group.add_argument('--vmd-script', - type=str, default=None, + type=str, help='Write VMD commands to define the Multi-map ' 'collective variable into this file.') @@ -1091,8 +1141,9 @@ def check_prep_args(args): "\"; values of this 2D map will be multiplied by " "the deformation amplitudes. " "The mean is automatically subtracted.") - args.deformation_map = np.loadtxt(args.deformation_map, dtype=np.float) - args.deformation_map_mask = np.ones(shape=args.deformation_map.shape) + args.deformation_map = np.loadtxt(args.deformation_map, dtype=np.float64) + if args.deformation_map_mask is None: + args.deformation_map_mask = np.ones(shape=args.deformation_map.shape) print("Subtracting the mean Z-coordinate of the map; before =", args.deformation_map.mean()) args.deformation_map -= args.deformation_map.mean() @@ -1106,7 +1157,7 @@ def check_prep_args(args): if type(args.deformation_map_mask) == str: print("Loading the mask from file \""+args.deformation_map_mask+"\"") args.deformation_map_mask = np.loadtxt(args.deformation_map_mask, - dtype=np.float) + dtype=np.float64) print("Mean value of the mask map =", args.deformation_map_mask.mean()) if args.mask_by_amplitude > 0.0: @@ -1115,7 +1166,7 @@ def check_prep_args(args): print("Masking the deformation map based on its magnitude") max_amplitude = np.abs(args.deformation_map).max() mask_threshold = args.mask_by_amplitude*max_amplitude - args.deformation_map_mask *= \ + args.deformation_map_mask = \ get_mask_from_deformation(deformation_map=args.deformation_map, threshold=mask_threshold, tolerance=0.2*mask_threshold) @@ -1123,8 +1174,13 @@ def check_prep_args(args): del mask_threshold, max_amplitude - args.n_deformations = max(len(args.input_map_files), - len(args.deformation_amplitudes)) + if args.system_dim == '3d': + args.n_deformations = len(args.input_map_files) + if args.system_dim == '2d': + args.n_deformations = len(args.deformation_amplitudes) // args.n_inputs + args.deformation_amplitudes = \ + np.array(args.deformation_amplitudes).reshape((args.n_inputs, + args.n_deformations)) if args.multimap_coefficients is None: args.multimap_coefficients = np.arange(1.0, args.n_deformations+1) @@ -1147,7 +1203,7 @@ def generate_maps_from_profiles(args): # TODO forbid triclinic cells # Object to hold the undeformed map for each input - map_input = create_grid_object(cell_sizes, args) + map_input = create_grid_object(cell_sizes, args.grid_npoints) print("Unit cell dimensions = ", format_array(cell_sizes)) print("Grid sizes (# points) = ", map_input.grid.shape) # Define cell_dx in a consistent manner, because the shape of Grid.delta @@ -1165,8 +1221,13 @@ def generate_maps_from_profiles(args): maps = {} # Loop over all input maps (e.g. upper and lower leaflet) - for z_profile, z_profile_com, input_label in \ - zip(z_profiles, z_profiles_coms, args.input_labels): + for i_input, z_profile, z_profile_com, input_label in \ + zip(range(args.n_inputs), z_profiles, z_profiles_coms, + args.input_labels): + + deformation_amplitudes = args.deformation_amplitudes[i_input,:] + + print("Amplitudes", deformation_amplitudes) maps[input_label] = [] @@ -1175,10 +1236,10 @@ def generate_maps_from_profiles(args): init_map_from_profile(z_profile=z_profile['p'], grid=map_input.grid) - for i, amplitude in enumerate(args.deformation_amplitudes): + for i, amplitude in enumerate(deformation_amplitudes): # Object to hold the deformed map - map_output = create_grid_object(cell_sizes, args) + map_output = create_grid_object(cell_sizes, args.grid_npoints) print(" Deformation amplitude =", format_array(amplitude)) amplitude_bins = np.divide(amplitude, cell_dx[2][2]) @@ -1205,7 +1266,7 @@ def generate_maps_from_profiles(args): if not args.deformation_map_mask is None: factor = args.deformation_map_mask.sum() / \ - np.float(args.deformation_map_mask.size) + np.float64(args.deformation_map_mask.size) print(" Masking the map (fraction = %.1f%%)" % (factor*100)) np.multiply(map_output.grid, args.deformation_map_mask[:,:,np.newaxis], @@ -1217,16 +1278,21 @@ def generate_maps_from_profiles(args): return maps -def write_namd_script(gridforces_script, map_labels, pdb_files, map_norms, - args): - with open(args.namd_script, 'w') as namd_script: +def write_namd_script(gridforces_script, map_labels, pdb_files, map_norms, map_files, args): + print("Writing NAMD script to file", args.namd_script) + pdb_file_cache = [] + with open(args.namd_script, 'w') as script: - print("Writing NAMD script to file", args.namd_script) + script.write("""# -*- tcl -*- - namd_script.write("""# -*- tcl -*- +## Use this script by setting "colvars on" and then "source %s" +""" % args.namd_script) + if not args.load_maps_internally: + script.write(""" -## Please ensure that "mGridForce on" and "colvars on" are defined before -## sourcing this file in a NAMD script. +## Please also ensure that "mGridForce on" is specified before sourcing this file in a NAMD script. +""") + script.write(""" # The width parameters sets the amplitude of the colvar's fluctuation, and # can be used to set the correct scale for harmonic restraint potentials @@ -1236,7 +1302,7 @@ def write_namd_script(gridforces_script, map_labels, pdb_files, map_norms, """) if args.system_dim == '3d': - namd_script.write(""" + script.write(""" # Define a default position for the center-of-mass restraint if { [info exists com_pos] == 0 } { set com_pos [list 0.0 0.0 0.0] @@ -1247,49 +1313,75 @@ def write_namd_script(gridforces_script, map_labels, pdb_files, map_norms, set com("z") [lindex ${com_pos} 2] """) - namd_script.write(""" -# Load gridForces maps -""") - namd_script.write(gridforces_script) + if not args.load_maps_internally: + script.write(""" + # Load gridForces maps + """) + script.write(gridforces_script) + write_colvars_script(script, map_labels, pdb_files, map_norms, map_files, pdb_file_cache, args) + + +def write_colvars_script(script, map_labels, pdb_files, map_norms, map_files, pdb_file_cache, args, indices=None): + if script is not None: + + # Disabled code path, using componentCoeff now rather than Tcl scripted_function = None + if not scripted_function is None: - namd_script.write(multimap_colvar_function_tcl(args)) - - mmcv_def = \ - multimap_colvar_def_tcl(name=args.colvar_name, - map_labels=map_labels, - map_norms=map_norms, - coefficients=\ - np.tile(args.multimap_coefficients, - args.n_inputs), - scripted_function=scripted_function, - use_pdb_weights=args.use_pdb_weights, - width='${multimap_cv_width}') - namd_script.write(mmcv_def) - - singles_def = \ - singlemap_colvars_def_tcl(map_labels=map_labels, - map_norms=map_norms, - use_pdb_weights=args.use_pdb_weights) - namd_script.write(singles_def) + script.write(multimap_colvar_function_tcl(args)) + + mmcv_def = multimap_colvar_def_tcl( + name=args.colvar_name, + map_labels=map_labels, + map_norms=map_norms, + indices=indices, + coefficients=np.tile( + args.multimap_coefficients, + args.n_inputs), + pdb_files=pdb_files, + use_pdb_weights=args.use_pdb_weights, + map_files=map_files, + scripted_function=scripted_function, + width='${multimap_cv_width}', + pdb_file_cache=pdb_file_cache) + + script.write(mmcv_def) + + if args.define_single_maps: + singles_def = singlemap_colvars_def_tcl( + map_labels=map_labels, + indices=indices, + map_norms=map_norms, + pdb_files=pdb_files, + use_pdb_weights=args.use_pdb_weights, + map_files=map_files, + pdb_file_cache=pdb_file_cache) + script.write(singles_def) + + # Center-of-mass and orientation restraint + + com_pdb_file = ori_pdb_file = args.com_restraint_pdb_file + if com_pdb_file is None: + com_pdb_file = ori_pdb_file = pdb_files[0] if args.system_dim == '3d': - # Center-of-mass restraint - namd_script.write(namd_com_restraint_def(pdb_files[0])) + if args.com_restraint: + script.write(namd_com_restraint_def(com_pdb_file)) if args.ori_restraint: - namd_script.write(namd_ori_restraint_def(pdb_files[0])) + script.write(namd_ori_restraint_def(ori_pdb_file)) - if args.system_dim == '2d': - namd_script.write(namd_com_z_restraint_def(pdb_files)) + if args.system_dim == '2d' and args.com_restraint: + script.write(namd_com_z_restraint_def(com_pdb_file)) -def write_vmd_script(vmd_load_map_cmds, map_labels, pdb_files, map_norms, args): +def write_vmd_script(vmd_load_map_cmds, map_labels, pdb_files, map_norms, map_files, args): + print("Writing VMD script to file", args.vmd_script) + pdb_file_cache = [] unique_files = list(set(pdb_files)) - with open(args.vmd_script, 'w') as vmd_script_file: - print("Writing VMD script to file", args.vmd_script) - vmd_script_file.write("""# -*- tcl -*- + with open(args.vmd_script, 'w') as script: + script.write("""# -*- tcl -*- ## VMD script to load map files onto the top molecule @@ -1301,7 +1393,7 @@ def write_vmd_script(vmd_load_map_cmds, map_labels, pdb_files, map_norms, args): """) if args.use_pdb_weights: - vmd_script_file.write(""" + script.write(""" # Precompute the per-atom weights by taking them from each PDB file foreach pdb_file { %s } { mol addfile $pdb_file type pdb waitfor all @@ -1312,43 +1404,22 @@ def write_vmd_script(vmd_load_map_cmds, map_labels, pdb_files, map_norms, args): } """ % ' '.join(set(unique_files))) - vmd_script_file.write(""" + script.write(""" # Attach Colvars to the top molecule cv molid top +""") + if not args.load_maps_internally: + script.write(""" # Load volumetric maps """) - vmd_script_file.write(vmd_load_map_cmds) - scripted_function = None - if not scripted_function is None: - vmd_script_file.write(multimap_colvar_function_tcl(args)) - - n = len(map_labels) - - mmcv_def = \ - multimap_colvar_def_tcl(name=args.colvar_name, - map_labels=map_labels, - map_norms=map_norms, - indices=range(n), - pdb_files=pdb_files, - coefficients=\ - np.tile(args.multimap_coefficients, - args.n_inputs), - scripted_function=scripted_function, - use_pdb_weights=args.use_pdb_weights, - width='${multimap_cv_width}') - vmd_script_file.write(mmcv_def) - - singles_def = \ - singlemap_colvars_def_tcl(map_labels=map_labels, - indices=range(n), - map_norms=map_norms, - pdb_files=pdb_files, - use_pdb_weights=args.use_pdb_weights) - vmd_script_file.write(singles_def) - - -def generate_multimap(args): + script.write(vmd_load_map_cmds) + + write_colvars_script(script, map_labels, pdb_files, map_norms, map_files, pdb_file_cache, args, + indices=range(len(map_labels))) + + +def gen_multimap(args): """ Process cmdline arguments, generate/load the maps, write input scripts """ @@ -1442,13 +1513,13 @@ def generate_multimap(args): else: map_file = map_files[input_label][i] - if not args.namd_script is None: + if not args.namd_script is None and not args.load_maps_internally: gridforces_script += \ get_gridforces_map_cmds(label=map_label, map_file=map_file, pdb_file=pdb_file) - if not args.vmd_script is None: + if not args.vmd_script is None and not args.load_maps_internally: vmd_load_map_cmds += \ get_vmd_load_map_cmds(map_file=map_file) @@ -1485,18 +1556,33 @@ def generate_multimap(args): all_pdb_files = [x for il in args.input_labels for x in pdb_files[il]] all_map_norms = [x for il in args.input_labels for x in map_norms[il]] - print("") + if args.load_maps_internally: + all_map_files = [x for il in args.input_labels for x in map_files[il]] + else: + all_map_files = None - if not args.namd_script is None: - write_namd_script(gridforces_script, all_map_labels, all_pdb_files, - all_map_norms, args) + print("") - if not args.vmd_script is None: - write_vmd_script(vmd_load_map_cmds, all_map_labels, all_pdb_files, - all_map_norms, args) + if args.namd_script: + write_namd_script( + gridforces_script, + all_map_labels, + all_pdb_files, + all_map_norms, + all_map_files, + args) + + if args.vmd_script: + write_vmd_script( + vmd_load_map_cmds, + all_map_labels, + all_pdb_files, + all_map_norms, + all_map_files, + args) if __name__ == '__main__': args = parse_cmdline_args() check_prep_args(args) - generate_multimap(args) + gen_multimap(args)