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.
47 #include "gromacs/fileio/futil.h"
48 #include "gromacs/commandline/pargs.h"
52 #include "gmx_fatal.h"
55 #include "gromacs/fileio/tpxio.h"
56 #include "gromacs/fileio/trxio.h"
61 static void gyro_eigen(double **gyr, double *eig, double **eigv, int *ord)
65 jacobi(gyr, DIM, eig, eigv, &nrot);
66 /* Order the eigenvalues */
69 for (d = 0; d < DIM; d++)
71 if (eig[d] > eig[ord[0]])
75 if (eig[d] < eig[ord[2]])
80 for (d = 0; d < DIM; d++)
82 if (ord[0] != d && ord[2] != d)
89 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
90 static void calc_int_dist(double *intd, rvec *x, int i0, int i1)
93 const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
94 int bd; /* Distance between beads */
97 for (bd = 1; bd < ml; bd++)
100 for (ii = i0; ii <= i1 - bd; ii++)
102 d += distance2(x[ii], x[ii+bd]);
109 int gmx_polystat(int argc, char *argv[])
111 const char *desc[] = {
112 "[THISMODULE] plots static properties of polymers as a function of time",
113 "and prints the average.[PAR]",
114 "By default it determines the average end-to-end distance and radii",
115 "of gyration of polymers. It asks for an index group and split this",
116 "into molecules. The end-to-end distance is then determined using",
117 "the first and the last atom in the index group for each molecules.",
118 "For the radius of gyration the total and the three principal components",
119 "for the average gyration tensor are written.",
120 "With option [TT]-v[tt] the eigenvectors are written.",
121 "With option [TT]-pc[tt] also the average eigenvalues of the individual",
122 "gyration tensors are written.",
123 "With option [TT]-i[tt] the mean square internal distances are",
125 "With option [TT]-p[tt] the persistence length is determined.",
126 "The chosen index group should consist of atoms that are",
127 "consecutively bonded in the polymer mainchains.",
128 "The persistence length is then determined from the cosine of",
129 "the angles between bonds with an index difference that is even,",
130 "the odd pairs are not used, because straight polymer backbones",
131 "are usually all trans and therefore only every second bond aligns.",
132 "The persistence length is defined as number of bonds where",
133 "the average cos reaches a value of 1/e. This point is determined",
134 "by a linear interpolation of log(<cos>)."
136 static gmx_bool bMW = TRUE, bPC = FALSE;
138 { "-mw", FALSE, etBOOL, {&bMW},
139 "Use the mass weighting for radii of gyration" },
140 { "-pc", FALSE, etBOOL, {&bPC},
141 "Plot average eigenvalues" }
145 { efTPX, NULL, NULL, ffREAD },
146 { efTRX, "-f", NULL, ffREAD },
147 { efNDX, NULL, NULL, ffOPTRD },
148 { efXVG, "-o", "polystat", ffWRITE },
149 { efXVG, "-v", "polyvec", ffOPTWR },
150 { efXVG, "-p", "persist", ffOPTWR },
151 { efXVG, "-i", "intdist", ffOPTWR }
153 #define NFILE asize(fnm)
158 int isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
162 rvec *x, *bond = NULL;
164 int natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = {0};
165 dvec cm, sum_eig = {0, 0, 0};
166 double **gyr, **gyr_all, eig[DIM], **eigv;
167 double sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
169 double *sum_inp = NULL, pers;
170 double *intd, ymax, ymin;
173 FILE *out, *outv, *outp, *outi;
174 const char *leg[8] = {
175 "end to end", "<R\\sg\\N>",
176 "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
177 "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>"
179 char **legp, buf[STRLEN];
180 gmx_rmpbc_t gpbc = NULL;
182 if (!parse_common_args(&argc, argv,
183 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
184 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
190 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
191 NULL, box, &natoms, NULL, NULL, NULL, top);
193 fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
194 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
195 1, &isize, &index, &grpname);
197 snew(molind, top->mols.nr+1);
200 for (i = 0; i < isize; i++)
202 if (i == 0 || index[i] >= top->mols.index[mol+1])
209 while (index[i] >= top->mols.index[mol+1]);
213 nat_min = top->atoms.nr;
215 for (mol = 0; mol < nmol; mol++)
217 nat_min = min(nat_min, molind[mol+1]-molind[mol]);
218 nat_max = max(nat_max, molind[mol+1]-molind[mol]);
220 fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
221 fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
224 sprintf(title, "Size of %d polymers", nmol);
225 out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
227 xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
229 if (opt2bSet("-v", NFILE, fnm))
231 outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
232 output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
234 for (d = 0; d < DIM; d++)
236 for (d2 = 0; d2 < DIM; d2++)
238 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
239 legp[d*DIM+d2] = strdup(buf);
242 xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
249 if (opt2bSet("-p", NFILE, fnm))
251 outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
252 output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
253 snew(bond, nat_max-1);
254 snew(sum_inp, nat_min/2);
255 snew(ninp, nat_min/2);
262 if (opt2bSet("-i", NFILE, fnm))
264 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
265 "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
266 i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
275 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
280 for (d = 0; d < DIM; d++)
283 snew(gyr_all[d], DIM);
292 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
296 gmx_rmpbc(gpbc, natoms, box, x);
299 for (d = 0; d < DIM; d++)
301 clear_dvec(gyr_all[d]);
311 for (i = 0; i < nat_min/2; i++)
318 for (mol = 0; mol < nmol; mol++)
321 ind1 = molind[mol+1];
323 /* Determine end to end distance */
324 sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
326 /* Determine internal distances */
329 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
332 /* Determine the radius of gyration */
334 for (d = 0; d < DIM; d++)
340 for (i = ind0; i < ind1; i++)
345 m = top->atoms.atom[a].m;
352 for (d = 0; d < DIM; d++)
355 for (d2 = 0; d2 < DIM; d2++)
357 gyr[d][d2] += m*x[a][d]*x[a][d2];
361 dsvmul(1/mmol, cm, cm);
362 for (d = 0; d < DIM; d++)
364 for (d2 = 0; d2 < DIM; d2++)
366 gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
367 gyr_all[d][d2] += gyr[d][d2];
372 gyro_eigen(gyr, eig, eigv, ord);
373 for (d = 0; d < DIM; d++)
375 sum_eig[d] += eig[ord[d]];
380 for (i = ind0; i < ind1-1; i++)
382 rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
383 unitv(bond[i-ind0], bond[i-ind0]);
385 for (i = ind0; i < ind1-1; i++)
387 for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
389 sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
398 for (d = 0; d < DIM; d++)
400 for (d2 = 0; d2 < DIM; d2++)
402 gyr_all[d][d2] /= nmol;
404 sum_gyro += gyr_all[d][d];
407 gyro_eigen(gyr_all, eig, eigv, ord);
409 fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
410 t*output_env_get_time_factor(oenv),
411 sqrt(sum_eed2), sqrt(sum_gyro),
412 sqrt(eig[ord[0]]), sqrt(eig[ord[1]]), sqrt(eig[ord[2]]));
415 for (d = 0; d < DIM; d++)
417 fprintf(out, " %8.4f", sqrt(sum_eig[d]/nmol));
424 fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
425 for (d = 0; d < DIM; d++)
427 for (d2 = 0; d2 < DIM; d2++)
429 fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
435 sum_eed2_tot += sum_eed2;
436 sum_gyro_tot += sum_gyro;
441 for (j = 0; j < nat_min/2; j += 2)
443 sum_inp[j] /= ninp[j];
444 if (i == -1 && sum_inp[j] <= exp(-1.0))
455 /* Do linear interpolation on a log scale */
457 + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
459 fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
460 sum_pers_tot += pers;
465 while (read_next_x(oenv, status, &t, x, box));
467 gmx_rmpbc_done(gpbc);
481 sum_eed2_tot /= frame;
482 sum_gyro_tot /= frame;
483 sum_pers_tot /= frame;
484 fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
486 fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n",
488 if (opt2bSet("-p", NFILE, fnm))
490 fprintf(stdout, "\nAverage persistence length: %.2f bonds\n",
494 /* Handle printing of internal distances. */
497 fprintf(outi, "@ xaxes scale Logarithmic\n");
500 j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
501 for (i = 0; i < j; i++)
503 intd[i] /= (i + 1) * frame * nmol;
513 xvgr_world(outi, 1, ymin, j, ymax, oenv);
514 for (i = 0; i < j; i++)
516 fprintf(outi, "%d %8.4f\n", i+1, intd[i]);
521 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
522 if (opt2bSet("-v", NFILE, fnm))
524 do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
526 if (opt2bSet("-p", NFILE, fnm))
528 do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");