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,2015,2016,2017,2018, 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.
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/commandline/viewit.h"
44 #include "gromacs/fileio/confio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/gmxana/gmx_ana.h"
48 #include "gromacs/gmxana/gstat.h"
49 #include "gromacs/math/functions.h"
50 #include "gromacs/math/utilities.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/pbcutil/rmpbc.h"
53 #include "gromacs/statistics/statistics.h"
54 #include "gromacs/topology/index.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/arraysize.h"
57 #include "gromacs/utility/fatalerror.h"
58 #include "gromacs/utility/futil.h"
59 #include "gromacs/utility/gmxassert.h"
60 #include "gromacs/utility/smalloc.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 *curr, int nx, const int index[], int nx0, rvec xc[],
100 const rvec dcom, gmx_bool bTen, matrix mat);
102 static real thistime(t_corr *curr)
104 return curr->time[curr->nframes];
107 static int in_data(t_corr *curr, int nx00)
109 return curr->nframes-curr->n_offs[nx00];
112 static t_corr *init_corr(int nrgrp, int type, int axis, real dim_factor,
113 int nmol, gmx_bool bTen, gmx_bool bMass, real dt, const t_topology *top,
114 real beginfit, real endfit)
120 curr->type = static_cast<msd_type>(type);
125 curr->beginfit = (1 - 2*GMX_REAL_EPS)*beginfit;
126 curr->endfit = (1 + 2*GMX_REAL_EPS)*endfit;
128 curr->n_offs = nullptr;
131 curr->dim_factor = dim_factor;
133 snew(curr->ndata, nrgrp);
134 snew(curr->data, nrgrp);
137 snew(curr->datam, nrgrp);
139 for (i = 0; (i < nrgrp); i++)
141 curr->ndata[i] = nullptr;
142 curr->data[i] = nullptr;
145 curr->datam[i] = nullptr;
148 curr->time = nullptr;
153 snew(curr->mass, curr->nmol);
154 for (i = 0; i < curr->nmol; i++)
163 const t_atoms *atoms = &top->atoms;
164 snew(curr->mass, atoms->nr);
165 for (i = 0; (i < atoms->nr); i++)
167 curr->mass[i] = atoms->atom[i].m;
175 static void corr_print(t_corr *curr, gmx_bool bTen, const char *fn, const char *title,
177 real msdtime, real beginfit, real endfit,
178 real *DD, real *SigmaD, char *grpname[],
179 const gmx_output_env_t *oenv)
184 out = xvgropen(fn, title, output_env_get_xvgr_tlabel(oenv), yaxis, oenv);
187 fprintf(out, "# MSD gathered over %g %s with %d restarts\n",
188 msdtime, output_env_get_time_unit(oenv).c_str(), curr->nrestart);
189 fprintf(out, "# Diffusion constants fitted from time %g to %g %s\n",
190 beginfit, endfit, output_env_get_time_unit(oenv).c_str());
191 for (i = 0; i < curr->ngrp; i++)
193 fprintf(out, "# D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)\n",
194 grpname[i], DD[i], SigmaD[i]);
197 for (i = 0; i < curr->nframes; i++)
199 fprintf(out, "%10g", output_env_conv_time(oenv, curr->time[i]));
200 for (j = 0; j < curr->ngrp; j++)
202 fprintf(out, " %10g", curr->data[j][i]);
205 fprintf(out, " %10g %10g %10g %10g %10g %10g",
206 curr->datam[j][i][XX][XX],
207 curr->datam[j][i][YY][YY],
208 curr->datam[j][i][ZZ][ZZ],
209 curr->datam[j][i][YY][XX],
210 curr->datam[j][i][ZZ][XX],
211 curr->datam[j][i][ZZ][YY]);
219 /* called from corr_loop, to do the main calculations */
220 static void calc_corr(t_corr *curr, int nr, int nx, int index[], rvec xc[],
221 gmx_bool bRmCOMM, rvec com, t_calc_func *calc1, gmx_bool bTen)
228 /* Check for new starting point */
229 if (curr->nlast < curr->nrestart)
231 if ((thistime(curr) >= (curr->nlast*curr->delta_t)) && (nr == 0))
233 std::memcpy(curr->x0[curr->nlast], xc, curr->ncoords*sizeof(xc[0]));
234 curr->n_offs[curr->nlast] = curr->nframes;
235 copy_rvec(com, curr->com[curr->nlast]);
240 /* nx0 appears to be the number of new starting points,
241 * so for all starting points, call calc1.
243 for (nx0 = 0; (nx0 < curr->nlast); nx0++)
247 rvec_sub(com, curr->com[nx0], dcom);
253 g = calc1(curr, nx, index, nx0, xc, dcom, bTen, mat);
255 printf("g[%d]=%g\n", nx0, g);
257 curr->data[nr][in_data(curr, nx0)] += g;
260 m_add(curr->datam[nr][in_data(curr, nx0)], mat,
261 curr->datam[nr][in_data(curr, nx0)]);
263 curr->ndata[nr][in_data(curr, nx0)]++;
267 /* the non-mass-weighted mean-squared displacement calcuation */
268 static real calc1_norm(t_corr *curr, int nx, const int index[], int nx0, rvec xc[],
269 const rvec dcom, gmx_bool bTen, matrix mat)
278 for (i = 0; (i < nx); i++)
285 for (m = 0; (m < DIM); m++)
287 rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
291 for (m2 = 0; m2 <= m; m2++)
293 mat[m][m2] += rv[m]*rv[m2];
301 r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
306 for (m = 0; (m < DIM); m++)
310 r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
316 gmx_fatal(FARGS, "Error: did not expect option value %d", curr->type);
321 msmul(mat, 1.0/nx, mat);
326 /* calculate the com of molecules in x and put it into xa */
327 static void calc_mol_com(int nmol, const int *molindex, const t_block *mols, const t_atoms *atoms,
334 for (m = 0; m < nmol; m++)
339 for (i = mols->index[mol]; i < mols->index[mol+1]; i++)
341 mass = atoms->atom[i].m;
342 for (d = 0; d < DIM; d++)
344 xm[d] += mass*x[i][d];
348 svmul(1/mtot, xm, xa[m]);
352 static real calc_one_mw(t_corr *curr, int ix, int nx0, rvec xc[], real *tm,
353 const rvec dcom, gmx_bool bTen, matrix mat)
369 for (m = 0; (m < DIM); m++)
371 rv[m] = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
372 r2 += mm*rv[m]*rv[m];
375 for (m2 = 0; m2 <= m; m2++)
377 mat[m][m2] += mm*rv[m]*rv[m2];
385 r = xc[ix][curr->type-X] - curr->x0[nx0][ix][curr->type-X] -
390 for (m = 0; (m < DIM); m++)
394 r = xc[ix][m] - curr->x0[nx0][ix][m] - dcom[m];
400 gmx_fatal(FARGS, "Options got screwed. Did not expect value %d\n", curr->type);
405 /* the normal, mass-weighted mean-squared displacement calcuation */
406 static real calc1_mw(t_corr *curr, int nx, const int index[], int nx0, rvec xc[],
407 const rvec dcom, gmx_bool bTen, matrix mat)
414 for (i = 0; (i < nx); i++)
416 g += calc_one_mw(curr, index[i], nx0, xc, &tm, dcom, bTen, mat);
422 msmul(mat, 1/tm, mat);
428 /* prepare the coordinates by removing periodic boundary crossings.
429 gnx = the number of atoms/molecules
431 xcur = the current coordinates
432 xprev = the previous coordinates
433 box = the box matrix */
434 static void prep_data(gmx_bool bMol, int gnx, const int index[],
435 rvec xcur[], rvec xprev[], matrix box)
440 /* Remove periodicity */
441 for (m = 0; (m < DIM); m++)
443 hbox[m] = 0.5*box[m][m];
446 for (i = 0; (i < gnx); i++)
457 for (m = DIM-1; m >= 0; m--)
463 while (xcur[ind][m]-xprev[ind][m] <= -hbox[m])
465 rvec_inc(xcur[ind], box[m]);
467 while (xcur[ind][m]-xprev[ind][m] > hbox[m])
469 rvec_dec(xcur[ind], box[m]);
475 /* calculate the center of mass for a group
476 gnx = the number of atoms/molecules
478 xcur = the current coordinates
479 xprev = the previous coordinates
481 atoms = atom data (for mass)
482 com(output) = center of mass */
483 static void calc_com(gmx_bool bMol, int gnx, int index[],
484 rvec xcur[], rvec xprev[], matrix box, const t_atoms *atoms,
495 prep_data(bMol, gnx, index, xcur, xprev, box);
496 for (i = 0; (i < gnx); i++)
508 mass = atoms->atom[ind].m;
509 for (m = 0; m < DIM; m++)
511 sx[m] += mass*xcur[ind][m];
515 for (m = 0; m < DIM; m++)
517 com[m] = sx[m]/tmass;
522 static real calc1_mol(t_corr *curr, int nx, const int gmx_unused index[], int nx0, rvec xc[],
523 const rvec dcom, gmx_bool bTen, matrix mat)
526 real g, tm, gtot, tt;
528 tt = curr->time[in_data(curr, nx0)];
532 for (i = 0; (i < nx); i++)
534 g = calc_one_mw(curr, i, nx0, xc, &tm, dcom, bTen, mat);
535 /* We don't need to normalize as the mass was set to 1 */
537 if (tt >= curr->beginfit && (curr->endfit < 0 || tt <= curr->endfit))
539 gmx_stats_add_point(curr->lsq[nx0][i], tt, g, 0, 0);
542 msmul(mat, 1.0/nx, mat);
547 static void printmol(t_corr *curr, const char *fn,
548 const char *fn_pdb, const int *molindex, const t_topology *top,
549 rvec *x, int ePBC, matrix box, const gmx_output_env_t *oenv)
555 real a, b, D, Dav, D2av, VarD, sqrtD, sqrtD_max, scale;
556 t_pdbinfo *pdbinfo = nullptr;
557 const int *mol2a = nullptr;
559 out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D", oenv);
563 pdbinfo = top->atoms.pdbinfo;
564 mol2a = top->mols.index;
569 for (i = 0; (i < curr->nmol); i++)
571 lsq1 = gmx_stats_init();
572 for (j = 0; (j < curr->nrestart); j++)
576 while (gmx_stats_get_point(curr->lsq[j][i], &xx, &yy, &dx, &dy, 0) == estatsOK)
578 gmx_stats_add_point(lsq1, xx, yy, dx, dy);
581 gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, nullptr, nullptr, nullptr, nullptr);
582 gmx_stats_free(lsq1);
583 D = a*FACTOR/curr->dim_factor;
589 D2av += gmx::square(D);
590 fprintf(out, "%10d %10g\n", i, D);
593 sqrtD = std::sqrt(D);
594 if (sqrtD > sqrtD_max)
598 for (j = mol2a[molindex[i]]; j < mol2a[molindex[i]+1]; j++)
600 pdbinfo[j].bfac = sqrtD;
605 do_view(oenv, fn, "-graphtype bar");
607 /* Compute variance, stddev and error */
610 VarD = D2av - gmx::square(Dav);
611 printf("<D> = %.4f Std. Dev. = %.4f Error = %.4f\n",
612 Dav, std::sqrt(VarD), std::sqrt(VarD/curr->nmol));
617 while (scale*sqrtD_max > 10)
621 while (scale*sqrtD_max < 0.1)
625 GMX_RELEASE_ASSERT(pdbinfo != nullptr, "Internal error - pdbinfo not set for PDB input");
626 for (i = 0; i < top->atoms.nr; i++)
628 pdbinfo[i].bfac *= scale;
630 write_sto_conf(fn_pdb, "molecular MSD", &top->atoms, x, nullptr, ePBC, box);
634 /* this is the main loop for the correlation type functions
635 * fx and nx are file pointers to things like read_first_x and
638 static int corr_loop(t_corr *curr, const char *fn, const t_topology *top, int ePBC,
639 gmx_bool bMol, int gnx[], int *index[],
640 t_calc_func *calc1, gmx_bool bTen, int *gnx_com, int *index_com[],
641 real dt, real t_pdb, rvec **x_pdb, matrix box_pdb,
642 const gmx_output_env_t *oenv)
644 rvec *x[2]; /* the coordinates to read */
645 rvec *xa[2]; /* the coordinates to calculate displacements for */
648 int natoms, i, j, cur = 0, maxframes = 0;
653 gmx_rmpbc_t gpbc = nullptr;
655 natoms = read_first_x(oenv, &status, fn, &curr->t0, &(x[cur]), box);
657 fprintf(stderr, "Read %d atoms for first frame\n", natoms);
659 if ((gnx_com != nullptr) && natoms < top->atoms.nr)
661 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);
664 snew(x[prev], natoms);
666 // if com is requested, the data structure needs to be large enough to do this
667 // to prevent overflow
668 if (bMol && !gnx_com)
670 curr->ncoords = curr->nmol;
671 snew(xa[0], curr->ncoords);
672 snew(xa[1], curr->ncoords);
676 curr->ncoords = natoms;
690 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
693 /* the loop over all frames */
696 if (x_pdb && ((bFirst && t_pdb < t) ||
698 t_pdb > t - 0.5*(t - t_prev) &&
699 t_pdb < t + 0.5*(t - t_prev))))
701 if (*x_pdb == nullptr)
703 snew(*x_pdb, natoms);
705 for (i = 0; i < natoms; i++)
707 copy_rvec(x[cur][i], (*x_pdb)[i]);
709 copy_mat(box, box_pdb);
713 /* check whether we've reached a restart point */
714 if (bRmod(t, curr->t0, dt))
718 srenew(curr->x0, curr->nrestart);
719 snew(curr->x0[curr->nrestart-1], curr->ncoords);
720 srenew(curr->com, curr->nrestart);
721 srenew(curr->n_offs, curr->nrestart);
722 srenew(curr->lsq, curr->nrestart);
723 snew(curr->lsq[curr->nrestart-1], curr->nmol);
724 for (i = 0; i < curr->nmol; i++)
726 curr->lsq[curr->nrestart-1][i] = gmx_stats_init();
731 fprintf(debug, "Extended data structures because of new restart %d\n",
735 /* create or extend the frame-based arrays */
736 if (curr->nframes >= maxframes-1)
740 for (i = 0; (i < curr->ngrp); i++)
742 curr->ndata[i] = nullptr;
743 curr->data[i] = nullptr;
746 curr->datam[i] = nullptr;
749 curr->time = nullptr;
752 for (i = 0; (i < curr->ngrp); i++)
754 srenew(curr->ndata[i], maxframes);
755 srenew(curr->data[i], maxframes);
758 srenew(curr->datam[i], maxframes);
760 for (j = maxframes-10; j < maxframes; j++)
762 curr->ndata[i][j] = 0;
763 curr->data[i][j] = 0;
766 clear_mat(curr->datam[i][j]);
770 srenew(curr->time, maxframes);
774 curr->time[curr->nframes] = t - curr->t0;
776 /* make the molecules whole */
779 gmx_rmpbc(gpbc, natoms, box, x[cur]);
782 /* calculate the molecules' centers of masses and put them into xa */
783 // NOTE and WARNING! If above both COM removal and individual molecules have been
784 // requested, x and xa point to the same memory, and the coordinate
785 // data becomes overwritten by the molecule data.
788 calc_mol_com(gnx[0], index[0], &top->mols, &top->atoms, x[cur], xa[cur]);
791 /* for the first frame, the previous frame is a copy of the first frame */
794 std::memcpy(xa[prev], xa[cur], curr->ncoords*sizeof(xa[prev][0]));
798 /* first remove the periodic boundary condition crossings */
799 for (i = 0; i < curr->ngrp; i++)
801 prep_data(bMol, gnx[i], index[i], xa[cur], xa[prev], box);
804 /* calculate the center of mass */
807 calc_com(bMol, gnx_com[0], index_com[0], xa[cur], xa[prev], box,
811 /* loop over all groups in index file */
812 for (i = 0; (i < curr->ngrp); i++)
814 /* calculate something useful, like mean square displacements */
815 calc_corr(curr, i, gnx[i], index[i], xa[cur], (gnx_com != nullptr), com,
823 while (read_next_x(oenv, status, &t, x[cur], box));
824 fprintf(stderr, "\nUsed %d restart points spaced %g %s over %g %s\n\n",
826 output_env_conv_time(oenv, dt), output_env_get_time_unit(oenv).c_str(),
827 output_env_conv_time(oenv, curr->time[curr->nframes-1]),
828 output_env_get_time_unit(oenv).c_str() );
832 gmx_rmpbc_done(gpbc);
840 static void index_atom2mol(int *n, int *index, const t_block *mols)
842 int nat, i, nmol, mol, j;
850 while (index[i] > mols->index[mol])
855 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
858 for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
860 if (i >= nat || index[i] != j)
862 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
869 fprintf(stderr, "Split group of %d atoms into %d molecules\n", nat, nmol);
874 static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file,
875 const char *mol_file, const char *pdb_file, real t_pdb,
876 int nrgrp, t_topology *top, int ePBC,
877 gmx_bool bTen, gmx_bool bMW, gmx_bool bRmCOMM,
878 int type, real dim_factor, int axis,
879 real dt, real beginfit, real endfit, const gmx_output_env_t *oenv)
882 int *gnx; /* the selected groups' sizes */
883 int **index; /* selected groups' indices */
885 int i, i0, i1, j, N, nat_trx;
886 real *DD, *SigmaD, a, a2, b, r, chi2;
889 int *gnx_com = nullptr; /* the COM removal group size */
890 int **index_com = nullptr; /* the COM removal group atom indices */
891 char **grpname_com = nullptr; /* the COM removal group name */
895 snew(grpname, nrgrp);
897 fprintf(stderr, "\nSelect a group to calculate mean squared displacement for:\n");
898 get_index(&top->atoms, ndx_file, nrgrp, gnx, index, grpname);
904 snew(grpname_com, 1);
906 fprintf(stderr, "\nNow select a group for center of mass removal:\n");
907 get_index(&top->atoms, ndx_file, 1, gnx_com, index_com, grpname_com);
912 index_atom2mol(&gnx[0], index[0], &top->mols);
915 msd = init_corr(nrgrp, type, axis, dim_factor,
916 mol_file == nullptr ? 0 : gnx[0], bTen, bMW, dt, top,
920 corr_loop(msd, trx_file, top, ePBC, mol_file ? gnx[0] : false, gnx, index,
921 (mol_file != nullptr) ? calc1_mol : (bMW ? calc1_mw : calc1_norm),
922 bTen, gnx_com, index_com, dt, t_pdb,
923 pdb_file ? &x : nullptr, box, oenv);
925 /* Correct for the number of points */
926 for (j = 0; (j < msd->ngrp); j++)
928 for (i = 0; (i < msd->nframes); i++)
930 msd->data[j][i] /= msd->ndata[j][i];
933 msmul(msd->datam[j][i], 1.0/msd->ndata[j][i], msd->datam[j][i]);
940 if (pdb_file && x == nullptr)
942 fprintf(stderr, "\nNo frame found need time tpdb = %g ps\n"
943 "Can not write %s\n\n", t_pdb, pdb_file);
946 top->atoms.nr = nat_trx;
947 if (pdb_file && top->atoms.pdbinfo == nullptr)
949 snew(top->atoms.pdbinfo, top->atoms.nr);
951 printmol(msd, mol_file, pdb_file, index[0], top, x, ePBC, box, oenv);
960 i0 = static_cast<int>(0.1*(msd->nframes - 1) + 0.5);
961 beginfit = msd->time[i0];
965 for (i0 = 0; i0 < msd->nframes && msd->time[i0] < beginfit; i0++)
973 i1 = static_cast<int>(0.9*(msd->nframes - 1) + 0.5) + 1;
974 endfit = msd->time[i1-1];
978 for (i1 = i0; i1 < msd->nframes && msd->time[i1] <= endfit; i1++)
983 fprintf(stdout, "Fitting from %g to %g %s\n\n", beginfit, endfit,
984 output_env_get_time_unit(oenv).c_str());
989 fprintf(stdout, "Not enough points for fitting (%d).\n"
990 "Can not determine the diffusion constant.\n", N);
995 snew(SigmaD, msd->ngrp);
996 for (j = 0; j < msd->ngrp; j++)
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] = std::abs(a-a2);
1008 lsq_y_ax_b(N, &(msd->time[i0]), &(msd->data[j][i0]), &(DD[j]), &b, &r, &chi2);
1009 DD[j] *= FACTOR/msd->dim_factor;
1010 SigmaD[j] *= FACTOR/msd->dim_factor;
1011 if (DD[j] > 0.01 && DD[j] < 1e4)
1013 fprintf(stdout, "D[%10s] %.4f (+/- %.4f) 1e-5 cm^2/s\n",
1014 grpname[j], DD[j], SigmaD[j]);
1018 fprintf(stdout, "D[%10s] %.4g (+/- %.4g) 1e-5 cm^2/s\n",
1019 grpname[j], DD[j], SigmaD[j]);
1023 /* Print mean square displacement */
1024 corr_print(msd, bTen, msd_file,
1025 "Mean Square Displacement",
1027 msd->time[msd->nframes-1], beginfit, endfit, DD, SigmaD, grpname, oenv);
1030 int gmx_msd(int argc, char *argv[])
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 ",
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] is -1, fitting starts at 10%",
1065 "and when [TT]-endfit[tt] is -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 [REF].pdb[ref] 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]."
1076 static const char *normtype[] = { nullptr, "no", "x", "y", "z", nullptr };
1077 static const char *axtitle[] = { nullptr, "no", "x", "y", "z", nullptr };
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 = FALSE;
1084 static gmx_bool bMW = TRUE;
1085 static gmx_bool bRmCOMM = FALSE;
1087 { "-type", FALSE, etENUM, {normtype},
1088 "Compute diffusion coefficient in one direction" },
1089 { "-lateral", FALSE, etENUM, {axtitle},
1090 "Calculate the lateral diffusion in a plane perpendicular to" },
1091 { "-ten", FALSE, etBOOL, {&bTen},
1092 "Calculate the full tensor" },
1093 { "-ngroup", FALSE, etINT, {&ngroup},
1094 "Number of groups to calculate MSD for" },
1095 { "-mw", FALSE, etBOOL, {&bMW},
1096 "Mass weighted MSD" },
1097 { "-rmcomm", FALSE, etBOOL, {&bRmCOMM},
1098 "Remove center of mass motion" },
1099 { "-tpdb", FALSE, etTIME, {&t_pdb},
1100 "The frame to use for option [TT]-pdb[tt] (%t)" },
1101 { "-trestart", FALSE, etTIME, {&dt},
1102 "Time between restarting points in trajectory (%t)" },
1103 { "-beginfit", FALSE, etTIME, {&beginfit},
1104 "Start time for fitting the MSD (%t), -1 is 10%" },
1105 { "-endfit", FALSE, etTIME, {&endfit},
1106 "End time for fitting the MSD (%t), -1 is 90%" }
1110 { efTRX, nullptr, nullptr, ffREAD },
1111 { efTPS, nullptr, nullptr, ffREAD },
1112 { efNDX, nullptr, nullptr, ffOPTRD },
1113 { efXVG, nullptr, "msd", ffWRITE },
1114 { efXVG, "-mol", "diff_mol", ffOPTWR },
1115 { efPDB, "-pdb", "diff_mol", ffOPTWR }
1117 #define NFILE asize(fnm)
1122 const char *trx_file, *tps_file, *ndx_file, *msd_file, *mol_file, *pdb_file;
1127 gmx_output_env_t *oenv;
1129 if (!parse_common_args(&argc, argv,
1130 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_TIME_UNIT,
1131 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
1135 trx_file = ftp2fn_null(efTRX, NFILE, fnm);
1136 tps_file = ftp2fn_null(efTPS, NFILE, fnm);
1137 ndx_file = ftp2fn_null(efNDX, NFILE, fnm);
1138 msd_file = ftp2fn_null(efXVG, NFILE, fnm);
1139 pdb_file = opt2fn_null("-pdb", NFILE, fnm);
1142 mol_file = opt2fn("-mol", NFILE, fnm);
1146 mol_file = opt2fn_null("-mol", NFILE, fnm);
1151 gmx_fatal(FARGS, "Must have at least 1 group (now %d)", ngroup);
1153 if (mol_file && ngroup > 1)
1155 gmx_fatal(FARGS, "With molecular msd can only have 1 group (now %d)",
1163 fprintf(stderr, "Calculating diffusion coefficients for molecules.\n");
1166 GMX_RELEASE_ASSERT(normtype[0] != nullptr, "Options inconsistency; normtype[0] is NULL");
1167 GMX_RELEASE_ASSERT(axtitle[0] != nullptr, "Options inconsistency; axtitle[0] is NULL");
1169 if (normtype[0][0] != 'n')
1171 type = normtype[0][0] - 'x' + X; /* See defines above */
1179 if ((type == NORMAL) && (axtitle[0][0] != 'n'))
1183 axis = (axtitle[0][0] - 'x'); /* See defines above */
1190 if (bTen && type != NORMAL)
1192 gmx_fatal(FARGS, "Can only calculate the full tensor for 3D msd");
1195 bTop = read_tps_conf(tps_file, &top, &ePBC, &xdum, nullptr, box, bMW || bRmCOMM);
1196 if (mol_file && !bTop)
1199 "Could not read a topology from %s. Try a tpr file instead.",
1203 do_corr(trx_file, ndx_file, msd_file, mol_file, pdb_file, t_pdb, ngroup,
1204 &top, ePBC, bTen, bMW, bRmCOMM, type, dim_factor, axis, dt, beginfit, endfit,
1207 view_all(oenv, NFILE, fnm);