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