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 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/utility/futil.h"
51 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/pbcutil/rmpbc.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/fileio/trxio.h"
61 #include "gromacs/utility/fatalerror.h"
63 real
calc_gyro(rvec x
[], int gnx
, atom_id index
[], t_atom atom
[], real tm
,
64 rvec gvec
, rvec d
, gmx_bool bQ
, gmx_bool bRot
, gmx_bool bMOI
, matrix trans
)
67 real gyro
, dx2
, m0
, Itot
;
72 principal_comp(gnx
, index
, atom
, x
, trans
, d
);
78 for (m
= 0; (m
< DIM
); m
++)
83 pr_rvecs(stderr
, 0, "trans", trans
, DIM
);
85 /* rotate_atoms(gnx,index,x,trans); */
88 for (i
= 0; (i
< gnx
); i
++)
93 m0
= fabs(atom
[ii
].q
);
99 for (m
= 0; (m
< DIM
); m
++)
101 dx2
= x
[ii
][m
]*x
[ii
][m
];
105 gyro
= comp
[XX
]+comp
[YY
]+comp
[ZZ
];
107 for (m
= 0; (m
< DIM
); m
++)
109 gvec
[m
] = sqrt((gyro
-comp
[m
])/tm
);
112 return sqrt(gyro
/tm
);
115 void calc_gyro_z(rvec x
[], matrix box
,
116 int gnx
, atom_id index
[], t_atom atom
[],
117 int nz
, real time
, FILE *out
)
119 static dvec
*inertia
= NULL
;
120 static double *tm
= NULL
;
122 real zf
, w
, sdet
, e1
, e2
;
130 for (i
= 0; i
< nz
; i
++)
132 clear_dvec(inertia
[i
]);
136 for (i
= 0; (i
< gnx
); i
++)
139 zf
= nz
*x
[ii
][ZZ
]/box
[ZZ
][ZZ
];
148 for (j
= 0; j
< 2; j
++)
155 w
= atom
[ii
].m
*(1 + cos(M_PI
*(zf
- zi
)));
156 inertia
[zi
][0] += w
*sqr(x
[ii
][YY
]);
157 inertia
[zi
][1] += w
*sqr(x
[ii
][XX
]);
158 inertia
[zi
][2] -= w
*x
[ii
][XX
]*x
[ii
][YY
];
162 fprintf(out
, "%10g", time
);
163 for (j
= 0; j
< nz
; j
++)
165 for (i
= 0; i
< 3; i
++)
167 inertia
[j
][i
] /= tm
[j
];
169 sdet
= sqrt(sqr(inertia
[j
][0] - inertia
[j
][1]) + 4*sqr(inertia
[j
][2]));
170 e1
= sqrt(0.5*(inertia
[j
][0] + inertia
[j
][1] + sdet
));
171 e2
= sqrt(0.5*(inertia
[j
][0] + inertia
[j
][1] - sdet
));
172 fprintf(out
, " %5.3f %5.3f", e1
, e2
);
177 int gmx_gyrate(int argc
, char *argv
[])
179 const char *desc
[] = {
180 "[THISMODULE] computes the radius of gyration of a molecule",
181 "and the radii of gyration about the [IT]x[it]-, [IT]y[it]- and [IT]z[it]-axes,",
182 "as a function of time. The atoms are explicitly mass weighted.[PAR]",
183 "With the [TT]-nmol[tt] option the radius of gyration will be calculated",
184 "for multiple molecules by splitting the analysis group in equally",
186 "With the option [TT]-nz[tt] 2D radii of gyration in the [IT]x-y[it] plane",
187 "of slices along the [IT]z[it]-axis are calculated."
189 static int nmol
= 1, nz
= 0;
190 static gmx_bool bQ
= FALSE
, bRot
= FALSE
, bMOI
= FALSE
;
192 { "-nmol", FALSE
, etINT
, {&nmol
},
193 "The number of molecules to analyze" },
194 { "-q", FALSE
, etBOOL
, {&bQ
},
195 "Use absolute value of the charge of an atom as weighting factor instead of mass" },
196 { "-p", FALSE
, etBOOL
, {&bRot
},
197 "Calculate the radii of gyration about the principal axes." },
198 { "-moi", FALSE
, etBOOL
, {&bMOI
},
199 "Calculate the moments of inertia (defined by the principal axes)." },
200 { "-nz", FALSE
, etINT
, {&nz
},
201 "Calculate the 2D radii of gyration of this number of slices along the z-axis" },
208 rvec xcm
, gvec
, gvec1
;
211 real
**moi_trans
= NULL
;
212 int max_moi
= 0, delta_moi
= 100;
213 rvec d
, d1
; /* eigenvalues of inertia tensor */
214 real t
, t0
, tm
, gyro
;
216 char *grpname
, title
[256];
217 int i
, j
, m
, gnx
, nam
, mol
;
220 gmx_rmpbc_t gpbc
= NULL
;
221 const char *leg
[] = { "Rg", "RgX", "RgY", "RgZ" };
222 const char *legI
[] = { "Itot", "I1", "I2", "I3" };
223 #define NLEG asize(leg)
225 { efTRX
, "-f", NULL
, ffREAD
},
226 { efTPS
, NULL
, NULL
, ffREAD
},
227 { efNDX
, NULL
, NULL
, ffOPTRD
},
228 { efXVG
, NULL
, "gyrate", ffWRITE
},
229 { efXVG
, "-acf", "moi-acf", ffOPTWR
},
231 #define NFILE asize(fnm)
236 ppa
= add_acf_pargs(&npargs
, pa
);
238 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
| PCA_BE_NICE
,
239 NFILE
, fnm
, npargs
, ppa
, asize(desc
), desc
, 0, NULL
, &oenv
))
243 bACF
= opt2bSet("-acf", NFILE
, fnm
);
244 if (bACF
&& nmol
!= 1)
246 gmx_fatal(FARGS
, "Can only do acf with nmol=1");
248 bRot
= bRot
|| bMOI
|| bACF
;
255 printf("Will rotate system along principal axes\n");
256 snew(moi_trans
, DIM
);
260 printf("Will print moments of inertia\n");
265 printf("Will print radius normalised by charge\n");
268 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), title
, &top
, &ePBC
, &x
, NULL
, box
, TRUE
);
269 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &gnx
, &index
, &grpname
);
271 if (nmol
> gnx
|| gnx
% nmol
!= 0)
273 gmx_fatal(FARGS
, "The number of atoms in the group (%d) is not a multiple of nmol (%d)", gnx
, nmol
);
277 natoms
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
284 out
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
),
285 "Radius of Charge", "Time (ps)", "Rg (nm)", oenv
);
289 out
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
),
290 "Moments of inertia", "Time (ps)", "I (a.m.u. nm\\S2\\N)", oenv
);
294 out
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
),
295 "Radius of gyration", "Time (ps)", "Rg (nm)", oenv
);
299 xvgr_legend(out
, NLEG
, legI
, oenv
);
305 if (output_env_get_print_xvgr_codes(oenv
))
307 fprintf(out
, "@ subtitle \"Axes are principal component axes\"\n");
310 xvgr_legend(out
, NLEG
, leg
, oenv
);
314 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natoms
);
320 gmx_rmpbc_copy(gpbc
, natoms
, box
, x
, x_s
);
325 for (mol
= 0; mol
< nmol
; mol
++)
327 tm
= sub_xcm(nz
== 0 ? x_s
: x
, nam
, index
+mol
*nam
, top
.atoms
.atom
, xcm
, bQ
);
330 gyro
+= calc_gyro(x_s
, nam
, index
+mol
*nam
, top
.atoms
.atom
,
331 tm
, gvec1
, d1
, bQ
, bRot
, bMOI
, trans
);
335 calc_gyro_z(x
, box
, nam
, index
+mol
*nam
, top
.atoms
.atom
, nz
, t
, out
);
337 rvec_inc(gvec
, gvec1
);
343 svmul(1.0/nmol
, gvec
, gvec
);
344 svmul(1.0/nmol
, d
, d
);
353 max_moi
+= delta_moi
;
354 for (m
= 0; (m
< DIM
); m
++)
356 srenew(moi_trans
[m
], max_moi
*DIM
);
359 for (m
= 0; (m
< DIM
); m
++)
361 copy_rvec(trans
[m
], moi_trans
[m
]+DIM
*j
);
363 fprintf(out
, "%10g %10g %10g %10g %10g\n",
364 t
, gyro
, d
[XX
], d
[YY
], d
[ZZ
]);
368 fprintf(out
, "%10g %10g %10g %10g %10g\n",
369 t
, gyro
, gvec
[XX
], gvec
[YY
], gvec
[ZZ
]);
374 while (read_next_x(oenv
, status
, &t
, x
, box
));
378 gmx_rmpbc_done(gpbc
);
385 int mode
= eacVector
;
387 do_autocorr(opt2fn("-acf", NFILE
, fnm
), oenv
,
388 "Moment of inertia vector ACF",
389 j
, 3, moi_trans
, (t
-t0
)/j
, mode
, FALSE
);
390 do_view(oenv
, opt2fn("-acf", NFILE
, fnm
), "-nxy");
393 do_view(oenv
, ftp2fn(efXVG
, NFILE
, fnm
), "-nxy");