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++) {
68 if (eig[d] > eig[ord[0]])
70 if (eig[d] < eig[ord[2]])
73 for(d=0; d<DIM; d++) {
74 if (ord[0] != d && ord[2] != d)
79 /* Calculate mean square internal distances (Auhl et al., JCP 119, 12718) */
80 static void calc_int_dist(double *intd, rvec *x, int i0, int i1)
83 const int ml = i1 - i0 + 1; /* Number of beads in molecule. */
84 int bd; /* Distance between beads */
87 for (bd = 1; bd < ml; bd++) {
89 for (ii = i0; ii <= i1 - bd; ii++)
90 d += distance2(x[ii], x[ii+bd]);
96 int gmx_polystat(int argc,char *argv[])
98 const char *desc[] = {
99 "[TT]g_polystat[tt] plots static properties of polymers as a function of time",
100 "and prints the average.[PAR]",
101 "By default it determines the average end-to-end distance and radii",
102 "of gyration of polymers. It asks for an index group and split this",
103 "into molecules. The end-to-end distance is then determined using",
104 "the first and the last atom in the index group for each molecules.",
105 "For the radius of gyration the total and the three principal components",
106 "for the average gyration tensor are written.",
107 "With option [TT]-v[tt] the eigenvectors are written.",
108 "With option [TT]-pc[tt] also the average eigenvalues of the individual",
109 "gyration tensors are written.",
110 "With option [TT]-i[tt] the mean square internal distances are",
112 "With option [TT]-p[tt] the persistence length is determined.",
113 "The chosen index group should consist of atoms that are",
114 "consecutively bonded in the polymer mainchains.",
115 "The persistence length is then determined from the cosine of",
116 "the angles between bonds with an index difference that is even,",
117 "the odd pairs are not used, because straight polymer backbones",
118 "are usually all trans and therefore only every second bond aligns.",
119 "The persistence length is defined as number of bonds where",
120 "the average cos reaches a value of 1/e. This point is determined",
121 "by a linear interpolation of log(<cos>)."
123 static gmx_bool bMW = TRUE, bPC = FALSE;
125 { "-mw", FALSE, etBOOL, {&bMW},
126 "Use the mass weighting for radii of gyration" },
127 { "-pc", FALSE, etBOOL, {&bPC},
128 "Plot average eigenvalues" }
132 { efTPX, NULL, NULL, ffREAD },
133 { efTRX, "-f", NULL, ffREAD },
134 { efNDX, NULL, NULL, ffOPTRD },
135 { efXVG, "-o", "polystat", ffWRITE },
136 { efXVG, "-v", "polyvec", ffOPTWR },
137 { efXVG, "-p", "persist", ffOPTWR },
138 { efXVG, "-i", "intdist", ffOPTWR }
140 #define NFILE asize(fnm)
145 int isize,*index,nmol,*molind,mol,nat_min=0,nat_max=0;
151 int natoms,i,j,frame,ind0,ind1,a,d,d2,ord[DIM]={0};
152 dvec cm,sum_eig={0,0,0};
153 double **gyr,**gyr_all,eig[DIM],**eigv;
154 double sum_eed2,sum_eed2_tot,sum_gyro,sum_gyro_tot,sum_pers_tot;
156 double *sum_inp=NULL,pers;
157 double *intd, ymax, ymin;
160 FILE *out,*outv,*outp, *outi;
161 const char *leg[8] = { "end to end", "<R\\sg\\N>",
162 "<R\\sg\\N> eig1", "<R\\sg\\N> eig2", "<R\\sg\\N> eig3",
163 "<R\\sg\\N eig1>", "<R\\sg\\N eig2>", "<R\\sg\\N eig3>" };
164 char **legp,buf[STRLEN];
165 gmx_rmpbc_t gpbc=NULL;
167 CopyRight(stderr,argv[0]);
168 parse_common_args(&argc,argv,
169 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
170 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
173 ePBC = read_tpx_top(ftp2fn(efTPX,NFILE,fnm),
174 NULL,box,&natoms,NULL,NULL,NULL,top);
176 fprintf(stderr,"Select a group of polymer mainchain atoms:\n");
177 get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),
178 1,&isize,&index,&grpname);
180 snew(molind,top->mols.nr+1);
183 for(i=0; i<isize; i++) {
184 if (i == 0 || index[i] >= top->mols.index[mol+1]) {
188 } while (index[i] >= top->mols.index[mol+1]);
192 nat_min = top->atoms.nr;
194 for(mol=0; mol<nmol; mol++) {
195 nat_min = min(nat_min,molind[mol+1]-molind[mol]);
196 nat_max = max(nat_max,molind[mol+1]-molind[mol]);
198 fprintf(stderr,"Group %s consists of %d molecules\n",grpname,nmol);
199 fprintf(stderr,"Group size per molecule, min: %d atoms, max %d atoms\n",
202 sprintf(title,"Size of %d polymers",nmol);
203 out = xvgropen(opt2fn("-o",NFILE,fnm),title,output_env_get_xvgr_tlabel(oenv),"(nm)",
205 xvgr_legend(out,bPC ? 8 : 5,leg,oenv);
207 if (opt2bSet("-v",NFILE,fnm)) {
208 outv = xvgropen(opt2fn("-v",NFILE,fnm),"Principal components",
209 output_env_get_xvgr_tlabel(oenv),"(nm)",oenv);
211 for(d=0; d<DIM; d++) {
212 for(d2=0; d2<DIM; d2++) {
213 sprintf(buf,"eig%d %c",d+1,'x'+d2);
214 legp[d*DIM+d2] = strdup(buf);
217 xvgr_legend(outv,DIM*DIM,(const char**)legp,oenv);
222 if (opt2bSet("-p",NFILE,fnm)) {
223 outp = xvgropen(opt2fn("-p",NFILE,fnm),"Persistence length",
224 output_env_get_xvgr_tlabel(oenv),"bonds",oenv);
225 snew(bond,nat_max-1);
226 snew(sum_inp,nat_min/2);
227 snew(ninp,nat_min/2);
232 if (opt2bSet("-i", NFILE, fnm)) {
233 outi = xvgropen(opt2fn("-i", NFILE, fnm), "Internal distances",
234 "n", "<R\\S2\\N(n)>/n (nm\\S2\\N)",oenv);
235 i = index[molind[1]-1] - index[molind[0]]; /* Length of polymer -1 */
242 natoms = read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
247 for(d=0; d<DIM; d++) {
249 snew(gyr_all[d],DIM);
258 gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
261 gmx_rmpbc(gpbc,natoms,box,x);
265 clear_dvec(gyr_all[d]);
271 for(i=0; i<nat_min/2; i++) {
277 for(mol=0; mol<nmol; mol++) {
279 ind1 = molind[mol+1];
281 /* Determine end to end distance */
282 sum_eed2 += distance2(x[index[ind0]],x[index[ind1-1]]);
284 /* Determine internal distances */
286 calc_int_dist(intd, x, index[ind0], index[ind1-1]);
289 /* Determine the radius of gyration */
295 for(i=ind0; i<ind1; i++) {
298 m = top->atoms.atom[a].m;
302 for(d=0; d<DIM; d++) {
304 for(d2=0; d2<DIM; d2++)
305 gyr[d][d2] += m*x[a][d]*x[a][d2];
308 dsvmul(1/mmol,cm,cm);
309 for(d=0; d<DIM; d++) {
310 for(d2=0; d2<DIM; d2++) {
311 gyr[d][d2] = gyr[d][d2]/mmol - cm[d]*cm[d2];
312 gyr_all[d][d2] += gyr[d][d2];
316 gyro_eigen(gyr,eig,eigv,ord);
318 sum_eig[d] += eig[ord[d]];
321 for(i=ind0; i<ind1-1; i++) {
322 rvec_sub(x[index[i+1]],x[index[i]],bond[i-ind0]);
323 unitv(bond[i-ind0],bond[i-ind0]);
325 for(i=ind0; i<ind1-1; i++) {
326 for(j=0; (i+j<ind1-1 && j<nat_min/2); j+=2) {
327 sum_inp[j] += iprod(bond[i-ind0],bond[i-ind0+j]);
336 for(d=0; d<DIM; d++) {
337 for(d2=0; d2<DIM; d2++)
338 gyr_all[d][d2] /= nmol;
339 sum_gyro += gyr_all[d][d];
342 gyro_eigen(gyr_all,eig,eigv,ord);
344 fprintf(out,"%10.3f %8.4f %8.4f %8.4f %8.4f %8.4f",
345 t*output_env_get_time_factor(oenv),
346 sqrt(sum_eed2),sqrt(sum_gyro),
347 sqrt(eig[ord[0]]),sqrt(eig[ord[1]]),sqrt(eig[ord[2]]));
350 fprintf(out," %8.4f",sqrt(sum_eig[d]/nmol));
355 fprintf(outv,"%10.3f",t*output_env_get_time_factor(oenv));
356 for(d=0; d<DIM; d++) {
357 for(d2=0; d2<DIM; d2++)
358 fprintf(outv," %6.3f",eigv[ord[d]][d2]);
363 sum_eed2_tot += sum_eed2;
364 sum_gyro_tot += sum_gyro;
368 for(j=0; j<nat_min/2; j+=2) {
369 sum_inp[j] /= ninp[j];
370 if (i == -1 && sum_inp[j] <= exp(-1.0))
376 /* Do linear interpolation on a log scale */
378 + 2*(log(sum_inp[i-2]) + 1)/(log(sum_inp[i-2]) - log(sum_inp[i]));
380 fprintf(outp,"%10.3f %8.4f\n",t*output_env_get_time_factor(oenv),pers);
381 sum_pers_tot += pers;
385 } while (read_next_x(oenv,status,&t,natoms,x,box));
387 gmx_rmpbc_done(gpbc);
397 sum_eed2_tot /= frame;
398 sum_gyro_tot /= frame;
399 sum_pers_tot /= frame;
400 fprintf(stdout,"\nAverage end to end distance: %.3f (nm)\n",
402 fprintf(stdout,"\nAverage radius of gyration: %.3f (nm)\n",
404 if (opt2bSet("-p",NFILE,fnm))
405 fprintf(stdout,"\nAverage persistence length: %.2f bonds\n",
408 /* Handle printing of internal distances. */
410 fprintf(outi, "@ xaxes scale Logarithmic\n");
413 j = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
414 for (i = 0; i < j; i++) {
415 intd[i] /= (i + 1) * frame * nmol;
421 xvgr_world(outi, 1, ymin, j, ymax,oenv);
422 for (i = 0; i < j; i++) {
423 fprintf(outi, "%d %8.4f\n", i+1, intd[i]);
428 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
429 if (opt2bSet("-v",NFILE,fnm))
430 do_view(oenv,opt2fn("-v",NFILE,fnm),"-nxy");
431 if (opt2bSet("-p",NFILE,fnm))
432 do_view(oenv,opt2fn("-p",NFILE,fnm),"-nxy");