src/gromacs/gmxana/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 "futil.h"
  49 #include "xvgr.h"
  50 #include "txtdump.h"
  51 #include "gmx_statistics.h"
  52 #include "gstat.h"
  53 #include "index.h"
  54 #include "random.h"
  55 #include "names.h"
  56 #include "physics.h"
  57 #include "calcmu.h"
  58 #include "enxio.h"
  59 #include "nrjac.h"
  60 #include "matio.h"
  61 #include "gmx_ana.h"
  62
  63
  64 #define e2d(x) ENM2DEBYE*(x)
  65 #define EANG2CM  E_CHARGE*1.0e-10       /* e Angstrom to Coulomb meter */
  66 #define CM2D  SPEED_OF_LIGHT*1.0e+24    /* Coulomb meter to Debye */
  67
  68 typedef struct {
  69     int      nelem;
  70     real     spacing, radius;
  71     real    *elem;
  72     int     *count;
  73     gmx_bool bPhi;
  74     int      nx, ny;
  75     real   **cmap;
  76 } t_gkrbin;
  77
  78 static t_gkrbin *mk_gkrbin(real radius, real rcmax, gmx_bool bPhi, int ndegrees)
  79 {
  80     t_gkrbin *gb;
  81     char     *ptr;
  82     int       i;
  83
  84     snew(gb, 1);
  85
  86     if ((ptr = getenv("GKRWIDTH")) != NULL)
  87     {
  88         double bw;
  89
  90         sscanf(ptr, "%lf", &bw);
  91         gb->spacing = bw;
  92     }
  93     else
  94     {
  95         gb->spacing = 0.01; /* nm */
  96     }
  97     gb->nelem   = 1 + radius/gb->spacing;
  98     if (rcmax == 0)
  99     {
 100         gb->nx = gb->nelem;
 101     }
 102     else
 103     {
 104         gb->nx = 1 + 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 = 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 = (M_PI+phi)*gb->ny/(2*M_PI);
 145         }
 146         else
 147         {
 148             cy = (alpha*gb->ny)/M_PI; /*((1+cosa)*0.5*(gb->ny));*/
 149         }
 150         cy = min(gb->ny-1, 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, index, grp0, grp1;
 175     real         qtot, q, r2, 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 || (cosa != 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     double hi, vol;
 269
 270     hi = 0;
 271     for (i = 0; (i < gb->nx); i++)
 272     {
 273         vol = 4*M_PI*sqr(gb->spacing*i)*gb->spacing;
 274         for (j = 0; (j < gb->ny); j++)
 275         {
 276             gb->cmap[i][j] /= vol;
 277             hi              = max(hi, gb->cmap[i][j]);
 278         }
 279     }
 280     if (hi <= 0)
 281     {
 282         gmx_fatal(FARGS, "No data in the cmap");
 283     }
 284     return hi;
 285 }
 286
 287 static void print_cmap(const char *cmap, t_gkrbin *gb, int *nlevels)
 288 {
 289     FILE   *out;
 290     int     i, j;
 291     real    hi;
 292
 293     real   *xaxis, *yaxis;
 294     t_rgb   rlo = { 1, 1, 1 };
 295     t_rgb   rhi = { 0, 0, 0 };
 296
 297     hi = normalize_cmap(gb);
 298     snew(xaxis, gb->nx+1);
 299     for (i = 0; (i < gb->nx+1); i++)
 300     {
 301         xaxis[i] = i*gb->spacing;
 302     }
 303     snew(yaxis, gb->ny);
 304     for (j = 0; (j < gb->ny); j++)
 305     {
 306         if (gb->bPhi)
 307         {
 308             yaxis[j] = (360.0*j)/(gb->ny-1.0)-180;
 309         }
 310         else
 311         {
 312             yaxis[j] = (180.0*j)/(gb->ny-1.0);
 313         }
 314         /*2.0*j/(gb->ny-1.0)-1.0;*/
 315     }
 316     out = ffopen(cmap, "w");
 317     write_xpm(out, 0,
 318               "Dipole Orientation Distribution", "Fraction", "r (nm)",
 319               gb->bPhi ? "Phi" : "Alpha",
 320               gb->nx, gb->ny, xaxis, yaxis,
 321               gb->cmap, 0, hi, rlo, rhi, nlevels);
 322     ffclose(out);
 323     sfree(xaxis);
 324     sfree(yaxis);
 325 }
 326
 327 static void print_gkrbin(const char *fn, t_gkrbin *gb,
 328                          int ngrp, int nframes, real volume,
 329                          const output_env_t oenv)
 330 {
 331     /* We compute Gk(r), gOO and hOO according to
 332      * Nymand & Linse, JCP 112 (2000) pp 6386-6395.
 333      * In this implementation the angle between dipoles is stored
 334      * rather than their inner product. This allows to take polarizible
 335      * models into account. The RDF is calculated as well, almost for free!
 336      */
 337     FILE       *fp;
 338     const char *leg[] = { "G\\sk\\N(r)", "< cos >", "h\\sOO\\N", "g\\sOO\\N", "Energy" };
 339     int         i, j, n, last;
 340     real        x0, x1, ggg, Gkr, vol_s, rho, gOO, hOO, cosav, ener;
 341     double      fac;
 342
 343     fp = xvgropen(fn, "Distance dependent Gk", "r (nm)", "G\\sk\\N(r)", oenv);
 344     xvgr_legend(fp, asize(leg), leg, oenv);
 345
 346     Gkr = 1;  /* Self-dipole inproduct = 1 */
 347     rho = ngrp/volume;
 348
 349     if (debug)
 350     {
 351         fprintf(debug, "Number density is %g molecules / nm^3\n", rho);
 352         fprintf(debug, "ngrp = %d, nframes = %d\n", ngrp, nframes);
 353     }
 354
 355     last = gb->nelem-1;
 356     while (last > 1 && gb->elem[last-1] == 0)
 357     {
 358         last--;
 359     }
 360
 361     /* Divide by dipole squared, by number of frames, by number of origins.
 362      * Multiply by 2 because we only take half the matrix of interactions
 363      * into account.
 364      */
 365     fac  = 2.0/((double) ngrp * (double) nframes);
 366
 367     x0 = 0;
 368     for (i = 0; i < last; i++)
 369     {
 370         /* Centre of the coordinate in the spherical layer */
 371         x1    = x0+gb->spacing;
 372
 373         /* Volume of the layer */
 374         vol_s = (4.0/3.0)*M_PI*(x1*x1*x1-x0*x0*x0);
 375
 376         /* gOO */
 377         gOO   = gb->count[i]*fac/(rho*vol_s);
 378
 379         /* Dipole correlation hOO, normalized by the relative number density, like
 380          * in a Radial distribution function.
 381          */
 382         ggg  = gb->elem[i]*fac;
 383         hOO  = 3.0*ggg/(rho*vol_s);
 384         Gkr += ggg;
 385         if (gb->count[i])
 386         {
 387             cosav = gb->elem[i]/gb->count[i];
 388         }
 389         else
 390         {
 391             cosav = 0;
 392         }
 393         ener = -0.5*cosav*ONE_4PI_EPS0/(x1*x1*x1);
 394
 395         fprintf(fp, "%10.5e %12.5e %12.5e %12.5e %12.5e  %12.5e\n",
 396                 x1, Gkr, cosav, hOO, gOO, ener);
 397
 398         /* Swap x0 and x1 */
 399         x0 = x1;
 400     }
 401     ffclose(fp);
 402 }
 403
 404 gmx_bool read_mu_from_enx(ener_file_t fmu, int Vol, ivec iMu, rvec mu, real *vol,
 405                           real *t, int nre, t_enxframe *fr)
 406 {
 407     int          i;
 408     gmx_bool     bCont;
 409     char         buf[22];
 410
 411     bCont = do_enx(fmu, fr);
 412     if (fr->nre != nre)
 413     {
 414         fprintf(stderr, "Something strange: expected %d entries in energy file at step %s\n(time %g) but found %d entries\n",
 415                 nre, gmx_step_str(fr->step, buf), fr->t, fr->nre);
 416     }
 417
 418     if (bCont)
 419     {
 420         if (Vol != -1)          /* we've got Volume in the energy file */
 421         {
 422             *vol = fr->ener[Vol].e;
 423         }
 424         for (i = 0; i < DIM; i++)
 425         {
 426             mu[i] = fr->ener[iMu[i]].e;
 427         }
 428         *t = fr->t;
 429     }
 430
 431     return bCont;
 432 }
 433
 434 static void neutralize_mols(int n, int *index, t_block *mols, t_atom *atom)
 435 {
 436     double mtot, qtot;
 437     int    ncharged, m, a0, a1, a;
 438
 439     ncharged = 0;
 440     for (m = 0; m < n; m++)
 441     {
 442         a0   = mols->index[index[m]];
 443         a1   = mols->index[index[m]+1];
 444         mtot = 0;
 445         qtot = 0;
 446         for (a = a0; a < a1; a++)
 447         {
 448             mtot += atom[a].m;
 449             qtot += atom[a].q;
 450         }
 451         /* This check is only for the count print */
 452         if (fabs(qtot) > 0.01)
 453         {
 454             ncharged++;
 455         }
 456         if (mtot > 0)
 457         {
 458             /* Remove the net charge at the center of mass */
 459             for (a = a0; a < a1; a++)
 460             {
 461                 atom[a].q -= qtot*atom[a].m/mtot;
 462             }
 463         }
 464     }
 465
 466     if (ncharged)
 467     {
 468         printf("There are %d charged molecules in the selection,\n"
 469                "will subtract their charge at their center of mass\n", ncharged);
 470     }
 471 }
 472
 473 static void mol_dip(int k0, int k1, rvec x[], t_atom atom[], rvec mu)
 474 {
 475     int  k, m;
 476     real q;
 477
 478     clear_rvec(mu);
 479     for (k = k0; (k < k1); k++)
 480     {
 481         q  = e2d(atom[k].q);
 482         for (m = 0; (m < DIM); m++)
 483         {
 484             mu[m] += q*x[k][m];
 485         }
 486     }
 487 }
 488
 489 static void mol_quad(int k0, int k1, rvec x[], t_atom atom[], rvec quad)
 490 {
 491     int      i, k, m, n, niter;
 492     real     q, r2, mass, masstot;
 493     rvec     com;        /* center of mass */
 494     rvec     r;          /* distance of atoms to center of mass */
 495     real     rcom_m, rcom_n;
 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, dc1, d, n5, 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, j, k, n, m, natom = 0, nmol, gnx_tot, teller, tel3;
 783     t_trxstatus  *status;
 784     int          *dipole_bin, ndipbin, ibin, iVol, step, idim = -1;
 785     unsigned long mode;
 786     char          buf[STRLEN];
 787     real          rcut = 0, t, t0, t1, dt, lambda, dd, rms_cos;
 788     rvec          dipaxis;
 789     matrix        box;
 790     gmx_bool      bCorr, bTotal, bCont;
 791     double        M_diff = 0, epsilon, invtel, vol_aver;
 792     double        mu_ave, mu_mol, M2_ave = 0, M_ave2 = 0, M_av[DIM], M_av2[DIM];
 793     double        M[3], M2[3], M4[3], Gk = 0, g_k = 0;
 794     gmx_stats_t   Mx, My, Mz, Msq, Vol, *Qlsq, mulsq, muframelsq = NULL;
 795     ivec          iMu;
 796     real        **muall        = NULL;
 797     rvec         *slab_dipoles = NULL;
 798     t_atom       *atom         = NULL;
 799     t_block      *mols         = NULL;
 800     gmx_rmpbc_t   gpbc         = NULL;
 801
 802     gnx_tot = gnx[0];
 803     if (ncos > 1)
 804     {
 805         gnx_tot += gnx[1];
 806     }
 807
 808     vol_aver = 0.0;
 809
 810     iVol = -1;
 811
 812     /* initialize to a negative value so we can see that it wasn't picked up */
 813     iMu[XX] = iMu[YY] = iMu[ZZ] = -1;
 814     if (bMU)
 815     {
 816         fmu = open_enx(mufn, "r");
 817         do_enxnms(fmu, &nre, &enm);
 818
 819         /* Determine the indexes of the energy grps we need */
 820         for (i = 0; (i < nre); i++)
 821         {
 822             if (strstr(enm[i].name, "Volume"))
 823             {
 824                 iVol = i;
 825             }
 826             else if (strstr(enm[i].name, "Mu-X"))
 827             {
 828                 iMu[XX] = i;
 829             }
 830             else if (strstr(enm[i].name, "Mu-Y"))
 831             {
 832                 iMu[YY] = i;
 833             }
 834             else if (strstr(enm[i].name, "Mu-Z"))
 835             {
 836                 iMu[ZZ] = i;
 837             }
 838         }
 839         if (iMu[XX] < 0 || iMu[YY] < 0 || iMu[ZZ] < 0)
 840         {
 841             gmx_fatal(FARGS, "No index for Mu-X, Mu-Y or Mu-Z energy group.");
 842         }
 843     }
 844     else
 845     {
 846         atom = top->atoms.atom;
 847         mols = &(top->mols);
 848     }
 849     if ((iVol == -1) && bMU)
 850     {
 851         printf("Using Volume from topology: %g nm^3\n", volume);
 852     }
 853
 854     /* Correlation stuff */
 855     bCorr  = (corrtype[0] != 'n');
 856     bTotal = (corrtype[0] == 't');
 857     if (bCorr)
 858     {
 859         if (bTotal)
 860         {
 861             snew(muall, 1);
 862             snew(muall[0], nframes*DIM);
 863         }
 864         else
 865         {
 866             snew(muall, gnx[0]);
 867             for (i = 0; (i < gnx[0]); i++)
 868             {
 869                 snew(muall[i], nframes*DIM);
 870             }
 871         }
 872     }
 873
 874     /* Allocate array which contains for every molecule in a frame the
 875      * dipole moment.
 876      */
 877     if (!bMU)
 878     {
 879         snew(dipole, gnx_tot);
 880     }
 881
 882     /* Statistics */
 883     snew(Qlsq, DIM);
 884     for (i = 0; (i < DIM); i++)
 885     {
 886         Qlsq[i] = gmx_stats_init();
 887     }
 888     mulsq = gmx_stats_init();
 889
 890     /* Open all the files */
 891     outmtot = xvgropen(out_mtot,
 892                        "Total dipole moment of the simulation box vs. time",
 893                        "Time (ps)", "Total Dipole Moment (Debye)", oenv);
 894     outeps  = xvgropen(out_eps, "Epsilon and Kirkwood factors",
 895                        "Time (ps)", "", oenv);
 896     outaver = xvgropen(out_aver, "Total dipole moment",
 897                        "Time (ps)", "D", oenv);
 898     if (bSlab)
 899     {
 900         idim = axtitle[0] - 'X';
 901         if ((idim < 0) || (idim >= DIM))
 902         {
 903             idim = axtitle[0] - 'x';
 904         }
 905         if ((idim < 0) || (idim >= DIM))
 906         {
 907             bSlab = FALSE;
 908         }
 909         if (nslices < 2)
 910         {
 911             bSlab = FALSE;
 912         }
 913         fprintf(stderr, "axtitle = %s, nslices = %d, idim = %d\n",
 914                 axtitle, nslices, idim);
 915         if (bSlab)
 916         {
 917             snew(slab_dipoles, nslices);
 918             fprintf(stderr, "Doing slab analysis\n");
 919         }
 920     }
 921
 922     if (fnadip)
 923     {
 924         adip = xvgropen(fnadip, "Average molecular dipole", "Dipole (D)", "", oenv);
 925         xvgr_legend(adip, NLEGADIP, leg_adip, oenv);
 926
 927     }
 928     if (cosaver)
 929     {
 930         caver = xvgropen(cosaver, bPairs ? "Average pair orientation" :
 931                          "Average absolute dipole orientation", "Time (ps)", "", oenv);
 932         xvgr_legend(caver, NLEGCOSAVER, bPairs ? leg_cosaver : &(leg_cosaver[1]),
 933                     oenv);
 934     }
 935
 936     if (fndip3d)
 937     {
 938         snew(dipsp, gnx_tot);
 939
 940         /* we need a dummy file for gnuplot */
 941         dip3d = (FILE *)ffopen("dummy.dat", "w");
 942         fprintf(dip3d, "%f %f %f", 0.0, 0.0, 0.0);
 943         ffclose(dip3d);
 944
 945         dip3d = (FILE *)ffopen(fndip3d, "w");
 946         fprintf(dip3d, "# This file was created by %s\n", Program());
 947         fprintf(dip3d, "# which is part of G R O M A C S:\n");
 948         fprintf(dip3d, "#\n");
 949     }
 950
 951     /* Write legends to all the files */
 952     xvgr_legend(outmtot, NLEGMTOT, leg_mtot, oenv);
 953     xvgr_legend(outaver, NLEGAVER, leg_aver, oenv);
 954
 955     if (bMU && (mu_aver == -1))
 956     {
 957         xvgr_legend(outeps, NLEGEPS-2, leg_eps, oenv);
 958     }
 959     else
 960     {
 961         xvgr_legend(outeps, NLEGEPS, leg_eps, oenv);
 962     }
 963
 964     snew(fr, 1);
 965     clear_rvec(mu_t);
 966     teller = 0;
 967     /* Read the first frame from energy or traj file */
 968     if (bMU)
 969     {
 970         do
 971         {
 972             bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
 973             if (bCont)
 974             {
 975                 timecheck = check_times(t);
 976                 if (timecheck < 0)
 977                 {
 978                     teller++;
 979                 }
 980                 if ((teller % 10) == 0)
 981                 {
 982                     fprintf(stderr, "\r Skipping Frame %6d, time: %8.3f", teller, t);
 983                 }
 984             }
 985             else
 986             {
 987                 printf("End of %s reached\n", mufn);
 988                 break;
 989             }
 990         }
 991         while (bCont && (timecheck < 0));
 992     }
 993     else
 994     {
 995         natom  = read_first_x(oenv, &status, fn, &t, &x, box);
 996     }
 997
 998     /* Calculate spacing for dipole bin (simple histogram) */
 999     ndipbin = 1+(mu_max/0.01);
1000     snew(dipole_bin, ndipbin);
1001     epsilon    = 0.0;
1002     mu_ave     = 0.0;
1003     for (m = 0; (m < DIM); m++)
1004     {
1005         M[m] = M2[m] = M4[m] = 0.0;
1006     }
1007
1008     if (bGkr)
1009     {
1010         /* Use 0.7 iso 0.5 to account for pressure scaling */
1011         /*  rcut   = 0.7*sqrt(max_cutoff2(box)); */
1012         rcut   = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
1013
1014         gkrbin = mk_gkrbin(rcut, rcmax, bPhi, ndegrees);
1015     }
1016     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natom);
1017
1018     /* Start while loop over frames */
1019     t1     = t0 = t;
1020     teller = 0;
1021     do
1022     {
1023         if (bCorr && (teller >= nframes))
1024         {
1025             nframes += 1000;
1026             if (bTotal)
1027             {
1028                 srenew(muall[0], nframes*DIM);
1029             }
1030             else
1031             {
1032                 for (i = 0; (i < gnx_tot); i++)
1033                 {
1034                     srenew(muall[i], nframes*DIM);
1035                 }
1036             }
1037         }
1038         t1 = t;
1039
1040         muframelsq = gmx_stats_init();
1041
1042         /* Initialise */
1043         for (m = 0; (m < DIM); m++)
1044         {
1045             M_av2[m] = 0;
1046         }
1047
1048         if (bMU)
1049         {
1050             /* Copy rvec into double precision local variable */
1051             for (m = 0; (m < DIM); m++)
1052             {
1053                 M_av[m]  = mu_t[m];
1054             }
1055         }
1056         else
1057         {
1058             /* Initialise */
1059             for (m = 0; (m < DIM); m++)
1060             {
1061                 M_av[m] = 0;
1062             }
1063
1064             gmx_rmpbc(gpbc, natom, box, x);
1065
1066             /* Begin loop of all molecules in frame */
1067             for (n = 0; (n < ncos); n++)
1068             {
1069                 for (i = 0; (i < gnx[n]); i++)
1070                 {
1071                     int gi, ind0, ind1;
1072
1073                     ind0  = mols->index[molindex[n][i]];
1074                     ind1  = mols->index[molindex[n][i]+1];
1075
1076                     mol_dip(ind0, ind1, x, atom, dipole[i]);
1077                     gmx_stats_add_point(mulsq, 0, norm(dipole[i]), 0, 0);
1078                     gmx_stats_add_point(muframelsq, 0, norm(dipole[i]), 0, 0);
1079                     if (bSlab)
1080                     {
1081                         update_slab_dipoles(ind0, ind1, x,
1082                                             dipole[i], idim, nslices, slab_dipoles, box);
1083                     }
1084                     if (bQuad)
1085                     {
1086                         mol_quad(ind0, ind1, x, atom, quad);
1087                         for (m = 0; (m < DIM); m++)
1088                         {
1089                             gmx_stats_add_point(Qlsq[m], 0, quad[m], 0, 0);
1090                         }
1091                     }
1092                     if (bCorr && !bTotal)
1093                     {
1094                         tel3              = DIM*teller;
1095                         muall[i][tel3+XX] = dipole[i][XX];
1096                         muall[i][tel3+YY] = dipole[i][YY];
1097                         muall[i][tel3+ZZ] = dipole[i][ZZ];
1098                     }
1099                     mu_mol = 0.0;
1100                     for (m = 0; (m < DIM); m++)
1101                     {
1102                         M_av[m]  += dipole[i][m];               /* M per frame */
1103                         mu_mol   += dipole[i][m]*dipole[i][m];  /* calc. mu for distribution */
1104                     }
1105                     mu_mol = sqrt(mu_mol);
1106
1107                     mu_ave += mu_mol;                         /* calc. the average mu */
1108
1109                     /* Update the dipole distribution */
1110                     ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
1111                     if (ibin < ndipbin)
1112                     {
1113                         dipole_bin[ibin]++;
1114                     }
1115
1116                     if (fndip3d)
1117                     {
1118                         rvec2sprvec(dipole[i], dipsp[i]);
1119
1120                         if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI)
1121                         {
1122                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1123                             {
1124                                 ncolour = 1;
1125                             }
1126                             else
1127                             {
1128                                 ncolour = 2;
1129                             }
1130                         }
1131                         else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI)
1132                         {
1133                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1134                             {
1135                                 ncolour = 3;
1136                             }
1137                             else
1138                             {
1139                                 ncolour = 4;
1140                             }
1141                         }
1142                         else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI)
1143                         {
1144                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1145                             {
1146                                 ncolour = 5;
1147                             }
1148                             else
1149                             {
1150                                 ncolour = 6;
1151                             }
1152                         }
1153                         else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI)
1154                         {
1155                             if (dipsp[i][ZZ] < 0.5 * M_PI)
1156                             {
1157                                 ncolour = 7;
1158                             }
1159                             else
1160                             {
1161                                 ncolour = 8;
1162                             }
1163                         }
1164                         if (dip3d)
1165                         {
1166                             fprintf(dip3d, "set arrow %d from %f, %f, %f to %f, %f, %f lt %d  # %d %d\n",
1167                                     i+1,
1168                                     x[ind0][XX],
1169                                     x[ind0][YY],
1170                                     x[ind0][ZZ],
1171                                     x[ind0][XX]+dipole[i][XX]/25,
1172                                     x[ind0][YY]+dipole[i][YY]/25,
1173                                     x[ind0][ZZ]+dipole[i][ZZ]/25,
1174                                     ncolour, ind0, i);
1175                         }
1176                     }
1177                 } /* End loop of all molecules in frame */
1178
1179                 if (dip3d)
1180                 {
1181                     fprintf(dip3d, "set title \"t = %4.3f\"\n", t);
1182                     fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1183                     fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1184                     fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1185                     fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1186                     fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1187                 }
1188             }
1189         }
1190         /* Compute square of total dipole */
1191         for (m = 0; (m < DIM); m++)
1192         {
1193             M_av2[m] = M_av[m]*M_av[m];
1194         }
1195
1196         if (cosaver)
1197         {
1198             compute_avercos(gnx_tot, dipole, &dd, dipaxis, bPairs);
1199             rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
1200                            sqr(dipaxis[YY]-0.5)+
1201                            sqr(dipaxis[ZZ]-0.5));
1202             if (bPairs)
1203             {
1204                 fprintf(caver, "%10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
1205                         t, dd, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1206             }
1207             else
1208             {
1209                 fprintf(caver, "%10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
1210                         t, rms_cos, dipaxis[XX], dipaxis[YY], dipaxis[ZZ]);
1211             }
1212         }
1213
1214         if (bGkr)
1215         {
1216             do_gkr(gkrbin, ncos, gnx, molindex, mols->index, x, dipole, ePBC, box,
1217                    atom, gkatom);
1218         }
1219
1220         if (bTotal)
1221         {
1222             tel3              = DIM*teller;
1223             muall[0][tel3+XX] = M_av[XX];
1224             muall[0][tel3+YY] = M_av[YY];
1225             muall[0][tel3+ZZ] = M_av[ZZ];
1226         }
1227
1228         /* Write to file the total dipole moment of the box, and its components
1229          * for this frame.
1230          */
1231         if ((skip == 0) || ((teller % skip) == 0))
1232         {
1233             fprintf(outmtot, "%10g  %12.8e %12.8e %12.8e %12.8e\n",
1234                     t, M_av[XX], M_av[YY], M_av[ZZ],
1235                     sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1236         }
1237
1238         for (m = 0; (m < DIM); m++)
1239         {
1240             M[m]  += M_av[m];
1241             M2[m] += M_av2[m];
1242             M4[m] += sqr(M_av2[m]);
1243         }
1244         /* Increment loop counter */
1245         teller++;
1246
1247         /* Calculate for output the running averages */
1248         invtel  = 1.0/teller;
1249         M2_ave  = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1250         M_ave2  = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1251         M_diff  = M2_ave - M_ave2;
1252
1253         /* Compute volume from box in traj, else we use the one from above */
1254         if (!bMU)
1255         {
1256             volume  = det(box);
1257         }
1258         vol_aver += volume;
1259
1260         epsilon = calc_eps(M_diff, (vol_aver/teller), epsilonRF, temp);
1261
1262         /* Calculate running average for dipole */
1263         if (mu_ave != 0)
1264         {
1265             mu_aver = (mu_ave/gnx_tot)*invtel;
1266         }
1267
1268         if ((skip == 0) || ((teller % skip) == 0))
1269         {
1270             /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1271              * the two. Here M is sum mu_i. Further write the finite system
1272              * Kirkwood G factor and epsilon.
1273              */
1274             fprintf(outaver, "%10g  %10.3e %10.3e %10.3e %10.3e\n",
1275                     t, M2_ave, M_ave2, M_diff, M_ave2/M2_ave);
1276
1277             if (fnadip)
1278             {
1279                 real aver;
1280                 gmx_stats_get_average(muframelsq, &aver);
1281                 fprintf(adip, "%10g %f \n", t, aver);
1282             }
1283             /*if (dipole)
1284                printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1285              */
1286             if (!bMU || (mu_aver != -1))
1287             {
1288                 /* Finite system Kirkwood G-factor */
1289                 Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1290                 /* Infinite system Kirkwood G-factor */
1291                 if (epsilonRF == 0.0)
1292                 {
1293                     g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1294                 }
1295                 else
1296                 {
1297                     g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1298                            Gk/(3*epsilon*(2*epsilonRF+1)));
1299                 }
1300
1301                 fprintf(outeps, "%10g  %10.3e %10.3e %10.3e\n", t, epsilon, Gk, g_k);
1302
1303             }
1304             else
1305             {
1306                 fprintf(outeps, "%10g  %12.8e\n", t, epsilon);
1307             }
1308         }
1309         gmx_stats_done(muframelsq);
1310
1311         if (bMU)
1312         {
1313             bCont = read_mu_from_enx(fmu, iVol, iMu, mu_t, &volume, &t, nre, fr);
1314         }
1315         else
1316         {
1317             bCont = read_next_x(oenv, status, &t, x, box);
1318         }
1319         timecheck = check_times(t);
1320     }
1321     while (bCont && (timecheck == 0) );
1322
1323     gmx_rmpbc_done(gpbc);
1324
1325     if (!bMU)
1326     {
1327         close_trj(status);
1328     }
1329
1330     ffclose(outmtot);
1331     ffclose(outaver);
1332     ffclose(outeps);
1333
1334     if (fnadip)
1335     {
1336         ffclose(adip);
1337     }
1338
1339     if (cosaver)
1340     {
1341         ffclose(caver);
1342     }
1343
1344     if (dip3d)
1345     {
1346         fprintf(dip3d, "set xrange [0.0:%4.2f]\n", box[XX][XX]);
1347         fprintf(dip3d, "set yrange [0.0:%4.2f]\n", box[YY][YY]);
1348         fprintf(dip3d, "set zrange [0.0:%4.2f]\n\n", box[ZZ][ZZ]);
1349         fprintf(dip3d, "splot 'dummy.dat' using 1:2:3 w vec\n");
1350         fprintf(dip3d, "pause -1 'Hit return to continue'\n");
1351         ffclose(dip3d);
1352     }
1353
1354     if (bSlab)
1355     {
1356         dump_slab_dipoles(slabfn, idim, nslices, slab_dipoles, box, teller, oenv);
1357         sfree(slab_dipoles);
1358     }
1359
1360     vol_aver /= teller;
1361     printf("Average volume over run is %g\n", vol_aver);
1362     if (bGkr)
1363     {
1364         print_gkrbin(gkrfn, gkrbin, gnx[0], teller, vol_aver, oenv);
1365         print_cmap(cmap, gkrbin, nlevels);
1366     }
1367     /* Autocorrelation function */
1368     if (bCorr)
1369     {
1370         if (teller < 2)
1371         {
1372             printf("Not enough frames for autocorrelation\n");
1373         }
1374         else
1375         {
1376             dt = (t1 - t0)/(teller-1);
1377             printf("t0 %g, t %g, teller %d\n", t0, t, teller);
1378
1379             mode = eacVector;
1380
1381             if (bTotal)
1382             {
1383                 do_autocorr(corf, oenv, "Autocorrelation Function of Total Dipole",
1384                             teller, 1, muall, dt, mode, TRUE);
1385             }
1386             else
1387             {
1388                 do_autocorr(corf, oenv, "Dipole Autocorrelation Function",
1389                             teller, gnx_tot, muall, dt,
1390                             mode, strcmp(corrtype, "molsep"));
1391             }
1392         }
1393     }
1394     if (!bMU)
1395     {
1396         real aver, sigma, error, lsq;
1397
1398         gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
1399         printf("\nDipole moment (Debye)\n");
1400         printf("---------------------\n");
1401         printf("Average  = %8.4f  Std. Dev. = %8.4f  Error = %8.4f\n",
1402                aver, sigma, error);
1403         if (bQuad)
1404         {
1405             rvec a, s, e;
1406             int  mm;
1407             for (m = 0; (m < DIM); m++)
1408             {
1409                 gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
1410             }
1411
1412             printf("\nQuadrupole moment (Debye-Ang)\n");
1413             printf("-----------------------------\n");
1414             printf("Averages  = %8.4f  %8.4f  %8.4f\n", a[XX], a[YY], a[ZZ]);
1415             printf("Std. Dev. = %8.4f  %8.4f  %8.4f\n", s[XX], s[YY], s[ZZ]);
1416             printf("Error     = %8.4f  %8.4f  %8.4f\n", e[XX], e[YY], e[ZZ]);
1417         }
1418         printf("\n");
1419     }
1420     printf("The following averages for the complete trajectory have been calculated:\n\n");
1421     printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1422     printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1423     printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1424
1425     printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1426     printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1427     printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1428
1429     printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1430     printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1431
1432     printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1433
1434     if (!bMU || (mu_aver != -1))
1435     {
1436         printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1437         printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1438     }
1439     printf("Epsilon = %g\n", epsilon);
1440
1441     if (!bMU)
1442     {
1443         /* Write to file the dipole moment distibution during the simulation.
1444          */
1445         outdd = xvgropen(dipdist, "Dipole Moment Distribution", "mu (Debye)", "", oenv);
1446         for (i = 0; (i < ndipbin); i++)
1447         {
1448             fprintf(outdd, "%10g  %10f\n",
1449                     (i*mu_max)/ndipbin, dipole_bin[i]/(double)teller);
1450         }
1451         ffclose(outdd);
1452         sfree(dipole_bin);
1453     }
1454     if (bGkr)
1455     {
1456         done_gkrbin(&gkrbin);
1457     }
1458 }
1459
1460 void dipole_atom2molindex(int *n, int *index, t_block *mols)
1461 {
1462     int nmol, i, j, m;
1463
1464     nmol = 0;
1465     i    = 0;
1466     while (i < *n)
1467     {
1468         m = 0;
1469         while (m < mols->nr && index[i] != mols->index[m])
1470         {
1471             m++;
1472         }
1473         if (m == mols->nr)
1474         {
1475             gmx_fatal(FARGS, "index[%d]=%d does not correspond to the first atom of a molecule", i+1, index[i]+1);
1476         }
1477         for (j = mols->index[m]; j < mols->index[m+1]; j++)
1478         {
1479             if (i >= *n || index[i] != j)
1480             {
1481                 gmx_fatal(FARGS, "The index group is not a set of whole molecules");
1482             }
1483             i++;
1484         }
1485         /* Modify the index in place */
1486         index[nmol++] = m;
1487     }
1488     printf("There are %d molecules in the selection\n", nmol);
1489
1490     *n = nmol;
1491 }
1492 int gmx_dipoles(int argc, char *argv[])
1493 {
1494     const char    *desc[] = {
1495         "[TT]g_dipoles[tt] computes the total dipole plus fluctuations of a simulation",
1496         "system. From this you can compute e.g. the dielectric constant for",
1497         "low-dielectric media.",
1498         "For molecules with a net charge, the net charge is subtracted at",
1499         "center of mass of the molecule.[PAR]",
1500         "The file [TT]Mtot.xvg[tt] contains the total dipole moment of a frame, the",
1501         "components as well as the norm of the vector.",
1502         "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",
1503         "simulation.",
1504         "The file [TT]dipdist.xvg[tt] contains the distribution of dipole moments during",
1505         "the simulation",
1506         "The value of [TT]-mumax[tt] is used as the highest value in the distribution graph.[PAR]",
1507         "Furthermore, the dipole autocorrelation function will be computed when",
1508         "option [TT]-corr[tt] is used. The output file name is given with the [TT]-c[tt]",
1509         "option.",
1510         "The correlation functions can be averaged over all molecules",
1511         "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1512         "or it can be computed over the total dipole moment of the simulation box",
1513         "([TT]total[tt]).[PAR]",
1514         "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1515         "G-factor, as well as the average cosine of the angle between the dipoles",
1516         "as a function of the distance. The plot also includes gOO and hOO",
1517         "according to Nymand & Linse, J. Chem. Phys. 112 (2000) pp 6386-6395. In the same plot, ",
1518         "we also include the energy per scale computed by taking the inner product of",
1519         "the dipoles divided by the distance to the third power.[PAR]",
1520         "[PAR]",
1521         "EXAMPLES[PAR]",
1522         "[TT]g_dipoles -corr mol -P 1 -o dip_sqr -mu 2.273 -mumax 5.0[tt][PAR]",
1523         "This will calculate the autocorrelation function of the molecular",
1524         "dipoles using a first order Legendre polynomial of the angle of the",
1525         "dipole vector and itself a time t later. For this calculation 1001",
1526         "frames will be used. Further, the dielectric constant will be calculated",
1527         "using an [TT]-epsilonRF[tt] of infinity (default), temperature of 300 K (default) and",
1528         "an average dipole moment of the molecule of 2.273 (SPC). For the",
1529         "distribution function a maximum of 5.0 will be used."
1530     };
1531     real           mu_max     = 5, mu_aver = -1, rcmax = 0;
1532     real           epsilonRF  = 0.0, temp = 300;
1533     gmx_bool       bAverCorr  = FALSE, bMolCorr = FALSE, bPairs = TRUE, bPhi = FALSE, bQuad = FALSE;
1534     const char    *corrtype[] = {NULL, "none", "mol", "molsep", "total", NULL};
1535     const char    *axtitle    = "Z";
1536     int            nslices    = 10; /* nr of slices defined       */
1537     int            skip       = 0, nFA = 0, nFB = 0, ncos = 1;
1538     int            nlevels    = 20, ndegrees = 90;
1539     output_env_t   oenv;
1540     t_pargs        pa[] = {
1541         { "-mu",       FALSE, etREAL, {&mu_aver},
1542           "dipole of a single molecule (in Debye)" },
1543         { "-mumax",    FALSE, etREAL, {&mu_max},
1544           "max dipole in Debye (for histogram)" },
1545         { "-epsilonRF", FALSE, etREAL, {&epsilonRF},
1546           "[GRK]epsilon[grk] of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1547         { "-skip",     FALSE, etINT, {&skip},
1548           "Skip steps in the output (but not in the computations)" },
1549         { "-temp",     FALSE, etREAL, {&temp},
1550           "Average temperature of the simulation (needed for dielectric constant calculation)" },
1551         { "-corr",     FALSE, etENUM, {corrtype},
1552           "Correlation function to calculate" },
1553         { "-pairs",    FALSE, etBOOL, {&bPairs},
1554           "Calculate [MAG][COS][GRK]theta[grk][cos][mag] between all pairs of molecules. May be slow" },
1555         { "-quad",     FALSE, etBOOL, {&bQuad},
1556           "Take quadrupole into account"},
1557         { "-ncos",     FALSE, etINT, {&ncos},
1558           "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." },
1559         { "-axis",     FALSE, etSTR, {&axtitle},
1560           "Take the normal on the computational box in direction X, Y or Z." },
1561         { "-sl",       FALSE, etINT, {&nslices},
1562           "Divide the box into this number of slices." },
1563         { "-gkratom",  FALSE, etINT, {&nFA},
1564           "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" },
1565         { "-gkratom2", FALSE, etINT, {&nFB},
1566           "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1567         { "-rcmax",    FALSE, etREAL, {&rcmax},
1568           "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterion based on the box length will be used." },
1569         { "-phi",      FALSE, etBOOL, {&bPhi},
1570           "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." },
1571         { "-nlevels",  FALSE, etINT, {&nlevels},
1572           "Number of colors in the cmap output" },
1573         { "-ndegrees", FALSE, etINT, {&ndegrees},
1574           "Number of divisions on the [IT]y[it]-axis in the cmap output (for 180 degrees)" }
1575     };
1576     int           *gnx;
1577     int            nFF[2];
1578     atom_id      **grpindex;
1579     char         **grpname = NULL;
1580     gmx_bool       bCorr, bGkr, bMU, bSlab;
1581     t_filenm       fnm[] = {
1582         { efEDR, "-en", NULL,         ffOPTRD },
1583         { efTRX, "-f", NULL,           ffREAD },
1584         { efTPX, NULL, NULL,           ffREAD },
1585         { efNDX, NULL, NULL,           ffOPTRD },
1586         { efXVG, "-o",   "Mtot",       ffWRITE },
1587         { efXVG, "-eps", "epsilon",    ffWRITE },
1588         { efXVG, "-a",   "aver",       ffWRITE },
1589         { efXVG, "-d",   "dipdist",    ffWRITE },
1590         { efXVG, "-c",   "dipcorr",    ffOPTWR },
1591         { efXVG, "-g",   "gkr",        ffOPTWR },
1592         { efXVG, "-adip", "adip",       ffOPTWR },
1593         { efXVG, "-dip3d", "dip3d",    ffOPTWR },
1594         { efXVG, "-cos", "cosaver",    ffOPTWR },
1595         { efXPM, "-cmap", "cmap",       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     bGkr  = opt2bSet("-g", NFILE, fnm);
1624     if (opt2parg_bSet("-ncos", asize(pa), pa))
1625     {
1626         if ((ncos != 1) && (ncos != 2))
1627         {
1628             gmx_fatal(FARGS, "ncos has to be either 1 or 2");
1629         }
1630         bGkr = TRUE;
1631     }
1632     bSlab = (opt2bSet("-slab", NFILE, fnm) || opt2parg_bSet("-sl", asize(pa), pa) ||
1633              opt2parg_bSet("-axis", asize(pa), pa));
1634     if (bMU)
1635     {
1636         bAverCorr = TRUE;
1637         if (bQuad)
1638         {
1639             printf("WARNING: Can not determine quadrupoles from energy file\n");
1640             bQuad = FALSE;
1641         }
1642         if (bGkr)
1643         {
1644             printf("WARNING: Can not determine Gk(r) from energy file\n");
1645             bGkr  = FALSE;
1646             ncos  = 1;
1647         }
1648         if (mu_aver == -1)
1649         {
1650             printf("WARNING: Can not calculate Gk and gk, since you did\n"
1651                    "         not enter a valid dipole for the molecules\n");
1652         }
1653     }
1654
1655     snew(top, 1);
1656     ePBC = read_tpx_top(ftp2fn(efTPX, NFILE, fnm), NULL, box,
1657                         &natoms, NULL, NULL, NULL, top);
1658
1659     snew(gnx, ncos);
1660     snew(grpname, ncos);
1661     snew(grpindex, ncos);
1662     get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),
1663               ncos, gnx, grpindex, grpname);
1664     for (k = 0; (k < ncos); k++)
1665     {
1666         dipole_atom2molindex(&gnx[k], grpindex[k], &(top->mols));
1667         neutralize_mols(gnx[k], grpindex[k], &(top->mols), top->atoms.atom);
1668     }
1669     nFF[0] = nFA;
1670     nFF[1] = nFB;
1671     do_dip(top, ePBC, det(box), ftp2fn(efTRX, NFILE, fnm),
1672            opt2fn("-o", NFILE, fnm), opt2fn("-eps", NFILE, fnm),
1673            opt2fn("-a", NFILE, fnm), opt2fn("-d", NFILE, fnm),
1674            opt2fn_null("-cos", NFILE, fnm),
1675            opt2fn_null("-dip3d", NFILE, fnm), opt2fn_null("-adip", NFILE, fnm),
1676            bPairs, corrtype[0],
1677            opt2fn("-c", NFILE, fnm),
1678            bGkr,    opt2fn("-g", NFILE, fnm),
1679            bPhi,    &nlevels,  ndegrees,
1680            ncos,
1681            opt2fn("-cmap", NFILE, fnm), rcmax,
1682            bQuad, bMU,     opt2fn("-en", NFILE, fnm),
1683            gnx, grpindex, mu_max, mu_aver, epsilonRF, temp, nFF, skip,
1684            bSlab, nslices, axtitle, opt2fn("-slab", NFILE, fnm), oenv);
1685
1686     do_view(oenv, opt2fn("-o", NFILE, fnm), "-autoscale xy -nxy");
1687     do_view(oenv, opt2fn("-eps", NFILE, fnm), "-autoscale xy -nxy");
1688     do_view(oenv, opt2fn("-a", NFILE, fnm), "-autoscale xy -nxy");
1689     do_view(oenv, opt2fn("-d", NFILE, fnm), "-autoscale xy");
1690     do_view(oenv, opt2fn("-c", NFILE, fnm), "-autoscale xy");
1691
1692     return 0;
1693 }