src/gromacs/gmxana/gmx_msd.c

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