gromacs/setup.py
def _setup_MD(dirname,
              deffnm='md', mdp=config.templates['md_OPLSAA.mdp'],
              struct=None,
              top='top/system.top', ndx=None,
              mainselection='"Protein"',
              qscript=config.qscript_template, qname=None, startdir=None, mdrun_opts="", budget=None, walltime=1/3.,
              dt=0.002, runtime=1e3, **mdp_kwargs):
    """Generic function to set up a ``mdrun`` MD simulation.

    See the user functions for usage.
    """

    if struct is None:
        raise ValueError('struct must be set to a input structure')
    structure = realpath(struct)
    topology = realpath(top)
    try:
        index = realpath(ndx)
    except AttributeError:  # (that's what realpath(None) throws...)
        index = None        # None is handled fine below

    qname = mdp_kwargs.pop('sgename', qname)    # compatibility for old scripts
    qscript = mdp_kwargs.pop('sge', qscript)    # compatibility for old scripts
    qscript_template = config.get_template(qscript)
    mdp_template = config.get_template(mdp)

    nsteps = int(float(runtime)/float(dt))

    mdp = deffnm + '.mdp'
    tpr = deffnm + '.tpr'
    mainindex = deffnm + '.ndx'
    final_structure = deffnm + '.gro'   # guess... really depends on templates,could also be DEFFNM.pdb

    # write the processed topology to the default output
    mdp_parameters = {'nsteps':nsteps, 'dt':dt, 'pp': 'processed.top'}
    mdp_parameters.update(mdp_kwargs)

    cbook.add_mdp_includes(topology, mdp_parameters)

    logger.info("[%(dirname)s] input mdp  = %(mdp_template)r", vars())
    with in_dir(dirname):
        if not (mdp_parameters.get('Tcoupl','').lower() == 'no' or mainselection is None):
            logger.info("[{dirname!s}] Automatic adjustment of T-coupling groups".format(**vars()))

            # make index file in almost all cases; with mainselection == None the user
            # takes FULL control and also has to provide the template or index
            groups = make_main_index(structure, selection=mainselection,
                                     oldndx=index, ndx=mainindex)
            natoms = {g['name']: float(g['natoms']) for g in groups}
            tc_group_names = ('__main__', '__environment__')   # defined in make_main_index()
            try:
                x = natoms['__main__']/natoms['__environment__']
            except KeyError:
                x = 0   # force using SYSTEM in code below
                wmsg = "Missing __main__ and/or __environment__ index group.\n" \
                       "This probably means that you have an atypical system. You can " \
                       "set mainselection=None and provide your own mdp and index files " \
                       "in order to set up temperature coupling.\n" \
                       "If no T-coupling is required then set Tcoupl='no'.\n" \
                       "For now we will just couple everything to 'System'."
                logger.warn(wmsg)
                warnings.warn(wmsg, category=AutoCorrectionWarning)
            if x < 0.1:
                # couple everything together
                tau_t = firstof(mdp_parameters.pop('tau_t', 0.1))
                ref_t = firstof(mdp_parameters.pop('ref_t', 300))
                # combine all in one T-coupling group
                mdp_parameters['tc-grps'] = 'System'
                mdp_parameters['tau_t'] = tau_t   # this overrides the commandline!
                mdp_parameters['ref_t'] = ref_t   # this overrides the commandline!
                mdp_parameters['gen-temp'] = mdp_parameters.pop('gen_temp', ref_t)
                wmsg = "Size of __main__ is only %.1f%% of __environment__ so " \
                       "we use 'System' for T-coupling and ref_t = %g K and " \
                       "tau_t = %g 1/ps (can be changed in mdp_parameters).\n" \
                       % (x * 100, ref_t, tau_t)
                logger.warn(wmsg)
                warnings.warn(wmsg, category=AutoCorrectionWarning)
            else:
                # couple protein and bath separately
                n_tc_groups = len(tc_group_names)
                tau_t = asiterable(mdp_parameters.pop('tau_t', 0.1))
                ref_t = asiterable(mdp_parameters.pop('ref_t', 300))

                if len(tau_t) != n_tc_groups:
                    tau_t = n_tc_groups * [tau_t[0]]
                    wmsg = "%d coupling constants should have been supplied for tau_t. "\
                        "Using %f 1/ps for all of them." % (n_tc_groups, tau_t[0])
                    logger.warn(wmsg)
                    warnings.warn(wmsg, category=AutoCorrectionWarning)
                if len(ref_t) != n_tc_groups:
                    ref_t = n_tc_groups * [ref_t[0]]
                    wmsg = "%d temperatures should have been supplied for ref_t. "\
                        "Using %g K for all of them." % (n_tc_groups, ref_t[0])
                    logger.warn(wmsg)
                    warnings.warn(wmsg, category=AutoCorrectionWarning)

                mdp_parameters['tc-grps'] = tc_group_names
                mdp_parameters['tau_t'] = tau_t
                mdp_parameters['ref_t'] = ref_t
                mdp_parameters['gen-temp'] = mdp_parameters.pop('gen_temp', ref_t[0])
            index = realpath(mainindex)
        if mdp_parameters.get('Tcoupl','').lower() == 'no':
            logger.info("Tcoupl == no: disabling all temperature coupling mdp options")
            mdp_parameters['tc-grps'] = ""
            mdp_parameters['tau_t'] = ""
            mdp_parameters['ref_t'] = ""
            mdp_parameters['gen-temp'] = ""
        if mdp_parameters.get('Pcoupl','').lower() == 'no':
            logger.info("Pcoupl == no: disabling all pressure coupling mdp options")
            mdp_parameters['tau_p'] = ""
            mdp_parameters['ref_p'] = ""
            mdp_parameters['compressibility'] = ""

        unprocessed = cbook.edit_mdp(mdp_template, new_mdp=mdp, **mdp_parameters)
        check_mdpargs(unprocessed)
        gromacs.grompp(f=mdp, p=topology, c=structure, n=index, o=tpr, **unprocessed)

        runscripts = qsub.generate_submit_scripts(
            qscript_template, deffnm=deffnm, jobname=qname, budget=budget,
            startdir=startdir, mdrun_opts=mdrun_opts, walltime=walltime)

    logger.info("[%(dirname)s] output mdp = %(mdp)r", vars())
    logger.info("[%(dirname)s] output ndx = %(ndx)r", vars())
    logger.info("[%(dirname)s] output tpr = %(tpr)r", vars())
    logger.info("[%(dirname)s] output runscripts = %(runscripts)r", vars())
    logger.info("[%(dirname)s] All files set up for a run time of %(runtime)g ps "
                "(dt=%(dt)g, nsteps=%(nsteps)g)" % vars())

    kwargs = {'struct': realpath(os.path.join(dirname, final_structure)),      # guess
              'top': topology,
              'ndx': index,            # possibly mainindex
              'qscript': runscripts,
              'mainselection': mainselection,
              'deffnm': deffnm,        # return deffnm (tpr = deffnm.tpr!)
              }
    kwargs.update(mdp_kwargs)  # return extra mdp args so that one can use them for prod run
    return kwargs
def _setup_MD(dirname,
              deffnm='md', mdp=config.templates['md_OPLSAA.mdp'],
              struct=None,
              top='top/system.top', ndx=None,
              mainselection='"Protein"',
              qscript=config.qscript_template, qname=None, startdir=None, mdrun_opts="", budget=None, walltime=1/3.,
              dt=0.002, runtime=1e3, **mdp_kwargs):
    """Generic function to set up a ``mdrun`` MD simulation.

    See the user functions for usage.
    """

    if struct is None:
        raise ValueError('struct must be set to a input structure')
    structure = realpath(struct)
    topology = realpath(top)
    try:
        index = realpath(ndx)
    except AttributeError:  # (that's what realpath(None) throws...)
        index = None        # None is handled fine below

    qname = mdp_kwargs.pop('sgename', qname)    # compatibility for old scripts
    qscript = mdp_kwargs.pop('sge', qscript)    # compatibility for old scripts
    qscript_template = config.get_template(qscript)
    mdp_template = config.get_template(mdp)

    nsteps = int(float(runtime)/float(dt))

    mdp = deffnm + '.mdp'
    tpr = deffnm + '.tpr'
    mainindex = deffnm + '.ndx'
    final_structure = deffnm + '.gro'   # guess... really depends on templates,could also be DEFFNM.pdb

    # write the processed topology to the default output
    mdp_parameters = {'nsteps':nsteps, 'dt':dt, 'pp': 'processed.top'}
    mdp_parameters.update(mdp_kwargs)

    cbook.add_mdp_includes(topology, mdp_parameters)

    logger.info("[%(dirname)s] input mdp  = %(mdp_template)r", vars())
    with in_dir(dirname):
        if not (mdp_parameters.get('Tcoupl','').lower() == 'no' or mainselection is None):
            logger.info("[{dirname!s}] Automatic adjustment of T-coupling groups".format(**vars()))

            # make index file in almost all cases; with mainselection == None the user
            # takes FULL control and also has to provide the template or index
            groups = make_main_index(structure, selection=mainselection,
                                     oldndx=index, ndx=mainindex)
            natoms = {g['name']: float(g['natoms']) for g in groups}
            tc_group_names = ('__main__', '__environment__')   # defined in make_main_index()
            try:
                x = natoms['__main__']/natoms['__environment__']
            except KeyError:
                x = 0   # force using SYSTEM in code below
                wmsg = "Missing __main__ and/or __environment__ index group.\n" \
                       "This probably means that you have an atypical system. You can " \
                       "set mainselection=None and provide your own mdp and index files " \
                       "in order to set up temperature coupling.\n" \
                       "If no T-coupling is required then set Tcoupl='no'.\n" \
                       "For now we will just couple everything to 'System'."
                logger.warn(wmsg)
                warnings.warn(wmsg, category=AutoCorrectionWarning)
            if x < 0.1:
                # couple everything together
                tau_t = firstof(mdp_parameters.pop('tau_t', 0.1))
                ref_t = firstof(mdp_parameters.pop('ref_t', 300))
                # combine all in one T-coupling group
                mdp_parameters['tc-grps'] = 'System'
                mdp_parameters['tau_t'] = tau_t   # this overrides the commandline!
                mdp_parameters['ref_t'] = ref_t   # this overrides the commandline!
                mdp_parameters['gen-temp'] = mdp_parameters.pop('gen_temp', ref_t)
                wmsg = "Size of __main__ is only %.1f%% of __environment__ so " \
                       "we use 'System' for T-coupling and ref_t = %g K and " \
                       "tau_t = %g 1/ps (can be changed in mdp_parameters).\n" \
                       % (x * 100, ref_t, tau_t)
                logger.warn(wmsg)
                warnings.warn(wmsg, category=AutoCorrectionWarning)
            else:
                # couple protein and bath separately
                n_tc_groups = len(tc_group_names)
                tau_t = asiterable(mdp_parameters.pop('tau_t', 0.1))
                ref_t = asiterable(mdp_parameters.pop('ref_t', 300))

                if len(tau_t) != n_tc_groups:
                    tau_t = n_tc_groups * [tau_t[0]]
                    wmsg = "%d coupling constants should have been supplied for tau_t. "\
                        "Using %f 1/ps for all of them." % (n_tc_groups, tau_t[0])
                    logger.warn(wmsg)
                    warnings.warn(wmsg, category=AutoCorrectionWarning)
                if len(ref_t) != n_tc_groups:
                    ref_t = n_tc_groups * [ref_t[0]]
                    wmsg = "%d temperatures should have been supplied for ref_t. "\
                        "Using %g K for all of them." % (n_tc_groups, ref_t[0])
                    logger.warn(wmsg)
                    warnings.warn(wmsg, category=AutoCorrectionWarning)

                mdp_parameters['tc-grps'] = tc_group_names
                mdp_parameters['tau_t'] = tau_t
                mdp_parameters['ref_t'] = ref_t
                mdp_parameters['gen-temp'] = mdp_parameters.pop('gen_temp', ref_t[0])
            index = realpath(mainindex)
        if mdp_parameters.get('Tcoupl','').lower() == 'no':
            logger.info("Tcoupl == no: disabling all temperature coupling mdp options")
            mdp_parameters['tc-grps'] = ""
            mdp_parameters['tau_t'] = ""
            mdp_parameters['ref_t'] = ""
            mdp_parameters['gen-temp'] = ""
        if mdp_parameters.get('Pcoupl','').lower() == 'no':
            logger.info("Pcoupl == no: disabling all pressure coupling mdp options")
            mdp_parameters['tau_p'] = ""
            mdp_parameters['ref_p'] = ""
            mdp_parameters['compressibility'] = ""

        unprocessed = cbook.edit_mdp(mdp_template, new_mdp=mdp, **mdp_parameters)
        check_mdpargs(unprocessed)
        gromacs.grompp(f=mdp, p=topology, c=structure, n=index, o=tpr, **unprocessed)

        runscripts = qsub.generate_submit_scripts(
            qscript_template, deffnm=deffnm, jobname=qname, budget=budget,
            startdir=startdir, mdrun_opts=mdrun_opts, walltime=walltime)

    logger.info("[%(dirname)s] output mdp = %(mdp)r", vars())
    logger.info("[%(dirname)s] output ndx = %(ndx)r", vars())
    logger.info("[%(dirname)s] output tpr = %(tpr)r", vars())
    logger.info("[%(dirname)s] output runscripts = %(runscripts)r", vars())
    logger.info("[%(dirname)s] All files set up for a run time of %(runtime)g ps "
                "(dt=%(dt)g, nsteps=%(nsteps)g)" % vars())

    kwargs = {'struct': realpath(os.path.join(dirname, final_structure)),      # guess
              'top': topology,
              'ndx': index,            # possibly mainindex
              'qscript': runscripts,
              'mainselection': mainselection,
              'deffnm': deffnm,        # return deffnm (tpr = deffnm.tpr!)
              }
    kwargs.update(mdp_kwargs)  # return extra mdp args so that one can use them for prod run
    return kwargs
Generic function to set up a ``mdrun`` MD simulation.

    See the user functions for usage.