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
50 #include "gmx_fatal.h"
58 static void gyro_eigen(double **gyr, double *eig, double **eigv, int *ord)
62 jacobi(gyr, DIM, eig, eigv, &nrot);
63 /* Order the eigenvalues */
66 for (d = 0; d < DIM; d++)
68 if (eig[d] > eig[ord[0]])
72 if (eig[d] < eig[ord[2]])
77 for (d = 0; d < DIM; d++)
79 if (ord[0] != d && ord[2] != d)
86 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
87 static void calc_int_dist(double *intd, rvec *x, int i0, int i1)
90 const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
91 int bd; /* Distance between beads */
94 for (bd = 1; bd < ml; bd++)
97 for (ii = i0; ii <= i1 - bd; ii++)
99 d += distance2(x[ii], x[ii+bd]);
106 int gmx_polystat(int argc, char *argv[])
108 const char *desc[] = {
109 "[TT]g_polystat[tt] plots static properties of polymers as a function of time",
110 "and prints the average.[PAR]",
111 "By default it determines the average end-to-end distance and radii",
112 "of gyration of polymers. It asks for an index group and split this",
113 "into molecules. The end-to-end distance is then determined using",
114 "the first and the last atom in the index group for each molecules.",
115 "For the radius of gyration the total and the three principal components",
116 "for the average gyration tensor are written.",
117 "With option [TT]-v[tt] the eigenvectors are written.",
118 "With option [TT]-pc[tt] also the average eigenvalues of the individual",
119 "gyration tensors are written.",
120 "With option [TT]-i[tt] the mean square internal distances are",
122 "With option [TT]-p[tt] the persistence length is determined.",
123 "The chosen index group should consist of atoms that are",
124 "consecutively bonded in the polymer mainchains.",
125 "The persistence length is then determined from the cosine of",
126 "the angles between bonds with an index difference that is even,",
127 "the odd pairs are not used, because straight polymer backbones",
128 "are usually all trans and therefore only every second bond aligns.",
129 "The persistence length is defined as number of bonds where",
130 "the average cos reaches a value of 1/e. This point is determined",
131 "by a linear interpolation of log(<cos>)."
133 static gmx_bool bMW = TRUE, bPC = FALSE;
135 { "-mw", FALSE, etBOOL, {&bMW},
136 "Use the mass weighting for radii of gyration" },
137 { "-pc", FALSE, etBOOL, {&bPC},
138 "Plot average eigenvalues" }
142 { efTPX, NULL, NULL, ffREAD },
143 { efTRX, "-f", NULL, ffREAD },
144 { efNDX, NULL, NULL, ffOPTRD },
145 { efXVG, "-o", "polystat", ffWRITE },
146 { efXVG, "-v", "polyvec", ffOPTWR },
147 { efXVG, "-p", "persist", ffOPTWR },
148 { efXVG, "-i", "intdist", ffOPTWR }
150 #define NFILE asize(fnm)
155 int isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
159 rvec *x, *bond = NULL;
161 int natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = {0};
162 dvec cm, sum_eig = {0, 0, 0};
163 double **gyr, **gyr_all, eig[DIM], **eigv;
164 double sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
166 double *sum_inp = NULL, pers;
167 double *intd, ymax, ymin;
170 FILE *out, *outv, *outp, *outi;
171 const char *leg[8] = {
172 "end to end", "<R\\sg\\N>",
173 "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
174 "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>"
176 char **legp, buf[STRLEN];
177 gmx_rmpbc_t gpbc = NULL;
179 if (!parse_common_args(&argc, argv,
180 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
181 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
187 ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm),
188 NULL, box, &natoms, NULL, NULL, NULL, top);
190 fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
191 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
192 1, &isize, &index, &grpname);
194 snew(molind, top->mols.nr+1);
197 for (i = 0; i < isize; i++)
199 if (i == 0 || index[i] >= top->mols.index[mol+1])
206 while (index[i] >= top->mols.index[mol+1]);
210 nat_min = top->atoms.nr;
212 for (mol = 0; mol < nmol; mol++)
214 nat_min = min(nat_min, molind[mol+1]-molind[mol]);
215 nat_max = max(nat_max, molind[mol+1]-molind[mol]);
217 fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
218 fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n",
221 sprintf(title, "Size of %d polymers", nmol);
222 out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",
224 xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
226 if (opt2bSet("-v", NFILE, fnm))
228 outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
229 output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
231 for (d = 0; d < DIM; d++)
233 for (d2 = 0; d2 < DIM; d2++)
235 sprintf(buf, "eig%d %c", d+1, 'x'+d2);
236 legp[d*DIM+d2] = strdup(buf);
239 xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);
246 if (opt2bSet("-p", NFILE, fnm))
248 outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
249 output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
250 snew(bond, nat_max-1);
251 snew(sum_inp, nat_min/2);
252 snew(ninp, nat_min/2);
259 if (opt2bSet("-i", NFILE, fnm))
261 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
262 "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
263 i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
272 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
277 for (d = 0; d < DIM; d++)
280 snew(gyr_all[d], DIM);
289 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
293 gmx_rmpbc(gpbc, natoms, box, x);
296 for (d = 0; d < DIM; d++)
298 clear_dvec(gyr_all[d]);
308 for (i = 0; i < nat_min/2; i++)
315 for (mol = 0; mol < nmol; mol++)
318 ind1 = molind[mol+1];
320 /* Determine end to end distance */
321 sum_eed2 += distance2(x[index[ind0]], x[index[ind1-1]]);
323 /* Determine internal distances */
326 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
329 /* Determine the radius of gyration */
331 for (d = 0; d < DIM; d++)
337 for (i = ind0; i < ind1; i++)
342 m = top->atoms.atom[a].m;
349 for (d = 0; d < DIM; d++)
352 for (d2 = 0; d2 < DIM; d2++)
354 gyr[d][d2] += m*x[a][d]*x[a][d2];
358 dsvmul(1/mmol, cm, cm);
359 for (d = 0; d < DIM; d++)
361 for (d2 = 0; d2 < DIM; d2++)
363 gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
364 gyr_all[d][d2] += gyr[d][d2];
369 gyro_eigen(gyr, eig, eigv, ord);
370 for (d = 0; d < DIM; d++)
372 sum_eig[d] += eig[ord[d]];
377 for (i = ind0; i < ind1-1; i++)
379 rvec_sub(x[index[i+1]], x[index[i]], bond[i-ind0]);
380 unitv(bond[i-ind0], bond[i-ind0]);
382 for (i = ind0; i < ind1-1; i++)
384 for (j = 0; (i+j < ind1-1 && j < nat_min/2); j += 2)
386 sum_inp[j] += iprod(bond[i-ind0], bond[i-ind0+j]);
395 for (d = 0; d < DIM; d++)
397 for (d2 = 0; d2 < DIM; d2++)
399 gyr_all[d][d2] /= nmol;
401 sum_gyro += gyr_all[d][d];
404 gyro_eigen(gyr_all, eig, eigv, ord);
406 fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
407 t*output_env_get_time_factor(oenv),
408 sqrt(sum_eed2), sqrt(sum_gyro),
409 sqrt(eig[ord[0]]), sqrt(eig[ord[1]]), sqrt(eig[ord[2]]));
412 for (d = 0; d < DIM; d++)
414 fprintf(out, " %8.4f", sqrt(sum_eig[d]/nmol));
421 fprintf(outv, "%10.3f", t*output_env_get_time_factor(oenv));
422 for (d = 0; d < DIM; d++)
424 for (d2 = 0; d2 < DIM; d2++)
426 fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
432 sum_eed2_tot += sum_eed2;
433 sum_gyro_tot += sum_gyro;
438 for (j = 0; j < nat_min/2; j += 2)
440 sum_inp[j] /= ninp[j];
441 if (i == -1 && sum_inp[j] <= exp(-1.0))
452 /* Do linear interpolation on a log scale */
454 + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
456 fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
457 sum_pers_tot += pers;
462 while (read_next_x(oenv, status, &t, x, box));
464 gmx_rmpbc_done(gpbc);
478 sum_eed2_tot /= frame;
479 sum_gyro_tot /= frame;
480 sum_pers_tot /= frame;
481 fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
483 fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n",
485 if (opt2bSet("-p", NFILE, fnm))
487 fprintf(stdout, "\nAverage persistence length: %.2f bonds\n",
491 /* Handle printing of internal distances. */
494 fprintf(outi, "@ xaxes scale Logarithmic\n");
497 j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
498 for (i = 0; i < j; i++)
500 intd[i] /= (i + 1) * frame * nmol;
510 xvgr_world(outi, 1, ymin, j, ymax, oenv);
511 for (i = 0; i < j; i++)
513 fprintf(outi, "%d %8.4f\n", i+1, intd[i]);
518 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
519 if (opt2bSet("-v", NFILE, fnm))
521 do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
523 if (opt2bSet("-p", NFILE, fnm))
525 do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");