2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
66 real calc_gyro(rvec x[],int gnx,atom_id index[],t_atom atom[],real tm,
67 rvec gvec,rvec d,gmx_bool bQ,gmx_bool bRot,gmx_bool bMOI,matrix trans)
70 real gyro,dx2,m0,Itot;
74 principal_comp(gnx,index,atom,x,trans,d);
78 for(m=0; (m<DIM); m++)
81 pr_rvecs(stderr,0,"trans",trans,DIM);
83 /* rotate_atoms(gnx,index,x,trans); */
86 for(i=0; (i<gnx); i++) {
92 for(m=0; (m<DIM); m++) {
93 dx2=x[ii][m]*x[ii][m];
97 gyro=comp[XX]+comp[YY]+comp[ZZ];
99 for(m=0; (m<DIM); m++)
100 gvec[m]=sqrt((gyro-comp[m])/tm);
102 return sqrt(gyro/tm);
105 void calc_gyro_z(rvec x[],matrix box,
106 int gnx,atom_id index[],t_atom atom[],
107 int nz,real time,FILE *out)
109 static dvec *inertia=NULL;
110 static double *tm=NULL;
112 real zf,w,sdet,e1,e2;
114 if (inertia == NULL) {
119 for(i=0; i<nz; i++) {
120 clear_dvec(inertia[i]);
124 for(i=0; (i<gnx); i++) {
126 zf = nz*x[ii][ZZ]/box[ZZ][ZZ];
135 w = atom[ii].m*(1 + cos(M_PI*(zf - zi)));
136 inertia[zi][0] += w*sqr(x[ii][YY]);
137 inertia[zi][1] += w*sqr(x[ii][XX]);
138 inertia[zi][2] -= w*x[ii][XX]*x[ii][YY];
142 fprintf(out,"%10g",time);
143 for(j=0; j<nz; j++) {
145 inertia[j][i] /= tm[j];
146 sdet = sqrt(sqr(inertia[j][0] - inertia[j][1]) + 4*sqr(inertia[j][2]));
147 e1 = sqrt(0.5*(inertia[j][0] + inertia[j][1] + sdet));
148 e2 = sqrt(0.5*(inertia[j][0] + inertia[j][1] - sdet));
149 fprintf(out," %5.3f %5.3f",e1,e2);
154 int gmx_gyrate(int argc,char *argv[])
156 const char *desc[] = {
157 "[TT]g_gyrate[tt] computes the radius of gyration of a group of atoms",
158 "and the radii of gyration about the [IT]x[it]-, [IT]y[it]- and [IT]z[it]-axes,",
159 "as a function of time. The atoms are explicitly mass weighted.[PAR]",
160 "With the [TT]-nmol[tt] option the radius of gyration will be calculated",
161 "for multiple molecules by splitting the analysis group in equally",
163 "With the option [TT]-nz[tt] 2D radii of gyration in the [IT]x-y[it] plane",
164 "of slices along the [IT]z[it]-axis are calculated."
166 static int nmol=1,nz=0;
167 static gmx_bool bQ=FALSE,bRot=FALSE,bMOI=FALSE;
169 { "-nmol", FALSE, etINT, {&nmol},
170 "The number of molecules to analyze" },
171 { "-q", FALSE, etBOOL, {&bQ},
172 "Use absolute value of the charge of an atom as weighting factor instead of mass" },
173 { "-p", FALSE, etBOOL, {&bRot},
174 "Calculate the radii of gyration about the principal axes." },
175 { "-moi", FALSE, etBOOL, {&bMOI},
176 "Calculate the moments of inertia (defined by the principal axes)." },
177 { "-nz", FALSE, etINT, {&nz},
178 "Calculate the 2D radii of gyration of this number of slices along the z-axis" },
188 real **moi_trans=NULL;
189 int max_moi=0,delta_moi=100;
190 rvec d,d1; /* eigenvalues of inertia tensor */
193 char *grpname,title[256];
194 int i,j,m,gnx,nam,mol;
197 gmx_rmpbc_t gpbc=NULL;
198 const char *leg[] = { "Rg", "RgX", "RgY", "RgZ" };
199 const char *legI[] = { "Itot", "I1", "I2", "I3" };
200 #define NLEG asize(leg)
202 { efTRX, "-f", NULL, ffREAD },
203 { efTPS, NULL, NULL, ffREAD },
204 { efNDX, NULL, NULL, ffOPTRD },
205 { efXVG, NULL, "gyrate", ffWRITE },
206 { efXVG, "-acf", "moi-acf", ffOPTWR },
208 #define NFILE asize(fnm)
212 CopyRight(stderr,argv[0]);
214 ppa = add_acf_pargs(&npargs,pa);
216 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
217 NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
218 bACF = opt2bSet("-acf",NFILE,fnm);
220 gmx_fatal(FARGS,"Can only do acf with nmol=1");
221 bRot = bRot || bMOI || bACF;
227 printf("Will rotate system along principal axes\n");
231 printf("Will print moments of inertia\n");
235 printf("Will print radius normalised by charge\n");
237 read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,TRUE);
238 get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
240 if (nmol > gnx || gnx % nmol != 0) {
241 gmx_fatal(FARGS,"The number of atoms in the group (%d) is not a multiple of nmol (%d)",gnx,nmol);
245 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
251 out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
252 "Radius of Charge","Time (ps)","Rg (nm)",oenv);
254 out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
255 "Moments of inertia","Time (ps)","I (a.m.u. nm\\S2\\N)",oenv);
257 out=xvgropen(ftp2fn(efXVG,NFILE,fnm),
258 "Radius of gyration","Time (ps)","Rg (nm)",oenv);
260 xvgr_legend(out,NLEG,legI,oenv);
263 if (output_env_get_print_xvgr_codes(oenv))
264 fprintf(out,"@ subtitle \"Axes are principal component axes\"\n");
265 xvgr_legend(out,NLEG,leg,oenv);
268 gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
271 gmx_rmpbc_copy(gpbc,natoms,box,x,x_s);
275 for(mol=0; mol<nmol; mol++) {
276 tm = sub_xcm(nz==0?x_s:x,nam,index+mol*nam,top.atoms.atom,xcm,bQ);
278 gyro += calc_gyro(x_s,nam,index+mol*nam,top.atoms.atom,
279 tm,gvec1,d1,bQ,bRot,bMOI,trans);
281 calc_gyro_z(x,box,nam,index+mol*nam,top.atoms.atom,nz,t,out);
282 rvec_inc(gvec,gvec1);
287 svmul(1.0/nmol,gvec,gvec);
294 max_moi += delta_moi;
295 for(m=0; (m<DIM); m++)
296 srenew(moi_trans[m],max_moi*DIM);
298 for(m=0; (m<DIM); m++)
299 copy_rvec(trans[m],moi_trans[m]+DIM*j);
300 fprintf(out,"%10g %10g %10g %10g %10g\n",
301 t,gyro,d[XX],d[YY],d[ZZ]); }
303 fprintf(out,"%10g %10g %10g %10g %10g\n",
304 t,gyro,gvec[XX],gvec[YY],gvec[ZZ]); }
307 } while(read_next_x(oenv,status,&t,natoms,x,box));
310 gmx_rmpbc_done(gpbc);
315 int mode = eacVector;
317 do_autocorr(opt2fn("-acf",NFILE,fnm),oenv,
318 "Moment of inertia vector ACF",
319 j,3,moi_trans,(t-t0)/j,mode,FALSE);
320 do_view(oenv,opt2fn("-acf",NFILE,fnm),"-nxy");
323 do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");