src/gromacs/gmxana/gmx_dipoles.cpp

   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 <string.h>
  40 #include <math.h>
  41
  42 #include <algorithm>
  43
  44 #include "gromacs/legacyheaders/macros.h"
  45 #include "gromacs/pbcutil/pbc.h"
  46 #include "gromacs/utility/futil.h"
  47 #include "gromacs/legacyheaders/viewit.h"
  48 #include "gromacs/legacyheaders/txtdump.h"
  49 #include "gromacs/statistics/statistics.h"
  50 #include "gstat.h"
  51 #include "gromacs/topology/index.h"
  52 #include "gromacs/random/random.h"
  53 #include "gromacs/legacyheaders/names.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/legacyheaders/calcmu.h"
  56 #include "gromacs/fileio/enxio.h"
  57 #include "gmx_ana.h"
  58 #include "gromacs/legacyheaders/copyrite.h"
  59 #include "gromacs/fileio/trxio.h"
  60
  61 #include "gromacs/bonded/bonded.h"
  62 #include "gromacs/commandline/pargs.h"
  63 #include "gromacs/fileio/matio.h"
  64 #include "gromacs/fileio/xvgr.h"
  65 #include "gromacs/linearalgebra/nrjac.h"
  66 #include "gromacs/math/vec.h"
  67 #include "gromacs/pbcutil/rmpbc.h"
  68 #include "gromacs/utility/exceptions.h"
  69 #include "gromacs/utility/fatalerror.h"
  70 #include "gromacs/utility/smalloc.h"
  71
  72 #define e2d(x) ENM2DEBYE*(x)
  73 #define EANG2CM  E_CHARGE*1.0e-10       /* e Angstrom to Coulomb meter */
  74 #define CM2D  SPEED_OF_LIGHT*1.0e+24    /* Coulomb meter to Debye */
  75
  76 typedef struct {
  77     int      nelem;
  78     real     spacing, radius;
  79     real    *elem;
  80     int     *count;
  81     gmx_bool bPhi;
  82     int      nx, ny;
  83     real   **cmap;
  84 } t_gkrbin;
  85
  86 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
  87 {
  88     t_gkrbin *gb;
  89     char     *ptr;
  90     int       i;
  91
  92     snew(gb, 1);
  93
  94     if ((ptr = getenv("GMX_DIPOLE_SPACING")) != NULL)
  95     {
  96         double bw = strtod(ptr, NULL);
  97         gb->spacing = bw;
  98     }
  99     else
 100     {
 101         gb->spacing = 0.01; /* nm */
 102     }
 103     gb->nelem   = 1 + static_cast<int>(radius/gb->spacing);
 104     if (rcmax == 0)
 105     {
 106         gb->nx = gb->nelem;
 107     }
 108     else
 109     {
 110         gb->nx = 1 + static_cast<int>(rcmax/gb->spacing);
 111     }
 112     gb->radius  = radius;
 113     snew(gb->elem, gb->nelem);
 114     snew(gb->count, gb->nelem);
 115
 116     snew(gb->cmap, gb->nx);
 117     gb->ny = std::max(2, ndegrees);
 118     for (i = 0; (i < gb->nx); i++)
 119     {
 120         snew(gb->cmap[i], gb->ny);
 121     }
 122     gb->bPhi = bPhi;
 123
 124     return gb;
 125 }
 126
 127 static void done_gkrbin(t_gkrbin **gb)
 128 {
 129     sfree((*gb)->elem);
 130     sfree((*gb)->count);
 131     sfree((*gb));
 132     *gb = NULL;
 133 }
 134
 135 static void add2gkr(t_gkrbin *gb, real r, real cosa, real phi)
 136 {
 137     int  cy, index = gmx_nint(r/gb->spacing);
 138     real alpha;
 139
 140     if (index < gb->nelem)
 141     {
 142         gb->elem[index]  += cosa;
 143         gb->count[index]++;
 144     }
 145     if (index < gb->nx)
 146     {
 147         alpha = acos(cosa);
 148         if (gb->bPhi)
 149         {
 150             cy = static_cast<int>((M_PI+phi)*gb->ny/(2*M_PI));
 151         }
 152         else
 153         {
 154             cy = static_cast<int>((alpha*gb->ny)/M_PI); /*((1+cosa)*0.5*(gb->ny));*/
 155         }
 156         cy = std::min(gb->ny-1, std::max(0, cy));
 157         if (debug)
 158         {
 159             fprintf(debug, "CY: %10f  %5d\n", alpha, cy);
 160         }
 161         gb->cmap[index][cy] += 1;
 162     }
 163 }
 164
 165 static void rvec2sprvec(rvec dipcart, rvec dipsp)
 166 {
 167     clear_rvec(dipsp);
 168     dipsp[0] = sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]+dipcart[ZZ]*dipcart[ZZ]); /* R */
 169     dipsp[1] = atan2(dipcart[YY], dipcart[XX]);                                               /* Theta */
 170     dipsp[2] = atan2(sqrt(dipcart[XX]*dipcart[XX]+dipcart[YY]*dipcart[YY]), dipcart[ZZ]);     /* Phi */
 171 }
 172
 173
 174
 175 void do_gkr(t_gkrbin *gb, int ncos, int *ngrp, int *molindex[],
 176             int mindex[], rvec x[], rvec mu[],
 177             int ePBC, matrix box, t_atom *atom, int *nAtom)
 178 {
 179     static rvec *xcm[2] = { NULL, NULL};
 180     int          gi, gj, j0, j1, i, j, k, n, grp0, grp1;
 181     real         qtot, q, cosa, rr, phi;
 182     rvec         dx;
 183     t_pbc        pbc;
 184
 185     for (n = 0; (n < ncos); n++)
 186     {
 187         if (!xcm[n])
 188         {
 189             snew(xcm[n], ngrp[n]);
 190         }
 191         for (i = 0; (i < ngrp[n]); i++)
 192         {
 193             /* Calculate center of mass of molecule */
 194             gi = molindex[n][i];
 195             j0 = mindex[gi];
 196
 197             if (nAtom[n] > 0)
 198             {
 199                 copy_rvec(x[j0+nAtom[n]-1], xcm[n][i]);
 200             }
 201             else
 202             {
 203                 j1 = mindex[gi+1];
 204                 clear_rvec(xcm[n][i]);
 205                 qtot = 0;
 206                 for (j = j0; j < j1; j++)
 207                 {
 208                     q     = fabs(atom[j].q);
 209                     qtot += q;
 210                     for (k = 0; k < DIM; k++)
 211                     {
 212                         xcm[n][i][k] += q*x[j][k];
 213                     }
 214                 }
 215                 svmul(1/qtot, xcm[n][i], xcm[n][i]);
 216             }
 217         }
 218     }
 219     set_pbc(&pbc, ePBC, box);
 220     grp0 = 0;
 221     grp1 = ncos-1;
 222     for (i = 0; i < ngrp[grp0]; i++)
 223     {
 224         gi = molindex[grp0][i];
 225         j0 = (ncos == 2) ? 0 : i+1;
 226         for (j = j0; j < ngrp[grp1]; j++)
 227         {
 228             gj   = molindex[grp1][j];
 229             if ((iprod(mu[gi], mu[gi]) > 0) && (iprod(mu[gj], mu[gj]) > 0))
 230             {
 231                 /* Compute distance between molecules including PBC */
 232                 pbc_dx(&pbc, xcm[grp0][i], xcm[grp1][j], dx);
 233                 rr = norm(dx);
 234
 235                 if (gb->bPhi)
 236                 {
 237                     rvec xi, xj, xk, xl;
 238                     rvec r_ij, r_kj, r_kl, mm, nn;
 239                     real sign;
 240                     int  t1, t2, t3;
 241
 242                     copy_rvec(xcm[grp0][i], xj);
 243                     copy_rvec(xcm[grp1][j], xk);
 244                     rvec_add(xj, mu[gi], xi);
 245                     rvec_add(xk, mu[gj], xl);
 246                     phi = dih_angle(xi, xj, xk, xl, &pbc,
 247                                     r_ij, r_kj, r_kl, mm, nn, /* out */
 248                                     &sign, &t1, &t2, &t3);
 249                     cosa = cos(phi);
 250                 }
 251                 else
 252                 {
 253                     cosa = cos_angle(mu[gi], mu[gj]);
 254                     phi  = 0;
 255                 }
 256                 if (debug || gmx_isnan(cosa))
 257                 {
 258                     fprintf(debug ? debug : stderr,
 259                             "mu[%d] = %5.2f %5.2f %5.2f |mi| = %5.2f, mu[%d] = %5.2f %5.2f %5.2f |mj| = %5.2f rr = %5.2f cosa = %5.2f\n",
 260                             gi, mu[gi][XX], mu[gi][YY], mu[gi][ZZ], norm(mu[gi]),
 261                             gj, mu[gj][XX], mu[gj][YY], mu[gj][ZZ], norm(mu[gj]),
 262                             rr, cosa);
 263                 }
 264
 265                 add2gkr(gb, rr, cosa, phi);
 266             }
 267         }
 268     }
 269 }
 270
 271 static real normalize_cmap(t_gkrbin *gb)
 272 {
 273     int    i, j;
 274     real   hi;
 275     double vol;
 276
 277     hi = 0;
 278     for (i = 0; (i < gb->nx); i++)
 279     {
 280         vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
 281         for (j = 0; (j < gb->ny); j++)
 282         {
 283             gb->cmap[i][j] /= vol;
 284             hi              = std::max(hi, gb->cmap[i][j]);
 285         }
 286     }
 287     if (hi <= 0)
 288     {
 289         gmx_fatal(FARGS, "No data in the cmap");
 290     }
 291     return hi;
 292 }
 293
 294 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
 295 {
 296     FILE   *out;
 297     int     i, j;
 298     real    hi;
 299
 300     real   *xaxis, *yaxis;
 301     t_rgb   rlo = { 1, 1, 1 };
 302     t_rgb   rhi = { 0, 0, 0 };
 303
 304     hi = normalize_cmap(gb);
 305     snew(xaxis, gb->nx+1);
 306     for (i = 0; (i < gb->nx+1); i++)
 307     {
 308         xaxis[i] = i*gb->spacing;
 309     }
 310     snew(yaxis, gb->ny);
 311     for (j = 0; (j < gb->ny); j++)
 312     {
 313         if (gb->bPhi)
 314         {
 315             yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
 316         }
 317         else
 318         {
 319             yaxis[j] = (180.0*j)/(gb->ny-1.0);
 320         }
 321         /*2.0*j/(gb->ny-1.0)-1.0;*/
 322     }
 323     out = gmx_ffopen(cmap, "w");
 324     write_xpm(out, 0,
 325               "Dipole Orientation Distribution", "Fraction", "r (nm)",
 326               gb->bPhi ? "Phi" : "Alpha",
 327               gb->nx, gb->ny, xaxis, yaxis,
 328               gb->cmap, 0, hi, rlo, rhi, nlevels);
 329     gmx_ffclose(out);
 330     sfree(xaxis);
 331     sfree(yaxis);
 332 }
 333
 334 static void print_gkrbin(const char *fn, t_gkrbin *gb,
 335                          int ngrp, int nframes, real volume,
 336                          const output_env_t oenv)
 337 {
 338     /* We compute Gk(r), gOO and hOO according to
 339      * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
 340      * In this implementation the angle between dipoles is stored
 341      * rather than their inner product. This allows to take polarizible
 342      * models into account. The RDF is calculated as well, almost for free!
 343      */
 344     FILE       *fp;
 345     const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
 346     int         i, last;
 347     real        x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
 348     double      fac;
 349
 350     fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
 351     xvgr_legend(fp, asize(leg), leg, oenv);
 352
 353     Gkr = 1;  /* Self-dipole inproduct = 1 */
 354     rho = ngrp/volume;
 355
 356     if (debug)
 357     {
 358         fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
 359         fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
 360     }
 361
 362     last = gb->nelem-1;
 363     while (last > 1 && gb->elem[last-1] == 0)
 364     {
 365         last--;
 366     }
 367
 368     /* Divide by dipole squared, by number of frames, by number of origins.
 369      * Multiply by 2 because we only take half the matrix of interactions
 370      * into account.
 371      */
 372     fac  = 2.0/((double) ngrp * (double) nframes);
 373
 374     x0 = 0;
 375     for (i = 0; i < last; i++)
 376     {
 377         /* Centre of the coordinate in the spherical layer */
 378         x1    = x0+gb->spacing;
 379
 380         /* Volume of the layer */
 381         vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
 382
 383         /* gOO */
 384         gOO   = gb->count[i]*fac/(rho*vol_s);
 385
 386         /* Dipole correlation hOO, normalized by the relative number density, like
 387          * in a Radial distribution function.
 388          */
 389         ggg  = gb->elem[i]*fac;
 390         hOO  = 3.0*ggg/(rho*vol_s);
 391         Gkr += ggg;
 392         if (gb->count[i])
 393         {
 394             cosav = gb->elem[i]/gb->count[i];
 395         }
 396         else
 397         {
 398             cosav = 0;
 399         }
 400         ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
 401
 402         fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e  %12.5e\n",
 403                 x1, Gkr, cosav, hOO, gOO, ener);
 404
 405         /* Swap x0 and x1 */
 406         x0 = x1;
 407     }
 408     gmx_ffclose(fp);
 409 }
 410
 411 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
 412                           real *t, int nre, t_enxframe *fr)
 413 {
 414     int          i;
 415     gmx_bool     bCont;
 416     char         buf[22];
 417
 418     bCont = do_enx(fmu, fr);
 419     if (fr->nre != nre)
 420     {
 421         fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
 422                 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
 423     }
 424
 425     if (bCont)
 426     {
 427         if (Vol != -1)          /* we've got Volume in the energy file */
 428         {
 429             *vol = fr->ener[Vol].e;
 430         }
 431         for (i = 0; i < DIM; i++)
 432         {
 433             mu[i] = fr->ener[iMu[i]].e;
 434         }
 435         *t = fr->t;
 436     }
 437
 438     return bCont;
 439 }
 440
 441 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
 442 {
 443     double mtot, qtot;
 444     int    ncharged, m, a0, a1, a;
 445
 446     ncharged = 0;
 447     for (m = 0; m < n; m++)
 448     {
 449         a0   = mols->index[index[m]];
 450         a1   = mols->index[index[m]+1];
 451         mtot = 0;
 452         qtot = 0;
 453         for (a = a0; a < a1; a++)
 454         {
 455             mtot += atom[a].m;
 456             qtot += atom[a].q;
 457         }
 458         /* This check is only for the count print */
 459         if (fabs(qtot) > 0.01)
 460         {
 461             ncharged++;
 462         }
 463         if (mtot > 0)
 464         {
 465             /* Remove the net charge at the center of mass */
 466             for (a = a0; a < a1; a++)
 467             {
 468                 atom[a].q -= qtot*atom[a].m/mtot;
 469             }
 470         }
 471     }
 472
 473     if (ncharged)
 474     {
 475         printf("There are %d charged molecules in the selection,\n"
 476                "will subtract their charge at their center of mass\n", ncharged);
 477     }
 478 }
 479
 480 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
 481 {
 482     int  k, m;
 483     real q;
 484
 485     clear_rvec(mu);
 486     for (k = k0; (k < k1); k++)
 487     {
 488         q  = e2d(atom[k].q);
 489         for (m = 0; (m < DIM); m++)
 490         {
 491             mu[m] += q*x[k][m];
 492         }
 493     }
 494 }
 495
 496 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
 497 {
 498     int      i, k, m, n, niter;
 499     real     q, r2, mass, masstot;
 500     rvec     com;        /* center of mass */
 501     rvec     r;          /* distance of atoms to center of mass */
 502     double **inten;
 503     double   dd[DIM], **ev, tmp;
 504
 505     snew(inten, DIM);
 506     snew(ev, DIM);
 507     for (i = 0; (i < DIM); i++)
 508     {
 509         snew(inten[i], DIM);
 510         snew(ev[i], DIM);
 511         dd[i] = 0.0;
 512     }
 513
 514     /* Compute center of mass */
 515     clear_rvec(com);
 516     masstot = 0.0;
 517     for (k = k0; (k < k1); k++)
 518     {
 519         mass     = atom[k].m;
 520         masstot += mass;
 521         for (i = 0; (i < DIM); i++)
 522         {
 523             com[i] += mass*x[k][i];
 524         }
 525     }
 526     svmul((1.0/masstot), com, com);
 527
 528     /* We want traceless quadrupole moments, so let us calculate the complete
 529      * quadrupole moment tensor and diagonalize this tensor to get
 530      * the individual components on the diagonal.
 531      */
 532
 533 #define delta(a, b) (( a == b ) ? 1.0 : 0.0)
 534
 535     for (m = 0; (m < DIM); m++)
 536     {
 537         for (n = 0; (n < DIM); n++)
 538         {
 539             inten[m][n] = 0;
 540         }
 541     }
 542     for (k = k0; (k < k1); k++)         /* loop over atoms in a molecule */
 543     {
 544         q  = (atom[k].q)*100.0;
 545         rvec_sub(x[k], com, r);
 546         r2 = iprod(r, r);
 547         for (m = 0; (m < DIM); m++)
 548         {
 549             for (n = 0; (n < DIM); n++)
 550             {
 551                 inten[m][n] += 0.5*q*(3.0*r[m]*r[n] - r2*delta(m, n))*EANG2CM*CM2D;
 552             }
 553         }
 554     }
 555     if (debug)
 556     {
 557         for (i = 0; (i < DIM); i++)
 558         {
 559             fprintf(debug, "Q[%d] = %8.3f  %8.3f  %8.3f\n",
 560                     i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
 561         }
 562     }
 563
 564     /* We've got the quadrupole tensor, now diagonalize the sucker */
 565     jacobi(inten, 3, dd, ev, &niter);
 566
 567     if (debug)
 568     {
 569         for (i = 0; (i < DIM); i++)
 570         {
 571             fprintf(debug, "ev[%d] = %8.3f  %8.3f  %8.3f\n",
 572                     i, ev[i][XX], ev[i][YY], ev[i][ZZ]);
 573         }
 574         for (i = 0; (i < DIM); i++)
 575         {
 576             fprintf(debug, "Q'[%d] = %8.3f  %8.3f  %8.3f\n",
 577                     i, inten[i][XX], inten[i][YY], inten[i][ZZ]);
 578         }
 579     }
 580     /* Sort the eigenvalues, for water we know that the order is as follows:
 581      *
 582      * Q_yy, Q_zz, Q_xx
 583      *
 584      * At the moment I have no idea how this will work out for other molecules...
 585      */
 586
 587 #define SWAP(i)                                 \
 588     if (dd[i+1] > dd[i]) {                      \
 589         tmp     = dd[i];                              \
 590         dd[i]   = dd[i+1];                          \
 591         dd[i+1] = tmp;                            \
 592     }
 593     SWAP(0);
 594     SWAP(1);
 595     SWAP(0);
 596
 597     for (m = 0; (m < DIM); m++)
 598     {
 599         quad[0] = dd[2];  /* yy */
 600         quad[1] = dd[0];  /* zz */
 601         quad[2] = dd[1];  /* xx */
 602     }
 603
 604     if (debug)
 605     {
 606         pr_rvec(debug, 0, "Quadrupole", quad, DIM, TRUE);
 607     }
 608
 609     /* clean-up */
 610     for (i = 0; (i < DIM); i++)
 611     {
 612         sfree(inten[i]);
 613         sfree(ev[i]);
 614     }
 615     sfree(inten);
 616     sfree(ev);
 617 }
 618
 619 /*
 620  * Calculates epsilon according to M. Neumann, Mol. Phys. 50, 841 (1983)
 621  */
 622 real calc_eps(double M_diff, double volume, double epsRF, double temp)
 623 {
 624     double eps, A, teller, noemer;
 625     double eps_0 = 8.854187817e-12;   /* epsilon_0 in C^2 J^-1 m^-1 */
 626     double fac   = 1.112650021e-59;   /* converts Debye^2 to C^2 m^2 */
 627
 628     A = M_diff*fac/(3*eps_0*volume*NANO*NANO*NANO*BOLTZMANN*temp);
 629
 630     if (epsRF == 0.0)
 631     {
 632         teller = 1 + A;
 633         noemer = 1;
 634     }
 635     else
 636     {
 637         teller = 1 + (A*2*epsRF/(2*epsRF+1));
 638         noemer = 1 - (A/(2*epsRF+1));
 639     }
 640     eps = teller / noemer;
 641
 642     return eps;
 643 }
 644
 645 static void update_slab_dipoles(int k0, int k1, rvec x[], rvec mu,
 646                                 int idim, int nslice, rvec slab_dipole[],
 647                                 matrix box)
 648 {
 649     int  k;
 650     real xdim = 0;
 651
 652     for (k = k0; (k < k1); k++)
 653     {
 654         xdim += x[k][idim];
 655     }
 656     xdim /= (k1-k0);
 657     k     = ((int)(xdim*nslice/box[idim][idim] + nslice)) % nslice;
 658     rvec_inc(slab_dipole[k], mu);
 659 }
 660
 661 static void dump_slab_dipoles(const char *fn, int idim, int nslice,
 662                               rvec slab_dipole[], matrix box, int nframes,
 663                               const output_env_t oenv)
 664 {
 665     FILE       *fp;
 666     char        buf[STRLEN];
 667     int         i;
 668     real        mutot;
 669     const char *leg_dim[4] = {
 670         "\\f{12}m\\f{4}\\sX\\N",
 671         "\\f{12}m\\f{4}\\sY\\N",
 672         "\\f{12}m\\f{4}\\sZ\\N",
 673         "\\f{12}m\\f{4}\\stot\\N"
 674     };
 675
 676     sprintf(buf, "Box-%c (nm)", 'X'+idim);
 677     fp = xvgropen(fn, "Average dipole moment per slab", buf, "\\f{12}m\\f{4} (D)",
 678                   oenv);
 679     xvgr_legend(fp, DIM, leg_dim, oenv);
 680     for (i = 0; (i < nslice); i++)
 681     {
 682         mutot = norm(slab_dipole[i])/nframes;
 683         fprintf(fp, "%10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
 684                 ((i+0.5)*box[idim][idim])/nslice,
 685                 slab_dipole[i][XX]/nframes,
 686                 slab_dipole[i][YY]/nframes,
 687                 slab_dipole[i][ZZ]/nframes,
 688                 mutot);
 689     }
 690     gmx_ffclose(fp);
 691     do_view(oenv, fn, "-autoscale xy -nxy");
 692 }
 693
 694 static void compute_avercos(int n, rvec dip[], real *dd, rvec axis, gmx_bool bPairs)
 695 {
 696     int    i, j, k;
 697     double dc, d, ddc1, ddc2, ddc3;
 698     rvec   xxx = { 1, 0, 0 };
 699     rvec   yyy = { 0, 1, 0 };
 700     rvec   zzz = { 0, 0, 1 };
 701
 702     d    = 0;
 703     ddc1 = ddc2 = ddc3 = 0;
 704     for (i = k = 0; (i < n); i++)
 705     {
 706         ddc1 += fabs(cos_angle(dip[i], xxx));
 707         ddc2 += fabs(cos_angle(dip[i], yyy));
 708         ddc3 += fabs(cos_angle(dip[i], zzz));
 709         if (bPairs)
 710         {
 711             for (j = i+1; (j < n); j++, k++)
 712             {
 713                 dc  = cos_angle(dip[i], dip[j]);
 714                 d  += fabs(dc);
 715             }
 716         }
 717     }
 718     *dd      = d/k;
 719     axis[XX] = ddc1/n;
 720     axis[YY] = ddc2/n;
 721     axis[ZZ] = ddc3/n;
 722 }
 723
 724 static void do_dip(t_topology *top, int ePBC, real volume,
 725                    const char *fn,
 726                    const char *out_mtot, const char *out_eps,
 727                    const char *out_aver, const char *dipdist,
 728                    const char *cosaver, const char *fndip3d,
 729                    const char *fnadip,  gmx_bool bPairs,
 730                    const char *corrtype, const char *corf,
 731                    gmx_bool bGkr,     const char *gkrfn,
 732                    gmx_bool bPhi,     int  *nlevels,  int ndegrees,
 733                    int  ncos,
 734                    const char *cmap,    real rcmax,
 735                    gmx_bool bQuad,
 736                    gmx_bool bMU,      const char *mufn,
 737                    int  *gnx,     int  *molindex[],
 738                    real mu_max,   real mu_aver,
 739                    real epsilonRF, real temp,
 740                    int  *gkatom,  int skip,
 741                    gmx_bool bSlab,    int nslices,
 742                    const char *axtitle, const char *slabfn,
 743                    const output_env_t oenv)
 744 {
 745     const char *leg_mtot[] = {
 746         "M\\sx \\N",
 747         "M\\sy \\N",
 748         "M\\sz \\N",
 749         "|M\\stot \\N|"
 750     };
 751 #define NLEGMTOT asize(leg_mtot)
 752     const char *leg_eps[] = {
 753         "epsilon",
 754         "G\\sk",
 755         "g\\sk"
 756     };
 757 #define NLEGEPS asize(leg_eps)
 758     const char *leg_aver[] = {
 759         "< |M|\\S2\\N >",
 760         "< |M| >\\S2\\N",
 761         "< |M|\\S2\\N > - < |M| >\\S2\\N",
 762         "< |M| >\\S2\\N / < |M|\\S2\\N >"
 763     };
 764 #define NLEGAVER asize(leg_aver)
 765     const char *leg_cosaver[] = {
 766         "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
 767         "RMSD cos",
 768         "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
 769         "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
 770         "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
 771     };
 772 #define NLEGCOSAVER asize(leg_cosaver)
 773     const char *leg_adip[] = {
 774         "<mu>",
 775         "Std. Dev.",
 776         "Error"
 777     };
 778 #define NLEGADIP asize(leg_adip)
 779
 780     FILE         *outdd, *outmtot, *outaver, *outeps, *caver = NULL;
 781     FILE         *dip3d = NULL, *adip = NULL;
 782     rvec         *x, *dipole = NULL, mu_t, quad, *dipsp = NULL;
 783     t_gkrbin     *gkrbin = NULL;
 784     gmx_enxnm_t  *enm    = NULL;
 785     t_enxframe   *fr;
 786     int           nframes = 1000, nre, timecheck = 0, ncolour = 0;
 787     ener_file_t   fmu     = NULL;
 788     int           i, n, m, natom = 0, gnx_tot, teller, tel3;
 789     t_trxstatus  *status;
 790     int          *dipole_bin, ndipbin, ibin, iVol, idim = -1;
 791     unsigned long mode;
 792     real          rcut = 0, t, t0, t1, dt, dd, rms_cos;
 793     rvec          dipaxis;
 794     matrix        box;
 795     gmx_bool      bCorr, bTotal, bCont;
 796     double        M_diff = 0, epsilon, invtel, vol_aver;
 797     double        mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
 798     double        M[3], M2[3], M4[3], Gk = 0, g_k = 0;
 799     gmx_stats_t  *Qlsq, mulsq, muframelsq = NULL;
 800     ivec          iMu;
 801     real        **muall        = NULL;
 802     rvec         *slab_dipoles = NULL;
 803     t_atom       *atom         = NULL;
 804     t_block      *mols         = NULL;
 805     gmx_rmpbc_t   gpbc         = NULL;
 806
 807     gnx_tot = gnx[0];
 808     if (ncos > 1)
 809     {
 810         gnx_tot += gnx[1];
 811     }
 812
 813     vol_aver = 0.0;
 814
 815     iVol = -1;
 816
 817     /* initialize to a negative value so we can see that it wasn't picked up */
 818     iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
 819     if (bMU)
 820     {
 821         fmu = open_enx(mufn, "r");
 822         do_enxnms(fmu, &nre, &enm);
 823
 824         /* Determine the indexes of the energy grps we need */
 825         for (i = 0; (i < nre); i++)
 826         {
 827             if (strstr(enm[i].name, "Volume"))
 828             {
 829                 iVol = i;
 830             }
 831             else if (strstr(enm[i].name, "Mu-X"))
 832             {
 833                 iMu[XX] = i;
 834             }
 835             else if (strstr(enm[i].name, "Mu-Y"))
 836             {
 837                 iMu[YY] = i;
 838             }
 839             else if (strstr(enm[i].name, "Mu-Z"))
 840             {
 841                 iMu[ZZ] = i;
 842             }
 843         }
 844         if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
 845         {
 846             gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
 847         }
 848     }
 849     else
 850     {
 851         atom = top->atoms.atom;
 852         mols = &(top->mols);
 853     }
 854     if ((iVol == -1) && bMU)
 855     {
 856         printf("Using Volume from topology: %g nm^3\n", volume);
 857     }
 858
 859     /* Correlation stuff */
 860     bCorr  = (corrtype[0] != 'n');
 861     bTotal = (corrtype[0] == 't');
 862     if (bCorr)
 863     {
 864         if (bTotal)
 865         {
 866             snew(muall, 1);
 867             snew(muall[0], nframes*DIM);
 868         }
 869         else
 870         {
 871             snew(muall, gnx[0]);
 872             for (i = 0; (i < gnx[0]); i++)
 873             {
 874                 snew(muall[i], nframes*DIM);
 875             }
 876         }
 877     }
 878
 879     /* Allocate array which contains for every molecule in a frame the
 880      * dipole moment.
 881      */
 882     if (!bMU)
 883     {
 884         snew(dipole, gnx_tot);
 885     }
 886
 887     /* Statistics */
 888     snew(Qlsq, DIM);
 889     for (i = 0; (i < DIM); i++)
 890     {
 891         Qlsq[i] = gmx_stats_init();
 892     }
 893     mulsq = gmx_stats_init();
 894
 895     /* Open all the files */
 896     outmtot = xvgropen(out_mtot,
 897                        "Total dipole moment of the simulation box vs. time",
 898                        "Time (ps)", "Total Dipole Moment (Debye)", oenv);
 899     outeps  = xvgropen(out_eps, "Epsilon and Kirkwood factors",
 900                        "Time (ps)", "", oenv);
 901     outaver = xvgropen(out_aver, "Total dipole moment",
 902                        "Time (ps)", "D", oenv);
 903     if (bSlab)
 904     {
 905         idim = axtitle[0] - 'X';
 906         if ((idim < 0) || (idim >= DIM))
 907         {
 908             idim = axtitle[0] - 'x';
 909         }
 910         if ((idim < 0) || (idim >= DIM))
 911         {
 912             bSlab = FALSE;
 913         }
 914         if (nslices < 2)
 915         {
 916             bSlab = FALSE;
 917         }
 918         fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
 919                 axtitle, nslices, idim);
 920         if (bSlab)
 921         {
 922             snew(slab_dipoles, nslices);
 923             fprintf(stderr, "Doing slab analysis\n");
 924         }
 925     }
 926
 927     if (fnadip)
 928     {
 929         adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
 930         xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
 931
 932     }
 933     if (cosaver)
 934     {
 935         caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
 936                          "Average absolute dipole orientation", "Time (ps)", "", oenv);
 937         xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
 938                     oenv);
 939     }
 940
 941     if (fndip3d)
 942     {
 943         snew(dipsp, gnx_tot);
 944
 945         /* we need a dummy file for gnuplot */
 946         dip3d = (FILE *)gmx_ffopen("dummy.dat", "w");
 947         fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
 948         gmx_ffclose(dip3d);
 949
 950         dip3d = (FILE *)gmx_ffopen(fndip3d, "w");
 951         try
 952         {
 953             gmx::BinaryInformationSettings settings;
 954             settings.generatedByHeader(true);
 955             settings.linePrefix("# ");
 956             gmx::printBinaryInformation(dip3d, output_env_get_program_context(oenv),
 957                                         settings);
 958         }
 959         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 960     }
 961
 962     /* Write legends to all the files */
 963     xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
 964     xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
 965
 966     if (bMU && (mu_aver == -1))
 967     {
 968         xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
 969     }
 970     else
 971     {
 972         xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
 973     }
 974
 975     snew(fr, 1);
 976     clear_rvec(mu_t);
 977     teller = 0;
 978     /* Read the first frame from energy or traj file */
 979     if (bMU)
 980     {
 981         do
 982         {
 983             bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
 984             if (bCont)
 985             {
 986                 timecheck = check_times(t);
 987                 if (timecheck < 0)
 988                 {
 989                     teller++;
 990                 }
 991                 if ((teller % 10) == 0)
 992                 {
 993                     fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
 994                 }
 995             }
 996             else
 997             {
 998                 printf("End of %s reached\n", mufn);
 999                 break;
1000             }
1001         }
1002         while (bCont && (timecheck < 0));
1003     }
1004     else
1005     {
1006         natom  = read_first_x(oenv, &status, fn, &t, &x, box);
1007     }
1008
1009     /* Calculate spacing for dipole bin (simple histogram) */
1010     ndipbin = 1 + static_cast<int>(mu_max/0.01);
1011     snew(dipole_bin, ndipbin);
1012     mu_ave     = 0.0;
1013     for (m = 0; (m < DIM); m++)
1014     {
1015         M[m] = M2[m] = M4[m] = 0.0;
1016     }
1017
1018     if (bGkr)
1019     {
1020         /* Use 0.7 iso 0.5 to account for pressure scaling */
1021         /*  rcut   = 0.7*sqrt(max_cutoff2(box)); */
1022         rcut   = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1023
1024         gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1025     }
1026     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1027
1028     /* Start while loop over frames */
1029     t0     = t;
1030     teller = 0;
1031     do
1032     {
1033         if (bCorr && (teller >= nframes))
1034         {
1035             nframes += 1000;
1036             if (bTotal)
1037             {
1038                 srenew(muall[0], nframes*DIM);
1039             }
1040             else
1041             {
1042                 for (i = 0; (i < gnx_tot); i++)
1043                 {
1044                     srenew(muall[i], nframes*DIM);
1045                 }
1046             }
1047         }
1048         t1 = t;
1049
1050         muframelsq = gmx_stats_init();
1051
1052         /* Initialise */
1053         for (m = 0; (m < DIM); m++)
1054         {
1055             M_av2[m] = 0;
1056         }
1057
1058         if (bMU)
1059         {
1060             /* Copy rvec into double precision local variable */
1061             for (m = 0; (m < DIM); m++)
1062             {
1063                 M_av[m]  = mu_t[m];
1064             }
1065         }
1066         else
1067         {
1068             /* Initialise */
1069             for (m = 0; (m < DIM); m++)
1070             {
1071                 M_av[m] = 0;
1072             }
1073
1074             gmx_rmpbc(gpbc, natom, box, x);
1075
1076             /* Begin loop of all molecules in frame */
1077             for (n = 0; (n < ncos); n++)
1078             {
1079                 for (i = 0; (i < gnx[n]); i++)
1080                 {
1081                     int ind0, ind1;
1082
1083                     ind0  = mols->index[molindex[n][i]];
1084                     ind1  = mols->index[molindex[n][i]+1];
1085
1086                     mol_dip(ind0, ind1, x, atom, dipole[i]);
1087                     gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1088                     gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1089                     if (bSlab)
1090                     {
1091                         update_slab_dipoles(ind0, ind1, x,
1092                                             dipole[i], idim, nslices, slab_dipoles, box);
1093                     }
1094                     if (bQuad)
1095                     {
1096                         mol_quad(ind0, ind1, x, atom, quad);
1097                         for (m = 0; (m < DIM); m++)
1098                         {
1099                             gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1100                         }
1101                     }
1102                     if (bCorr && !bTotal)
1103                     {
1104                         tel3              = DIM*teller;
1105                         muall[i][tel3+XX] = dipole[i][XX];
1106                         muall[i][tel3+YY] = dipole[i][YY];
1107                         muall[i][tel3+ZZ] = dipole[i][ZZ];
1108                     }
1109                     mu_mol = 0.0;
1110                     for (m = 0; (m < DIM); m++)
1111                     {
1112                         M_av[m]  += dipole[i][m];               /* M per frame */
1113                         mu_mol   += dipole[i][m]*dipole[i][m];  /* calc. mu for distribution */
1114                     }
1115                     mu_mol = sqrt(mu_mol);
1116
1117                     mu_ave += mu_mol;                         /* calc. the average mu */
1118
1119                     /* Update the dipole distribution */
1120                     ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1121                     if (ibin < ndipbin)
1122                     {
1123                         dipole_bin[ibin]++;
1124                     }
1125
1126                     if (fndip3d)
1127                     {
1128                         rvec2sprvec(dipole[i], dipsp[i]);
1129
1130                         if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1131                         {
1132                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1133                             {
1134                                 ncolour = 1;
1135                             }
1136                             else
1137                             {
1138                                 ncolour = 2;
1139                             }
1140                         }
1141                         else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1142                         {
1143                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1144                             {
1145                                 ncolour = 3;
1146                             }
1147                             else
1148                             {
1149                                 ncolour = 4;
1150                             }
1151                         }
1152                         else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1153                         {
1154                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1155                             {
1156                                 ncolour = 5;
1157                             }
1158                             else
1159                             {
1160                                 ncolour = 6;
1161                             }
1162                         }
1163                         else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1164                         {
1165                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1166                             {
1167                                 ncolour = 7;
1168                             }
1169                             else
1170                             {
1171                                 ncolour = 8;
1172                             }
1173                         }
1174                         if (dip3d)
1175                         {
1176                             fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d  # %d %d\n",
1177                                     i+1,
1178                                     x[ind0][XX],
1179                                     x[ind0][YY],
1180                                     x[ind0][ZZ],
1181                                     x[ind0][XX]+dipole[i][XX]/25,
1182                                     x[ind0][YY]+dipole[i][YY]/25,
1183                                     x[ind0][ZZ]+dipole[i][ZZ]/25,
1184                                     ncolour, ind0, i);
1185                         }
1186                     }
1187                 } /* End loop of all molecules in frame */
1188
1189                 if (dip3d)
1190                 {
1191                     fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1192                     fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1193                     fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1194                     fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1195                     fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1196                     fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1197                 }
1198             }
1199         }
1200         /* Compute square of total dipole */
1201         for (m = 0; (m < DIM); m++)
1202         {
1203             M_av2[m] = M_av[m]*M_av[m];
1204         }
1205
1206         if (cosaver)
1207         {
1208             compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1209             rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1210                            sqr(dipaxis[YY]-0.5)+
1211                            sqr(dipaxis[ZZ]-0.5));
1212             if (bPairs)
1213             {
1214                 fprintf(caver, "%10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
1215                         t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1216             }
1217             else
1218             {
1219                 fprintf(caver, "%10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
1220                         t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1221             }
1222         }
1223
1224         if (bGkr)
1225         {
1226             do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1227                    atom, gkatom);
1228         }
1229
1230         if (bTotal)
1231         {
1232             tel3              = DIM*teller;
1233             muall[0][tel3+XX] = M_av[XX];
1234             muall[0][tel3+YY] = M_av[YY];
1235             muall[0][tel3+ZZ] = M_av[ZZ];
1236         }
1237
1238         /* Write to file the total dipole moment of the box, and its components
1239          * for this frame.
1240          */
1241         if ((skip == 0) || ((teller % skip) == 0))
1242         {
1243             fprintf(outmtot, "%10g  %12.8e %12.8e %12.8e %12.8e\n",
1244                     t, M_av[XX], M_av[YY], M_av[ZZ],
1245                     sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1246         }
1247
1248         for (m = 0; (m < DIM); m++)
1249         {
1250             M[m]  += M_av[m];
1251             M2[m] += M_av2[m];
1252             M4[m] += sqr(M_av2[m]);
1253         }
1254         /* Increment loop counter */
1255         teller++;
1256
1257         /* Calculate for output the running averages */
1258         invtel  = 1.0/teller;
1259         M2_ave  = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1260         M_ave2  = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1261         M_diff  = M2_ave - M_ave2;
1262
1263         /* Compute volume from box in traj, else we use the one from above */
1264         if (!bMU)
1265         {
1266             volume  = det(box);
1267         }
1268         vol_aver += volume;
1269
1270         epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1271
1272         /* Calculate running average for dipole */
1273         if (mu_ave != 0)
1274         {
1275             mu_aver = (mu_ave/gnx_tot)*invtel;
1276         }
1277
1278         if ((skip == 0) || ((teller % skip) == 0))
1279         {
1280             /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1281              * the two. Here M is sum mu_i. Further write the finite system
1282              * Kirkwood G factor and epsilon.
1283              */
1284             fprintf(outaver, "%10g  %10.3e %10.3e %10.3e %10.3e\n",
1285                     t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1286
1287             if (fnadip)
1288             {
1289                 real aver;
1290                 gmx_stats_get_average(muframelsq, &aver);
1291                 fprintf(adip, "%10g %f \n", t, aver);
1292             }
1293             /*if (dipole)
1294                printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1295              */
1296             if (!bMU || (mu_aver != -1))
1297             {
1298                 /* Finite system Kirkwood G-factor */
1299                 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1300                 /* Infinite system Kirkwood G-factor */
1301                 if (epsilonRF == 0.0)
1302                 {
1303                     g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1304                 }
1305                 else
1306                 {
1307                     g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1308                            Gk/(3*epsilon*(2*epsilonRF+1)));
1309                 }
1310
1311                 fprintf(outeps, "%10g  %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1312
1313             }
1314             else
1315             {
1316                 fprintf(outeps, "%10g  %12.8e\n", t, epsilon);
1317             }
1318         }
1319         gmx_stats_done(muframelsq);
1320
1321         if (bMU)
1322         {
1323             bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1324         }
1325         else
1326         {
1327             bCont = read_next_x(oenv, status, &t, x, box);
1328         }
1329         timecheck = check_times(t);
1330     }
1331     while (bCont && (timecheck == 0) );
1332
1333     gmx_rmpbc_done(gpbc);
1334
1335     if (!bMU)
1336     {
1337         close_trj(status);
1338     }
1339
1340     gmx_ffclose(outmtot);
1341     gmx_ffclose(outaver);
1342     gmx_ffclose(outeps);
1343
1344     if (fnadip)
1345     {
1346         gmx_ffclose(adip);
1347     }
1348
1349     if (cosaver)
1350     {
1351         gmx_ffclose(caver);
1352     }
1353
1354     if (dip3d)
1355     {
1356         fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1357         fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1358         fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1359         fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1360         fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1361         gmx_ffclose(dip3d);
1362     }
1363
1364     if (bSlab)
1365     {
1366         dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1367         sfree(slab_dipoles);
1368     }
1369
1370     vol_aver /= teller;
1371     printf("Average volume over run is %g\n", vol_aver);
1372     if (bGkr)
1373     {
1374         print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1375         print_cmap(cmap, gkrbin, nlevels);
1376     }
1377     /* Autocorrelation function */
1378     if (bCorr)
1379     {
1380         if (teller < 2)
1381         {
1382             printf("Not enough frames for autocorrelation\n");
1383         }
1384         else
1385         {
1386             dt = (t1 - t0)/(teller-1);
1387             printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1388
1389             mode = eacVector;
1390
1391             if (bTotal)
1392             {
1393                 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1394                             teller, 1, muall, dt, mode, TRUE);
1395             }
1396             else
1397             {
1398                 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1399                             teller, gnx_tot, muall, dt,
1400                             mode, strcmp(corrtype, "molsep"));
1401             }
1402         }
1403     }
1404     if (!bMU)
1405     {
1406         real aver, sigma, error;
1407
1408         gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1409         printf("\nDipole moment (Debye)\n");
1410         printf("---------------------\n");
1411         printf("Average  = %8.4f  Std. Dev. = %8.4f  Error = %8.4f\n",
1412                aver, sigma, error);
1413         if (bQuad)
1414         {
1415             rvec a, s, e;
1416             for (m = 0; (m < DIM); m++)
1417             {
1418                 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1419             }
1420
1421             printf("\nQuadrupole moment (Debye-Ang)\n");
1422             printf("-----------------------------\n");
1423             printf("Averages  = %8.4f  %8.4f  %8.4f\n", a[XX], a[YY], a[ZZ]);
1424             printf("Std. Dev. = %8.4f  %8.4f  %8.4f\n", s[XX], s[YY], s[ZZ]);
1425             printf("Error     = %8.4f  %8.4f  %8.4f\n", e[XX], e[YY], e[ZZ]);
1426         }
1427         printf("\n");
1428     }
1429     printf("The following averages for the complete trajectory have been calculated:\n\n");
1430     printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1431     printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1432     printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1433
1434     printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1435     printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1436     printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1437
1438     printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1439     printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1440
1441     printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1442
1443     if (!bMU || (mu_aver != -1))
1444     {
1445         printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1446         printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1447     }
1448     printf("Epsilon = %g\n", epsilon);
1449
1450     if (!bMU)
1451     {
1452         /* Write to file the dipole moment distibution during the simulation.
1453          */
1454         outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1455         for (i = 0; (i < ndipbin); i++)
1456         {
1457             fprintf(outdd, "%10g  %10f\n",
1458                     (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1459         }
1460         gmx_ffclose(outdd);
1461         sfree(dipole_bin);
1462     }
1463     if (bGkr)
1464     {
1465         done_gkrbin(&gkrbin);
1466     }
1467 }
1468
1469 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1470 {
1471     int nmol, i, j, m;
1472
1473     nmol = 0;
1474     i    = 0;
1475     while (i < *n)
1476     {
1477         m = 0;
1478         while (m < mols->nr && index[i] != mols->index[m])
1479         {
1480             m++;
1481         }
1482         if (m == mols->nr)
1483         {
1484             gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1485         }
1486         for (j = mols->index[m]; j < mols->index[m+1]; j++)
1487         {
1488             if (i >= *n || index[i] != j)
1489             {
1490                 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1491             }
1492             i++;
1493         }
1494         /* Modify the index in place */
1495         index[nmol++] = m;
1496     }
1497     printf("There are %d molecules in the selection\n", nmol);
1498
1499     *n = nmol;
1500 }
1501 int gmx_dipoles(int argc, char *argv[])
1502 {
1503     const char    *desc[] = {
1504         "[THISMODULE] computes the total dipole plus fluctuations of a simulation",
1505         "system. From this you can compute e.g. the dielectric constant for",
1506         "low-dielectric media.",
1507         "For molecules with a net charge, the net charge is subtracted at",
1508         "center of mass of the molecule.[PAR]",
1509         "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1510         "components as well as the norm of the vector.",
1511         "The file [TT]aver.xvg[tt] contains [CHEVRON][MAG][GRK]mu[grk][mag]^2[chevron] and [MAG][CHEVRON][GRK]mu[grk][chevron][mag]^2 during the",
1512         "simulation.",
1513         "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1514         "the simulation",
1515         "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1516         "Furthermore, the dipole autocorrelation function will be computed when",
1517         "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1518         "option.",
1519         "The correlation functions can be averaged over all molecules",
1520         "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1521         "or it can be computed over the total dipole moment of the simulation box",
1522         "([TT]total[tt]).[PAR]",
1523         "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1524         "G-factor, as well as the average cosine of the angle between the dipoles",
1525         "as a function of the distance. The plot also includes gOO and hOO",
1526         "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1527         "we also include the energy per scale computed by taking the inner product of",
1528         "the dipoles divided by the distance to the third power.[PAR]",
1529         "[PAR]",
1530         "EXAMPLES[PAR]",
1531         "[TT]gmx dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1532         "This will calculate the autocorrelation function of the molecular",
1533         "dipoles using a first order Legendre polynomial of the angle of the",
1534         "dipole vector and itself a time t later. For this calculation 1001",
1535         "frames will be used. Further, the dielectric constant will be calculated",
1536         "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1537         "an average dipole moment of the molecule of 2.273 (SPC). For the",
1538         "distribution function a maximum of 5.0 will be used."
1539     };
1540     real           mu_max     = 5, mu_aver = -1, rcmax = 0;
1541     real           epsilonRF  = 0.0, temp = 300;
1542     gmx_bool       bPairs     = TRUE, bPhi = FALSE, bQuad = FALSE;
1543     const char    *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1544     const char    *axtitle    = "Z";
1545     int            nslices    = 10; /* nr of slices defined       */
1546     int            skip       = 0, nFA = 0, nFB = 0, ncos = 1;
1547     int            nlevels    = 20, ndegrees = 90;
1548     output_env_t   oenv;
1549     t_pargs        pa[] = {
1550         { "-mu",       FALSE, etREAL, {&mu_aver},
1551           "dipole of a single molecule (in Debye)" },
1552         { "-mumax",    FALSE, etREAL, {&mu_max},
1553           "max dipole in Debye (for histogram)" },
1554         { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1555           "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1556         { "-skip",     FALSE, etINT, {&skip},
1557           "Skip steps in the output (but not in the computations)" },
1558         { "-temp",     FALSE, etREAL, {&temp},
1559           "Average temperature of the simulation (needed for dielectric constant calculation)" },
1560         { "-corr",     FALSE, etENUM, {corrtype},
1561           "Correlation function to calculate" },
1562         { "-pairs",    FALSE, etBOOL, {&bPairs},
1563           "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1564         { "-quad",     FALSE, etBOOL, {&bQuad},
1565           "Take quadrupole into account"},
1566         { "-ncos",     FALSE, etINT, {&ncos},
1567           "Must be 1 or 2. Determines whether the [CHEVRON][COS][GRK]theta[grk][cos][chevron] is computed between all molecules in one group, or between molecules in two different groups. This turns on the [TT]-g[tt] flag." },
1568         { "-axis",     FALSE, etSTR, {&axtitle},
1569           "Take the normal on the computational box in direction X, Y or Z." },
1570         { "-sl",       FALSE, etINT, {&nslices},
1571           "Divide the box into this number of slices." },
1572         { "-gkratom",  FALSE, etINT, {&nFA},
1573           "Use the n-th atom of a molecule (starting from 1) to calculate the distance between molecules rather than the center of charge (when 0) in the calculation of distance dependent Kirkwood factors" },
1574         { "-gkratom2", FALSE, etINT, {&nFB},
1575           "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1576         { "-rcmax",    FALSE, etREAL, {&rcmax},
1577           "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1578         { "-phi",      FALSE, etBOOL, {&bPhi},
1579           "Plot the 'torsion angle' defined as the rotation of the two dipole vectors around the distance vector between the two molecules in the [TT].xpm[tt] file from the [TT]-cmap[tt] option. By default the cosine of the angle between the dipoles is plotted." },
1580         { "-nlevels",  FALSE, etINT, {&nlevels},
1581           "Number of colors in the cmap output" },
1582         { "-ndegrees", FALSE, etINT, {&ndegrees},
1583           "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1584     };
1585     int           *gnx;
1586     int            nFF[2];
1587     atom_id      **grpindex;
1588     char         **grpname = NULL;
1589     gmx_bool       bGkr, bMU, bSlab;
1590     t_filenm       fnm[] = {
1591         { efEDR, "-en", NULL,         ffOPTRD },
1592         { efTRX, "-f", NULL,           ffREAD },
1593         { efTPX, NULL, NULL,           ffREAD },
1594         { efNDX, NULL, NULL,           ffOPTRD },
1595         { efXVG, "-o",   "Mtot",       ffWRITE },
1596         { efXVG, "-eps", "epsilon",    ffWRITE },
1597         { efXVG, "-a",   "aver",       ffWRITE },
1598         { efXVG, "-d",   "dipdist",    ffWRITE },
1599         { efXVG, "-c",   "dipcorr",    ffOPTWR },
1600         { efXVG, "-g",   "gkr",        ffOPTWR },
1601         { efXVG, "-adip", "adip",       ffOPTWR },
1602         { efXVG, "-dip3d", "dip3d",    ffOPTWR },
1603         { efXVG, "-cos", "cosaver",    ffOPTWR },
1604         { efXPM, "-cmap", "cmap",       ffOPTWR },
1605         { efXVG, "-slab", "slab",       ffOPTWR }
1606     };
1607 #define NFILE asize(fnm)
1608     int            npargs;
1609     t_pargs       *ppa;
1610     t_topology    *top;
1611     int            ePBC;
1612     int            k, natoms;
1613     matrix         box;
1614
1615     npargs = asize(pa);
1616     ppa    = add_acf_pargs(&npargs, pa);
1617     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
1618                            NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
1619     {
1620         return 0;
1621     }
1622
1623     printf("Using %g as mu_max and %g as the dipole moment.\n",
1624            mu_max, mu_aver);
1625     if (epsilonRF == 0.0)
1626     {
1627         printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1628     }
1629
1630     bMU   = opt2bSet("-en", NFILE, fnm);
1631     if (bMU)
1632     {
1633         gmx_fatal(FARGS, "Due to new ways of treating molecules in GROMACS the total dipole in the energy file may be incorrect, because molecules can be split over periodic boundary conditions before computing the dipole. Please use your trajectory file.");
1634     }
1635     bGkr  = opt2bSet("-g", NFILE, fnm);
1636     if (opt2parg_bSet("-ncos", asize(pa), pa))
1637     {
1638         if ((ncos != 1) && (ncos != 2))
1639         {
1640             gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1641         }
1642         bGkr = TRUE;
1643     }
1644     bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1645              opt2parg_bSet("-axis", asize(pa), pa));
1646     if (bMU)
1647     {
1648         if (bQuad)
1649         {
1650             printf("WARNING: Can not determine quadrupoles from energy file\n");
1651             bQuad = FALSE;
1652         }
1653         if (bGkr)
1654         {
1655             printf("WARNING: Can not determine Gk(r) from energy file\n");
1656             bGkr  = FALSE;
1657             ncos  = 1;
1658         }
1659         if (mu_aver == -1)
1660         {
1661             printf("WARNING: Can not calculate Gk and gk, since you did\n"
1662                    "         not enter a valid dipole for the molecules\n");
1663         }
1664     }
1665
1666     snew(top, 1);
1667     ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1668                         &natoms, NULL, NULL, NULL, top);
1669
1670     snew(gnx, ncos);
1671     snew(grpname, ncos);
1672     snew(grpindex, ncos);
1673     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1674               ncos, gnx, grpindex, grpname);
1675     for (k = 0; (k < ncos); k++)
1676     {
1677         dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1678         neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1679     }
1680     nFF[0] = nFA;
1681     nFF[1] = nFB;
1682     do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1683            opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1684            opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1685            opt2fn_null("-cos", NFILE, fnm),
1686            opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1687            bPairs, corrtype[0],
1688            opt2fn("-c", NFILE, fnm),
1689            bGkr,    opt2fn("-g", NFILE, fnm),
1690            bPhi,    &nlevels,  ndegrees,
1691            ncos,
1692            opt2fn("-cmap", NFILE, fnm), rcmax,
1693            bQuad, bMU,     opt2fn("-en", NFILE, fnm),
1694            gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1695            bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1696
1697     do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1698     do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1699     do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1700     do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1701     do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");
1702
1703     return 0;
1704 }