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