Gromacs

Tabulated Potentials

    Non-bonded tabulated interactions

    Theory

    It is possible to specify intermolecular interaction functions that are different from the normal Lennard-Jones or Buckingham type.  To implement such functions, one can make use of tabulated potentials or otherwise modify the source code, which is much more laborious.  For information on tabulated potentials, see the manual, section 6.7.

    The generic form of the nonbonded potential function used by GROMACS is:

     

    Where the functions (for a normal Lennard-Jones potential) are defined as:

     

     

     

    For a Buckingham potential, f(r) and g(r) are the same, but h(r) is defined as:

     

    The parameters A, B, and C are read from the topology.  The type of interaction, Lennard-Jones or Buckingham, is determined by the nbfunc parameter in the [ defaults ] directive of the topology.

    To utilize a user-defined potential, several elements are necessary:

    1. Proper settings in the .mdp file for vdwtype and coulombtype.  The setting for vdwtype should be User and coulombtype can be any one of: User, PME-User, or PME-User-Switch. See the manual for details.

    2. Appropriate table.xvg files containing the values pertaining to the functions described above.  More on this below.

    Constructing the Table(s)

    The GROMACS manual, section 6.7, describes the requirements for the construction of the tabulated potentials.  The format of the lines in the table must be:

    The spacing of the table (i.e., the interval between values of r) must be uniform.  For single-precision GROMACS, a spacing of 0.002 nm is adequate.  For double-precision, use 0.0005 nm.  Values between those given in the table.xvg file will be interpolated by GROMACS, using a cubic spline method.  The values of r must go up to a value of rc + 1, where rc is the longest cutoff value defined in the .mdp file.  Example tables can be found in the /share/gromacs/top subdirectory of your GROMACS installation.

    The following Fortran code generates a 9-6 table (i.e., table6-9.xvg in /share/gromacs/top):

    program gen_table
    implicit none
    real,parameter :: delr=0.002,rcut=1.0
    real :: r
    integer :: nbins,j
    
    nbins=int( (rcut+1)/delr) + 1
    
    do j=0,nbins
        r=delr*j
        write(6,*)r, 1/r, 1/(r*r),−1/(r**6),−6/(r**7), 1/(r**9), 9/(r**10)
    end do

    The following C code generates a 6-12 table (table6-12.xvg):

    #include <stdio.h>
    #include <math.h>
    
    main()
    {
        FILE    *fout;
        double  r;
    
        fout = fopen("table_example.xvg", "w");
        fprintf(fout, "#\n# Example LJ 6-12 Potential\n#\n");
    
        for (r=0; r<=3; r+=0.002) {
    
            double f = 1/r;
            double fprime = 1/(pow(r,2));
    
            double g = -1/(pow(r,6));
            double gprime = -6/(pow(r,7));
    
            double h = 1/(pow(r,12));
            double hprime = 12/(pow(r,13));
    
            /* print output */
            if (r<0.04) {
                fprintf(fout, "%12.10e   %12.10e %12.10e   %12.10e %12.10e   %12.10e %12.10e\n", r,0.0,0.0,0.0,0.0,0.0,0.0);
            } else {
                fprintf(fout, "%12.10e   %12.10e %12.10e   %12.10e %12.10e   %12.10e %12.10e\n", r,f,fprime,g,gprime,h,hprime);
            }
        }
    
        fclose(fout);
        return(0);
    }

    Note that when r=0 you will have either infinite or undefined values.  The C code above assigns zero values for all functions at any value of r<0.04 nm, as in the example table (/share/gromacs/top/table6-12.xvg).  The only value for which this is absolutely required is r=0.

     

    Putting It All Together

    Here is an example of a system for which there are two particle types, called A and B.  Nonbonded parameters (A and C in the case of an LJ potential; A, B, and C for Buckingham) for both particle types are specified in the topology.  To tell grompp (and ultimately mdrun) that tabulated potentials should be used, set the following:

    vdwtype = User
    coulombtype = User
    energygrps = A B
    energygrp_table = A A B B
    

    Now mdrun will expect to find three files: table.xvg (for A-B interactions), table_A_A.xvg (for A-A interactions), and table_B_B.xvg (for B-B interactions).  If necessary, you must build the appropriate index groups for A and B before running grompp.

     

    For more detailed instructions and a few examples, see the attached file, provided by Gareth Tribello.

    Bonded tabulated interactions

    Nobody's written anything here yet, but you should start with reading manual section 4.2.13

    Page last modified 21:21, 9 Feb 2014 by JLemkul?