src/tools/gmx_dipoles.c

   1 /*
   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     ddc1 += fabs(cos_angle(dip[i],xxx));
 598     ddc2 += fabs(cos_angle(dip[i],yyy));
 599     ddc3 += fabs(cos_angle(dip[i],zzz));
 600     if (bPairs)
 601       for(j=i+1; (j<n); j++,k++) {
 602         dc  = cos_angle(dip[i],dip[j]);
 603         d  += fabs(dc);
 604       }
 605   }
 606   *dd  = d/k;
 607   axis[XX] = ddc1/n;
 608   axis[YY] = ddc2/n;
 609   axis[ZZ] = ddc3/n;
 610 }
 611
 612 static void do_dip(t_topology *top,int ePBC,real volume,
 613                    const char *fn,
 614                    const char *out_mtot,const char *out_eps,
 615                    const char *out_aver, const char *dipdist,
 616                    const char *cosaver, const char *fndip3d,
 617                    const char *fnadip,  gmx_bool bPairs,
 618                    const char *corrtype,const char *corf,
 619                    gmx_bool bGkr,     const char *gkrfn,
 620                    gmx_bool bPhi,     int  *nlevels,  int ndegrees,
 621                    int  ncos,
 622                    const char *cmap,    real rcmax,
 623                    gmx_bool bQuad,    const char *quadfn,
 624                    gmx_bool bMU,      const char *mufn,
 625                    int  *gnx,     int  *molindex[],
 626                    real mu_max,   real mu_aver,
 627                    real epsilonRF,real temp,
 628                    int  *gkatom,  int skip,
 629                    gmx_bool bSlab,    int nslices,
 630                    const char *axtitle, const char *slabfn,
 631                    const output_env_t oenv)
 632 {
 633   const char *leg_mtot[] = {
 634     "M\\sx \\N",
 635     "M\\sy \\N",
 636     "M\\sz \\N",
 637     "|M\\stot \\N|"
 638   };
 639 #define NLEGMTOT asize(leg_mtot)
 640   const char *leg_eps[] = {
 641     "epsilon",
 642     "G\\sk",
 643     "g\\sk"
 644   };
 645 #define NLEGEPS asize(leg_eps)
 646   const char *leg_aver[] = {
 647     "< |M|\\S2\\N >",
 648     "< |M| >\\S2\\N",
 649     "< |M|\\S2\\N > - < |M| >\\S2\\N",
 650     "< |M| >\\S2\\N / < |M|\\S2\\N >"
 651   };
 652 #define NLEGAVER asize(leg_aver)
 653   const char *leg_cosaver[] = {
 654     "\\f{4}<|cos\\f{12}q\\f{4}\\sij\\N|>",
 655     "RMSD cos",
 656     "\\f{4}<|cos\\f{12}q\\f{4}\\siX\\N|>",
 657     "\\f{4}<|cos\\f{12}q\\f{4}\\siY\\N|>",
 658     "\\f{4}<|cos\\f{12}q\\f{4}\\siZ\\N|>"
 659   };
 660 #define NLEGCOSAVER asize(leg_cosaver)
 661   const char *leg_adip[] = {
 662     "<mu>",
 663     "Std. Dev.",
 664     "Error"
 665   };
 666 #define NLEGADIP asize(leg_adip)
 667
 668   FILE       *outdd,*outmtot,*outaver,*outeps,*caver=NULL;
 669   FILE       *dip3d=NULL,*adip=NULL;
 670   rvec       *x,*dipole=NULL,mu_t,quad,*dipsp=NULL;
 671   t_gkrbin   *gkrbin = NULL;
 672   gmx_enxnm_t *enm=NULL;
 673   t_enxframe *fr;
 674   int        nframes=1000,nre,timecheck=0,ncolour=0;
 675   ener_file_t fmu=NULL;
 676   int        i,j,k,n,m,natom=0,nmol,gnx_tot,teller,tel3;
 677   t_trxstatus *status;
 678   int        *dipole_bin,ndipbin,ibin,iVol,step,idim=-1;
 679   unsigned long mode;
 680   char       buf[STRLEN];
 681   real       rcut=0,t,t0,t1,dt,lambda,dd,rms_cos;
 682   rvec       dipaxis;
 683   matrix     box;
 684   gmx_bool       bCorr,bTotal,bCont;
 685   double     M_diff=0,epsilon,invtel,vol_aver;
 686   double     mu_ave,mu_mol,M2_ave=0,M_ave2=0,M_av[DIM],M_av2[DIM];
 687   double     M[3],M2[3],M4[3],Gk=0,g_k=0;
 688   gmx_stats_t Mx,My,Mz,Msq,Vol,*Qlsq,mulsq,muframelsq=NULL;
 689   ivec       iMu;
 690   real       **muall=NULL;
 691   rvec       *slab_dipoles=NULL;
 692   t_atom     *atom=NULL;
 693   t_block    *mols=NULL;
 694   gmx_rmpbc_t gpbc=NULL;
 695
 696   gnx_tot = gnx[0];
 697   if (ncos > 1) {
 698     gnx_tot += gnx[1];
 699   }
 700
 701   vol_aver = 0.0;
 702
 703   iVol=-1;
 704   if (bMU) {
 705     fmu = open_enx(mufn,"r");
 706     do_enxnms(fmu,&nre,&enm);
 707
 708     /* Determine the indexes of the energy grps we need */
 709     for (i=0; (i<nre); i++) {
 710       if (strstr(enm[i].name,"Volume"))
 711         iVol=i;
 712       else if (strstr(enm[i].name,"Mu-X"))
 713         iMu[XX]=i;
 714       else if (strstr(enm[i].name,"Mu-Y"))
 715         iMu[YY]=i;
 716       else if (strstr(enm[i].name,"Mu-Z"))
 717         iMu[ZZ]=i;
 718     }
 719   }
 720   else {
 721     atom = top->atoms.atom;
 722     mols = &(top->mols);
 723   }
 724
 725   if (iVol == -1)
 726     printf("Using Volume from topology: %g nm^3\n",volume);
 727
 728   /* Correlation stuff */
 729   bCorr  = (corrtype[0] != 'n');
 730   bTotal = (corrtype[0] == 't');
 731   if (bCorr) {
 732     if (bTotal) {
 733       snew(muall,1);
 734       snew(muall[0],nframes*DIM);
 735     }
 736     else {
 737       snew(muall,gnx[0]);
 738       for(i=0; (i<gnx[0]); i++)
 739         snew(muall[i],nframes*DIM);
 740     }
 741   }
 742
 743   /* Allocate array which contains for every molecule in a frame the
 744    * dipole moment.
 745    */
 746   if (!bMU)
 747     snew(dipole,gnx_tot);
 748
 749   /* Statistics */
 750   snew(Qlsq,DIM);
 751   for(i=0; (i<DIM); i++)
 752     Qlsq[i] = gmx_stats_init();
 753   mulsq = gmx_stats_init();
 754
 755   /* Open all the files */
 756   outmtot = xvgropen(out_mtot,
 757                      "Total dipole moment of the simulation box vs. time",
 758                      "Time (ps)","Total Dipole Moment (Debye)",oenv);
 759   outeps  = xvgropen(out_eps,"Epsilon and Kirkwood factors",
 760                      "Time (ps)","",oenv);
 761   outaver = xvgropen(out_aver,"Total dipole moment",
 762                      "Time (ps)","D",oenv);
 763   if (bSlab) {
 764     idim = axtitle[0] - 'X';
 765     if ((idim < 0) || (idim >= DIM))
 766       idim = axtitle[0] - 'x';
 767     if ((idim < 0) || (idim >= DIM))
 768       bSlab = FALSE;
 769     if (nslices < 2)
 770       bSlab = FALSE;
 771     fprintf(stderr,"axtitle = %s, nslices = %d, idim = %d\n",
 772             axtitle,nslices,idim);
 773     if (bSlab) {
 774       snew(slab_dipoles,nslices);
 775       fprintf(stderr,"Doing slab analysis\n");
 776     }
 777   }
 778
 779   if (fnadip) {
 780     adip = xvgropen(fnadip, "Average molecular dipole","Dipole (D)","",oenv);
 781     xvgr_legend(adip,NLEGADIP,leg_adip, oenv);
 782
 783   }
 784   if (cosaver) {
 785     caver = xvgropen(cosaver,bPairs ? "Average pair orientation" :
 786                      "Average absolute dipole orientation","Time (ps)","",oenv);
 787     xvgr_legend(caver,NLEGCOSAVER,bPairs ? leg_cosaver : &(leg_cosaver[1]),
 788                 oenv);
 789   }
 790
 791   if (fndip3d) {
 792     snew(dipsp,gnx_tot);
 793
 794     /* we need a dummy file for gnuplot */
 795     dip3d = (FILE *)ffopen("dummy.dat","w");
 796     fprintf(dip3d,"%f %f %f", 0.0,0.0,0.0);
 797     ffclose(dip3d);
 798
 799     dip3d = (FILE *)ffopen(fndip3d,"w");
 800     fprintf(dip3d,"# This file was created by %s\n",Program());
 801     fprintf(dip3d,"# which is part of G R O M A C S:\n");
 802     fprintf(dip3d,"#\n");
 803   }
 804
 805   /* Write legends to all the files */
 806   xvgr_legend(outmtot,NLEGMTOT,leg_mtot,oenv);
 807   xvgr_legend(outaver,NLEGAVER,leg_aver,oenv);
 808
 809   if (bMU && (mu_aver == -1))
 810     xvgr_legend(outeps,NLEGEPS-2,leg_eps,oenv);
 811   else
 812     xvgr_legend(outeps,NLEGEPS,leg_eps,oenv);
 813
 814   snew(fr,1);
 815   clear_rvec(mu_t);
 816   teller = 0;
 817   /* Read the first frame from energy or traj file */
 818   if (bMU)
 819     do {
 820       bCont = read_mu_from_enx(fmu,iVol,iMu,mu_t,&volume,&t,nre,fr);
 821       if (bCont) {
 822         timecheck=check_times(t);
 823         if (timecheck < 0)
 824           teller++;
 825         if ((teller % 10) == 0)
 826           fprintf(stderr,"\r Skipping Frame %6d, time: %8.3f", teller, t);
 827       }
 828       else {
 829         printf("End of %s reached\n",mufn);
 830         break;
 831       }
 832     } while (bCont && (timecheck < 0));
 833   else
 834     natom  = read_first_x(oenv,&status,fn,&t,&x,box);
 835
 836   /* Calculate spacing for dipole bin (simple histogram) */
 837   ndipbin = 1+(mu_max/0.01);
 838   snew(dipole_bin, ndipbin);
 839   epsilon    = 0.0;
 840   mu_ave     = 0.0;
 841   for(m=0; (m<DIM); m++) {
 842     M[m] = M2[m] = M4[m] = 0.0;
 843   }
 844
 845   if (bGkr) {
 846     /* Use 0.7 iso 0.5 to account for pressure scaling */
 847     /*  rcut   = 0.7*sqrt(max_cutoff2(box)); */
 848     rcut   = 0.7*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])+sqr(box[ZZ][ZZ]));
 849
 850     gkrbin = mk_gkrbin(rcut,rcmax,bPhi,ndegrees);
 851   }
 852   gpbc = gmx_rmpbc_init(&top->idef,ePBC,natom,box);
 853
 854   /* Start while loop over frames */
 855   t1 = t0 = t;
 856   teller = 0;
 857   do {
 858     if (bCorr && (teller >= nframes)) {
 859       nframes += 1000;
 860       if (bTotal) {
 861         srenew(muall[0],nframes*DIM);
 862       }
 863       else {
 864         for(i=0; (i<gnx_tot); i++)
 865           srenew(muall[i],nframes*DIM);
 866       }
 867     }
 868     t1 = t;
 869
 870     if (bMU) {
 871       /* Copy rvec into double precision local variable */
 872       for(m=0; (m<DIM); m++)
 873         M_av[m]  = mu_t[m];
 874     }
 875     else {
 876       /* Initialise */
 877       for(m=0; (m<DIM); m++) {
 878         M_av[m] = 0;
 879         M_av2[m] = 0;
 880       }
 881       gmx_rmpbc(gpbc,natom,box,x);
 882
 883       muframelsq = gmx_stats_init();
 884       /* Begin loop of all molecules in frame */
 885       for(n=0; (n<ncos); n++) {
 886         for(i=0; (i<gnx[n]); i++) {
 887           int gi,ind0,ind1;
 888
 889           ind0  = mols->index[molindex[n][i]];
 890           ind1  = mols->index[molindex[n][i]+1];
 891
 892           mol_dip(ind0,ind1,x,atom,dipole[i]);
 893           gmx_stats_add_point(mulsq,0,norm(dipole[i]),0,0);
 894           gmx_stats_add_point(muframelsq,0,norm(dipole[i]),0,0);
 895           if (bSlab)
 896             update_slab_dipoles(ind0,ind1,x,
 897                                 dipole[i],idim,nslices,slab_dipoles,box);
 898           if (bQuad) {
 899             mol_quad(ind0,ind1,x,atom,quad);
 900             for(m=0; (m<DIM); m++)
 901               gmx_stats_add_point(Qlsq[m],0,quad[m],0,0);
 902           }
 903           if (bCorr && !bTotal) {
 904             tel3=DIM*teller;
 905             muall[i][tel3+XX] = dipole[i][XX];
 906             muall[i][tel3+YY] = dipole[i][YY];
 907             muall[i][tel3+ZZ] = dipole[i][ZZ];
 908           }
 909           mu_mol = 0.0;
 910           for(m=0; (m<DIM); m++) {
 911             M_av[m]  += dipole[i][m];               /* M per frame */
 912             mu_mol   += dipole[i][m]*dipole[i][m];  /* calc. mu for distribution */
 913           }
 914           mu_mol = sqrt(mu_mol);
 915
 916           mu_ave += mu_mol;                         /* calc. the average mu */
 917
 918           /* Update the dipole distribution */
 919           ibin = (int)(ndipbin*mu_mol/mu_max + 0.5);
 920           if (ibin < ndipbin)
 921             dipole_bin[ibin]++;
 922
 923           if (fndip3d) {
 924             rvec2sprvec(dipole[i],dipsp[i]);
 925
 926             if (dipsp[i][YY] > -M_PI && dipsp[i][YY] < -0.5*M_PI) {
 927               if (dipsp[i][ZZ] < 0.5 * M_PI) {
 928                 ncolour = 1;
 929               } else {
 930                 ncolour = 2;
 931               }
 932             }else if (dipsp[i][YY] > -0.5*M_PI && dipsp[i][YY] < 0.0*M_PI) {
 933               if (dipsp[i][ZZ] < 0.5 * M_PI) {
 934                 ncolour = 3;
 935               } else {
 936                 ncolour = 4;
 937               }
 938             }else if (dipsp[i][YY] > 0.0 && dipsp[i][YY] < 0.5*M_PI) {
 939               if (dipsp[i][ZZ] < 0.5 * M_PI) {
 940                 ncolour = 5;
 941               } else {
 942                 ncolour = 6;
 943               }
 944             }else if (dipsp[i][YY] > 0.5*M_PI && dipsp[i][YY] < M_PI) {
 945               if (dipsp[i][ZZ] < 0.5 * M_PI) {
 946                 ncolour = 7;
 947               } else {
 948                 ncolour = 8;
 949               }
 950             }
 951             if (dip3d)
 952               fprintf(dip3d,"set arrow %d from %f, %f, %f to %f, %f, %f lt %d  # %d %d\n",
 953                       i+1,
 954                       x[ind0][XX],
 955                       x[ind0][YY],
 956                       x[ind0][ZZ],
 957                       x[ind0][XX]+dipole[i][XX]/25,
 958                       x[ind0][YY]+dipole[i][YY]/25,
 959                       x[ind0][ZZ]+dipole[i][ZZ]/25,
 960                       ncolour, ind0, i);
 961           }
 962         } /* End loop of all molecules in frame */
 963
 964         if (dip3d) {
 965           fprintf(dip3d,"set title \"t = %4.3f\"\n",t);
 966           fprintf(dip3d,"set xrange [0.0:%4.2f]\n",box[XX][XX]);
 967           fprintf(dip3d,"set yrange [0.0:%4.2f]\n",box[YY][YY]);
 968           fprintf(dip3d,"set zrange [0.0:%4.2f]\n\n",box[ZZ][ZZ]);
 969           fprintf(dip3d,"splot 'dummy.dat' using 1:2:3 w vec\n");
 970           fprintf(dip3d,"pause -1 'Hit return to continue'\n");
 971         }
 972       }
 973     }
 974     /* Compute square of total dipole */
 975     for(m=0; (m<DIM); m++)
 976       M_av2[m] = M_av[m]*M_av[m];
 977
 978     if (cosaver) {
 979       compute_avercos(gnx_tot,dipole,&dd,dipaxis,bPairs);
 980       rms_cos = sqrt(sqr(dipaxis[XX]-0.5)+
 981                      sqr(dipaxis[YY]-0.5)+
 982                      sqr(dipaxis[ZZ]-0.5));
 983       if (bPairs)
 984         fprintf(caver,"%10.3e  %10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
 985                 t,dd,rms_cos,dipaxis[XX],dipaxis[YY],dipaxis[ZZ]);
 986       else
 987         fprintf(caver,"%10.3e  %10.3e  %10.3e  %10.3e  %10.3e\n",
 988                 t,rms_cos,dipaxis[XX],dipaxis[YY],dipaxis[ZZ]);
 989     }
 990
 991     if (bGkr) {
 992       do_gkr(gkrbin,ncos,gnx,molindex,mols->index,x,dipole,ePBC,box,
 993              atom,gkatom);
 994     }
 995
 996     if (bTotal) {
 997       tel3 = DIM*teller;
 998       muall[0][tel3+XX] = M_av[XX];
 999       muall[0][tel3+YY] = M_av[YY];
1000       muall[0][tel3+ZZ] = M_av[ZZ];
1001     }
1002
1003     /* Write to file the total dipole moment of the box, and its components
1004      * for this frame.
1005      */
1006     if ((skip == 0) || ((teller % skip) == 0))
1007       fprintf(outmtot,"%10g  %12.8e %12.8e %12.8e %12.8e\n",
1008               t,M_av[XX],M_av[YY],M_av[ZZ],
1009               sqrt(M_av2[XX]+M_av2[YY]+M_av2[ZZ]));
1010
1011     for(m=0; (m<DIM); m++) {
1012       M[m]  += M_av[m];
1013       M2[m] += M_av2[m];
1014       M4[m] += sqr(M_av2[m]);
1015     }
1016     /* Increment loop counter */
1017     teller++;
1018
1019     /* Calculate for output the running averages */
1020     invtel  = 1.0/teller;
1021     M2_ave  = (M2[XX]+M2[YY]+M2[ZZ])*invtel;
1022     M_ave2  = invtel*(invtel*(M[XX]*M[XX] + M[YY]*M[YY] + M[ZZ]*M[ZZ]));
1023     M_diff  = M2_ave - M_ave2;
1024
1025     /* Compute volume from box in traj, else we use the one from above */
1026     if (!bMU)
1027       volume  = det(box);
1028     vol_aver += volume;
1029
1030     epsilon = calc_eps(M_diff,(vol_aver/teller),epsilonRF,temp);
1031
1032     /* Calculate running average for dipole */
1033     if (mu_ave != 0)
1034       mu_aver = (mu_ave/gnx_tot)*invtel;
1035
1036     if ((skip == 0) || ((teller % skip) == 0)) {
1037       /* Write to file < |M|^2 >, |< M >|^2. And the difference between
1038        * the two. Here M is sum mu_i. Further write the finite system
1039        * Kirkwood G factor and epsilon.
1040        */
1041       fprintf(outaver,"%10g  %10.3e %10.3e %10.3e %10.3e\n",
1042               t,M2_ave,M_ave2,M_diff,M_ave2/M2_ave);
1043
1044       if (fnadip) {
1045         real aver;
1046         gmx_stats_get_average(muframelsq,&aver);
1047         fprintf(adip, "%10g %f \n", t,aver);
1048       }
1049       /*if (dipole)
1050         printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
1051       */
1052       if (!bMU || (mu_aver != -1)) {
1053         /* Finite system Kirkwood G-factor */
1054         Gk = M_diff/(gnx_tot*mu_aver*mu_aver);
1055         /* Infinite system Kirkwood G-factor */
1056         if (epsilonRF == 0.0)
1057           g_k = ((2*epsilon+1)*Gk/(3*epsilon));
1058         else
1059           g_k = ((2*epsilonRF+epsilon)*(2*epsilon+1)*
1060                  Gk/(3*epsilon*(2*epsilonRF+1)));
1061
1062         fprintf(outeps,"%10g  %10.3e %10.3e %10.3e\n",t,epsilon,Gk,g_k);
1063
1064       }
1065       else
1066         fprintf(outeps,"%10g  %12.8e\n",t,epsilon);
1067     }
1068
1069     if (bMU)
1070       bCont = read_mu_from_enx(fmu,iVol,iMu,mu_t,&volume,&t,nre,fr);
1071     else
1072       bCont = read_next_x(oenv,status,&t,natom,x,box);
1073   } while (bCont);
1074
1075   gmx_rmpbc_done(gpbc);
1076
1077   if (!bMU)
1078     close_trj(status);
1079
1080   ffclose(outmtot);
1081   ffclose(outaver);
1082   ffclose(outeps);
1083
1084   if (fnadip)
1085     ffclose(adip);
1086
1087   if (cosaver)
1088     ffclose(caver);
1089
1090   if (dip3d) {
1091     fprintf(dip3d,"set xrange [0.0:%4.2f]\n",box[XX][XX]);
1092     fprintf(dip3d,"set yrange [0.0:%4.2f]\n",box[YY][YY]);
1093     fprintf(dip3d,"set zrange [0.0:%4.2f]\n\n",box[ZZ][ZZ]);
1094     fprintf(dip3d,"splot 'dummy.dat' using 1:2:3 w vec\n");
1095     fprintf(dip3d,"pause -1 'Hit return to continue'\n");
1096     ffclose(dip3d);
1097   }
1098
1099   if (bSlab) {
1100     dump_slab_dipoles(slabfn,idim,nslices,slab_dipoles,box,teller,oenv);
1101     sfree(slab_dipoles);
1102   }
1103
1104   vol_aver /= teller;
1105   printf("Average volume over run is %g\n",vol_aver);
1106   if (bGkr) {
1107     print_gkrbin(gkrfn,gkrbin,gnx[0],teller,vol_aver,oenv);
1108     print_cmap(cmap,gkrbin,nlevels);
1109   }
1110   /* Autocorrelation function */
1111   if (bCorr) {
1112     if (teller < 2) {
1113       printf("Not enough frames for autocorrelation\n");
1114     }
1115     else {
1116       dt=(t1 - t0)/(teller-1);
1117       printf("t0 %g, t %g, teller %d\n", t0,t,teller);
1118
1119       mode = eacVector;
1120
1121       if (bTotal)
1122         do_autocorr(corf,oenv,"Autocorrelation Function of Total Dipole",
1123                     teller,1,muall,dt,mode,TRUE);
1124       else
1125         do_autocorr(corf,oenv,"Dipole Autocorrelation Function",
1126                     teller,gnx_tot,muall,dt,
1127                     mode,strcmp(corrtype,"molsep"));
1128     }
1129   }
1130   if (!bMU) {
1131     real aver,sigma,error,lsq;
1132
1133     gmx_stats_get_ase(mulsq,&aver,&sigma,&error);
1134     printf("\nDipole moment (Debye)\n");
1135     printf("---------------------\n");
1136     printf("Average  = %8.4f  Std. Dev. = %8.4f  Error = %8.4f\n",
1137            aver,sigma,error);
1138     if (bQuad) {
1139       rvec a,s,e;
1140       int mm;
1141       for(m=0; (m<DIM); m++)
1142         gmx_stats_get_ase(mulsq,&(a[m]),&(s[m]),&(e[m]));
1143
1144       printf("\nQuadrupole moment (Debye-Ang)\n");
1145       printf("-----------------------------\n");
1146       printf("Averages  = %8.4f  %8.4f  %8.4f\n",a[XX],a[YY],a[ZZ]);
1147       printf("Std. Dev. = %8.4f  %8.4f  %8.4f\n",s[XX],s[YY],s[ZZ]);
1148       printf("Error     = %8.4f  %8.4f  %8.4f\n",e[XX],e[YY],e[ZZ]);
1149     }
1150     printf("\n");
1151   }
1152   printf("The following averages for the complete trajectory have been calculated:\n\n");
1153   printf(" Total < M_x > = %g Debye\n", M[XX]/teller);
1154   printf(" Total < M_y > = %g Debye\n", M[YY]/teller);
1155   printf(" Total < M_z > = %g Debye\n\n", M[ZZ]/teller);
1156
1157   printf(" Total < M_x^2 > = %g Debye^2\n", M2[XX]/teller);
1158   printf(" Total < M_y^2 > = %g Debye^2\n", M2[YY]/teller);
1159   printf(" Total < M_z^2 > = %g Debye^2\n\n", M2[ZZ]/teller);
1160
1161   printf(" Total < |M|^2 > = %g Debye^2\n", M2_ave);
1162   printf(" Total |< M >|^2 = %g Debye^2\n\n", M_ave2);
1163
1164   printf(" < |M|^2 > - |< M >|^2 = %g Debye^2\n\n", M_diff);
1165
1166   if (!bMU || (mu_aver != -1)) {
1167     printf("Finite system Kirkwood g factor G_k = %g\n", Gk);
1168     printf("Infinite system Kirkwood g factor g_k = %g\n\n", g_k);
1169   }
1170   printf("Epsilon = %g\n", epsilon);
1171
1172   if (!bMU) {
1173     /* Write to file the dipole moment distibution during the simulation.
1174      */
1175     outdd=xvgropen(dipdist,"Dipole Moment Distribution","mu (Debye)","",oenv);
1176     for(i=0; (i<ndipbin); i++)
1177       fprintf(outdd,"%10g  %10f\n",
1178               (i*mu_max)/ndipbin,dipole_bin[i]/(double)teller);
1179     ffclose(outdd);
1180     sfree(dipole_bin);
1181   }
1182   if (bGkr)
1183     done_gkrbin(&gkrbin);
1184 }
1185
1186 void dipole_atom2molindex(int *n,int *index,t_block *mols)
1187 {
1188   int nmol,i,j,m;
1189
1190   nmol = 0;
1191   i=0;
1192   while (i < *n) {
1193     m=0;
1194     while(m < mols->nr && index[i] != mols->index[m])
1195       m++;
1196     if (m == mols->nr)
1197       gmx_fatal(FARGS,"index[%d]=%d does not correspond to the first atom of a molecule",i+1,index[i]+1);
1198     for(j=mols->index[m]; j<mols->index[m+1]; j++) {
1199       if (i >= *n || index[i] != j)
1200         gmx_fatal(FARGS,"The index group is not a set of whole molecules");
1201       i++;
1202     }
1203     /* Modify the index in place */
1204     index[nmol++] = m;
1205   }
1206   printf("There are %d molecules in the selection\n",nmol);
1207
1208   *n = nmol;
1209 }
1210 int gmx_dipoles(int argc,char *argv[])
1211 {
1212   const char *desc[] = {
1213     "g_dipoles computes the total dipole plus fluctuations of a simulation",
1214     "system. From this you can compute e.g. the dielectric constant for",
1215     "low dielectric media.",
1216     "For molecules with a net charge, the net charge is subtracted at",
1217     "center of mass of the molecule.[PAR]",
1218     "The file Mtot.xvg contains the total dipole moment of a frame, the",
1219     "components as well as the norm of the vector.",
1220     "The file aver.xvg contains < |Mu|^2 > and |< Mu >|^2 during the",
1221     "simulation.",
1222     "The file dipdist.xvg contains the distribution of dipole moments during",
1223     "the simulation",
1224     "The mu_max is used as the highest value in the distribution graph.[PAR]",
1225     "Furthermore the dipole autocorrelation function will be computed when",
1226     "option -corr is used. The output file name is given with the [TT]-c[tt]",
1227     "option.",
1228     "The correlation functions can be averaged over all molecules",
1229     "([TT]mol[tt]), plotted per molecule separately ([TT]molsep[tt])",
1230     "or it can be computed over the total dipole moment of the simulation box",
1231     "([TT]total[tt]).[PAR]",
1232     "Option [TT]-g[tt] produces a plot of the distance dependent Kirkwood",
1233     "G-factor, as well as the average cosine of the angle between the dipoles",
1234     "as a function of the distance. The plot also includes gOO and hOO",
1235     "according to Nymand & Linse, JCP 112 (2000) pp 6386-6395. In the same plot",
1236     "we also include the energy per scale computed by taking the inner product of",
1237     "the dipoles divided by the distance to the third power.[PAR]",
1238     "[PAR]",
1239     "EXAMPLES[PAR]",
1240     "g_dipoles -corr mol -P1 -o dip_sqr -mu 2.273 -mumax 5.0 -nofft[PAR]",
1241     "This will calculate the autocorrelation function of the molecular",
1242     "dipoles using a first order Legendre polynomial of the angle of the",
1243     "dipole vector and itself a time t later. For this calculation 1001",
1244     "frames will be used. Further the dielectric constant will be calculated",
1245     "using an epsilonRF of infinity (default), temperature of 300 K (default) and",
1246     "an average dipole moment of the molecule of 2.273 (SPC). For the",
1247     "distribution function a maximum of 5.0 will be used."
1248   };
1249   real mu_max=5, mu_aver=-1,rcmax=0;
1250   real epsilonRF=0.0, temp=300;
1251   gmx_bool bAverCorr=FALSE,bMolCorr=FALSE,bPairs=TRUE,bPhi=FALSE;
1252   const char *corrtype[]={NULL, "none", "mol", "molsep", "total", NULL};
1253   const char *axtitle="Z";
1254   int  nslices = 10;      /* nr of slices defined       */
1255   int  skip=0,nFA=0,nFB=0,ncos=1;
1256   int  nlevels=20,ndegrees=90;
1257   output_env_t oenv;
1258   t_pargs pa[] = {
1259     { "-mu",       FALSE, etREAL, {&mu_aver},
1260       "dipole of a single molecule (in Debye)" },
1261     { "-mumax",    FALSE, etREAL, {&mu_max},
1262       "max dipole in Debye (for histrogram)" },
1263     { "-epsilonRF",FALSE, etREAL, {&epsilonRF},
1264       "epsilon of the reaction field used during the simulation, needed for dielectric constant calculation. WARNING: 0.0 means infinity (default)" },
1265     { "-skip",     FALSE, etINT, {&skip},
1266       "Skip steps in the output (but not in the computations)" },
1267     { "-temp",     FALSE, etREAL, {&temp},
1268       "Average temperature of the simulation (needed for dielectric constant calculation)" },
1269     { "-corr",     FALSE, etENUM, {corrtype},
1270       "Correlation function to calculate" },
1271     { "-pairs",    FALSE, etBOOL, {&bPairs},
1272       "Calculate |cos theta| between all pairs of molecules. May be slow" },
1273     { "-ncos",     FALSE, etINT, {&ncos},
1274       "Must be 1 or 2. Determines whether the <cos> is computed between all mole cules in one group, or between molecules in two different groups. This turns on the -gkr flag." },
1275     { "-axis",     FALSE, etSTR, {&axtitle},
1276       "Take the normal on the computational box in direction X, Y or Z." },
1277     { "-sl",       FALSE, etINT, {&nslices},
1278       "Divide the box in #nr slices." },
1279     { "-gkratom",  FALSE, etINT, {&nFA},
1280       "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" },
1281     { "-gkratom2", FALSE, etINT, {&nFB},
1282       "Same as previous option in case ncos = 2, i.e. dipole interaction between two groups of molecules" },
1283     { "-rcmax",    FALSE, etREAL, {&rcmax},
1284       "Maximum distance to use in the dipole orientation distribution (with ncos == 2). If zero, a criterium based on the box length will be used." },
1285     { "-phi",      FALSE, etBOOL, {&bPhi},
1286       "Plot the 'torsion angle' defined as the rotation of the two dipole vectors around the distance vector between the two molecules in the xpm file from the -cmap option. By default the cosine of the angle between the dipoles is plotted." },
1287     { "-nlevels",  FALSE, etINT, {&nlevels},
1288       "Number of colors in the cmap output" },
1289     { "-ndegrees", FALSE, etINT, {&ndegrees},
1290       "Number of divisions on the y-axis in the camp output (for 180 degrees)" }
1291   };
1292   int          *gnx;
1293   int          nFF[2];
1294   atom_id      **grpindex;
1295   char         **grpname=NULL;
1296   gmx_bool         bCorr,bQuad,bGkr,bMU,bSlab;
1297   t_filenm fnm[] = {
1298     { efEDR, "-en", NULL,         ffOPTRD },
1299     { efTRX, "-f", NULL,           ffREAD },
1300     { efTPX, NULL, NULL,           ffREAD },
1301     { efNDX, NULL, NULL,           ffOPTRD },
1302     { efXVG, "-o",   "Mtot",       ffWRITE },
1303     { efXVG, "-eps", "epsilon",    ffWRITE },
1304     { efXVG, "-a",   "aver",       ffWRITE },
1305     { efXVG, "-d",   "dipdist",    ffWRITE },
1306     { efXVG, "-c",   "dipcorr",    ffOPTWR },
1307     { efXVG, "-g",   "gkr",        ffOPTWR },
1308     { efXVG, "-adip","adip",       ffOPTWR },
1309     { efXVG, "-dip3d", "dip3d",    ffOPTWR },
1310     { efXVG, "-cos", "cosaver",    ffOPTWR },
1311     { efXPM, "-cmap","cmap",       ffOPTWR },
1312     { efXVG, "-q",   "quadrupole", ffOPTWR },
1313     { efXVG, "-slab","slab",       ffOPTWR }
1314   };
1315 #define NFILE asize(fnm)
1316   int     npargs;
1317   t_pargs *ppa;
1318   t_topology *top;
1319   int     ePBC;
1320   int     k,natoms;
1321   matrix  box;
1322
1323   CopyRight(stderr,argv[0]);
1324   npargs = asize(pa);
1325   ppa    = add_acf_pargs(&npargs,pa);
1326   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
1327                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1328
1329   printf("Using %g as mu_max and %g as the dipole moment.\n",
1330           mu_max,mu_aver);
1331   if (epsilonRF == 0.0)
1332     printf("WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity\n");
1333
1334   bMU   = opt2bSet("-en",NFILE,fnm);
1335   bQuad = opt2bSet("-q",NFILE,fnm);
1336   bGkr  = opt2bSet("-g",NFILE,fnm);
1337   if (opt2parg_bSet("-ncos",asize(pa),pa)) {
1338     if ((ncos != 1) && (ncos != 2))
1339       gmx_fatal(FARGS,"ncos has to be either 1 or 2");
1340     bGkr = TRUE;
1341   }
1342   bSlab = (opt2bSet("-slab",NFILE,fnm) || opt2parg_bSet("-sl",asize(pa),pa) ||
1343            opt2parg_bSet("-axis",asize(pa),pa));
1344   if (bMU) {
1345     bAverCorr = TRUE;
1346     if (bQuad) {
1347       printf("WARNING: Can not determine quadrupoles from energy file\n");
1348       bQuad = FALSE;
1349     }
1350     if (bGkr) {
1351       printf("WARNING: Can not determine Gk(r) from energy file\n");
1352       bGkr  = FALSE;
1353       ncos = 1;
1354     }
1355     if (mu_aver == -1)
1356       printf("WARNING: Can not calculate Gk and gk, since you did\n"
1357              "         not enter a valid dipole for the molecules\n");
1358   }
1359
1360   snew(top,1);
1361   ePBC = read_tpx_top(ftp2fn(efTPX,NFILE,fnm),NULL,box,
1362                       &natoms,NULL,NULL,NULL,top);
1363
1364   snew(gnx,ncos);
1365   snew(grpname,ncos);
1366   snew(grpindex,ncos);
1367   get_index(&top->atoms,ftp2fn_null(efNDX,NFILE,fnm),
1368             ncos,gnx,grpindex,grpname);
1369   for(k=0; (k<ncos); k++) {
1370     dipole_atom2molindex(&gnx[k],grpindex[k],&(top->mols));
1371     neutralize_mols(gnx[k],grpindex[k],&(top->mols),top->atoms.atom);
1372   }
1373   nFF[0] = nFA;
1374   nFF[1] = nFB;
1375   do_dip(top,ePBC,det(box),ftp2fn(efTRX,NFILE,fnm),
1376          opt2fn("-o",NFILE,fnm),opt2fn("-eps",NFILE,fnm),
1377          opt2fn("-a",NFILE,fnm),opt2fn("-d",NFILE,fnm),
1378          opt2fn_null("-cos",NFILE,fnm),
1379          opt2fn_null("-dip3d",NFILE,fnm),opt2fn_null("-adip",NFILE,fnm),
1380          bPairs,corrtype[0],
1381          opt2fn("-c",NFILE,fnm),
1382          bGkr,    opt2fn("-g",NFILE,fnm),
1383          bPhi,    &nlevels,  ndegrees,
1384          ncos,
1385          opt2fn("-cmap",NFILE,fnm),rcmax,
1386          bQuad,   opt2fn("-q",NFILE,fnm),
1387          bMU,     opt2fn("-en",NFILE,fnm),
1388          gnx,grpindex,mu_max,mu_aver,epsilonRF,temp,nFF,skip,
1389          bSlab,nslices,axtitle,opt2fn("-slab",NFILE,fnm),oenv);
1390
1391   do_view(oenv,opt2fn("-o",NFILE,fnm),"-autoscale xy -nxy");
1392   do_view(oenv,opt2fn("-eps",NFILE,fnm),"-autoscale xy -nxy");
1393   do_view(oenv,opt2fn("-a",NFILE,fnm),"-autoscale xy -nxy");
1394   do_view(oenv,opt2fn("-d",NFILE,fnm),"-autoscale xy");
1395   do_view(oenv,opt2fn("-c",NFILE,fnm),"-autoscale xy");
1396
1397   thanx(stderr);
1398
1399   return 0;
1400 }