src/gromacs/gmxana/gmx_dipoles.cpp

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