PotsExport2Gromacs

PotsExport2Gromacs(pots, npoints=2500, Umax=6000,
Rmaxtable=2.5, PHImaxtable=180, filename='', noplot=False, hardcopy=True, figsize=(14, 7.5), dpi=120, sigma=0.5,
zeroforce=True, interpol=True, ofilename='')

This procedure exports a set of potentials into GROMACS's .xvg format. Such an export is a rather advanced and multistage process, here a detailed description of the procedure is given.

The first thing to take into account is that Gromacs requires typically a higher resolution of the grid for tabulated potentials than typical resolution of the RDF inversion which is used in MagiC. Furthermore, the potential and force should be smooth functions of the distance, and do not have discontinuities at the end of the intervals where they are determined. Thus the original range where the potentials and RDFs were determined during the inversion procedure, should be extended in a smooth fashion to nearby range of distances.

Thus in general the question is the following: we start with a tabulated potential which is defined on a r-grid $[r_{min}:r_{max}]$ with a given density (figure 4-A, red circles); and as the result we need to obtain a tabulated potential, which is defined on a r-grid $[0:r_{vdw}+r_{table-extension}]$ with larger density (figure 4-D, cyan line).

In order to do that, we need to introduce left-side and right-side extensions of the potential which should be smoothly connected to the original potential (figure 4-B, blue circles). The left side extension should represent strongly repulsive core so it is approximated by

$$U^{left}(r)=ar^2+br+U_{max}$$ (7)

and coefficients $a$ and $b$ are chosen to provide continuity of $\frac{d}{dr}U(r)$ at $r=r_{min}$ . For the right side extension we should take into account the origin of the potential: Short-range intermolecular potentials should decay to zero when $r>r_{max}$ :

$$U^{right}_{short\,range}(r)=U(r_{max})\cdot\exp[{-10\frac{(r-r_{max})}{r_{vdw+table-extension}-r_{max}}}]$$ (8)

Here 10 is just some pre-defined coefficient, $r_{vdw+table-extension}$ is a range of the table required by GROMACS. Angle bending potentials should also decay to zero at $\phi=180^{\circ}$ :

$$U^{right}_{angle}(\phi)=U(\phi_{max})\cdot\exp[{-100\frac{(\phi-\phi_{max})}{180^{\circ}-\phi_{max}}}]$$ (9)

In contrast, pairwise bond potentials should have an attractive wall at the right side, which is approximated by harmonic wall in the same way as the repulsive wall at the left side:

$$U^{right}_{pair\,bond}(r)=ar^2+br+U_{max}$$ (10)

coefficients $a$ and $b$ are chosen to provide continuity of $\frac{d}{dr}U(r)$ at $r=r_{max}$ .

Once the extension has been made, we can interpolate all the points to a denser grid (figure 4-C). Number of nodes in the grid is defined by npoints. If interpol=False, a grid of original density is used. The interpolation is made by Gaussian smoothing, e.g. every point of the resulting potential is obtained as

$$U^{new}(r)=\frac{1}{Z(r)}\sum_{r_i=r_{min}}^{r_{max}}{U^{orig}(r_i)\cdot\exp{\frac{-(r-r_i)^2}{2\sigma^2}}}$$ (11)

$$Z(r)=\sum_{r_i=r_{min}}^{r_{max}}{\exp{\frac{-(r-r_i)^2}{2\sigma^2}}}$$ (12)

where sigma defines how broad is the averaging. By default $\sigma=0.5 \Delta r_i$

Once the interpolation is done, each resulting potential and force based on it are written to a xvg-file. The xvg-file is named by the name and type of the corresponding original potential. Forces can be suppressed by setting noforce=True, then zeros will be written to the file. In such case GROMACS should automatically calculate forces from a given potential.

In order to control the results, each original potential, the extrapolated part and the interpolated potential are plotted together. The plotting is controlled by parameters noplot, hardcopy, figsize, dpi, which are explained below.

Figure 4: Original potential - red circles

pots*

- set of potentials to export into GROMACS .xvg format (mandatory argument)

ofilename

- Prefix used for naming of the output .xvg files (optional argument).

zeroforce

- do not write forces into .xvg-file, but write zeros instead (optional argument). In such case GROMACS should automatically calculate forces from potentials. Default: False - forces are to be written. NB: Even if zeroforce=True force values are plotted.

npoints

- number of points in the resulting table in .xvg file (optional argument). Default value - 2500 points.

Umax

- height of potential wall at r=0 in case of non-bonding potential and at $r=r_{min}$ , $r=r_{max}$ in case of bonding potential (optional argument). Default value 6000 kJ/mol

Rmaxtable

- cutoff range of the resulting potentials (both intermolecular and pair bonds) in nanometers (optional argument). Default value is 2.5 nm. Note that GROMACS requires tabulated potentials to be defined up to r=rvdw+table-extension (in terms of GROMACS parameters).

sigma

- Gaussian smoothing parameter measured in resolution of the original potential (optional parameter). Default value 0.5

interpol

- if the potentials should be interpolated, otherwise original resolution of the table will be kept and npoints value will be ignored (optional argument). Default - True, potentials are interpolated.

noplot

- if interpolated potentials potential/forces should be plotted (optional argument). Default: False - the graphs are plotted.

hardcopy

- if the plots should be saved as *.eps files (optional argument). Default value - True.

figsize

- size of the plot in inches (x,y) (optional argument). Default size is (14,8).

dpi

- resolution of the plot in dpi (optional argument). Default value: 120 dpi.

Examples:
MagicTools.PotsExport2Gromacs(Pots) - simplest call with default parameters.
MagicTools.PotsExport2Gromacs(Pots[0:10], sigma=0.1, zeroforce=True, interpol=False ) - export only first 10 potentials from the set; do not use dense net for interpolation; use averaging with sigma=0.1 (very narrow gaussian); do not export forces, write zeros instead.