3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
51 #include "gmx_fatal.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 "[TT]g_polystat[tt] 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 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);
185 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
186 NULL, box, &natoms, NULL, NULL, NULL, top);
188 fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
189 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
190 1, &isize, &index, &grpname);
192 snew(molind, top->mols.nr+1);
195 for (i = 0; i < isize; i++)
197 if (i == 0 || index[i] >= top->mols.index[mol+1])
204 while (index[i] >= top->mols.index[mol+1]);
208 nat_min = top->atoms.nr;
210 for (mol = 0; mol < nmol; mol++)
212 nat_min = min(nat_min, molind[mol+1]-molind[mol]);
213 nat_max = max(nat_max, molind[mol+1]-molind[mol]);
215 fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
216 fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
219 sprintf(title, "Size of %d polymers", nmol);
220 out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
222 xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
224 if (opt2bSet("-v", NFILE, fnm))
226 outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
227 output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
229 for (d = 0; d < DIM; d++)
231 for (d2 = 0; d2 < DIM; d2++)
233 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
234 legp[d*DIM+d2] = strdup(buf);
237 xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
244 if (opt2bSet("-p", NFILE, fnm))
246 outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
247 output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
248 snew(bond, nat_max-1);
249 snew(sum_inp, nat_min/2);
250 snew(ninp, nat_min/2);
257 if (opt2bSet("-i", NFILE, fnm))
259 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
260 "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
261 i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
270 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
275 for (d = 0; d < DIM; d++)
278 snew(gyr_all[d], DIM);
287 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
291 gmx_rmpbc(gpbc, natoms, box, x);
294 for (d = 0; d < DIM; d++)
296 clear_dvec(gyr_all[d]);
306 for (i = 0; i < nat_min/2; i++)
313 for (mol = 0; mol < nmol; mol++)
316 ind1 = molind[mol+1];
318 /* Determine end to end distance */
319 sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
321 /* Determine internal distances */
324 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
327 /* Determine the radius of gyration */
329 for (d = 0; d < DIM; d++)
335 for (i = ind0; i < ind1; i++)
340 m = top->atoms.atom[a].m;
347 for (d = 0; d < DIM; d++)
350 for (d2 = 0; d2 < DIM; d2++)
352 gyr[d][d2] += m*x[a][d]*x[a][d2];
356 dsvmul(1/mmol, cm, cm);
357 for (d = 0; d < DIM; d++)
359 for (d2 = 0; d2 < DIM; d2++)
361 gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
362 gyr_all[d][d2] += gyr[d][d2];
367 gyro_eigen(gyr, eig, eigv, ord);
368 for (d = 0; d < DIM; d++)
370 sum_eig[d] += eig[ord[d]];
375 for (i = ind0; i < ind1-1; i++)
377 rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
378 unitv(bond[i-ind0], bond[i-ind0]);
380 for (i = ind0; i < ind1-1; i++)
382 for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
384 sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
393 for (d = 0; d < DIM; d++)
395 for (d2 = 0; d2 < DIM; d2++)
397 gyr_all[d][d2] /= nmol;
399 sum_gyro += gyr_all[d][d];
402 gyro_eigen(gyr_all, eig, eigv, ord);
404 fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
405 t*output_env_get_time_factor(oenv),
406 sqrt(sum_eed2), sqrt(sum_gyro),
407 sqrt(eig[ord[0]]), sqrt(eig[ord[1]]), sqrt(eig[ord[2]]));
410 for (d = 0; d < DIM; d++)
412 fprintf(out, " %8.4f", sqrt(sum_eig[d]/nmol));
419 fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
420 for (d = 0; d < DIM; d++)
422 for (d2 = 0; d2 < DIM; d2++)
424 fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
430 sum_eed2_tot += sum_eed2;
431 sum_gyro_tot += sum_gyro;
436 for (j = 0; j < nat_min/2; j += 2)
438 sum_inp[j] /= ninp[j];
439 if (i == -1 && sum_inp[j] <= exp(-1.0))
450 /* Do linear interpolation on a log scale */
452 + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
454 fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
455 sum_pers_tot += pers;
460 while (read_next_x(oenv, status, &t, natoms, x, box));
462 gmx_rmpbc_done(gpbc);
476 sum_eed2_tot /= frame;
477 sum_gyro_tot /= frame;
478 sum_pers_tot /= frame;
479 fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
481 fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n",
483 if (opt2bSet("-p", NFILE, fnm))
485 fprintf(stdout, "\nAverage persistence length: %.2f bonds\n",
489 /* Handle printing of internal distances. */
492 fprintf(outi, "@ xaxes scale Logarithmic\n");
495 j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
496 for (i = 0; i < j; i++)
498 intd[i] /= (i + 1) * frame * nmol;
508 xvgr_world(outi, 1, ymin, j, ymax, oenv);
509 for (i = 0; i < j; i++)
511 fprintf(outi, "%d %8.4f\n", i+1, intd[i]);
516 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
517 if (opt2bSet("-v", NFILE, fnm))
519 do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
521 if (opt2bSet("-p", NFILE, fnm))
523 do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");