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