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