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
54 #include "gmx_statistics.h"
62 #define FACTOR 1000.0 /* Convert nm^2/ps to 10e-5 cm^2/s */
63 /* NORMAL = total diffusion coefficient (default). X,Y,Z is diffusion
64 coefficient in X,Y,Z direction. LATERAL is diffusion coefficient in
65 plane perpendicular to axis
68 NOT_USED, NORMAL, X, Y, Z, LATERAL
72 real t0; /* start time and time increment between */
73 real delta_t; /* time between restart points */
74 real beginfit, /* the begin/end time for fits as reals between */
76 real dim_factor; /* the dimensionality factor for the diffusion
78 real **data; /* the displacement data. First index is the group
79 number, second is frame number */
80 real *time; /* frame time */
81 real *mass; /* masses for mass-weighted msd */
83 rvec **x0; /* original positions */
84 rvec *com; /* center of mass correction for each frame */
85 gmx_stats_t **lsq; /* fitting stats for individual molecule msds */
86 msd_type type; /* the type of msd to calculate (lateral, etc.)*/
87 int axis; /* the axis along which to calculate */
89 int nrestart; /* number of restart points */
90 int nmol; /* number of molecules (for bMol) */
91 int nframes; /* number of frames */
93 int ngrp; /* number of groups to use for msd calculation */
95 int **ndata; /* the number of msds (particles/mols) per data
99 typedef real t_calc_func (t_corr *, int, atom_id[], int, rvec[], rvec, gmx_bool, matrix,
100 const output_env_t oenv);
102 static real thistime(t_corr *curr)
104 return curr->time[curr->nframes];
107 static gmx_bool in_data(t_corr *curr, int nx00)
109 return curr->nframes-curr->n_offs[nx00];
112 t_corr *init_corr(int nrgrp, int type, int axis, real dim_factor,
113 int nmol, gmx_bool bTen, gmx_bool bMass, real dt, t_topology *top,
114 real beginfit, real endfit)
121 curr->type = (msd_type)type;
126 curr->beginfit = (1 - 2*GMX_REAL_EPS)*beginfit;
127 curr->endfit = (1 + 2*GMX_REAL_EPS)*endfit;
132 curr->dim_factor = dim_factor;
134 snew(curr->ndata, nrgrp);
135 snew(curr->data, nrgrp);
138 snew(curr->datam, nrgrp);
140 for (i = 0; (i < nrgrp); i++)
142 curr->ndata[i] = NULL;
143 curr->data[i] = NULL;
146 curr->datam[i] = NULL;
154 snew(curr->mass, curr->nmol);
155 for (i = 0; i < curr->nmol; i++)
165 snew(curr->mass, atoms->nr);
166 for (i = 0; (i < atoms->nr); i++)
168 curr->mass[i] = atoms->atom[i].m;
176 static void corr_print(t_corr *curr, gmx_bool bTen, const char *fn, const char *title,
178 real msdtime, real beginfit, real endfit,
179 real *DD, real *SigmaD, char *grpname[],
180 const output_env_t oenv)
185 out = xvgropen(fn, title, output_env_get_xvgr_tlabel(oenv), yaxis, oenv);
188 fprintf(out, "# MSD gathered over %g %s with %d restarts\n",
189 msdtime, output_env_get_time_unit(oenv), curr->nrestart);
190 fprintf(out, "# Diffusion constants fitted from time %g to %g %s\n",
191 beginfit, endfit, output_env_get_time_unit(oenv));
192 for (i = 0; i < curr->ngrp; i++)
194 fprintf(out, "# D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)\n",
195 grpname[i], DD[i], SigmaD[i]);
198 for (i = 0; i < curr->nframes; i++)
200 fprintf(out, "%10g", output_env_conv_time(oenv, curr->time[i]));
201 for (j = 0; j < curr->ngrp; j++)
203 fprintf(out, " %10g", curr->data[j][i]);
206 fprintf(out, " %10g %10g %10g %10g %10g %10g",
207 curr->datam[j][i][XX][XX],
208 curr->datam[j][i][YY][YY],
209 curr->datam[j][i][ZZ][ZZ],
210 curr->datam[j][i][YY][XX],
211 curr->datam[j][i][ZZ][XX],
212 curr->datam[j][i][ZZ][YY]);
220 /* called from corr_loop, to do the main calculations */
221 static void calc_corr(t_corr *curr, int nr, int nx, atom_id index[], rvec xc[],
222 gmx_bool bRmCOMM, rvec com, t_calc_func *calc1, gmx_bool bTen,
223 const output_env_t oenv)
230 /* Check for new starting point */
231 if (curr->nlast < curr->nrestart)
233 if ((thistime(curr) >= (curr->nlast*curr->delta_t)) && (nr == 0))
235 memcpy(curr->x0[curr->nlast], xc, curr->ncoords*sizeof(xc[0]));
236 curr->n_offs[curr->nlast] = curr->nframes;
237 copy_rvec(com, curr->com[curr->nlast]);
242 /* nx0 appears to be the number of new starting points,
243 * so for all starting points, call calc1.
245 for (nx0 = 0; (nx0 < curr->nlast); nx0++)
249 rvec_sub(com, curr->com[nx0], dcom);
255 g = calc1(curr, nx, index, nx0, xc, dcom, bTen, mat, oenv);
257 printf("g[%d]=%g\n", nx0, g);
259 curr->data[nr][in_data(curr, nx0)] += g;
262 m_add(curr->datam[nr][in_data(curr, nx0)], mat,
263 curr->datam[nr][in_data(curr, nx0)]);
265 curr->ndata[nr][in_data(curr, nx0)]++;
269 /* the non-mass-weighted mean-squared displacement calcuation */
270 static real calc1_norm(t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
271 rvec dcom, gmx_bool bTen, matrix mat, const output_env_t oenv)
280 for (i = 0; (i < nx); i++)
287 for (m = 0; (m < DIM); m++)
289 rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
293 for (m2 = 0; m2 <= m; m2++)
295 mat[m][m2] += rv[m]*rv[m2];
303 r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
308 for (m = 0; (m < DIM); m++)
312 r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
318 gmx_fatal(FARGS, "Error: did not expect option value %d", curr->type);
323 msmul(mat, 1.0/nx, mat);
328 /* calculate the com of molecules in x and put it into xa */
329 static void calc_mol_com(int nmol, int *molindex, t_block *mols, t_atoms *atoms,
336 for (m = 0; m < nmol; m++)
341 for (i = mols->index[mol]; i < mols->index[mol+1]; i++)
343 mass = atoms->atom[i].m;
344 for (d = 0; d < DIM; d++)
346 xm[d] += mass*x[i][d];
350 svmul(1/mtot, xm, xa[m]);
354 static real calc_one_mw(t_corr *curr, int ix, int nx0, rvec xc[], real *tm,
355 rvec dcom, gmx_bool bTen, matrix mat)
371 for (m = 0; (m < DIM); m++)
373 rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
374 r2 += mm*rv[m]*rv[m];
377 for (m2 = 0; m2 <= m; m2++)
379 mat[m][m2] += mm*rv[m]*rv[m2];
387 r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
392 for (m = 0; (m < DIM); m++)
396 r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
402 gmx_fatal(FARGS, "Options got screwed. Did not expect value %d\n", curr->type);
407 /* the normal, mass-weighted mean-squared displacement calcuation */
408 static real calc1_mw(t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
409 rvec dcom, gmx_bool bTen, matrix mat, const output_env_t oenv)
416 for (i = 0; (i < nx); i++)
418 g += calc_one_mw(curr, index[i], nx0, xc, &tm, dcom, bTen, mat);
424 msmul(mat, 1/tm, mat);
430 /* prepare the coordinates by removing periodic boundary crossings.
431 gnx = the number of atoms/molecules
433 xcur = the current coordinates
434 xprev = the previous coordinates
435 box = the box matrix */
436 static void prep_data(gmx_bool bMol, int gnx, atom_id index[],
437 rvec xcur[], rvec xprev[], matrix box)
442 /* Remove periodicity */
443 for (m = 0; (m < DIM); m++)
445 hbox[m] = 0.5*box[m][m];
448 for (i = 0; (i < gnx); i++)
459 for (m = DIM-1; m >= 0; m--)
465 while (xcur[ind][m]-xprev[ind][m] <= -hbox[m])
467 rvec_inc(xcur[ind], box[m]);
469 while (xcur[ind][m]-xprev[ind][m] > hbox[m])
471 rvec_dec(xcur[ind], box[m]);
477 /* calculate the center of mass for a group
478 gnx = the number of atoms/molecules
480 xcur = the current coordinates
481 xprev = the previous coordinates
483 atoms = atom data (for mass)
484 com(output) = center of mass */
485 static void calc_com(gmx_bool bMol, int gnx, atom_id index[],
486 rvec xcur[], rvec xprev[], matrix box, t_atoms *atoms,
498 prep_data(bMol, gnx, index, xcur, xprev, box);
499 for (i = 0; (i < gnx); i++)
511 mass = atoms->atom[ind].m;
512 for (m = 0; m < DIM; m++)
514 sx[m] += mass*xcur[ind][m];
518 for (m = 0; m < DIM; m++)
520 com[m] = sx[m]/tmass;
525 static real calc1_mol(t_corr *curr, int nx, atom_id index[], int nx0, rvec xc[],
526 rvec dcom, gmx_bool bTen, matrix mat, const output_env_t oenv)
529 real g, tm, gtot, tt;
531 tt = curr->time[in_data(curr, nx0)];
535 for (i = 0; (i < nx); i++)
537 g = calc_one_mw(curr, i, nx0, xc, &tm, dcom, bTen, mat);
538 /* We don't need to normalize as the mass was set to 1 */
540 if (tt >= curr->beginfit && (curr->endfit < 0 || tt <= curr->endfit))
542 gmx_stats_add_point(curr->lsq[nx0][i], tt, g, 0, 0);
545 msmul(mat, 1.0/nx, mat);
550 void printmol(t_corr *curr, const char *fn,
551 const char *fn_pdb, int *molindex, t_topology *top,
552 rvec *x, int ePBC, matrix box, const output_env_t oenv)
558 real a, b, D, Dav, D2av, VarD, sqrtD, sqrtD_max, scale;
559 t_pdbinfo *pdbinfo = NULL;
562 out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D", oenv);
566 if (top->atoms.pdbinfo == NULL)
568 snew(top->atoms.pdbinfo, top->atoms.nr);
570 pdbinfo = top->atoms.pdbinfo;
571 mol2a = top->mols.index;
576 for (i = 0; (i < curr->nmol); i++)
578 lsq1 = gmx_stats_init();
579 for (j = 0; (j < curr->nrestart); j++)
583 while (gmx_stats_get_point(curr->lsq[j][i], &xx, &yy, &dx, &dy, 0) == estatsOK)
585 gmx_stats_add_point(lsq1, xx, yy, dx, dy);
588 gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, NULL, NULL, NULL, NULL);
589 gmx_stats_done(lsq1);
591 D = a*FACTOR/curr->dim_factor;
598 fprintf(out, "%10d %10g\n", i, D);
602 if (sqrtD > sqrtD_max)
606 for (j = mol2a[molindex[i]]; j < mol2a[molindex[i]+1]; j++)
608 pdbinfo[j].bfac = sqrtD;
613 do_view(oenv, fn, "-graphtype bar");
615 /* Compute variance, stddev and error */
618 VarD = D2av - sqr(Dav);
619 printf("<D> = %.4f Std. Dev. = %.4f Error = %.4f\n",
620 Dav, sqrt(VarD), sqrt(VarD/curr->nmol));
625 while (scale*sqrtD_max > 10)
629 while (scale*sqrtD_max < 0.1)
633 for (i = 0; i < top->atoms.nr; i++)
635 pdbinfo[i].bfac *= scale;
637 write_sto_conf(fn_pdb, "molecular MSD", &top->atoms, x, NULL, ePBC, box);
641 /* this is the main loop for the correlation type functions
642 * fx and nx are file pointers to things like read_first_x and
645 int corr_loop(t_corr *curr, const char *fn, t_topology *top, int ePBC,
646 gmx_bool bMol, int gnx[], atom_id *index[],
647 t_calc_func *calc1, gmx_bool bTen, int *gnx_com, atom_id *index_com[],
648 real dt, real t_pdb, rvec **x_pdb, matrix box_pdb,
649 const output_env_t oenv)
651 rvec *x[2]; /* the coordinates to read */
652 rvec *xa[2]; /* the coordinates to calculate displacements for */
655 int natoms, i, j, cur = 0, maxframes = 0;
660 gmx_rmpbc_t gpbc = NULL;
662 natoms = read_first_x(oenv, &status, fn, &curr->t0, &(x[cur]), box);
664 fprintf(stderr, "Read %d atoms for first frame\n", natoms);
666 if ((gnx_com != NULL) && natoms < top->atoms.nr)
668 fprintf(stderr, "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);
671 snew(x[prev], natoms);
675 curr->ncoords = curr->nmol;
676 snew(xa[0], curr->ncoords);
677 snew(xa[1], curr->ncoords);
681 curr->ncoords = natoms;
695 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms, box);
698 /* the loop over all frames */
701 if (x_pdb && ((bFirst && t_pdb < t) ||
703 t_pdb > t - 0.5*(t - t_prev) &&
704 t_pdb < t + 0.5*(t - t_prev))))
708 snew(*x_pdb, natoms);
710 for (i = 0; i < natoms; i++)
712 copy_rvec(x[cur][i], (*x_pdb)[i]);
714 copy_mat(box, box_pdb);
718 /* check whether we've reached a restart point */
719 if (bRmod(t, curr->t0, dt))
723 srenew(curr->x0, curr->nrestart);
724 snew(curr->x0[curr->nrestart-1], curr->ncoords);
725 srenew(curr->com, curr->nrestart);
726 srenew(curr->n_offs, curr->nrestart);
727 srenew(curr->lsq, curr->nrestart);
728 snew(curr->lsq[curr->nrestart-1], curr->nmol);
729 for (i = 0; i < curr->nmol; i++)
731 curr->lsq[curr->nrestart-1][i] = gmx_stats_init();
736 fprintf(debug, "Extended data structures because of new restart %d\n",
740 /* create or extend the frame-based arrays */
741 if (curr->nframes >= maxframes-1)
745 for (i = 0; (i < curr->ngrp); i++)
747 curr->ndata[i] = NULL;
748 curr->data[i] = NULL;
751 curr->datam[i] = NULL;
757 for (i = 0; (i < curr->ngrp); i++)
759 srenew(curr->ndata[i], maxframes);
760 srenew(curr->data[i], maxframes);
763 srenew(curr->datam[i], maxframes);
765 for (j = maxframes-10; j < maxframes; j++)
767 curr->ndata[i][j] = 0;
768 curr->data[i][j] = 0;
771 clear_mat(curr->datam[i][j]);
775 srenew(curr->time, maxframes);
779 curr->time[curr->nframes] = t - curr->t0;
781 /* for the first frame, the previous frame is a copy of the first frame */
784 memcpy(xa[prev], xa[cur], curr->ncoords*sizeof(xa[prev][0]));
788 /* make the molecules whole */
791 gmx_rmpbc(gpbc, natoms, box, x[cur]);
794 /* calculate the molecules' centers of masses and put them into xa */
797 calc_mol_com(gnx[0], index[0], &top->mols, &top->atoms, x[cur], xa[cur]);
800 /* first remove the periodic boundary condition crossings */
801 for (i = 0; i < curr->ngrp; i++)
803 prep_data(bMol, gnx[i], index[i], xa[cur], xa[prev], box);
806 /* calculate the center of mass */
809 prep_data(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev], box);
810 calc_com(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev], box,
814 /* loop over all groups in index file */
815 for (i = 0; (i < curr->ngrp); i++)
817 /* calculate something useful, like mean square displacements */
818 calc_corr(curr, i, gnx[i], index[i], xa[cur], (gnx_com != NULL), com,
826 while (read_next_x(oenv, status, &t, natoms, x[cur], box));
827 fprintf(stderr, "\nUsed %d restart points spaced %g %s over %g %s\n\n",
829 output_env_conv_time(oenv, dt), output_env_get_time_unit(oenv),
830 output_env_conv_time(oenv, curr->time[curr->nframes-1]),
831 output_env_get_time_unit(oenv) );
835 gmx_rmpbc_done(gpbc);
843 static void index_atom2mol(int *n, int *index, t_block *mols)
845 int nat, i, nmol, mol, j;
853 while (index[i] > mols->index[mol])
858 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
861 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
863 if (i >= nat || index[i] != j)
865 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
872 fprintf(stderr, "Split group of %d atoms into %d molecules\n", nat, nmol);
877 void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file,
878 const char *mol_file, const char *pdb_file, real t_pdb,
879 int nrgrp, t_topology *top, int ePBC,
880 gmx_bool bTen, gmx_bool bMW, gmx_bool bRmCOMM,
881 int type, real dim_factor, int axis,
882 real dt, real beginfit, real endfit, const output_env_t oenv)
885 int *gnx; /* the selected groups' sizes */
886 atom_id **index; /* selected groups' indices */
888 int i, i0, i1, j, N, nat_trx;
889 real *DD, *SigmaD, a, a2, b, r, chi2;
892 int *gnx_com = NULL; /* the COM removal group size */
893 atom_id **index_com = NULL; /* the COM removal group atom indices */
894 char **grpname_com = NULL; /* the COM removal group name */
898 snew(grpname, nrgrp);
900 fprintf(stderr, "\nSelect a group to calculate mean squared displacement for:\n");
901 get_index(&top->atoms, ndx_file, nrgrp, gnx, index, grpname);
907 snew(grpname_com, 1);
909 fprintf(stderr, "\nNow select a group for center of mass removal:\n");
910 get_index(&top->atoms, ndx_file, 1, gnx_com, index_com, grpname_com);
915 index_atom2mol(&gnx[0], index[0], &top->mols);
918 msd = init_corr(nrgrp, type, axis, dim_factor,
919 mol_file == NULL ? 0 : gnx[0], bTen, bMW, dt, top,
923 corr_loop(msd, trx_file, top, ePBC, mol_file ? gnx[0] : 0, gnx, index,
924 (mol_file != NULL) ? calc1_mol : (bMW ? calc1_mw : calc1_norm),
925 bTen, gnx_com, index_com, dt, t_pdb,
926 pdb_file ? &x : NULL, box, oenv);
928 /* Correct for the number of points */
929 for (j = 0; (j < msd->ngrp); j++)
931 for (i = 0; (i < msd->nframes); i++)
933 msd->data[j][i] /= msd->ndata[j][i];
936 msmul(msd->datam[j][i], 1.0/msd->ndata[j][i], msd->datam[j][i]);
943 if (pdb_file && x == NULL)
945 fprintf(stderr, "\nNo frame found need time tpdb = %g ps\n"
946 "Can not write %s\n\n", t_pdb, pdb_file);
949 top->atoms.nr = nat_trx;
950 printmol(msd, mol_file, pdb_file, index[0], top, x, ePBC, box, oenv);
959 i0 = (int)(0.1*(msd->nframes - 1) + 0.5);
960 beginfit = msd->time[i0];
964 for (i0 = 0; i0 < msd->nframes && msd->time[i0] < beginfit; i0++)
972 i1 = (int)(0.9*(msd->nframes - 1) + 0.5) + 1;
973 endfit = msd->time[i1-1];
977 for (i1 = i0; i1 < msd->nframes && msd->time[i1] <= endfit; i1++)
982 fprintf(stdout, "Fitting from %g to %g %s\n\n", beginfit, endfit,
983 output_env_get_time_unit(oenv));
988 fprintf(stdout, "Not enough points for fitting (%d).\n"
989 "Can not determine the diffusion constant.\n", N);
994 snew(SigmaD, msd->ngrp);
995 for (j = 0; j < msd->ngrp; j++)
999 lsq_y_ax_b(N/2, &(msd->time[i0]), &(msd->data[j][i0]), &a, &b, &r, &chi2);
1000 lsq_y_ax_b(N/2, &(msd->time[i0+N/2]), &(msd->data[j][i0+N/2]), &a2, &b, &r, &chi2);
1001 SigmaD[j] = fabs(a-a2);
1007 lsq_y_ax_b(N, &(msd->time[i0]), &(msd->data[j][i0]), &(DD[j]), &b, &r, &chi2);
1008 DD[j] *= FACTOR/msd->dim_factor;
1009 SigmaD[j] *= FACTOR/msd->dim_factor;
1010 if (DD[j] > 0.01 && DD[j] < 1e4)
1012 fprintf(stdout, "D[%10s] %.4f (+/- %.4f) 1e-5 cm^2/s\n",
1013 grpname[j], DD[j], SigmaD[j]);
1017 fprintf(stdout, "D[%10s] %.4g (+/- %.4g) 1e-5 cm^2/s\n",
1018 grpname[j], DD[j], SigmaD[j]);
1022 /* Print mean square displacement */
1023 corr_print(msd, bTen, msd_file,
1024 "Mean Square Displacement",
1026 msd->time[msd->nframes-1], beginfit, endfit, DD, SigmaD, grpname, oenv);
1029 int gmx_msd(int argc, char *argv[])
1031 const char *desc[] = {
1032 "[TT]g_msd[tt] computes the mean square displacement (MSD) of atoms from",
1033 "a set of initial positions. This provides an easy way to compute",
1034 "the diffusion constant using the Einstein relation.",
1035 "The time between the reference points for the MSD calculation",
1036 "is set with [TT]-trestart[tt].",
1037 "The diffusion constant is calculated by least squares fitting a",
1038 "straight line (D*t + c) through the MSD(t) from [TT]-beginfit[tt] to",
1039 "[TT]-endfit[tt] (note that t is time from the reference positions,",
1040 "not simulation time). An error estimate given, which is the difference",
1041 "of the diffusion coefficients obtained from fits over the two halves",
1042 "of the fit interval.[PAR]",
1043 "There are three, mutually exclusive, options to determine different",
1044 "types of mean square displacement: [TT]-type[tt], [TT]-lateral[tt]",
1045 "and [TT]-ten[tt]. Option [TT]-ten[tt] writes the full MSD tensor for",
1046 "each group, the order in the output is: trace xx yy zz yx zx zy.[PAR]",
1047 "If [TT]-mol[tt] is set, [TT]g_msd[tt] plots the MSD for individual molecules",
1048 "(including making molecules whole across periodic boundaries): ",
1049 "for each individual molecule a diffusion constant is computed for ",
1050 "its center of mass. The chosen index group will be split into ",
1052 "The default way to calculate a MSD is by using mass-weighted averages.",
1053 "This can be turned off with [TT]-nomw[tt].[PAR]",
1054 "With the option [TT]-rmcomm[tt], the center of mass motion of a ",
1055 "specific group can be removed. For trajectories produced with ",
1056 "GROMACS this is usually not necessary, ",
1057 "as [TT]mdrun[tt] usually already removes the center of mass motion.",
1058 "When you use this option be sure that the whole system is stored",
1059 "in the trajectory file.[PAR]",
1060 "The diffusion coefficient is determined by linear regression of the MSD,",
1061 "where, unlike for the normal output of D, the times are weighted",
1062 "according to the number of reference points, i.e. short times have",
1063 "a higher weight. Also when [TT]-beginfit[tt]=-1,fitting starts at 10%",
1064 "and when [TT]-endfit[tt]=-1, fitting goes to 90%.",
1065 "Using this option one also gets an accurate error estimate",
1066 "based on the statistics between individual molecules.",
1067 "Note that this diffusion coefficient and error estimate are only",
1068 "accurate when the MSD is completely linear between",
1069 "[TT]-beginfit[tt] and [TT]-endfit[tt].[PAR]",
1070 "Option [TT]-pdb[tt] writes a [TT].pdb[tt] file with the coordinates of the frame",
1071 "at time [TT]-tpdb[tt] with in the B-factor field the square root of",
1072 "the diffusion coefficient of the molecule.",
1073 "This option implies option [TT]-mol[tt]."
1075 static const char *normtype[] = { NULL, "no", "x", "y", "z", NULL };
1076 static const char *axtitle[] = { NULL, "no", "x", "y", "z", NULL };
1077 static int ngroup = 1;
1078 static real dt = 10;
1079 static real t_pdb = 0;
1080 static real beginfit = -1;
1081 static real endfit = -1;
1082 static gmx_bool bTen = FALSE;
1083 static gmx_bool bMW = TRUE;
1084 static gmx_bool bRmCOMM = FALSE;
1086 { "-type", FALSE, etENUM, {normtype},
1087 "Compute diffusion coefficient in one direction" },
1088 { "-lateral", FALSE, etENUM, {axtitle},
1089 "Calculate the lateral diffusion in a plane perpendicular to" },
1090 { "-ten", FALSE, etBOOL, {&bTen},
1091 "Calculate the full tensor" },
1092 { "-ngroup", FALSE, etINT, {&ngroup},
1093 "Number of groups to calculate MSD for" },
1094 { "-mw", FALSE, etBOOL, {&bMW},
1095 "Mass weighted MSD" },
1096 { "-rmcomm", FALSE, etBOOL, {&bRmCOMM},
1097 "Remove center of mass motion" },
1098 { "-tpdb", FALSE, etTIME, {&t_pdb},
1099 "The frame to use for option [TT]-pdb[tt] (%t)" },
1100 { "-trestart", FALSE, etTIME, {&dt},
1101 "Time between restarting points in trajectory (%t)" },
1102 { "-beginfit", FALSE, etTIME, {&beginfit},
1103 "Start time for fitting the MSD (%t), -1 is 10%" },
1104 { "-endfit", FALSE, etTIME, {&endfit},
1105 "End time for fitting the MSD (%t), -1 is 90%" }
1109 { efTRX, NULL, NULL, ffREAD },
1110 { efTPS, NULL, NULL, ffREAD },
1111 { efNDX, NULL, NULL, ffOPTRD },
1112 { efXVG, NULL, "msd", ffWRITE },
1113 { efXVG, "-mol", "diff_mol", ffOPTWR },
1114 { efPDB, "-pdb", "diff_mol", ffOPTWR }
1116 #define NFILE asize(fnm)
1122 const char *trx_file, *tps_file, *ndx_file, *msd_file, *mol_file, *pdb_file;
1129 parse_common_args(&argc, argv,
1130 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_TIME_UNIT | PCA_BE_NICE,
1131 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
1132 trx_file = ftp2fn_null(efTRX, NFILE, fnm);
1133 tps_file = ftp2fn_null(efTPS, NFILE, fnm);
1134 ndx_file = ftp2fn_null(efNDX, NFILE, fnm);
1135 msd_file = ftp2fn_null(efXVG, NFILE, fnm);
1136 pdb_file = opt2fn_null("-pdb", NFILE, fnm);
1139 mol_file = opt2fn("-mol", NFILE, fnm);
1143 mol_file = opt2fn_null("-mol", NFILE, fnm);
1148 gmx_fatal(FARGS, "Must have at least 1 group (now %d)", ngroup);
1150 if (mol_file && ngroup > 1)
1152 gmx_fatal(FARGS, "With molecular msd can only have 1 group (now %d)",
1160 fprintf(stderr, "Calculating diffusion coefficients for molecules.\n");
1163 if (normtype[0][0] != 'n')
1165 type = normtype[0][0] - 'x' + X; /* See defines above */
1173 if ((type == NORMAL) && (axtitle[0][0] != 'n'))
1177 axis = (axtitle[0][0] - 'x'); /* See defines above */
1184 if (bTen && type != NORMAL)
1186 gmx_fatal(FARGS, "Can only calculate the full tensor for 3D msd");
1189 bTop = read_tps_conf(tps_file, title, &top, &ePBC, &xdum, NULL, box, bMW || bRmCOMM);
1190 if (mol_file && !bTop)
1193 "Could not read a topology from %s. Try a tpr file instead.",
1197 do_corr(trx_file, ndx_file, msd_file, mol_file, pdb_file, t_pdb, ngroup,
1198 &top, ePBC, bTen, bMW, bRmCOMM, type, dim_factor, axis, dt, beginfit, endfit,
1201 view_all(oenv, NFILE, fnm);