src/gromacs/gmxana/gmx_dipoles.cpp

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