[MD-010] Polymer Trajectory Analysis (structure and hydration) using Gromacs

Gromacs embedded module for polymer trajectory analysis

In a previous post, the general procedure of MD trajectory analysis has been discussed [here]. Here more detailed analysis methods are demonstrated for polymers trajectory analysis using the embedded modules in Gromacs, which is more than a simulation engine.

Structures and hydration analysis

Command	Description
gmx polystat	This calculates the polymer statistics data as a function of time, including the end-to-end distance, radius of gyration and its three principal components. With -p flag, it also calculates the persistence length of a polymer when only the backbone atoms are provided
gmx gyrate	This calculates the polymer statistics data as a function of time, including the radius of gyration and its the mass-weighted root-mean-square of the radii components orthogonal to each axis. With -p flag, it calculates based on the principal axes, which is helpful in obtaining the principal components of a polymer, see below
gmx hbond	This calculates the hydrogen bonds between two groups of atoms
gmx sasa	This calculates the solvent accessible surface area of a polymer
gmx distance	This calculates the distance between two groups
gmx angle	This calculates the angle between three groups

In calculating the persistence length, Gromacs uses the backbone atoms and calculates when the bond orientation correlation decays to exp(-1) using every other two atoms as a step and return how many bonds (e.g., c-c bonds rather than repeat units) reaching that point.

Characterization of shape

In polymer science, the shape of polymer is an indication for polymer conformation and thus can be used to infer the morphology of a polymer, e.g., a globule, coil or a rod. Very often the following shape factors are calculated based on the radius of gyration (R_g) and/or the end-to-end distance (R_end) of a polymer:

aspect ratio: R_end/R_g
asphericity: $R^{2}_{g1} - \frac{1}{2} (R^{2}_{g2} + R^{2}_{g3})$
anisotropy: $\frac{3}{2} \frac{R^{4}_{g1}+R^{4}_{g2}+R^{4}_{g3}}{(R^{2}_{g1}+R^{2}_{g2}+R^{2}_{g3})^2} - \frac{1}{2}$

where the components are defined by:

$R^{2}_{g} = R^{2}_{g1} + R^{2}_{g2} + R^{2}_{g3} \text{ and } R^{2}_{g1} \geq R^{2}_{g2} \geq R^{2}_{g3}$

gmx polystat vs gmx gyrate

Note that there are similarities between gmx polystat and gmx gyrate, however, this is not redundant and can be used when one fails in some calculations. For example, for polymer micelles, if one wants to calculate the radius of gyration of the micelle, rather than each individual chain, then one can use gmx gyrate since gmx polystat will divide the group into many individual chains and calculate the averages among them.

For the same reason, for polymer-grafted nanostructures, if one wants to calculate the chain statistics of each graft, then gmx gyrate can be used with -nmol flag which divides the groups of atoms into nmol chains. The results would be the average among all grafts.

Magic with -p flag in gmx gyrate: without -p, each row of gmx gyrate is:

$t\;R_{g}\;\sqrt{R^{2}_{gy}+R^{2}_{gz}}\;\sqrt{R^{2}_{gx}+R^{2}_{gz}}\;\sqrt{R^{2}_{gx}+R^{2}_{gy}}$

the components are around coordinate system axes, instead of the principal axes of the gyration tensor.

With the -p flag, each row of output becomes:

$t\;R_{g}\;\sqrt{R^{2}_{g2}+R^{2}_{g3}}\;\sqrt{R^{2}_{g1}+R^{2}_{g3}}\;\sqrt{R^{2}_{g1}+R^{2}_{g2}}$

one easily sees that upon one further step, the three principal components can be obtained.