Gromacs

Patching mdrun

    Table of contents
    No headers

    Editing the mdrun program can be a daunting task, and one that scares off many potential developers. Here we treat an example of a simple extension to mdrun that implements a complicated electric field, which could be useful for instance in membrane simulations. For the example we will consider a box of water (you can use the one from the tutor directory) and add an electric field that has a triangular shape.

    The field is then given by the following pseudo-code:

     if (x < box/2) 
       E = 2*x/box;
     else
       E = 1-2*x/box;
    

    Now, go to the gmx/src/mdlib directory and open the file sim_util.c in a text editor. Scroll down until you find the function calc_f_el which already calculatues forces due to an electric field, but with different shape. Since this is a hacking exercise, we will relentlessly remove the code within the routine but not the variables. Just delete the code, and replace it by:

    for(i=start; (i<start+homenr); i++) {
      if (x[i][ZZ] < box[ZZ][ZZ]/2) {
       Ext[ZZ] = FIELDFAC*Ex[ZZ].a[0]*2*x[i][ZZ]/box[ZZ][ZZ];
      }
      else {
       Ext[ZZ] = FIELDFAC*Ex[ZZ].a[0]*(1- 2*x[i][ZZ]/box[ZZ][ZZ]);
      }
      f[i][ZZ] += charge[i]*Ext[ZZ];
    } 
    

    Now, if you were to save this file and try to compile it, you would find that the box variable is not declared inside this routine, and we have to add it and pass it on from the calling routine. Go up to the function declaration and add, after real t a declaration for matrix box (put a comma before it!). Now, finally, go down to the function do_force and look up the call to the calc_f_el function, and, as the last argument in the function call, add the box variable. Now try compiling the code, in the mdlib directory, type make (this assumes you have previously configured the source tree) and if there are errors, correct them and recompile. If not, you go to the gmx/src/kernel directory and type make install. Now you are ready to test the changes.

    In the above example we have forced the field to lie along the Z variable in the box, there is no real need for this but this is just for simplicity. The variable Ex[ZZ].a is actually set in the mdp file (field E-z). You may want to add the following to your mdp file:

     E-z = 1 5 0  
    

    This would give you a field with a strength of 5 V/nm. Now you can play with the field strength in the .mdp file and doing short, e.g. 100 ps simulatios on a box of 216 water molecules.

    Page last modified 00:32, 13 Oct 2009 by JLemkul?