src/tools/gmx_dipoles.c

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