Skip to content

Commit

Permalink
Add support for mapFile keyword
Browse files Browse the repository at this point in the history
  • Loading branch information
@giacomofiorin
giacomofiorin committed Oct 26, 2024
1 parent 8ab255f commit bae0383
Showing 1 changed file with 90 additions and 46 deletions.
136 changes: 90 additions & 46 deletions colvartools/gen_multimap.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -482,6 +482,7 @@ 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,
map_files=None,
scripted_function=None, width=None, pdb_file_cache=[]):
"""
Write a Tcl script to define the collective variable in VMD or NAMD.
Expand Down Expand Up @@ -512,13 +513,20 @@ 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
width : string
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:
Expand Down Expand Up @@ -566,14 +574,20 @@ 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
""" % label
else:
conf += """\
mapID %d
""" % indices[i]
else:
conf += """\
mapID %d
""" % indices[i]
mapFile %s
""" % map_files[i]


if not pdb_files is None:
if not pdb_files[i] in pdb_file_cache:
Expand Down Expand Up @@ -609,11 +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, pdb_file_cache=[]):
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):
Expand All @@ -629,16 +643,20 @@ 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
""" % label

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 += """\
Expand Down Expand Up @@ -677,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
Expand Down Expand Up @@ -705,6 +724,7 @@ 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
Expand Down Expand Up @@ -736,6 +756,7 @@ def namd_ori_restraint_def(pdb_file, name='ori', force_constant=5.0):


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
Expand All @@ -751,7 +772,7 @@ def namd_com_z_restraint_def(pdb_file, name='com_dist', force_constant=5.0):
ref { dummyAtom (0.0, 0.0, 0.0) }
}
}
""" % name, pdb_file
""" % (name, pdb_file)

conf += """\
Expand Down Expand Up @@ -994,7 +1015,14 @@ 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',
Expand All @@ -1007,7 +1035,6 @@ def parse_cmdline_args(namespace=None):
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 '
Expand All @@ -1023,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.')

Expand Down Expand Up @@ -1251,12 +1278,11 @@ def generate_maps_from_profiles(args):
return maps


def write_namd_script(gridforces_script, map_labels, pdb_files, map_norms, args):
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.script)

script.write("""# -*- tcl -*-
## Please ensure that "mGridForce on" and "colvars on" are defined before
Expand Down Expand Up @@ -1286,12 +1312,11 @@ def write_namd_script(gridforces_script, map_labels, pdb_files, map_norms, args)
""")
script.write(gridforces_script)

write_colvars_script(script, map_labels, pdb_files, map_norms, args)
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, args):
if script:
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
Expand All @@ -1303,49 +1328,52 @@ def write_colvars_script(script, map_labels, pdb_files, map_norms, args):
name=args.colvar_name,
map_labels=map_labels,
map_norms=map_norms,
indices=None,
pdb_files=pdb_files,
indices=indices,
coefficients=np.tile(
args.multimap_coefficients,
args.n_inputs),
scripted_function=scripted_function,
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=None,
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]
com_pdb_file = ori_pdb_file = pdb_files[0]

if args.system_dim == '3d':
if args.com_restraint:
script.write(namd_com_restraint_def(pdb_files[0]))
script.write(namd_com_restraint_def(com_pdb_file))
if args.ori_restraint:
script.write(namd_ori_restraint_def(pdb_files[0]))
script.write(namd_ori_restraint_def(ori_pdb_file))

if args.system_dim == '2d' and args.com_restraint:
script.write(namd_com_z_restraint_def(pdb_files))
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 script:
print("Writing VMD script to file", args.script)
script.write("""# -*- tcl -*-
## VMD script to load map files onto the top molecule
Expand Down Expand Up @@ -1377,7 +1405,8 @@ def write_vmd_script(vmd_load_map_cmds, map_labels, pdb_files, map_norms, args):
""")
script.write(vmd_load_map_cmds)

write_colvars_script(script, map_labels, pdb_files, map_norms, args):
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):
Expand Down Expand Up @@ -1517,15 +1546,30 @@ def gen_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__':
Expand Down

0 comments on commit bae0383

Please sign in to comment.