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