Gromacs games


Additional energy terms in custom Hamiltonians, exchanges

The standard energy function in looks something like this:

E = E_bonds + E_angles + E_dihedrals + E_LJ + E_Culombic

This energy term, at some temperature T then becomes:

U = exp(-E/kb*T)

Let's say one has a n such Hamiltonians, U1...Un. These can then exchange, if the energies overlap. The details of the exchange scheme are skipped here, for gromacs implementation if the exchange see src/kernel/repl_ex.c

In temperature replica exchange one uses a ladder of T values, hoping that U1...Un overlap. In solute tampering replica exchange, one modifies the energy function E, keeping the T constant, hoping again to see exchange.

In deciding if a pair of replicas on a ladder will exchange, one compares their potential energy (Epot).

When including extra energy terms into the hamiltonian (custom biases, for example) one has to make sure they will be included in the Epot.

For instance, let's see if the gromacs pull code adds its energy to epot (such that it can be exchanged)

pull_potential_wrapper defined in src/mdlib/sim_util.c will call the pull_potential function and and save the energy to an enerd object:

enerd->term[F_COM_PULL] += pull_potential(...);

pull_potential_wrapper is called both in do_force_cutsGROUP and do_force_cutsVERLET (all inside of src/mdlib/sim_util.c).

Immediately after computing the pull potential, the code calls sum_epot function, let's have a look at that:

void sum_epot(t_grpopts *opts, gmx_grppairener_t *grpp, real *epot)
{
    int i;

    /* Accumulate energies */
    epot[F_COUL_SR]  = sum_v(grpp->nener, grpp->ener[egCOULSR]);
    epot[F_LJ]       = sum_v(grpp->nener, grpp->ener[egLJSR]);
    epot[F_LJ14]     = sum_v(grpp->nener, grpp->ener[egLJ14]);
    epot[F_COUL14]   = sum_v(grpp->nener, grpp->ener[egCOUL14]);
    epot[F_COUL_LR]  = sum_v(grpp->nener, grpp->ener[egCOULLR]);
    epot[F_LJ_LR]    = sum_v(grpp->nener, grpp->ener[egLJLR]);
    /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
    epot[F_GBPOL]   += sum_v(grpp->nener, grpp->ener[egGB]);

/* lattice part of LR doesnt belong to any group
 * and has been added earlier
 */
    epot[F_BHAM]     = sum_v(grpp->nener, grpp->ener[egBHAMSR]);
    epot[F_BHAM_LR]  = sum_v(grpp->nener, grpp->ener[egBHAMLR]);

    /* Sum all the values in the epot (list of floats) to a value at F_EPOT index */
    epot[F_EPOT] = 0;
    for (i = 0; (i < F_EPOT); i++)
    {
        if (i != F_DISRESVIOL && i != F_ORIRESDEV)
        {
            epot[F_EPOT] += epot[i];
        }
    }
}

In the first half of this function call, nothing useful happens. The second half where epot is iterated over (epot is the same as enerd->term earlier on), and accumulated.

Because F_COM_PULL comes before F_EPOT in include/types/idef.h, it will be included in the epot summation.

Multiple input files

This is a minor, utility hack. Normally, mdrun uses a single input file - for example:

mdrun -s topol.tpr -o traj.xtc

For replica exchange simulations however, these become mulitple files - the command

mdrun -multi 2 -s topol.tpr -o traj.xtc


Will initialize from 2 .tpr files topol0.tpr and topol1.tpr, it will also write two traj output files.

How to add additional files that will be read in with replica exchange simulations? Modify the following file

# src/gmxlib/main.c
     else if (bParFn)
    {
        /* Patch output and tpx, cpt and rerun input file names */
        for (i = 0; (i < nfile); i++)
        {
            /* Because of possible multiple extensions per type we must look
             * at the actual file name
             */
            if (is_output(&fnm[i]) ||
                fnm[i].ftp == efTPX || fnm[i].ftp == efCPT ||
                strcmp(fnm[i].opt, "-rerun") == 0)
            {
                ftp = fn2ftp(fnm[i].fns[0]);
                par_fn(fnm[i].fns[0], ftp, cr, TRUE, FALSE, buf, 255);
                sfree(fnm[i].fns[0]);
                fnm[i].fns[0] = gmx_strdup(buf);
            }
        }
    }

Which should become:

fnm[i].ftp == efTPX || fnm[i].ftp == efCPT || fnm[i].ftp == efCUSTOM ||

Comments

Post a Comment

Popular posts from this blog

React.js – edit and delete comments in CommentBox

It appears customary for high-intelligence, high-skills groups of engineers to "roll their own" whenever they have a chance to. Counter-intuitively to most managerial/MBA types out there, re-inventing the wheel can give you... a better wheel. Success is not guaranteed though – but these people may leave sooner rather than later if you don't let them... Having looked at number of JS frameworks (ember, backbone, angular) there was no clear winner. They were all a mess, at least to me. So I approached React.js with no hopes – how much better could it possible be? It turns out, a lot! This post will be an extension of the official React.js tutorial . The tutorial stops at the stage where you can add comments to a dynamic lists, that will grow as you add comments. This falls short of the functionality displayed on facebook: once you add a comment you can edit or delete it, so that's what're going to do. Full disclosure, I only started using React this weekendo so m...
Read more

Example slurm cluster on your laptop (multiple VMs via vagrant)

To anybody who's considering a switch to modern queuing systems, like slurm, this will be an useful guide. Rather than doing a roll out on production nodes, we'll try things out in vagrant. This text assumes you have vagrant installed on your local machine (it's very easy). The ultimate aim will be to setup 2 servers (worker nodes) and one master/controller, and run a simple job on them. There are multiple great guides out there: most of them are out of date, and won't work without some tweaking on Ubuntu 16.04 Xenial. The entirety of the code is hosted at https://github.com/jandom/gromacs-slurm-openmpi-vagrant The README.md file will show you how to set it up, with checks on every step. References https://mussolblog.wordpress.com/2013/07/17/setting-up-a-testing-slurm-cluster/ https://github.com/gabrieleiannetti/slurm_cluster_wiki/wiki/Installing-a-Slurm-Cluster http://philipwfowler.me/2016/04/14/how-to-setup-a-gramble/ https://github.com/dakl/slurm-clus...
Read more

Wrapping openbabel in python - using cython

The previous post showed a contrived example of how one could access the openbabel functionality in python using the boost libraries. There is alternative to boost to wrap around c/c++ code, it's called cython. Here is an example openbabel wrapper, exposing some very simple functionality. Just before diving in, you should know that there already exists a wrapper around openbabel in python (it is called pybel). What I'm showing here doesn't even come close to pybel in terms of usefulness, the idea here is to show a demo of how easy wrapping things in python is. #! /usr/bin/python from ez_setup import use_setuptools use_setuptools() from setuptools import setup, Extension, find_packages from Cython.Distutils import build_ext import sys, os import glob import subprocess as sub def get_include(name="openbabel"): p = sub.Popen('locate %s' % name ,stdout=sub.PIPE,stderr=sub.PIPE, shell=True) output, errors = p.communicate() ...
Read more