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