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