Bug Summary

File:gromacs/gmxana/gmx_msd.c
Location:line 497, column 5
Description:Value stored to 'mass' is never read

Annotated Source Code

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, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <string.h>
43
44#include "gromacs/utility/smalloc.h"
45#include "macros.h"
46#include "gromacs/commandline/pargs.h"
47#include "gromacs/math/utilities.h"
48#include "gromacs/utility/futil.h"
49#include "index.h"
50#include "typedefs.h"
51#include "gromacs/fileio/xvgr.h"
52#include "viewit.h"
53#include "gstat.h"
54#include "gromacs/statistics/statistics.h"
55#include "gromacs/fileio/tpxio.h"
56#include "gromacs/fileio/trxio.h"
57#include "pbc.h"
58#include "gromacs/math/vec.h"
59#include "gromacs/fileio/confio.h"
60#include "gmx_ana.h"
61
62#include "gromacs/utility/fatalerror.h"
63
64#define FACTOR1000.0 1000.0 /* Convert nm^2/ps to 10e-5 cm^2/s */
65/* NORMAL = total diffusion coefficient (default). X,Y,Z is diffusion
66 coefficient in X,Y,Z direction. LATERAL is diffusion coefficient in
67 plane perpendicular to axis
68 */
69typedef enum {
70 NOT_USED, NORMAL, X, Y, Z, LATERAL
71} msd_type;
72
73typedef struct {
74 real t0; /* start time and time increment between */
75 real delta_t; /* time between restart points */
76 real beginfit, /* the begin/end time for fits as reals between */
77 endfit; /* 0 and 1 */
78 real dim_factor; /* the dimensionality factor for the diffusion
79 constant */
80 real **data; /* the displacement data. First index is the group
81 number, second is frame number */
82 real *time; /* frame time */
83 real *mass; /* masses for mass-weighted msd */
84 matrix **datam;
85 rvec **x0; /* original positions */
86 rvec *com; /* center of mass correction for each frame */
87 gmx_stats_t **lsq; /* fitting stats for individual molecule msds */
88 msd_type type; /* the type of msd to calculate (lateral, etc.)*/
89 int axis; /* the axis along which to calculate */
90 int ncoords;
91 int nrestart; /* number of restart points */
92 int nmol; /* number of molecules (for bMol) */
93 int nframes; /* number of frames */
94 int nlast;
95 int ngrp; /* number of groups to use for msd calculation */
96 int *n_offs;
97 int **ndata; /* the number of msds (particles/mols) per data
98 point. */
99} t_corr;
100
101typedef real t_calc_func (t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
102 rvec dcom, gmx_bool bTen, matrix mat);
103
104static real thistime(t_corr *curr)
105{
106 return curr->time[curr->nframes];
107}
108
109static gmx_bool in_data(t_corr *curr, int nx00)
110{
111 return curr->nframes-curr->n_offs[nx00];
112}
113
114t_corr *init_corr(int nrgrp, int type, int axis, real dim_factor,
115 int nmol, gmx_bool bTen, gmx_bool bMass, real dt, t_topology *top,
116 real beginfit, real endfit)
117{
118 t_corr *curr;
119 t_atoms *atoms;
120 int i;
121
122 snew(curr, 1)(curr) = save_calloc("curr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 122, (1), sizeof(*(curr)));
123 curr->type = (msd_type)type;
124 curr->axis = axis;
125 curr->ngrp = nrgrp;
126 curr->nrestart = 0;
127 curr->delta_t = dt;
128 curr->beginfit = (1 - 2*GMX_REAL_EPS5.96046448E-08)*beginfit;
129 curr->endfit = (1 + 2*GMX_REAL_EPS5.96046448E-08)*endfit;
130 curr->x0 = NULL((void*)0);
131 curr->n_offs = NULL((void*)0);
132 curr->nframes = 0;
133 curr->nlast = 0;
134 curr->dim_factor = dim_factor;
135
136 snew(curr->ndata, nrgrp)(curr->ndata) = save_calloc("curr->ndata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 136, (nrgrp), sizeof(*(curr->ndata)));
137 snew(curr->data, nrgrp)(curr->data) = save_calloc("curr->data", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 137, (nrgrp), sizeof(*(curr->data)));
138 if (bTen)
139 {
140 snew(curr->datam, nrgrp)(curr->datam) = save_calloc("curr->datam", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 140, (nrgrp), sizeof(*(curr->datam)));
141 }
142 for (i = 0; (i < nrgrp); i++)
143 {
144 curr->ndata[i] = NULL((void*)0);
145 curr->data[i] = NULL((void*)0);
146 if (bTen)
147 {
148 curr->datam[i] = NULL((void*)0);
149 }
150 }
151 curr->time = NULL((void*)0);
152 curr->lsq = NULL((void*)0);
153 curr->nmol = nmol;
154 if (curr->nmol > 0)
155 {
156 snew(curr->mass, curr->nmol)(curr->mass) = save_calloc("curr->mass", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 156, (curr->nmol), sizeof(*(curr->mass)));
157 for (i = 0; i < curr->nmol; i++)
158 {
159 curr->mass[i] = 1;
160 }
161 }
162 else
163 {
164 if (bMass)
165 {
166 atoms = &top->atoms;
167 snew(curr->mass, atoms->nr)(curr->mass) = save_calloc("curr->mass", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 167, (atoms->nr), sizeof(*(curr->mass)));
168 for (i = 0; (i < atoms->nr); i++)
169 {
170 curr->mass[i] = atoms->atom[i].m;
171 }
172 }
173 }
174
175 return curr;
176}
177
178static void corr_print(t_corr *curr, gmx_bool bTen, const char *fn, const char *title,
179 const char *yaxis,
180 real msdtime, real beginfit, real endfit,
181 real *DD, real *SigmaD, char *grpname[],
182 const output_env_t oenv)
183{
184 FILE *out;
185 int i, j;
186
187 out = xvgropen(fn, title, output_env_get_xvgr_tlabel(oenv), yaxis, oenv);
188 if (DD)
189 {
190 fprintf(out, "# MSD gathered over %g %s with %d restarts\n",
191 msdtime, output_env_get_time_unit(oenv), curr->nrestart);
192 fprintf(out, "# Diffusion constants fitted from time %g to %g %s\n",
193 beginfit, endfit, output_env_get_time_unit(oenv));
194 for (i = 0; i < curr->ngrp; i++)
195 {
196 fprintf(out, "# D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)\n",
197 grpname[i], DD[i], SigmaD[i]);
198 }
199 }
200 for (i = 0; i < curr->nframes; i++)
201 {
202 fprintf(out, "%10g", output_env_conv_time(oenv, curr->time[i]));
203 for (j = 0; j < curr->ngrp; j++)
204 {
205 fprintf(out, " %10g", curr->data[j][i]);
206 if (bTen)
207 {
208 fprintf(out, " %10g %10g %10g %10g %10g %10g",
209 curr->datam[j][i][XX0][XX0],
210 curr->datam[j][i][YY1][YY1],
211 curr->datam[j][i][ZZ2][ZZ2],
212 curr->datam[j][i][YY1][XX0],
213 curr->datam[j][i][ZZ2][XX0],
214 curr->datam[j][i][ZZ2][YY1]);
215 }
216 }
217 fprintf(out, "\n");
218 }
219 gmx_ffclose(out);
220}
221
222/* called from corr_loop, to do the main calculations */
223static void calc_corr(t_corr *curr, int nr, int nx, atom_id index[], rvec xc[],
224 gmx_bool bRmCOMM, rvec com, t_calc_func *calc1, gmx_bool bTen)
225{
226 int nx0;
227 real g;
228 matrix mat;
229 rvec dcom;
230
231 /* Check for new starting point */
232 if (curr->nlast < curr->nrestart)
233 {
234 if ((thistime(curr) >= (curr->nlast*curr->delta_t)) && (nr == 0))
235 {
236 memcpy(curr->x0[curr->nlast], xc, curr->ncoords*sizeof(xc[0]));
237 curr->n_offs[curr->nlast] = curr->nframes;
238 copy_rvec(com, curr->com[curr->nlast]);
239 curr->nlast++;
240 }
241 }
242
243 /* nx0 appears to be the number of new starting points,
244 * so for all starting points, call calc1.
245 */
246 for (nx0 = 0; (nx0 < curr->nlast); nx0++)
247 {
248 if (bRmCOMM)
249 {
250 rvec_sub(com, curr->com[nx0], dcom);
251 }
252 else
253 {
254 clear_rvec(dcom);
255 }
256 g = calc1(curr, nx, index, nx0, xc, dcom, bTen, mat);
257#ifdef DEBUG2
258 printf("g[%d]=%g\n", nx0, g);
259#endif
260 curr->data[nr][in_data(curr, nx0)] += g;
261 if (bTen)
262 {
263 m_add(curr->datam[nr][in_data(curr, nx0)], mat,
264 curr->datam[nr][in_data(curr, nx0)]);
265 }
266 curr->ndata[nr][in_data(curr, nx0)]++;
267 }
268}
269
270/* the non-mass-weighted mean-squared displacement calcuation */
271static real calc1_norm(t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
272 rvec dcom, gmx_bool bTen, matrix mat)
273{
274 int i, ix, m, m2;
275 real g, r, r2;
276 rvec rv;
277
278 g = 0.0;
279 clear_mat(mat);
280
281 for (i = 0; (i < nx); i++)
282 {
283 ix = index[i];
284 r2 = 0.0;
285 switch (curr->type)
286 {
287 case NORMAL:
288 for (m = 0; (m < DIM3); m++)
289 {
290 rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
291 r2 += rv[m]*rv[m];
292 if (bTen)
293 {
294 for (m2 = 0; m2 <= m; m2++)
295 {
296 mat[m][m2] += rv[m]*rv[m2];
297 }
298 }
299 }
300 break;
301 case X:
302 case Y:
303 case Z:
304 r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
305 dcom[curr->type-X];
306 r2 += r*r;
307 break;
308 case LATERAL:
309 for (m = 0; (m < DIM3); m++)
310 {
311 if (m != curr->axis)
312 {
313 r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
314 r2 += r*r;
315 }
316 }
317 break;
318 default:
319 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 319, "Error: did not expect option value %d", curr->type);
320 }
321 g += r2;
322 }
323 g /= nx;
324 msmul(mat, 1.0/nx, mat);
325
326 return g;
327}
328
329/* calculate the com of molecules in x and put it into xa */
330static void calc_mol_com(int nmol, int *molindex, t_block *mols, t_atoms *atoms,
331 rvec *x, rvec *xa)
332{
333 int m, mol, i, d;
334 rvec xm;
335 real mass, mtot;
336
337 for (m = 0; m < nmol; m++)
338 {
339 mol = molindex[m];
340 clear_rvec(xm);
341 mtot = 0;
342 for (i = mols->index[mol]; i < mols->index[mol+1]; i++)
343 {
344 mass = atoms->atom[i].m;
345 for (d = 0; d < DIM3; d++)
346 {
347 xm[d] += mass*x[i][d];
348 }
349 mtot += mass;
350 }
351 svmul(1/mtot, xm, xa[m]);
352 }
353}
354
355static real calc_one_mw(t_corr *curr, int ix, int nx0, rvec xc[], real *tm,
356 rvec dcom, gmx_bool bTen, matrix mat)
357{
358 real r2, r, mm;
359 rvec rv;
360 int m, m2;
361
362 mm = curr->mass[ix];
363 if (mm == 0)
364 {
365 return 0;
366 }
367 (*tm) += mm;
368 r2 = 0.0;
369 switch (curr->type)
370 {
371 case NORMAL:
372 for (m = 0; (m < DIM3); m++)
373 {
374 rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
375 r2 += mm*rv[m]*rv[m];
376 if (bTen)
377 {
378 for (m2 = 0; m2 <= m; m2++)
379 {
380 mat[m][m2] += mm*rv[m]*rv[m2];
381 }
382 }
383 }
384 break;
385 case X:
386 case Y:
387 case Z:
388 r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
389 dcom[curr->type-X];
390 r2 = mm*r*r;
391 break;
392 case LATERAL:
393 for (m = 0; (m < DIM3); m++)
394 {
395 if (m != curr->axis)
396 {
397 r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
398 r2 += mm*r*r;
399 }
400 }
401 break;
402 default:
403 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 403, "Options got screwed. Did not expect value %d\n", curr->type);
404 } /* end switch */
405 return r2;
406}
407
408/* the normal, mass-weighted mean-squared displacement calcuation */
409static real calc1_mw(t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
410 rvec dcom, gmx_bool bTen, matrix mat)
411{
412 int i;
413 real g, tm;
414
415 g = tm = 0.0;
416 clear_mat(mat);
417 for (i = 0; (i < nx); i++)
418 {
419 g += calc_one_mw(curr, index[i], nx0, xc, &tm, dcom, bTen, mat);
420 }
421
422 g /= tm;
423 if (bTen)
424 {
425 msmul(mat, 1/tm, mat);
426 }
427
428 return g;
429}
430
431/* prepare the coordinates by removing periodic boundary crossings.
432 gnx = the number of atoms/molecules
433 index = the indices
434 xcur = the current coordinates
435 xprev = the previous coordinates
436 box = the box matrix */
437static void prep_data(gmx_bool bMol, int gnx, atom_id index[],
438 rvec xcur[], rvec xprev[], matrix box)
439{
440 int i, m, ind;
441 rvec hbox;
442
443 /* Remove periodicity */
444 for (m = 0; (m < DIM3); m++)
445 {
446 hbox[m] = 0.5*box[m][m];
447 }
448
449 for (i = 0; (i < gnx); i++)
450 {
451 if (bMol)
452 {
453 ind = i;
454 }
455 else
456 {
457 ind = index[i];
458 }
459
460 for (m = DIM3-1; m >= 0; m--)
461 {
462 if (hbox[m] == 0)
463 {
464 continue;
465 }
466 while (xcur[ind][m]-xprev[ind][m] <= -hbox[m])
467 {
468 rvec_inc(xcur[ind], box[m]);
469 }
470 while (xcur[ind][m]-xprev[ind][m] > hbox[m])
471 {
472 rvec_dec(xcur[ind], box[m]);
473 }
474 }
475 }
476}
477
478/* calculate the center of mass for a group
479 gnx = the number of atoms/molecules
480 index = the indices
481 xcur = the current coordinates
482 xprev = the previous coordinates
483 box = the box matrix
484 atoms = atom data (for mass)
485 com(output) = center of mass */
486static void calc_com(gmx_bool bMol, int gnx, atom_id index[],
487 rvec xcur[], rvec xprev[], matrix box, t_atoms *atoms,
488 rvec com)
489{
490 int i, m, ind;
491 real mass;
492 double tmass;
493 dvec sx;
494
495 clear_dvec(sx);
496 tmass = 0;
497 mass = 1;
Value stored to 'mass' is never read
498
499 prep_data(bMol, gnx, index, xcur, xprev, box);
500 for (i = 0; (i < gnx); i++)
501 {
502 if (bMol)
503 {
504 ind = i;
505 }
506 else
507 {
508 ind = index[i];
509 }
510
511
512 mass = atoms->atom[ind].m;
513 for (m = 0; m < DIM3; m++)
514 {
515 sx[m] += mass*xcur[ind][m];
516 }
517 tmass += mass;
518 }
519 for (m = 0; m < DIM3; m++)
520 {
521 com[m] = sx[m]/tmass;
522 }
523}
524
525
526static real calc1_mol(t_corr *curr, int nx, atom_id gmx_unused__attribute__ ((unused)) index[], int nx0, rvec xc[],
527 rvec dcom, gmx_bool bTen, matrix mat)
528{
529 int i;
530 real g, tm, gtot, tt;
531
532 tt = curr->time[in_data(curr, nx0)];
533 gtot = 0;
534 tm = 0;
535 clear_mat(mat);
536 for (i = 0; (i < nx); i++)
537 {
538 g = calc_one_mw(curr, i, nx0, xc, &tm, dcom, bTen, mat);
539 /* We don't need to normalize as the mass was set to 1 */
540 gtot += g;
541 if (tt >= curr->beginfit && (curr->endfit < 0 || tt <= curr->endfit))
542 {
543 gmx_stats_add_point(curr->lsq[nx0][i], tt, g, 0, 0);
544 }
545 }
546 msmul(mat, 1.0/nx, mat);
547
548 return gtot/nx;
549}
550
551void printmol(t_corr *curr, const char *fn,
552 const char *fn_pdb, int *molindex, t_topology *top,
553 rvec *x, int ePBC, matrix box, const output_env_t oenv)
554{
555#define NDIST100 100
556 FILE *out;
557 gmx_stats_t lsq1;
558 int i, j;
559 real a, b, D, Dav, D2av, VarD, sqrtD, sqrtD_max, scale;
560 t_pdbinfo *pdbinfo = NULL((void*)0);
561 int *mol2a = NULL((void*)0);
562
563 out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D", oenv);
564
565 if (fn_pdb)
566 {
567 if (top->atoms.pdbinfo == NULL((void*)0))
568 {
569 snew(top->atoms.pdbinfo, top->atoms.nr)(top->atoms.pdbinfo) = save_calloc("top->atoms.pdbinfo"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 569, (top->atoms.nr), sizeof(*(top->atoms.pdbinfo)));
570 }
571 pdbinfo = top->atoms.pdbinfo;
572 mol2a = top->mols.index;
573 }
574
575 Dav = D2av = 0;
576 sqrtD_max = 0;
577 for (i = 0; (i < curr->nmol); i++)
578 {
579 lsq1 = gmx_stats_init();
580 for (j = 0; (j < curr->nrestart); j++)
581 {
582 real xx, yy, dx, dy;
583
584 while (gmx_stats_get_point(curr->lsq[j][i], &xx, &yy, &dx, &dy, 0) == estatsOK)
585 {
586 gmx_stats_add_point(lsq1, xx, yy, dx, dy);
587 }
588 }
589 gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0));
590 gmx_stats_done(lsq1);
591 sfree(lsq1)save_free("lsq1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 591, (lsq1));
592 D = a*FACTOR1000.0/curr->dim_factor;
593 if (D < 0)
594 {
595 D = 0;
596 }
597 Dav += D;
598 D2av += sqr(D);
599 fprintf(out, "%10d %10g\n", i, D);
600 if (pdbinfo)
601 {
602 sqrtD = sqrt(D);
603 if (sqrtD > sqrtD_max)
604 {
605 sqrtD_max = sqrtD;
606 }
607 for (j = mol2a[molindex[i]]; j < mol2a[molindex[i]+1]; j++)
608 {
609 pdbinfo[j].bfac = sqrtD;
610 }
611 }
612 }
613 gmx_ffclose(out);
614 do_view(oenv, fn, "-graphtype bar");
615
616 /* Compute variance, stddev and error */
617 Dav /= curr->nmol;
618 D2av /= curr->nmol;
619 VarD = D2av - sqr(Dav);
620 printf("<D> = %.4f Std. Dev. = %.4f Error = %.4f\n",
621 Dav, sqrt(VarD), sqrt(VarD/curr->nmol));
622
623 if (fn_pdb && x)
624 {
625 scale = 1;
626 while (scale*sqrtD_max > 10)
627 {
628 scale *= 0.1;
629 }
630 while (scale*sqrtD_max < 0.1)
631 {
632 scale *= 10;
633 }
634 for (i = 0; i < top->atoms.nr; i++)
635 {
636 pdbinfo[i].bfac *= scale;
637 }
638 write_sto_conf(fn_pdb, "molecular MSD", &top->atoms, x, NULL((void*)0), ePBC, box);
639 }
640}
641
642/* this is the main loop for the correlation type functions
643 * fx and nx are file pointers to things like read_first_x and
644 * read_next_x
645 */
646int corr_loop(t_corr *curr, const char *fn, t_topology *top, int ePBC,
647 gmx_bool bMol, int gnx[], atom_id *index[],
648 t_calc_func *calc1, gmx_bool bTen, int *gnx_com, atom_id *index_com[],
649 real dt, real t_pdb, rvec **x_pdb, matrix box_pdb,
650 const output_env_t oenv)
651{
652 rvec *x[2]; /* the coordinates to read */
653 rvec *xa[2]; /* the coordinates to calculate displacements for */
654 rvec com = {0};
655 real t, t_prev = 0;
656 int natoms, i, j, cur = 0, maxframes = 0;
657 t_trxstatus *status;
658#define prev(1-cur) (1-cur)
659 matrix box;
660 gmx_bool bFirst;
661 gmx_rmpbc_t gpbc = NULL((void*)0);
662
663 natoms = read_first_x(oenv, &status, fn, &curr->t0, &(x[cur]), box);
664#ifdef DEBUG
665 fprintf(stderrstderr, "Read %d atoms for first frame\n", natoms);
666#endif
667 if ((gnx_com != NULL((void*)0)) && natoms < top->atoms.nr)
668 {
669 fprintf(stderrstderr, "WARNING: The trajectory only contains part of the system (%d of %d atoms) and therefore the COM motion of only this part of the system will be removed\n", natoms, top->atoms.nr);
670 }
671
672 snew(x[prev], natoms)(x[(1-cur)]) = save_calloc("x[prev]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 672, (natoms), sizeof(*(x[(1-cur)])));
673
674 if (bMol)
675 {
676 curr->ncoords = curr->nmol;
677 snew(xa[0], curr->ncoords)(xa[0]) = save_calloc("xa[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 677, (curr->ncoords), sizeof(*(xa[0])));
678 snew(xa[1], curr->ncoords)(xa[1]) = save_calloc("xa[1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 678, (curr->ncoords), sizeof(*(xa[1])));
679 }
680 else
681 {
682 curr->ncoords = natoms;
683 xa[0] = x[0];
684 xa[1] = x[1];
685 }
686
687 bFirst = TRUE1;
688 t = curr->t0;
689 if (x_pdb)
690 {
691 *x_pdb = NULL((void*)0);
692 }
693
694 if (bMol)
695 {
696 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
697 }
698
699 /* the loop over all frames */
700 do
701 {
702 if (x_pdb && ((bFirst && t_pdb < t) ||
703 (!bFirst &&
704 t_pdb > t - 0.5*(t - t_prev) &&
705 t_pdb < t + 0.5*(t - t_prev))))
706 {
707 if (*x_pdb == NULL((void*)0))
708 {
709 snew(*x_pdb, natoms)(*x_pdb) = save_calloc("*x_pdb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 709, (natoms), sizeof(*(*x_pdb)));
710 }
711 for (i = 0; i < natoms; i++)
712 {
713 copy_rvec(x[cur][i], (*x_pdb)[i]);
714 }
715 copy_mat(box, box_pdb);
716 }
717
718
719 /* check whether we've reached a restart point */
720 if (bRmod(t, curr->t0, dt)bRmod_fd(t, curr->t0, dt, 0))
721 {
722 curr->nrestart++;
723
724 srenew(curr->x0, curr->nrestart)(curr->x0) = save_realloc("curr->x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 724, (curr->x0), (curr->nrestart), sizeof(*(curr->
x0)));
725 snew(curr->x0[curr->nrestart-1], curr->ncoords)(curr->x0[curr->nrestart-1]) = save_calloc("curr->x0[curr->nrestart-1]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 725, (curr->ncoords), sizeof(*(curr->x0[curr->nrestart
-1])));
726 srenew(curr->com, curr->nrestart)(curr->com) = save_realloc("curr->com", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 726, (curr->com), (curr->nrestart), sizeof(*(curr->
com)));
727 srenew(curr->n_offs, curr->nrestart)(curr->n_offs) = save_realloc("curr->n_offs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 727, (curr->n_offs), (curr->nrestart), sizeof(*(curr->
n_offs)));
728 srenew(curr->lsq, curr->nrestart)(curr->lsq) = save_realloc("curr->lsq", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 728, (curr->lsq), (curr->nrestart), sizeof(*(curr->
lsq)));
729 snew(curr->lsq[curr->nrestart-1], curr->nmol)(curr->lsq[curr->nrestart-1]) = save_calloc("curr->lsq[curr->nrestart-1]"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 729, (curr->nmol), sizeof(*(curr->lsq[curr->nrestart
-1])));
730 for (i = 0; i < curr->nmol; i++)
731 {
732 curr->lsq[curr->nrestart-1][i] = gmx_stats_init();
733 }
734
735 if (debug)
736 {
737 fprintf(debug, "Extended data structures because of new restart %d\n",
738 curr->nrestart);
739 }
740 }
741 /* create or extend the frame-based arrays */
742 if (curr->nframes >= maxframes-1)
743 {
744 if (maxframes == 0)
745 {
746 for (i = 0; (i < curr->ngrp); i++)
747 {
748 curr->ndata[i] = NULL((void*)0);
749 curr->data[i] = NULL((void*)0);
750 if (bTen)
751 {
752 curr->datam[i] = NULL((void*)0);
753 }
754 }
755 curr->time = NULL((void*)0);
756 }
757 maxframes += 10;
758 for (i = 0; (i < curr->ngrp); i++)
759 {
760 srenew(curr->ndata[i], maxframes)(curr->ndata[i]) = save_realloc("curr->ndata[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 760, (curr->ndata[i]), (maxframes), sizeof(*(curr->ndata
[i])));
761 srenew(curr->data[i], maxframes)(curr->data[i]) = save_realloc("curr->data[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 761, (curr->data[i]), (maxframes), sizeof(*(curr->data
[i])));
762 if (bTen)
763 {
764 srenew(curr->datam[i], maxframes)(curr->datam[i]) = save_realloc("curr->datam[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 764, (curr->datam[i]), (maxframes), sizeof(*(curr->datam
[i])));
765 }
766 for (j = maxframes-10; j < maxframes; j++)
767 {
768 curr->ndata[i][j] = 0;
769 curr->data[i][j] = 0;
770 if (bTen)
771 {
772 clear_mat(curr->datam[i][j]);
773 }
774 }
775 }
776 srenew(curr->time, maxframes)(curr->time) = save_realloc("curr->time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 776, (curr->time), (maxframes), sizeof(*(curr->time))
);
777 }
778
779 /* set the time */
780 curr->time[curr->nframes] = t - curr->t0;
781
782 /* for the first frame, the previous frame is a copy of the first frame */
783 if (bFirst)
784 {
785 memcpy(xa[prev(1-cur)], xa[cur], curr->ncoords*sizeof(xa[prev(1-cur)][0]));
786 bFirst = FALSE0;
787 }
788
789 /* make the molecules whole */
790 if (bMol)
791 {
792 gmx_rmpbc(gpbc, natoms, box, x[cur]);
793 }
794
795 /* calculate the molecules' centers of masses and put them into xa */
796 if (bMol)
797 {
798 calc_mol_com(gnx[0], index[0], &top->mols, &top->atoms, x[cur], xa[cur]);
799 }
800
801 /* first remove the periodic boundary condition crossings */
802 for (i = 0; i < curr->ngrp; i++)
803 {
804 prep_data(bMol, gnx[i], index[i], xa[cur], xa[prev(1-cur)], box);
805 }
806
807 /* calculate the center of mass */
808 if (gnx_com)
809 {
810 prep_data(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev(1-cur)], box);
811 calc_com(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev(1-cur)], box,
812 &top->atoms, com);
813 }
814
815 /* loop over all groups in index file */
816 for (i = 0; (i < curr->ngrp); i++)
817 {
818 /* calculate something useful, like mean square displacements */
819 calc_corr(curr, i, gnx[i], index[i], xa[cur], (gnx_com != NULL((void*)0)), com,
820 calc1, bTen);
821 }
822 cur = prev(1-cur);
823 t_prev = t;
824
825 curr->nframes++;
826 }
827 while (read_next_x(oenv, status, &t, x[cur], box));
828 fprintf(stderrstderr, "\nUsed %d restart points spaced %g %s over %g %s\n\n",
829 curr->nrestart,
830 output_env_conv_time(oenv, dt), output_env_get_time_unit(oenv),
831 output_env_conv_time(oenv, curr->time[curr->nframes-1]),
832 output_env_get_time_unit(oenv) );
833
834 if (bMol)
835 {
836 gmx_rmpbc_done(gpbc);
837 }
838
839 close_trj(status);
840
841 return natoms;
842}
843
844static void index_atom2mol(int *n, int *index, t_block *mols)
845{
846 int nat, i, nmol, mol, j;
847
848 nat = *n;
849 i = 0;
850 nmol = 0;
851 mol = 0;
852 while (i < nat)
853 {
854 while (index[i] > mols->index[mol])
855 {
856 mol++;
857 if (mol >= mols->nr)
858 {
859 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 859, "Atom index out of range: %d", index[i]+1);
860 }
861 }
862 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
863 {
864 if (i >= nat || index[i] != j)
865 {
866 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 866, "The index group does not consist of whole molecules");
867 }
868 i++;
869 }
870 index[nmol++] = mol;
871 }
872
873 fprintf(stderrstderr, "Split group of %d atoms into %d molecules\n", nat, nmol);
874
875 *n = nmol;
876}
877
878void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file,
879 const char *mol_file, const char *pdb_file, real t_pdb,
880 int nrgrp, t_topology *top, int ePBC,
881 gmx_bool bTen, gmx_bool bMW, gmx_bool bRmCOMM,
882 int type, real dim_factor, int axis,
883 real dt, real beginfit, real endfit, const output_env_t oenv)
884{
885 t_corr *msd;
886 int *gnx; /* the selected groups' sizes */
887 atom_id **index; /* selected groups' indices */
888 char **grpname;
889 int i, i0, i1, j, N, nat_trx;
890 real *DD, *SigmaD, a, a2, b, r, chi2;
891 rvec *x;
892 matrix box;
893 int *gnx_com = NULL((void*)0); /* the COM removal group size */
894 atom_id **index_com = NULL((void*)0); /* the COM removal group atom indices */
895 char **grpname_com = NULL((void*)0); /* the COM removal group name */
896
897 snew(gnx, nrgrp)(gnx) = save_calloc("gnx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 897, (nrgrp), sizeof(*(gnx)));
898 snew(index, nrgrp)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 898, (nrgrp), sizeof(*(index)));
899 snew(grpname, nrgrp)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 899, (nrgrp), sizeof(*(grpname)));
900
901 fprintf(stderrstderr, "\nSelect a group to calculate mean squared displacement for:\n");
902 get_index(&top->atoms, ndx_file, nrgrp, gnx, index, grpname);
903
904 if (bRmCOMM)
905 {
906 snew(gnx_com, 1)(gnx_com) = save_calloc("gnx_com", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 906, (1), sizeof(*(gnx_com)));
907 snew(index_com, 1)(index_com) = save_calloc("index_com", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 907, (1), sizeof(*(index_com)));
908 snew(grpname_com, 1)(grpname_com) = save_calloc("grpname_com", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 908, (1), sizeof(*(grpname_com)));
909
910 fprintf(stderrstderr, "\nNow select a group for center of mass removal:\n");
911 get_index(&top->atoms, ndx_file, 1, gnx_com, index_com, grpname_com);
912 }
913
914 if (mol_file)
915 {
916 index_atom2mol(&gnx[0], index[0], &top->mols);
917 }
918
919 msd = init_corr(nrgrp, type, axis, dim_factor,
920 mol_file == NULL((void*)0) ? 0 : gnx[0], bTen, bMW, dt, top,
921 beginfit, endfit);
922
923 nat_trx =
924 corr_loop(msd, trx_file, top, ePBC, mol_file ? gnx[0] : 0, gnx, index,
925 (mol_file != NULL((void*)0)) ? calc1_mol : (bMW ? calc1_mw : calc1_norm),
926 bTen, gnx_com, index_com, dt, t_pdb,
927 pdb_file ? &x : NULL((void*)0), box, oenv);
928
929 /* Correct for the number of points */
930 for (j = 0; (j < msd->ngrp); j++)
931 {
932 for (i = 0; (i < msd->nframes); i++)
933 {
934 msd->data[j][i] /= msd->ndata[j][i];
935 if (bTen)
936 {
937 msmul(msd->datam[j][i], 1.0/msd->ndata[j][i], msd->datam[j][i]);
938 }
939 }
940 }
941
942 if (mol_file)
943 {
944 if (pdb_file && x == NULL((void*)0))
945 {
946 fprintf(stderrstderr, "\nNo frame found need time tpdb = %g ps\n"
947 "Can not write %s\n\n", t_pdb, pdb_file);
948 }
949 i = top->atoms.nr;
950 top->atoms.nr = nat_trx;
951 printmol(msd, mol_file, pdb_file, index[0], top, x, ePBC, box, oenv);
952 top->atoms.nr = i;
953 }
954
955 DD = NULL((void*)0);
956 SigmaD = NULL((void*)0);
957
958 if (beginfit == -1)
959 {
960 i0 = (int)(0.1*(msd->nframes - 1) + 0.5);
961 beginfit = msd->time[i0];
962 }
963 else
964 {
965 for (i0 = 0; i0 < msd->nframes && msd->time[i0] < beginfit; i0++)
966 {
967 ;
968 }
969 }
970
971 if (endfit == -1)
972 {
973 i1 = (int)(0.9*(msd->nframes - 1) + 0.5) + 1;
974 endfit = msd->time[i1-1];
975 }
976 else
977 {
978 for (i1 = i0; i1 < msd->nframes && msd->time[i1] <= endfit; i1++)
979 {
980 ;
981 }
982 }
983 fprintf(stdoutstdout, "Fitting from %g to %g %s\n\n", beginfit, endfit,
984 output_env_get_time_unit(oenv));
985
986 N = i1-i0;
987 if (N <= 2)
988 {
989 fprintf(stdoutstdout, "Not enough points for fitting (%d).\n"
990 "Can not determine the diffusion constant.\n", N);
991 }
992 else
993 {
994 snew(DD, msd->ngrp)(DD) = save_calloc("DD", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 994, (msd->ngrp), sizeof(*(DD)));
995 snew(SigmaD, msd->ngrp)(SigmaD) = save_calloc("SigmaD", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 995, (msd->ngrp), sizeof(*(SigmaD)));
996 for (j = 0; j < msd->ngrp; j++)
997 {
998 if (N >= 4)
999 {
1000 lsq_y_ax_b(N/2, &(msd->time[i0]), &(msd->data[j][i0]), &a, &b, &r, &chi2);
1001 lsq_y_ax_b(N/2, &(msd->time[i0+N/2]), &(msd->data[j][i0+N/2]), &a2, &b, &r, &chi2);
1002 SigmaD[j] = fabs(a-a2);
1003 }
1004 else
1005 {
1006 SigmaD[j] = 0;
1007 }
1008 lsq_y_ax_b(N, &(msd->time[i0]), &(msd->data[j][i0]), &(DD[j]), &b, &r, &chi2);
1009 DD[j] *= FACTOR1000.0/msd->dim_factor;
1010 SigmaD[j] *= FACTOR1000.0/msd->dim_factor;
1011 if (DD[j] > 0.01 && DD[j] < 1e4)
1012 {
1013 fprintf(stdoutstdout, "D[%10s] %.4f (+/- %.4f) 1e-5 cm^2/s\n",
1014 grpname[j], DD[j], SigmaD[j]);
1015 }
1016 else
1017 {
1018 fprintf(stdoutstdout, "D[%10s] %.4g (+/- %.4g) 1e-5 cm^2/s\n",
1019 grpname[j], DD[j], SigmaD[j]);
1020 }
1021 }
1022 }
1023 /* Print mean square displacement */
1024 corr_print(msd, bTen, msd_file,
1025 "Mean Square Displacement",
1026 "MSD (nm\\S2\\N)",
1027 msd->time[msd->nframes-1], beginfit, endfit, DD, SigmaD, grpname, oenv);
1028}
1029
1030int gmx_msd(int argc, char *argv[])
1031{
1032 const char *desc[] = {
1033 "[THISMODULE] computes the mean square displacement (MSD) of atoms from",
1034 "a set of initial positions. This provides an easy way to compute",
1035 "the diffusion constant using the Einstein relation.",
1036 "The time between the reference points for the MSD calculation",
1037 "is set with [TT]-trestart[tt].",
1038 "The diffusion constant is calculated by least squares fitting a",
1039 "straight line (D*t + c) through the MSD(t) from [TT]-beginfit[tt] to",
1040 "[TT]-endfit[tt] (note that t is time from the reference positions,",
1041 "not simulation time). An error estimate given, which is the difference",
1042 "of the diffusion coefficients obtained from fits over the two halves",
1043 "of the fit interval.[PAR]",
1044 "There are three, mutually exclusive, options to determine different",
1045 "types of mean square displacement: [TT]-type[tt], [TT]-lateral[tt]",
1046 "and [TT]-ten[tt]. Option [TT]-ten[tt] writes the full MSD tensor for",
1047 "each group, the order in the output is: trace xx yy zz yx zx zy.[PAR]",
1048 "If [TT]-mol[tt] is set, [THISMODULE] plots the MSD for individual molecules",
1049 "(including making molecules whole across periodic boundaries): ",
1050 "for each individual molecule a diffusion constant is computed for ",
1051 "its center of mass. The chosen index group will be split into ",
1052 "molecules.[PAR]",
1053 "The default way to calculate a MSD is by using mass-weighted averages.",
1054 "This can be turned off with [TT]-nomw[tt].[PAR]",
1055 "With the option [TT]-rmcomm[tt], the center of mass motion of a ",
1056 "specific group can be removed. For trajectories produced with ",
1057 "GROMACS this is usually not necessary, ",
1058 "as [gmx-mdrun] usually already removes the center of mass motion.",
1059 "When you use this option be sure that the whole system is stored",
1060 "in the trajectory file.[PAR]",
1061 "The diffusion coefficient is determined by linear regression of the MSD,",
1062 "where, unlike for the normal output of D, the times are weighted",
1063 "according to the number of reference points, i.e. short times have",
1064 "a higher weight. Also when [TT]-beginfit[tt]=-1,fitting starts at 10%",
1065 "and when [TT]-endfit[tt]=-1, fitting goes to 90%.",
1066 "Using this option one also gets an accurate error estimate",
1067 "based on the statistics between individual molecules.",
1068 "Note that this diffusion coefficient and error estimate are only",
1069 "accurate when the MSD is completely linear between",
1070 "[TT]-beginfit[tt] and [TT]-endfit[tt].[PAR]",
1071 "Option [TT]-pdb[tt] writes a [TT].pdb[tt] file with the coordinates of the frame",
1072 "at time [TT]-tpdb[tt] with in the B-factor field the square root of",
1073 "the diffusion coefficient of the molecule.",
1074 "This option implies option [TT]-mol[tt]."
1075 };
1076 static const char *normtype[] = { NULL((void*)0), "no", "x", "y", "z", NULL((void*)0) };
1077 static const char *axtitle[] = { NULL((void*)0), "no", "x", "y", "z", NULL((void*)0) };
1078 static int ngroup = 1;
1079 static real dt = 10;
1080 static real t_pdb = 0;
1081 static real beginfit = -1;
1082 static real endfit = -1;
1083 static gmx_bool bTen = FALSE0;
1084 static gmx_bool bMW = TRUE1;
1085 static gmx_bool bRmCOMM = FALSE0;
1086 t_pargs pa[] = {
1087 { "-type", FALSE0, etENUM, {normtype},
1088 "Compute diffusion coefficient in one direction" },
1089 { "-lateral", FALSE0, etENUM, {axtitle},
1090 "Calculate the lateral diffusion in a plane perpendicular to" },
1091 { "-ten", FALSE0, etBOOL, {&bTen},
1092 "Calculate the full tensor" },
1093 { "-ngroup", FALSE0, etINT, {&ngroup},
1094 "Number of groups to calculate MSD for" },
1095 { "-mw", FALSE0, etBOOL, {&bMW},
1096 "Mass weighted MSD" },
1097 { "-rmcomm", FALSE0, etBOOL, {&bRmCOMM},
1098 "Remove center of mass motion" },
1099 { "-tpdb", FALSE0, etTIME, {&t_pdb},
1100 "The frame to use for option [TT]-pdb[tt] (%t)" },
1101 { "-trestart", FALSE0, etTIME, {&dt},
1102 "Time between restarting points in trajectory (%t)" },
1103 { "-beginfit", FALSE0, etTIME, {&beginfit},
1104 "Start time for fitting the MSD (%t), -1 is 10%" },
1105 { "-endfit", FALSE0, etTIME, {&endfit},
1106 "End time for fitting the MSD (%t), -1 is 90%" }
1107 };
1108
1109 t_filenm fnm[] = {
1110 { efTRX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
1111 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
1112 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
1113 { efXVG, NULL((void*)0), "msd", ffWRITE1<<2 },
1114 { efXVG, "-mol", "diff_mol", ffOPTWR(1<<2| 1<<3) },
1115 { efPDB, "-pdb", "diff_mol", ffOPTWR(1<<2| 1<<3) }
1116 };
1117#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
1118
1119 t_topology top;
1120 int ePBC;
1121 matrix box;
1122 char title[256];
1123 const char *trx_file, *tps_file, *ndx_file, *msd_file, *mol_file, *pdb_file;
1124 rvec *xdum;
1125 gmx_bool bTop;
1126 int axis, type;
1127 real dim_factor;
1128 output_env_t oenv;
1129
1130 if (!parse_common_args(&argc, argv,
1131 PCA_CAN_VIEW(1<<5) | PCA_CAN_BEGIN(1<<6) | PCA_CAN_END(1<<7) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13),
1132 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
1133 {
1134 return 0;
1135 }
1136 trx_file = ftp2fn_null(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1137 tps_file = ftp2fn_null(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1138 ndx_file = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1139 msd_file = ftp2fn_null(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1140 pdb_file = opt2fn_null("-pdb", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1141 if (pdb_file)
1142 {
1143 mol_file = opt2fn("-mol", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1144 }
1145 else
1146 {
1147 mol_file = opt2fn_null("-mol", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1148 }
1149
1150 if (ngroup < 1)
1151 {
1152 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 1152, "Must have at least 1 group (now %d)", ngroup);
1153 }
1154 if (mol_file && ngroup > 1)
1155 {
1156 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 1156, "With molecular msd can only have 1 group (now %d)",
1157 ngroup);
1158 }
1159
1160
1161 if (mol_file)
1162 {
1163 bMW = TRUE1;
1164 fprintf(stderrstderr, "Calculating diffusion coefficients for molecules.\n");
1165 }
1166
1167 if (normtype[0][0] != 'n')
1168 {
1169 type = normtype[0][0] - 'x' + X; /* See defines above */
1170 dim_factor = 2.0;
1171 }
1172 else
1173 {
1174 type = NORMAL;
1175 dim_factor = 6.0;
1176 }
1177 if ((type == NORMAL) && (axtitle[0][0] != 'n'))
1178 {
1179 type = LATERAL;
1180 dim_factor = 4.0;
1181 axis = (axtitle[0][0] - 'x'); /* See defines above */
1182 }
1183 else
1184 {
1185 axis = 0;
1186 }
1187
1188 if (bTen && type != NORMAL)
1189 {
1190 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 1190, "Can only calculate the full tensor for 3D msd");
1191 }
1192
1193 bTop = read_tps_conf(tps_file, title, &top, &ePBC, &xdum, NULL((void*)0), box, bMW || bRmCOMM);
1194 if (mol_file && !bTop)
1195 {
1196 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_msd.c"
, 1196,
1197 "Could not read a topology from %s. Try a tpr file instead.",
1198 tps_file);
1199 }
1200
1201 do_corr(trx_file, ndx_file, msd_file, mol_file, pdb_file, t_pdb, ngroup,
1202 &top, ePBC, bTen, bMW, bRmCOMM, type, dim_factor, axis, dt, beginfit, endfit,
1203 oenv);
1204
1205 view_all(oenv, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1206
1207 return 0;
1208}