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