2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/tpxio.h"
44 #include "gromacs/fileio/trxio.h"
45 #include "gromacs/fileio/xvgr.h"
46 #include "gromacs/gmxana/gmx_ana.h"
47 #include "gromacs/legacyheaders/macros.h"
48 #include "gromacs/legacyheaders/typedefs.h"
49 #include "gromacs/legacyheaders/viewit.h"
50 #include "gromacs/linearalgebra/nrjac.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/pbcutil/rmpbc.h"
53 #include "gromacs/topology/index.h"
54 #include "gromacs/utility/futil.h"
55 #include "gromacs/utility/smalloc.h"
57 static void gyro_eigen(double **gyr, double *eig, double **eigv, int *ord)
61 jacobi(gyr, DIM, eig, eigv, &nrot);
62 /* Order the eigenvalues */
65 for (d = 0; d < DIM; d++)
67 if (eig[d] > eig[ord[0]])
71 if (eig[d] < eig[ord[2]])
76 for (d = 0; d < DIM; d++)
78 if (ord[0] != d && ord[2] != d)
85 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
86 static void calc_int_dist(double *intd, rvec *x, int i0, int i1)
89 const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
90 int bd; /* Distance between beads */
93 for (bd = 1; bd < ml; bd++)
96 for (ii = i0; ii <= i1 - bd; ii++)
98 d += distance2(x[ii], x[ii+bd]);
105 int gmx_polystat(int argc, char *argv[])
107 const char *desc[] = {
108 "[THISMODULE] plots static properties of polymers as a function of time",
109 "and prints the average.[PAR]",
110 "By default it determines the average end-to-end distance and radii",
111 "of gyration of polymers. It asks for an index group and split this",
112 "into molecules. The end-to-end distance is then determined using",
113 "the first and the last atom in the index group for each molecules.",
114 "For the radius of gyration the total and the three principal components",
115 "for the average gyration tensor are written.",
116 "With option [TT]-v[tt] the eigenvectors are written.",
117 "With option [TT]-pc[tt] also the average eigenvalues of the individual",
118 "gyration tensors are written.",
119 "With option [TT]-i[tt] the mean square internal distances are",
121 "With option [TT]-p[tt] the persistence length is determined.",
122 "The chosen index group should consist of atoms that are",
123 "consecutively bonded in the polymer mainchains.",
124 "The persistence length is then determined from the cosine of",
125 "the angles between bonds with an index difference that is even,",
126 "the odd pairs are not used, because straight polymer backbones",
127 "are usually all trans and therefore only every second bond aligns.",
128 "The persistence length is defined as number of bonds where",
129 "the average cos reaches a value of 1/e. This point is determined",
130 "by a linear interpolation of log(<cos>)."
132 static gmx_bool bMW = TRUE, bPC = FALSE;
134 { "-mw", FALSE, etBOOL, {&bMW},
135 "Use the mass weighting for radii of gyration" },
136 { "-pc", FALSE, etBOOL, {&bPC},
137 "Plot average eigenvalues" }
141 { efTPR, NULL, NULL, ffREAD },
142 { efTRX, "-f", NULL, ffREAD },
143 { efNDX, NULL, NULL, ffOPTRD },
144 { efXVG, "-o", "polystat", ffWRITE },
145 { efXVG, "-v", "polyvec", ffOPTWR },
146 { efXVG, "-p", "persist", ffOPTWR },
147 { efXVG, "-i", "intdist", ffOPTWR }
149 #define NFILE asize(fnm)
154 int isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
158 rvec *x, *bond = NULL;
160 int natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = {0};
161 dvec cm, sum_eig = {0, 0, 0};
162 double **gyr, **gyr_all, eig[DIM], **eigv;
163 double sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
165 double *sum_inp = NULL, pers;
166 double *intd, ymax, ymin;
169 FILE *out, *outv, *outp, *outi;
170 const char *leg[8] = {
171 "end to end", "<R\\sg\\N>",
172 "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
173 "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>"
175 char **legp, buf[STRLEN];
176 gmx_rmpbc_t gpbc = NULL;
178 if (!parse_common_args(&argc, argv,
179 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
180 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
186 ePBC = read_tpx_top(ftp2fn(efTPR, NFILE, fnm),
187 NULL, box, &natoms, NULL, NULL, NULL, top);
189 fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
190 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
191 1, &isize, &index, &grpname);
193 snew(molind, top->mols.nr+1);
196 for (i = 0; i < isize; i++)
198 if (i == 0 || index[i] >= top->mols.index[mol+1])
205 while (index[i] >= top->mols.index[mol+1]);
209 nat_min = top->atoms.nr;
211 for (mol = 0; mol < nmol; mol++)
213 nat_min = min(nat_min, molind[mol+1]-molind[mol]);
214 nat_max = max(nat_max, molind[mol+1]-molind[mol]);
216 fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
217 fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
220 sprintf(title, "Size of %d polymers", nmol);
221 out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
223 xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
225 if (opt2bSet("-v", NFILE, fnm))
227 outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
228 output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
230 for (d = 0; d < DIM; d++)
232 for (d2 = 0; d2 < DIM; d2++)
234 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
235 legp[d*DIM+d2] = gmx_strdup(buf);
238 xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
245 if (opt2bSet("-p", NFILE, fnm))
247 outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
248 output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
249 snew(bond, nat_max-1);
250 snew(sum_inp, nat_min/2);
251 snew(ninp, nat_min/2);
258 if (opt2bSet("-i", NFILE, fnm))
260 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
261 "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
262 i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
271 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
276 for (d = 0; d < DIM; d++)
279 snew(gyr_all[d], DIM);
288 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
292 gmx_rmpbc(gpbc, natoms, box, x);
295 for (d = 0; d < DIM; d++)
297 clear_dvec(gyr_all[d]);
307 for (i = 0; i < nat_min/2; i++)
314 for (mol = 0; mol < nmol; mol++)
317 ind1 = molind[mol+1];
319 /* Determine end to end distance */
320 sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
322 /* Determine internal distances */
325 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
328 /* Determine the radius of gyration */
330 for (d = 0; d < DIM; d++)
336 for (i = ind0; i < ind1; i++)
341 m = top->atoms.atom[a].m;
348 for (d = 0; d < DIM; d++)
351 for (d2 = 0; d2 < DIM; d2++)
353 gyr[d][d2] += m*x[a][d]*x[a][d2];
357 dsvmul(1/mmol, cm, cm);
358 for (d = 0; d < DIM; d++)
360 for (d2 = 0; d2 < DIM; d2++)
362 gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
363 gyr_all[d][d2] += gyr[d][d2];
368 gyro_eigen(gyr, eig, eigv, ord);
369 for (d = 0; d < DIM; d++)
371 sum_eig[d] += eig[ord[d]];
376 for (i = ind0; i < ind1-1; i++)
378 rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
379 unitv(bond[i-ind0], bond[i-ind0]);
381 for (i = ind0; i < ind1-1; i++)
383 for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
385 sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
394 for (d = 0; d < DIM; d++)
396 for (d2 = 0; d2 < DIM; d2++)
398 gyr_all[d][d2] /= nmol;
400 sum_gyro += gyr_all[d][d];
403 gyro_eigen(gyr_all, eig, eigv, ord);
405 fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
406 t*output_env_get_time_factor(oenv),
407 sqrt(sum_eed2), sqrt(sum_gyro),
408 sqrt(eig[ord[0]]), sqrt(eig[ord[1]]), sqrt(eig[ord[2]]));
411 for (d = 0; d < DIM; d++)
413 fprintf(out, " %8.4f", sqrt(sum_eig[d]/nmol));
420 fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
421 for (d = 0; d < DIM; d++)
423 for (d2 = 0; d2 < DIM; d2++)
425 fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
431 sum_eed2_tot += sum_eed2;
432 sum_gyro_tot += sum_gyro;
437 for (j = 0; j < nat_min/2; j += 2)
439 sum_inp[j] /= ninp[j];
440 if (i == -1 && sum_inp[j] <= exp(-1.0))
451 /* Do linear interpolation on a log scale */
453 + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
455 fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
456 sum_pers_tot += pers;
461 while (read_next_x(oenv, status, &t, x, box));
463 gmx_rmpbc_done(gpbc);
477 sum_eed2_tot /= frame;
478 sum_gyro_tot /= frame;
479 sum_pers_tot /= frame;
480 fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
482 fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n",
484 if (opt2bSet("-p", NFILE, fnm))
486 fprintf(stdout, "\nAverage persistence length: %.2f bonds\n",
490 /* Handle printing of internal distances. */
493 if (output_env_get_print_xvgr_codes(oenv))
495 fprintf(outi, "@ xaxes scale Logarithmic\n");
499 j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
500 for (i = 0; i < j; i++)
502 intd[i] /= (i + 1) * frame * nmol;
512 xvgr_world(outi, 1, ymin, j, ymax, oenv);
513 for (i = 0; i < j; i++)
515 fprintf(outi, "%d %8.4f\n", i+1, intd[i]);
520 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
521 if (opt2bSet("-v", NFILE, fnm))
523 do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
525 if (opt2bSet("-p", NFILE, fnm))
527 do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");