
Tabulated PotentialsTable of contentsNonbonded tabulated interactionsTheoryIt is possible to specify intermolecular interaction functions that are different from the normal LennardJones 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 LennardJones 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, LennardJones or Buckingham, is determined by the To utilize a userdefined potential, several elements are necessary:
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 singleprecision GROMACS, a spacing of 0.002 nm is adequate. For doubleprecision, 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 r_{c} + 1, where r_{c} 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 96 table (i.e., table69.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 612 table (table612.xvg): #include <stdio.h> #include <math.h> main() { FILE *fout; double r; fout = fopen("table_example.xvg", "w"); fprintf(fout, "#\n# Example LJ 612 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/table612.xvg). The only value for which this is absolutely required is r=0.
Putting It All TogetherHere 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 AB interactions), table_A_A.xvg (for AA interactions), and table_B_B.xvg (for BB 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 interactionsNobody's written anything here yet, but you should start with reading manual section 4.2.13 