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