src/tools/gmx_dipoles.c

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