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
63 real calc_gyro(rvec x[],int gnx,atom_id index[],t_atom atom[],real tm,
64 rvec gvec,rvec d,bool bQ,bool bRot,bool bMOI,matrix trans)
67 real gyro,dx2,m0,Itot;
71 principal_comp(gnx,index,atom,x,trans,d);
75 for(m=0; (m<DIM); m++)
78 pr_rvecs(stderr,0,"trans",trans,DIM);
80 /* rotate_atoms(gnx,index,x,trans); */
83 for(i=0; (i<gnx); i++) {
89 for(m=0; (m<DIM); m++) {
90 dx2=x[ii][m]*x[ii][m];
94 gyro=comp[XX]+comp[YY]+comp[ZZ];
96 for(m=0; (m<DIM); m++)
97 gvec[m]=sqrt((gyro-comp[m])/tm);
102 void calc_gyro_z(rvec x[],matrix box,
103 int gnx,atom_id index[],t_atom atom[],
104 int nz,real time,FILE *out)
106 static dvec *inertia=NULL;
107 static double *tm=NULL;
109 real zf,w,sdet,e1,e2;
111 if (inertia == NULL) {
116 for(i=0; i<nz; i++) {
117 clear_dvec(inertia[i]);
121 for(i=0; (i<gnx); i++) {
123 zf = nz*x[ii][ZZ]/box[ZZ][ZZ];
132 w = atom[ii].m*(1 + cos(M_PI*(zf - zi)));
133 inertia[zi][0] += w*sqr(x[ii][YY]);
134 inertia[zi][1] += w*sqr(x[ii][XX]);
135 inertia[zi][2] -= w*x[ii][XX]*x[ii][YY];
139 fprintf(out,"%10g",time);
140 for(j=0; j<nz; j++) {
142 inertia[j][i] /= tm[j];
143 sdet = sqrt(sqr(inertia[j][0] - inertia[j][1]) + 4*sqr(inertia[j][2]));
144 e1 = sqrt(0.5*(inertia[j][0] + inertia[j][1] + sdet));
145 e2 = sqrt(0.5*(inertia[j][0] + inertia[j][1] - sdet));
146 fprintf(out," %5.3f %5.3f",e1,e2);
151 int gmx_gyrate(int argc,char *argv[])
153 const char *desc[] = {
154 "g_gyrate computes the radius of gyration of a group of atoms",
155 "and the radii of gyration about the x, y and z axes,",
156 "as a function of time. The atoms are explicitly mass weighted.[PAR]",
157 "With the [TT]-nmol[tt] option the radius of gyration will be calculated",
158 "for multiple molecules by splitting the analysis group in equally",
160 "With the option [TT]-nz[tt] 2D radii of gyration in the x-y plane",
161 "of slices along the z-axis are calculated."
163 static int nmol=1,nz=0;
164 static bool bQ=FALSE,bRot=FALSE,bMOI=FALSE;
166 { "-nmol", FALSE, etINT, {&nmol},
167 "The number of molecules to analyze" },
168 { "-q", FALSE, etBOOL, {&bQ},
169 "Use absolute value of the charge of an atom as weighting factor instead of mass" },
170 { "-p", FALSE, etBOOL, {&bRot},
171 "Calculate the radii of gyration about the principal axes." },
172 { "-moi", FALSE, etBOOL, {&bMOI},
173 "Calculate the moments of inertia (defined by the principal axes)." },
174 { "-nz", FALSE, etINT, {&nz},
175 "Calculate the 2D radii of gyration of # slices along the z-axis" },
185 real **moi_trans=NULL;
186 int max_moi=0,delta_moi=100;
187 rvec d,d1; /* eigenvalues of inertia tensor */
190 char *grpname,title[256];
191 int i,j,m,gnx,nam,mol;
194 gmx_rmpbc_t gpbc=NULL;
195 const char *leg[] = { "Rg", "RgX", "RgY", "RgZ" };
196 const char *legI[] = { "Itot", "I1", "I2", "I3" };
197 #define NLEG asize(leg)
199 { efTRX, "-f", NULL, ffREAD },
200 { efTPS, NULL, NULL, ffREAD },
201 { efNDX, NULL, NULL, ffOPTRD },
202 { efXVG, NULL, "gyrate", ffWRITE },
203 { efXVG, "-acf", "moi-acf", ffOPTWR },
205 #define NFILE asize(fnm)
209 CopyRight(stderr,argv[0]);
211 ppa = add_acf_pargs(&npargs,pa);
213 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
214 NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
215 bACF = opt2bSet("-acf",NFILE,fnm);
217 gmx_fatal(FARGS,"Can only do acf with nmol=1");
218 bRot = bRot || bMOI || bACF;
224 printf("Will rotate system along principal axes\n");
228 printf("Will print moments of inertia\n");
232 printf("Will print radius normalised by charge\n");
234 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
235 get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
237 if (nmol > gnx || gnx % nmol != 0) {
238 gmx_fatal(FARGS,"The number of atoms in the group (%d) is not a multiple of nmol (%d)",gnx,nmol);
242 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
248 out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
249 "Radius of Charge","Time (ps)","Rg (nm)",oenv);
251 out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
252 "Moments of inertia","Time (ps)","I (a.m.u. nm\\S2\\N)",oenv);
254 out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
255 "Radius of gyration","Time (ps)","Rg (nm)",oenv);
257 xvgr_legend(out,NLEG,legI,oenv);
260 if (output_env_get_print_xvgr_codes(oenv))
261 fprintf(out,"@ subtitle \"Axes are principal component axes\"\n");
262 xvgr_legend(out,NLEG,leg,oenv);
265 gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
268 gmx_rmpbc_copy(gpbc,natoms,box,x,x_s);
272 for(mol=0; mol<nmol; mol++) {
273 tm = sub_xcm(nz==0?x_s:x,nam,index+mol*nam,top.atoms.atom,xcm,bQ);
275 gyro += calc_gyro(x_s,nam,index+mol*nam,top.atoms.atom,
276 tm,gvec1,d1,bQ,bRot,bMOI,trans);
278 calc_gyro_z(x,box,nam,index+mol*nam,top.atoms.atom,nz,t,out);
279 rvec_inc(gvec,gvec1);
284 svmul(1.0/nmol,gvec,gvec);
291 max_moi += delta_moi;
292 for(m=0; (m<DIM); m++)
293 srenew(moi_trans[m],max_moi*DIM);
295 for(m=0; (m<DIM); m++)
296 copy_rvec(trans[m],moi_trans[m]+DIM*j);
297 fprintf(out,"%10g %10g %10g %10g %10g\n",
298 t,gyro,d[XX],d[YY],d[ZZ]); }
300 fprintf(out,"%10g %10g %10g %10g %10g\n",
301 t,gyro,gvec[XX],gvec[YY],gvec[ZZ]); }
304 } while(read_next_x(oenv,status,&t,natoms,x,box));
307 gmx_rmpbc_done(gpbc);
312 int mode = eacVector;
314 do_autocorr(opt2fn("-acf",NFILE,fnm),oenv,
315 "Moment of inertia vector ACF",
316 j,3,moi_trans,(t-t0)/j,mode,FALSE);
317 do_view(oenv,opt2fn("-acf",NFILE,fnm),"-nxy");
320 do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");