1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2017,2018 by the GROMACS development team.
7 * Copyright (c) 2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #include "gmxpre.h"
39
40 #include <cmath>
41 #include <cstring>
42
43 #include <algorithm>
44
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/fileio/tpxio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/linearalgebra/nrjac.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/pbcutil/rmpbc.h"
54 #include "gromacs/topology/index.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/arraysize.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/smalloc.h"
60
gyro_eigen(double ** gyr,double * eig,double ** eigv,int * ord)61 static void gyro_eigen(double** gyr, double* eig, double** eigv, int* ord)
62 {
63 int nrot, d;
64
65 jacobi(gyr, DIM, eig, eigv, &nrot);
66 /* Order the eigenvalues */
67 ord[0] = 0;
68 ord[2] = 2;
69 for (d = 0; d < DIM; d++)
70 {
71 if (eig[d] > eig[ord[0]])
72 {
73 ord[0] = d;
74 }
75 if (eig[d] < eig[ord[2]])
76 {
77 ord[2] = d;
78 }
79 }
80 for (d = 0; d < DIM; d++)
81 {
82 if (ord[0] != d && ord[2] != d)
83 {
84 ord[1] = d;
85 }
86 }
87 }
88
89 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
calc_int_dist(double * intd,rvec * x,int i0,int i1)90 static void calc_int_dist(double* intd, rvec* x, int i0, int i1)
91 {
92 int ii;
93 const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
94 int bd; /* Distance between beads */
95 double d;
96
97 for (bd = 1; bd < ml; bd++)
98 {
99 d = 0.;
100 for (ii = i0; ii <= i1 - bd; ii++)
101 {
102 d += distance2(x[ii], x[ii + bd]);
103 }
104 d /= ml - bd;
105 intd[bd - 1] += d;
106 }
107 }
108
gmx_polystat(int argc,char * argv[])109 int gmx_polystat(int argc, char* argv[])
110 {
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",
124 "written.[PAR]",
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>[log]."
135 };
136 static gmx_bool bMW = TRUE, bPC = FALSE;
137 t_pargs pa[] = {
138 { "-mw", FALSE, etBOOL, { &bMW }, "Use the mass weighting for radii of gyration" },
139 { "-pc", FALSE, etBOOL, { &bPC }, "Plot average eigenvalues" }
140 };
141
142 t_filenm fnm[] = { { efTPR, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
143 { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, "-o", "polystat", ffWRITE },
144 { efXVG, "-v", "polyvec", ffOPTWR }, { efXVG, "-p", "persist", ffOPTWR },
145 { efXVG, "-i", "intdist", ffOPTWR } };
146 #define NFILE asize(fnm)
147
148 t_topology* top;
149 gmx_output_env_t* oenv;
150 PbcType pbcType;
151 int isize, *index, nmol, *molind, mol, nat_min = 0, nat_max = 0;
152 char* grpname;
153 t_trxstatus* status;
154 real t;
155 rvec * x, *bond = nullptr;
156 matrix box;
157 int natoms, i, j, frame, ind0, ind1, a, d, d2, ord[DIM] = { 0 };
158 dvec cm, sum_eig = { 0, 0, 0 };
159 double ** gyr, **gyr_all, eig[DIM], **eigv;
160 double sum_eed2, sum_eed2_tot, sum_gyro, sum_gyro_tot, sum_pers_tot;
161 int* ninp = nullptr;
162 double * sum_inp = nullptr, pers;
163 double * intd, ymax, ymin;
164 double mmol, m;
165 char title[STRLEN];
166 FILE * out, *outv, *outp, *outi;
167 const char* leg[8] = { "end to end", "<R\\sg\\N>", "<R\\sg\\N> eig1",
168 "<R\\sg\\N> eig2", "<R\\sg\\N> eig3", "<R\\sg\\N eig1>",
169 "<R\\sg\\N eig2>", "<R\\sg\\N eig3>" };
170 char ** legp, buf[STRLEN];
171 gmx_rmpbc_t gpbc = nullptr;
172
173 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm,
174 asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
175 {
176 return 0;
177 }
178
179 snew(top, 1);
180 pbcType = read_tpx_top(ftp2fn(efTPR, NFILE, fnm), nullptr, box, &natoms, nullptr, nullptr, top);
181
182 fprintf(stderr, "Select a group of polymer mainchain atoms:\n");
183 get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
184
185 snew(molind, top->mols.nr + 1);
186 nmol = 0;
187 mol = -1;
188 for (i = 0; i < isize; i++)
189 {
190 if (i == 0 || index[i] >= top->mols.index[mol + 1])
191 {
192 molind[nmol++] = i;
193 do
194 {
195 mol++;
196 } while (index[i] >= top->mols.index[mol + 1]);
197 }
198 }
199 molind[nmol] = i;
200 nat_min = top->atoms.nr;
201 nat_max = 0;
202 for (mol = 0; mol < nmol; mol++)
203 {
204 nat_min = std::min(nat_min, molind[mol + 1] - molind[mol]);
205 nat_max = std::max(nat_max, molind[mol + 1] - molind[mol]);
206 }
207 fprintf(stderr, "Group %s consists of %d molecules\n", grpname, nmol);
208 fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms\n", nat_min, nat_max);
209
210 sprintf(title, "Size of %d polymers", nmol);
211 out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
212 xvgr_legend(out, bPC ? 8 : 5, leg, oenv);
213
214 if (opt2bSet("-v", NFILE, fnm))
215 {
216 outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",
217 output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);
218 snew(legp, DIM * DIM);
219 for (d = 0; d < DIM; d++)
220 {
221 for (d2 = 0; d2 < DIM; d2++)
222 {
223 sprintf(buf, "eig%d %c", d + 1, 'x' + d2);
224 legp[d * DIM + d2] = gmx_strdup(buf);
225 }
226 }
227 xvgr_legend(outv, DIM * DIM, legp, oenv);
228 }
229 else
230 {
231 outv = nullptr;
232 }
233
234 if (opt2bSet("-p", NFILE, fnm))
235 {
236 outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",
237 output_env_get_xvgr_tlabel(oenv), "bonds", oenv);
238 snew(bond, nat_max - 1);
239 snew(sum_inp, nat_min / 2);
240 snew(ninp, nat_min / 2);
241 }
242 else
243 {
244 outp = nullptr;
245 }
246
247 if (opt2bSet("-i", NFILE, fnm))
248 {
249 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances", "n",
250 "<R\\S2\\N(n)>/n (nm\\S2\\N)", oenv);
251 i = index[molind[1] - 1] - index[molind[0]]; /* Length of polymer -1 */
252 snew(intd, i);
253 }
254 else
255 {
256 intd = nullptr;
257 outi = nullptr;
258 }
259
260 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
261
262 snew(gyr, DIM);
263 snew(gyr_all, DIM);
264 snew(eigv, DIM);
265 for (d = 0; d < DIM; d++)
266 {
267 snew(gyr[d], DIM);
268 snew(gyr_all[d], DIM);
269 snew(eigv[d], DIM);
270 }
271
272 frame = 0;
273 sum_eed2_tot = 0;
274 sum_gyro_tot = 0;
275 sum_pers_tot = 0;
276
277 gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
278
279 do
280 {
281 gmx_rmpbc(gpbc, natoms, box, x);
282
283 sum_eed2 = 0;
284 for (d = 0; d < DIM; d++)
285 {
286 clear_dvec(gyr_all[d]);
287 }
288
289 if (bPC)
290 {
291 clear_dvec(sum_eig);
292 }
293
294 if (outp)
295 {
296 for (i = 0; i < nat_min / 2; i++)
297 {
298 sum_inp[i] = 0;
299 ninp[i] = 0;
300 }
301 }
302
303 for (mol = 0; mol < nmol; mol++)
304 {
305 ind0 = molind[mol];
306 ind1 = molind[mol + 1];
307
308 /* Determine end to end distance */
309 sum_eed2 += distance2(x[index[ind0]], x[index[ind1 - 1]]);
310
311 /* Determine internal distances */
312 if (outi)
313 {
314 calc_int_dist(intd, x, index[ind0], index[ind1 - 1]);
315 }
316
317 /* Determine the radius of gyration */
318 clear_dvec(cm);
319 for (d = 0; d < DIM; d++)
320 {
321 clear_dvec(gyr[d]);
322 }
323 mmol = 0;
324
325 for (i = ind0; i < ind1; i++)
326 {
327 a = index[i];
328 if (bMW)
329 {
330 m = top->atoms.atom[a].m;
331 }
332 else
333 {
334 m = 1;
335 }
336 mmol += m;
337 for (d = 0; d < DIM; d++)
338 {
339 cm[d] += m * x[a][d];
340 for (d2 = 0; d2 < DIM; d2++)
341 {
342 gyr[d][d2] += m * x[a][d] * x[a][d2];
343 }
344 }
345 }
346 dsvmul(1 / mmol, cm, cm);
347 for (d = 0; d < DIM; d++)
348 {
349 for (d2 = 0; d2 < DIM; d2++)
350 {
351 gyr[d][d2] = gyr[d][d2] / mmol - cm[d] * cm[d2];
352 gyr_all[d][d2] += gyr[d][d2];
353 }
354 }
355 if (bPC)
356 {
357 gyro_eigen(gyr, eig, eigv, ord);
358 for (d = 0; d < DIM; d++)
359 {
360 sum_eig[d] += eig[ord[d]];
361 }
362 }
363 if (outp)
364 {
365 for (i = ind0; i < ind1 - 1; i++)
366 {
367 rvec_sub(x[index[i + 1]], x[index[i]], bond[i - ind0]);
368 unitv(bond[i - ind0], bond[i - ind0]);
369 }
370 for (i = ind0; i < ind1 - 1; i++)
371 {
372 for (j = 0; (i + j < ind1 - 1 && j < nat_min / 2); j += 2)
373 {
374 sum_inp[j] += iprod(bond[i - ind0], bond[i - ind0 + j]);
375 ninp[j]++;
376 }
377 }
378 }
379 }
380 sum_eed2 /= nmol;
381
382 sum_gyro = 0;
383 for (d = 0; d < DIM; d++)
384 {
385 for (d2 = 0; d2 < DIM; d2++)
386 {
387 gyr_all[d][d2] /= nmol;
388 }
389 sum_gyro += gyr_all[d][d];
390 }
391
392 gyro_eigen(gyr_all, eig, eigv, ord);
393
394 fprintf(out, "%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f", t * output_env_get_time_factor(oenv),
395 std::sqrt(sum_eed2), sqrt(sum_gyro), std::sqrt(eig[ord[0]]), std::sqrt(eig[ord[1]]),
396 std::sqrt(eig[ord[2]]));
397 if (bPC)
398 {
399 for (d = 0; d < DIM; d++)
400 {
401 fprintf(out, " %8.4f", std::sqrt(sum_eig[d] / nmol));
402 }
403 }
404 fprintf(out, "\n");
405
406 if (outv)
407 {
408 fprintf(outv, "%10.3f", t * output_env_get_time_factor(oenv));
409 for (d = 0; d < DIM; d++)
410 {
411 for (d2 = 0; d2 < DIM; d2++)
412 {
413 fprintf(outv, " %6.3f", eigv[ord[d]][d2]);
414 }
415 }
416 fprintf(outv, "\n");
417 }
418
419 sum_eed2_tot += sum_eed2;
420 sum_gyro_tot += sum_gyro;
421
422 if (outp)
423 {
424 i = -1;
425 for (j = 0; j < nat_min / 2; j += 2)
426 {
427 sum_inp[j] /= ninp[j];
428 if (i == -1 && sum_inp[j] <= std::exp(-1.0))
429 {
430 i = j;
431 }
432 }
433 if (i == -1)
434 {
435 pers = j;
436 }
437 else
438 {
439 /* Do linear interpolation on a log scale */
440 pers = i - 2.0
441 + 2.0 * (std::log(sum_inp[i - 2]) + 1.0)
442 / (std::log(sum_inp[i - 2]) - std::log(sum_inp[i]));
443 }
444 fprintf(outp, "%10.3f %8.4f\n", t * output_env_get_time_factor(oenv), pers);
445 sum_pers_tot += pers;
446 }
447
448 frame++;
449 } while (read_next_x(oenv, status, &t, x, box));
450
451 gmx_rmpbc_done(gpbc);
452
453 close_trx(status);
454
455 xvgrclose(out);
456 if (outv)
457 {
458 xvgrclose(outv);
459 }
460 if (outp)
461 {
462 xvgrclose(outp);
463 }
464
465 sum_eed2_tot /= frame;
466 sum_gyro_tot /= frame;
467 sum_pers_tot /= frame;
468 fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n", std::sqrt(sum_eed2_tot));
469 fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n", std::sqrt(sum_gyro_tot));
470 if (opt2bSet("-p", NFILE, fnm))
471 {
472 fprintf(stdout, "\nAverage persistence length: %.2f bonds\n", sum_pers_tot);
473 }
474
475 /* Handle printing of internal distances. */
476 if (outi)
477 {
478 if (output_env_get_print_xvgr_codes(oenv))
479 {
480 fprintf(outi, "@ xaxes scale Logarithmic\n");
481 }
482 ymax = -1;
483 ymin = 1e300;
484 j = index[molind[1] - 1] - index[molind[0]]; /* Polymer length -1. */
485 for (i = 0; i < j; i++)
486 {
487 intd[i] /= (i + 1) * frame * nmol;
488 if (intd[i] > ymax)
489 {
490 ymax = intd[i];
491 }
492 if (intd[i] < ymin)
493 {
494 ymin = intd[i];
495 }
496 }
497 xvgr_world(outi, 1, ymin, j, ymax, oenv);
498 for (i = 0; i < j; i++)
499 {
500 fprintf(outi, "%d %8.4f\n", i + 1, intd[i]);
501 }
502 xvgrclose(outi);
503 }
504
505 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
506 if (opt2bSet("-v", NFILE, fnm))
507 {
508 do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
509 }
510 if (opt2bSet("-p", NFILE, fnm))
511 {
512 do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");
513 }
514
515 return 0;
516 }
517