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.
45 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/utility/futil.h"
47 #include "gromacs/commandline/pargs.h"
53 #include "gromacs/fileio/tpxio.h"
54 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/linearalgebra/nrjac.h"
59 static void gyro_eigen(double **gyr, double *eig, double **eigv, int *ord)
63 jacobi(gyr, DIM, eig, eigv, &nrot);
64 /* Order the eigenvalues */
67 for (d = 0; d < DIM; d++)
69 if (eig[d] > eig[ord[0]])
73 if (eig[d] < eig[ord[2]])
78 for (d = 0; d < DIM; d++)
80 if (ord[0] != d && ord[2] != d)
87 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
88 static void calc_int_dist(double *intd, rvec *x, int i0, int i1)
91 const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
92 int bd; /* Distance between beads */
95 for (bd = 1; bd < ml; bd++)
98 for (ii = i0; ii <= i1 - bd; ii++)
100 d += distance2(x[ii], x[ii+bd]);
107 int gmx_polystat(int argc, char *argv[])
109 const char *desc[] = {
110 "[THISMODULE] plots static properties of polymers as a function of time",
111 "and prints the average.[PAR]",
112 "By default it determines the average end-to-end distance and radii",
113 "of gyration of polymers. It asks for an index group and split this",
114 "into molecules. The end-to-end distance is then determined using",
115 "the first and the last atom in the index group for each molecules.",
116 "For the radius of gyration the total and the three principal components",
117 "for the average gyration tensor are written.",
118 "With option [TT]-v[tt] the eigenvectors are written.",
119 "With option [TT]-pc[tt] also the average eigenvalues of the individual",
120 "gyration tensors are written.",
121 "With option [TT]-i[tt] the mean square internal distances are",
123 "With option [TT]-p[tt] the persistence length is determined.",
124 "The chosen index group should consist of atoms that are",
125 "consecutively bonded in the polymer mainchains.",
126 "The persistence length is then determined from the cosine of",
127 "the angles between bonds with an index difference that is even,",
128 "the odd pairs are not used, because straight polymer backbones",
129 "are usually all trans and therefore only every second bond aligns.",
130 "The persistence length is defined as number of bonds where",
131 "the average cos reaches a value of 1/e. This point is determined",
132 "by a linear interpolation of log(<cos>)."
134 static gmx_bool bMW = TRUE, bPC = FALSE;
136 { "-mw", FALSE, etBOOL, {&bMW},
137 "Use the mass weighting for radii of gyration" },
138 { "-pc", FALSE, etBOOL, {&bPC},
139 "Plot average eigenvalues" }
143 { efTPX, NULL, NULL, ffREAD },
144 { efTRX, "-f", NULL, ffREAD },
145 { efNDX, NULL, NULL, ffOPTRD },
146 { efXVG, "-o", "polystat", ffWRITE },
147 { efXVG, "-v", "polyvec", ffOPTWR },
148 { efXVG, "-p", "persist", ffOPTWR },
149 { efXVG, "-i", "intdist", ffOPTWR }
151 #define NFILE asize(fnm)
156 int isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
160 rvec *x, *bond = NULL;
162 int natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = {0};
163 dvec cm, sum_eig = {0, 0, 0};
164 double **gyr, **gyr_all, eig[DIM], **eigv;
165 double sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
167 double *sum_inp = NULL, pers;
168 double *intd, ymax, ymin;
171 FILE *out, *outv, *outp, *outi;
172 const char *leg[8] = {
173 "end to end", "<R\\sg\\N>",
174 "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
175 "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>"
177 char **legp, buf[STRLEN];
178 gmx_rmpbc_t gpbc = NULL;
180 if (!parse_common_args(&argc, argv,
181 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
182 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
188 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
189 NULL, box, &natoms, NULL, NULL, NULL, top);
191 fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
192 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
193 1, &isize, &index, &grpname);
195 snew(molind, top->mols.nr+1);
198 for (i = 0; i < isize; i++)
200 if (i == 0 || index[i] >= top->mols.index[mol+1])
207 while (index[i] >= top->mols.index[mol+1]);
211 nat_min = top->atoms.nr;
213 for (mol = 0; mol < nmol; mol++)
215 nat_min = min(nat_min, molind[mol+1]-molind[mol]);
216 nat_max = max(nat_max, molind[mol+1]-molind[mol]);
218 fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
219 fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
222 sprintf(title, "Size of %d polymers", nmol);
223 out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
225 xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
227 if (opt2bSet("-v", NFILE, fnm))
229 outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
230 output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
232 for (d = 0; d < DIM; d++)
234 for (d2 = 0; d2 < DIM; d2++)
236 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
237 legp[d*DIM+d2] = strdup(buf);
240 xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
247 if (opt2bSet("-p", NFILE, fnm))
249 outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
250 output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
251 snew(bond, nat_max-1);
252 snew(sum_inp, nat_min/2);
253 snew(ninp, nat_min/2);
260 if (opt2bSet("-i", NFILE, fnm))
262 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
263 "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
264 i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
273 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
278 for (d = 0; d < DIM; d++)
281 snew(gyr_all[d], DIM);
290 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
294 gmx_rmpbc(gpbc, natoms, box, x);
297 for (d = 0; d < DIM; d++)
299 clear_dvec(gyr_all[d]);
309 for (i = 0; i < nat_min/2; i++)
316 for (mol = 0; mol < nmol; mol++)
319 ind1 = molind[mol+1];
321 /* Determine end to end distance */
322 sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
324 /* Determine internal distances */
327 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
330 /* Determine the radius of gyration */
332 for (d = 0; d < DIM; d++)
338 for (i = ind0; i < ind1; i++)
343 m = top->atoms.atom[a].m;
350 for (d = 0; d < DIM; d++)
353 for (d2 = 0; d2 < DIM; d2++)
355 gyr[d][d2] += m*x[a][d]*x[a][d2];
359 dsvmul(1/mmol, cm, cm);
360 for (d = 0; d < DIM; d++)
362 for (d2 = 0; d2 < DIM; d2++)
364 gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
365 gyr_all[d][d2] += gyr[d][d2];
370 gyro_eigen(gyr, eig, eigv, ord);
371 for (d = 0; d < DIM; d++)
373 sum_eig[d] += eig[ord[d]];
378 for (i = ind0; i < ind1-1; i++)
380 rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
381 unitv(bond[i-ind0], bond[i-ind0]);
383 for (i = ind0; i < ind1-1; i++)
385 for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
387 sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
396 for (d = 0; d < DIM; d++)
398 for (d2 = 0; d2 < DIM; d2++)
400 gyr_all[d][d2] /= nmol;
402 sum_gyro += gyr_all[d][d];
405 gyro_eigen(gyr_all, eig, eigv, ord);
407 fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
408 t*output_env_get_time_factor(oenv),
409 sqrt(sum_eed2), sqrt(sum_gyro),
410 sqrt(eig[ord[0]]), sqrt(eig[ord[1]]), sqrt(eig[ord[2]]));
413 for (d = 0; d < DIM; d++)
415 fprintf(out, " %8.4f", sqrt(sum_eig[d]/nmol));
422 fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
423 for (d = 0; d < DIM; d++)
425 for (d2 = 0; d2 < DIM; d2++)
427 fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
433 sum_eed2_tot += sum_eed2;
434 sum_gyro_tot += sum_gyro;
439 for (j = 0; j < nat_min/2; j += 2)
441 sum_inp[j] /= ninp[j];
442 if (i == -1 && sum_inp[j] <= exp(-1.0))
453 /* Do linear interpolation on a log scale */
455 + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
457 fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
458 sum_pers_tot += pers;
463 while (read_next_x(oenv, status, &t, x, box));
465 gmx_rmpbc_done(gpbc);
479 sum_eed2_tot /= frame;
480 sum_gyro_tot /= frame;
481 sum_pers_tot /= frame;
482 fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
484 fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n",
486 if (opt2bSet("-p", NFILE, fnm))
488 fprintf(stdout, "\nAverage persistence length: %.2f bonds\n",
492 /* Handle printing of internal distances. */
495 fprintf(outi, "@ xaxes scale Logarithmic\n");
498 j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
499 for (i = 0; i < j; i++)
501 intd[i] /= (i + 1) * frame * nmol;
511 xvgr_world(outi, 1, ymin, j, ymax, oenv);
512 for (i = 0; i < j; i++)
514 fprintf(outi, "%d %8.4f\n", i+1, intd[i]);
519 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
520 if (opt2bSet("-v", NFILE, fnm))
522 do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
524 if (opt2bSet("-p", NFILE, fnm))
526 do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");