src/tools/gmx_current.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  *                        VERSION 3.0
   9  *
  10  * Copyright (c) 1991-2001
  11  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  12  * University of Groningen, The Netherlands
  13  *
  14  * This program is free software; you can redistribute it and/or
  15  *
  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  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  31  *
  32  * And Hey:
  33  * Gyas ROwers Mature At Cryogenic Speed
  34  *
  35  * finished FD 09/07/08
  36  *
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42
  43 #include "statutil.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "vec.h"
  47 #include "copyrite.h"
  48 #include "statutil.h"
  49 #include "tpxio.h"
  50 #include "xvgr.h"
  51 #include "rmpbc.h"
  52 #include "pbc.h"
  53 #include "physics.h"
  54 #include "index.h"
  55 #include "gmx_statistics.h"
  56 #include "gmx_ana.h"
  57
  58 #define SQR(x) (pow(x,2.0))
  59 #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */
  60
  61
  62 static void index_atom2mol(int *n,int *index,t_block *mols)
  63 {
  64   int nat,i,nmol,mol,j;
  65
  66   nat = *n;
  67   i = 0;
  68   nmol = 0;
  69   mol = 0;
  70   while (i < nat) {
  71     while (index[i] > mols->index[mol]) {
  72       mol++;
  73       if (mol >= mols->nr)
  74         gmx_fatal(FARGS,"Atom index out of range: %d",index[i]+1);
  75     }
  76     for(j=mols->index[mol]; j<mols->index[mol+1]; j++) {
  77       if (i >= nat || index[i] != j)
  78         gmx_fatal(FARGS,"The index group does not consist of whole molecules");
  79       i++;
  80     }
  81     index[nmol++] = mol;
  82   }
  83
  84   fprintf(stderr,"\nSplit group of %d atoms into %d molecules\n",nat,nmol);
  85
  86   *n = nmol;
  87 }
  88
  89 static gmx_bool precalc(t_topology top,real mass2[],real qmol[]){
  90
  91   real mtot;
  92   real qtot;
  93   real qall;
  94   int i,j,k,l;
  95   int ai,ci;
  96   gmx_bool bNEU;
  97   ai=0;
  98   ci=0;
  99   qall=0.0;
 100
 101
 102
 103   for(i=0;i<top.mols.nr;i++){
 104     k=top.mols.index[i];
 105     l=top.mols.index[i+1];
 106     mtot=0.0;
 107     qtot=0.0;
 108
 109     for(j=k;j<l;j++){
 110       mtot+=top.atoms.atom[j].m;
 111       qtot+=top.atoms.atom[j].q;
 112     }
 113
 114     for(j=k;j<l;j++){
 115       top.atoms.atom[j].q-=top.atoms.atom[j].m*qtot/mtot;
 116       mass2[j]=top.atoms.atom[j].m/mtot;
 117       qmol[j]=qtot;
 118 }
 119
 120
 121     qall+=qtot;
 122
 123     if(qtot<0.0)
 124                 ai++;
 125     if(qtot>0.0)
 126                 ci++;
 127   }
 128
 129   if(abs(qall)>0.01){
 130           printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n",qall);
 131           bNEU=FALSE;
 132   }else
 133           bNEU=TRUE;
 134
 135   return bNEU;
 136
 137 }
 138
 139 static void remove_jump(matrix box,int natoms,rvec xp[],rvec x[]){
 140
 141   rvec hbox;
 142   int d,i,m;
 143
 144   for(d=0; d<DIM; d++)
 145     hbox[d] = 0.5*box[d][d];
 146   for(i=0; i<natoms; i++)
 147     for(m=DIM-1; m>=0; m--) {
 148       while (x[i][m]-xp[i][m] <= -hbox[m])
 149         for(d=0; d<=m; d++)
 150           x[i][d] += box[m][d];
 151       while (x[i][m]-xp[i][m] > hbox[m])
 152         for(d=0; d<=m; d++)
 153           x[i][d] -= box[m][d];
 154     }
 155 }
 156
 157 static void calc_mj(t_topology top,int ePBC,matrix box,gmx_bool bNoJump,int isize,int index0[],\
 158                 rvec fr[],rvec mj,real mass2[],real qmol[]){
 159
 160   int   i,j,k,l;
 161   rvec  tmp;
 162   rvec  center;
 163   rvec  mt1,mt2;
 164   t_pbc pbc;
 165
 166
 167   calc_box_center(ecenterRECT,box,center);
 168
 169   if(!bNoJump)
 170                 set_pbc(&pbc,ePBC,box);
 171
 172   clear_rvec(tmp);
 173
 174
 175   for(i=0;i<isize;i++){
 176         clear_rvec(mt1);
 177                 clear_rvec(mt2);
 178                 k=top.mols.index[index0[i]];
 179                 l=top.mols.index[index0[i+1]];
 180                 for(j=k;j<l;j++){
 181                         svmul(mass2[j],fr[j],tmp);
 182                         rvec_inc(mt1,tmp);
 183                 }
 184
 185                 if(bNoJump)
 186                         svmul(qmol[k],mt1,mt1);
 187                 else{
 188                         pbc_dx(&pbc,mt1,center,mt2);
 189                         svmul(qmol[k],mt2,mt1);
 190                 }
 191
 192                 rvec_inc(mj,mt1);
 193
 194   }
 195
 196 }
 197
 198
 199 static real calceps(real prefactor,real md2,real mj2,real cor,real eps_rf,gmx_bool bCOR){
 200
 201         /* bCOR determines if the correlation is computed via
 202          * static properties (FALSE) or the correlation integral (TRUE)
 203          */
 204
 205         real epsilon=0.0;
 206
 207
 208         if (bCOR) epsilon=md2-2.0*cor+mj2;
 209         else epsilon=md2+mj2+2.0*cor;
 210
 211         if (eps_rf==0.0){
 212                 epsilon=1.0+prefactor*epsilon;
 213
 214       }
 215         else{
 216                 epsilon=2.0*eps_rf+1.0+2.0*eps_rf*prefactor*epsilon;
 217                 epsilon/=2.0*eps_rf+1.0-prefactor*epsilon;
 218                 }
 219
 220
 221         return epsilon;
 222
 223     }
 224
 225
 226 static real calc_cacf(FILE *fcacf,real prefactor,real cacf[],real time[],int nfr,int vfr[],int ei,int nshift){
 227
 228         int i;
 229         real corint;
 230         real deltat=0.0;
 231         real rfr;
 232         real sigma=0.0;
 233         real sigma_ret=0.0;
 234         corint=0.0;
 235
 236         if(nfr>1){
 237   i=0;
 238
 239                 while(i<nfr){
 240
 241                                 rfr=(real) (nfr/nshift-i/nshift);
 242                                 cacf[i]/=rfr;
 243
 244                                 if(time[vfr[i]]<=time[vfr[ei]])
 245                                         sigma_ret=sigma;
 246
 247                                 fprintf(fcacf,"%.3f\t%10.6g\t%10.6g\n",time[vfr[i]],cacf[i],sigma);
 248
 249                                 if((i+1)<nfr){
 250                                         deltat=(time[vfr[i+1]]-time[vfr[i]]);
 251                                 }
 252                                         corint=2.0*deltat*cacf[i]*prefactor;
 253                                         if(i==0 || (i+1)==nfr)
 254                                                 corint*=0.5;
 255                                         sigma+=corint;
 256
 257       i++;
 258     }
 259
 260         }else
 261                 printf("Too less points.\n");
 262
 263         return sigma_ret;
 264
 265 }
 266
 267 static void calc_mjdsp(FILE *fmjdsp,real vol,real temp,real prefactor,rvec mjdsp[],real dsp2[],real time[],int nfr,real refr[]){
 268
 269         int             i;
 270         real    rtmp;
 271         real    rfr;
 272
 273
 274         fprintf(fmjdsp,"#Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n",prefactor);
 275
 276
 277
 278         for(i=0;i<nfr;i++){
 279
 280                 if(refr[i]!=0.0){
 281                         dsp2[i]*=prefactor/refr[i];
 282                         fprintf(fmjdsp,"%.3f\t%10.6g\n",time[i],dsp2[i]);
 283   }
 284
 285
 286 }
 287
 288 }
 289
 290
 291 static void dielectric(FILE *fmj,FILE *fmd,FILE *outf,FILE *fcur,FILE *mcor,
 292                        FILE *fmjdsp, gmx_bool bNoJump,gmx_bool bACF,gmx_bool bINT,
 293                        int ePBC,t_topology top, t_trxframe fr,real temp,
 294                        real trust,real bfit,real efit,real bvit,real evit,
 295                        t_trxstatus *status,int isize,int nmols, int nshift,
 296                        atom_id *index0,int indexm[],real mass2[],
 297                        real qmol[], real eps_rf, const output_env_t oenv)
 298 {
 299   int   i,j,k,l,f;
 300   int           valloc,nalloc,nfr,nvfr,m,itrust=0;
 301   int           vshfr;
 302   real  *xshfr=NULL;
 303   int   *vfr=NULL;
 304   real  refr=0.0;
 305   real  rfr=0.0;
 306   real  *cacf=NULL;
 307   real  *time=NULL;
 308   real  *djc=NULL;
 309   real  corint=0.0;
 310   real  prefactorav=0.0;
 311   real  prefactor=0.0;
 312   real  volume;
 313   real  volume_av=0.0;
 314   real  dk_s,dk_d,dk_f;
 315   real  dm_s,dm_d;
 316   real  mj=0.0;
 317   real  mj2=0.0;
 318   real  mjd=0.0;
 319   real  mjdav=0.0;
 320   real  md2=0.0;
 321   real  mdav2=0.0;
 322   real  sgk;
 323   rvec  mja_tmp;
 324   rvec  mjd_tmp;
 325   rvec  mdvec;
 326   rvec  *mu=NULL;
 327   rvec  *xp=NULL;
 328   rvec  *v0=NULL;
 329   rvec  *mjdsp=NULL;
 330   real  *dsp2=NULL;
 331   real  t0;
 332   real  rtmp;
 333   real  qtmp;
 334
 335
 336
 337   rvec  tmp;
 338   rvec  *mtrans=NULL;
 339
 340   /*
 341    * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo
 342    */
 343
 344   int bi,ei,ie,ii;
 345   real rest=0.0;
 346   real sigma=0.0;
 347   real malt=0.0;
 348   real err=0.0;
 349   real *xfit;
 350   real *yfit;
 351   gmx_rmpbc_t  gpbc=NULL;
 352
 353   /*
 354    * indices for EH
 355    */
 356
 357   ei=0;
 358   bi=0;
 359
 360   /*
 361    * indices for GK
 362    */
 363
 364   ii=0;
 365   ie=0;
 366   t0=0;
 367   sgk=0.0;
 368
 369
 370   /* This is the main loop over frames */
 371
 372
 373   nfr = 0;
 374
 375
 376   nvfr = 0;
 377   vshfr = 0;
 378   nalloc = 0;
 379   valloc = 0;
 380
 381   clear_rvec(mja_tmp);
 382   clear_rvec(mjd_tmp);
 383   clear_rvec(mdvec);
 384   clear_rvec(tmp);
 385   gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box);
 386
 387   do{
 388
 389     refr=(real)(nfr+1);
 390
 391     if(nfr >= nalloc){
 392       nalloc+=100;
 393       srenew(time,nalloc);
 394       srenew(mu,nalloc);
 395       srenew(mjdsp,nalloc);
 396       srenew(dsp2,nalloc);
 397       srenew(mtrans,nalloc);
 398       srenew(xshfr,nalloc);
 399
 400       for(i=nfr;i<nalloc;i++){
 401         clear_rvec(mjdsp[i]);
 402         clear_rvec(mu[i]);
 403         clear_rvec(mtrans[i]);
 404         dsp2[i]=0.0;
 405         xshfr[i]=0.0;
 406       }
 407     }
 408
 409
 410
 411     if (nfr==0){
 412       t0=fr.time;
 413
 414     }
 415
 416     time[nfr]=fr.time-t0;
 417
 418     if(time[nfr]<=bfit)
 419         bi=nfr;
 420                 if(time[nfr]<=efit)
 421                         ei=nfr;
 422
 423                 if(bNoJump){
 424
 425     if (xp)
 426       remove_jump(fr.box,fr.natoms,xp,fr.x);
 427                         else
 428       snew(xp,fr.natoms);
 429
 430                         for(i=0; i<fr.natoms;i++)
 431       copy_rvec(fr.x[i],xp[i]);
 432
 433     }
 434
 435                 gmx_rmpbc_trxfr(gpbc,&fr);
 436
 437                 calc_mj(top,ePBC,fr.box,bNoJump,nmols,indexm,fr.x,mtrans[nfr],mass2,qmol);
 438
 439     for(i=0;i<isize;i++){
 440       j=index0[i];
 441       svmul(top.atoms.atom[j].q,fr.x[j],fr.x[j]);
 442       rvec_inc(mu[nfr],fr.x[j]);
 443     }
 444
 445     /*if(mod(nfr,nshift)==0){*/
 446                 if(nfr%nshift==0){
 447                         for(j=nfr;j>=0;j--){
 448                                         rvec_sub(mtrans[nfr],mtrans[j],tmp);
 449                                         dsp2[nfr-j]+=norm2(tmp);
 450                                         xshfr[nfr-j]+=1.0;
 451                         }
 452                 }
 453
 454                 if (fr.bV){
 455         if(nvfr >= valloc){
 456                 valloc+=100;
 457                 srenew(vfr,valloc);
 458                 if(bINT)
 459                         srenew(djc,valloc);
 460                 srenew(v0,valloc);
 461                 if(bACF)
 462                         srenew(cacf,valloc);
 463         }
 464         if(time[nfr]<=bvit)
 465                 ii=nvfr;
 466         if(time[nfr]<=evit)
 467                                 ie=nvfr;
 468         vfr[nvfr]=nfr;
 469         clear_rvec(v0[nvfr]);
 470         if(bACF)
 471                 cacf[nvfr]=0.0;
 472         if(bINT)
 473                 djc[nvfr]=0.0;
 474         for(i=0;i<isize;i++){
 475                 j=index0[i];
 476                 svmul(mass2[j],fr.v[j],fr.v[j]);
 477                 svmul(qmol[j],fr.v[j],fr.v[j]);
 478                 rvec_inc(v0[nvfr],fr.v[j]);
 479       }
 480
 481         fprintf(fcur,"%.3f\t%.6f\t%.6f\t%.6f\n",time[nfr],v0[nfr][XX],v0[nfr][YY],v0[nfr][ZZ]);
 482         if(bACF || bINT)
 483                 /*if(mod(nvfr,nshift)==0){*/
 484                                 if(nvfr%nshift==0){
 485                         for(j=nvfr;j>=0;j--){
 486                                 if(bACF)
 487                                         cacf[nvfr-j]+=iprod(v0[nvfr],v0[j]);
 488                                 if(bINT)
 489                                         djc[nvfr-j]+=iprod(mu[vfr[j]],v0[nvfr]);
 490                         }
 491                         vshfr++;
 492                 }
 493                         nvfr++;
 494     }
 495
 496     volume = det(fr.box);
 497     volume_av += volume;
 498
 499     rvec_inc(mja_tmp,mtrans[nfr]);
 500     mjd+=iprod(mu[nfr],mtrans[nfr]);
 501     rvec_inc(mdvec,mu[nfr]);
 502
 503     mj2+=iprod(mtrans[nfr],mtrans[nfr]);
 504     md2+=iprod(mu[nfr],mu[nfr]);
 505
 506     fprintf(fmj,"%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n",time[nfr],mtrans[nfr][XX],mtrans[nfr][YY],mtrans[nfr][ZZ],mj2/refr,norm(mja_tmp)/refr);
 507     fprintf(fmd,"%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n",    \
 508             time[nfr],mu[nfr][XX],mu[nfr][YY],mu[nfr][ZZ],md2/refr,norm(mdvec)/refr);
 509
 510     nfr++;
 511
 512   }while(read_next_frame(oenv,status,&fr));
 513
 514   gmx_rmpbc_done(gpbc);
 515
 516   volume_av/=refr;
 517
 518         prefactor=1.0;
 519         prefactor/=3.0*EPSILON0*volume_av*BOLTZ*temp;
 520
 521
 522         prefactorav=E_CHARGE*E_CHARGE;
 523         prefactorav/=volume_av*BOLTZMANN*temp*NANO*6.0;
 524
 525         fprintf(stderr,"Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n",prefactorav);
 526
 527         calc_mjdsp(fmjdsp,volume_av,temp,prefactorav,mjdsp,dsp2,time,nfr,xshfr);
 528
 529   /*
 530    * Now we can average and calculate the correlation functions
 531    */
 532
 533
 534   mj2/=refr;
 535   mjd/=refr;
 536   md2/=refr;
 537
 538   svmul(1.0/refr,mdvec,mdvec);
 539   svmul(1.0/refr,mja_tmp,mja_tmp);
 540
 541   mdav2=norm2(mdvec);
 542   mj=norm2(mja_tmp);
 543   mjdav=iprod(mdvec,mja_tmp);
 544
 545
 546   printf("\n\nAverage translational dipole moment M_J [enm] after %d frames (|M|^2): %f %f %f (%f)\n",nfr,mja_tmp[XX],mja_tmp[YY],mja_tmp[ZZ],mj2);
 547   printf("\n\nAverage molecular dipole moment M_D [enm] after %d frames (|M|^2): %f %f %f (%f)\n",nfr,mdvec[XX],mdvec[YY],mdvec[ZZ],md2);
 548
 549   if (v0!=NULL){
 550         trust*=(real) nvfr;
 551                 itrust=(int) trust;
 552
 553         if(bINT){
 554
 555                 printf("\nCalculating M_D - current correlation integral ... \n");
 556                 corint=calc_cacf(mcor,prefactorav/EPSI0,djc,time,nvfr,vfr,ie,nshift);
 557
 558         }
 559
 560                 if(bACF){
 561
 562     printf("\nCalculating current autocorrelation ... \n");
 563                         sgk=calc_cacf(outf,prefactorav/PICO,cacf,time,nvfr,vfr,ie,nshift);
 564
 565                         if (ie>ii){
 566
 567                                 snew(xfit,ie-ii+1);
 568                                 snew(yfit,ie-ii+1);
 569
 570                                 for(i=ii;i<=ie;i++){
 571
 572                                         xfit[i-ii]=log(time[vfr[i]]);
 573                                         rtmp=fabs(cacf[i]);
 574                                         yfit[i-ii]=log(rtmp);
 575
 576   }
 577
 578                                 lsq_y_ax_b(ie-ii,xfit,yfit,&sigma,&malt,&err,&rest);
 579
 580                                 malt=exp(malt);
 581
 582                                 sigma+=1.0;
 583
 584                                 malt*=prefactorav*2.0e12/sigma;
 585
 586                                 sfree(xfit);
 587                                 sfree(yfit);
 588
 589                         }
 590                 }
 591   }
 592
 593
 594   /* Calculation of the dielectric constant */
 595
 596   fprintf(stderr,"\n********************************************\n");
 597         dk_s=calceps(prefactor,md2,mj2,mjd,eps_rf,FALSE);
 598         fprintf(stderr,"\nAbsolute values:\n epsilon=%f\n",dk_s);
 599         fprintf(stderr," <M_D^2> , <M_J^2>, <(M_J*M_D)^2>:  (%f, %f, %f)\n\n",md2,mj2,mjd);
 600         fprintf(stderr,"********************************************\n");
 601
 602
 603         dk_f=calceps(prefactor,md2-mdav2,mj2-mj,mjd-mjdav,eps_rf,FALSE);
 604         fprintf(stderr,"\n\nFluctuations:\n epsilon=%f\n\n",dk_f);
 605         fprintf(stderr,"\n deltaM_D , deltaM_J, deltaM_JD:  (%f, %f, %f)\n",md2-mdav2,mj2-mj,mjd-mjdav);
 606         fprintf(stderr,"\n********************************************\n");
 607         if (bINT){
 608                 dk_d=calceps(prefactor,md2-mdav2,mj2-mj,corint,eps_rf,TRUE);
 609                 fprintf(stderr,"\nStatic dielectric constant using integral and fluctuations: %f\n",dk_d);
 610                 fprintf(stderr,"\n < M_JM_D > via integral:  %.3f\n",-1.0*corint);
 611         }
 612
 613         fprintf(stderr,"\n***************************************************");
 614         fprintf(stderr,"\n\nAverage volume V=%f nm^3 at T=%f K\n",volume_av,temp);
 615         fprintf(stderr,"and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n",prefactor);
 616
 617
 618
 619         if(bACF){
 620                 fprintf(stderr,"Integral and integrated fit to the current acf yields at t=%f:\n",time[vfr[ii]]);
 621                 fprintf(stderr,"sigma=%8.3f (pure integral: %.3f)\n",sgk-malt*pow(time[vfr[ii]],sigma),sgk);
 622         }
 623
 624         if (ei>bi){
 625                         fprintf(stderr,"\nStart fit at %f ps (%f).\n",time[bi],bfit);
 626                         fprintf(stderr,"End fit at %f ps (%f).\n\n",time[ei],efit);
 627
 628                         snew(xfit,ei-bi+1);
 629                         snew(yfit,ei-bi+1);
 630
 631                         for(i=bi;i<=ei;i++){
 632                                 xfit[i-bi]=time[i];
 633                                 yfit[i-bi]=dsp2[i];
 634                         }
 635
 636                         lsq_y_ax_b(ei-bi,xfit,yfit,&sigma,&malt,&err,&rest);
 637
 638                         sigma*=1e12;
 639                         dk_d=calceps(prefactor,md2,0.5*malt/prefactorav,corint,eps_rf,TRUE);
 640
 641                         fprintf(stderr,"Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
 642                         fprintf(stderr,"sigma=%.4f\n",sigma);
 643                         fprintf(stderr,"translational fraction of M^2: %.4f\n",0.5*malt/prefactorav);
 644                         fprintf(stderr,"Dielectric constant using EH: %.4f\n",dk_d);
 645
 646                         sfree(xfit);
 647                         sfree(yfit);
 648   }
 649   else{
 650                         fprintf(stderr,"Too less points for a fit.\n");
 651   }
 652
 653
 654         if (v0!=NULL)
 655     sfree(v0);
 656         if(bACF)
 657                 sfree(cacf);
 658         if(bINT)
 659                 sfree(djc);
 660
 661   sfree(time);
 662
 663
 664         sfree(mjdsp);
 665   sfree(mu);
 666 }
 667
 668
 669
 670 int gmx_current(int argc,char *argv[])
 671 {
 672
 673   static int  nshift=1000;
 674   static real temp=300.0;
 675   static real trust=0.25;
 676   static real eps_rf=0.0;
 677   static gmx_bool bNoJump=TRUE;
 678   static real bfit=100.0;
 679   static real bvit=0.5;
 680   static real efit=400.0;
 681   static real evit=5.0;
 682   t_pargs pa[] = {
 683     { "-sh", FALSE, etINT, {&nshift},
 684       "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
 685     { "-nojump", FALSE, etBOOL, {&bNoJump},
 686       "Removes jumps of atoms across the box."},
 687     { "-eps", FALSE, etREAL, {&eps_rf},
 688                 "Dielectric constant of the surrounding medium. The value zero corresponds to infinity (tin-foil boundary conditions)."},
 689         { "-bfit", FALSE, etREAL, {&bfit},
 690                 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
 691         { "-efit", FALSE, etREAL, {&efit},
 692                 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
 693         { "-bvit", FALSE, etREAL, {&bvit},
 694                                 "Begin of the fit of the current autocorrelation function to a*t^b."},
 695         { "-evit", FALSE, etREAL, {&evit},
 696                                 "End of the fit of the current autocorrelation function to a*t^b."},
 697     { "-tr", FALSE, etREAL, {&trust},
 698                 "Fraction of the trajectory taken into account for the integral."},
 699     { "-temp", FALSE, etREAL, {&temp},
 700                 "Temperature for calculating epsilon."
 701     }
 702   };
 703
 704   output_env_t oenv;
 705   t_topology top;
 706   char       title[STRLEN];
 707   char       **grpname=NULL;
 708   const char *indexfn;
 709   t_trxframe fr;
 710   real       *mass2=NULL;
 711   rvec       *xtop,*vtop;
 712   matrix     box;
 713   atom_id    *index0;
 714   int                                   *indexm=NULL;
 715   int        isize;
 716   t_trxstatus *status;
 717   int        flags = 0;
 718   gmx_bool           bTop;
 719   gmx_bool               bNEU;
 720   gmx_bool               bACF;
 721   gmx_bool               bINT;
 722   int        ePBC=-1;
 723   int        natoms;
 724   int                   nmols;
 725   int        i,j,k=0,l;
 726   int        step;
 727   real       t;
 728   real       lambda;
 729   real                   *qmol;
 730   FILE       *outf=NULL;
 731   FILE       *outi=NULL;
 732   FILE       *tfile=NULL;
 733   FILE       *mcor=NULL;
 734   FILE       *fmj=NULL;
 735   FILE       *fmd=NULL;
 736   FILE           *fmjdsp=NULL;
 737   FILE          *fcur=NULL;
 738   t_filenm fnm[] = {
 739     { efTPS,  NULL,  NULL, ffREAD },   /* this is for the topology */
 740     { efNDX, NULL, NULL, ffOPTRD },
 741     { efTRX, "-f", NULL, ffREAD },     /* and this for the trajectory */
 742     { efXVG, "-o", "current.xvg", ffWRITE },
 743     { efXVG, "-caf", "caf.xvg", ffOPTWR },
 744     { efXVG, "-dsp", "dsp.xvg", ffWRITE },
 745     { efXVG, "-md", "md.xvg", ffWRITE },
 746     { efXVG, "-mj", "mj.xvg", ffWRITE},
 747     { efXVG, "-mc", "mc.xvg", ffOPTWR }
 748   };
 749
 750 #define NFILE asize(fnm)
 751
 752
 753     const char *desc[] = {
 754     "[TT]g_current[tt] is a tool for calculating the current autocorrelation function, the correlation",
 755     "of the rotational and translational dipole moment of the system, and the resulting static",
 756     "dielectric constant. To obtain a reasonable result, the index group has to be neutral.",
 757     "Furthermore, the routine is capable of extracting the static conductivity from the current ",
 758     "autocorrelation function, if velocities are given. Additionally, an Einstein-Helfand fit ",
 759     "can be used to obtain the static conductivity."
 760     "[PAR]",
 761     "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
 762     "correlation of the rotational and translational part of the dipole moment in the corresponding",
 763     "file. However, this option is only available for trajectories containing velocities.",
 764     "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
 765     "autocorrelation functions. Since averaging proceeds by shifting the starting point",
 766     "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
 767     "starting points. Towards the end, statistical inaccuracy grows and integrating the",
 768     "correlation function only yields reliable values until a certain point, depending on",
 769     "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
 770     "for calculating the static dielectric constant.",
 771     "[PAR]",
 772     "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
 773     "[PAR]",
 774     "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
 775     "a Reaction Field or dipole corrections of the Ewald summation ([TT]-eps[tt]=0 corresponds to",
 776     "tin-foil boundary conditions).",
 777     "[PAR]",
 778     "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
 779     "translational dipole moment, required for the Einstein-Helfand fit. The results from the fit allow",
 780     "the determination of the dielectric constant for system of charged molecules. However, it is also possible to extract",
 781     "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
 782     "option has to be used with care, since only very short time spans fulfill the approximation that the density",
 783     "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
 784     "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
 785     "the translational part of the dielectric constant."
 786   };
 787
 788
 789   /* At first the arguments will be parsed and the system information processed */
 790
 791   CopyRight(stderr,argv[0]);
 792
 793   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
 794                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 795
 796   bACF = opt2bSet("-caf",NFILE,fnm);
 797   bINT = opt2bSet("-mc",NFILE,fnm);
 798
 799   bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,&vtop,box,TRUE);
 800
 801
 802
 803   sfree(xtop);
 804   sfree(vtop);
 805   indexfn = ftp2fn_null(efNDX,NFILE,fnm);
 806   snew(grpname,1);
 807
 808   get_index(&(top.atoms),indexfn,1,&isize,&index0,grpname);
 809
 810   flags = flags | TRX_READ_X | TRX_READ_V;
 811
 812   read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
 813
 814   snew(mass2,top.atoms.nr);
 815   snew(qmol,top.atoms.nr);
 816
 817   bNEU=precalc(top,mass2,qmol);
 818
 819
 820   snew(indexm,isize);
 821
 822   for(i=0;i<isize;i++)
 823         indexm[i]=index0[i];
 824
 825   nmols=isize;
 826
 827
 828   index_atom2mol(&nmols,indexm,&top.mols);
 829
 830   if (fr.bV){
 831       if(bACF){
 832           outf =xvgropen(opt2fn("-caf",NFILE,fnm),
 833                   "Current autocorrelation function",output_env_get_xvgr_tlabel(oenv),
 834                   "ACF (e nm/ps)\\S2",oenv);
 835           fprintf(outf,"# time\t acf\t average \t std.dev\n");
 836       }
 837       fcur =xvgropen(opt2fn("-o",NFILE,fnm),
 838               "Current",output_env_get_xvgr_tlabel(oenv),"J(t) (e nm/ps)",oenv);
 839       fprintf(fcur,"# time\t Jx\t Jy \t J_z \n");
 840       if(bINT){
 841           mcor = xvgropen(opt2fn("-mc",NFILE,fnm),
 842                   "M\\sD\\N - current  autocorrelation function",
 843                   output_env_get_xvgr_tlabel(oenv),
 844                   "< M\\sD\\N (0)\\c7\\CJ(t) >  (e nm/ps)\\S2",oenv);
 845           fprintf(mcor,"# time\t M_D(0) J(t) acf \t Integral acf\n");
 846     }
 847   }
 848
 849   fmj = xvgropen(opt2fn("-mj",NFILE,fnm),
 850                  "Averaged translational part of M",output_env_get_xvgr_tlabel(oenv),
 851                  "< M\\sJ\\N > (enm)",oenv);
 852   fprintf(fmj,"# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
 853   fmd = xvgropen(opt2fn("-md",NFILE,fnm),
 854                  "Averaged rotational part of M",output_env_get_xvgr_tlabel(oenv),
 855                  "< M\\sD\\N > (enm)",oenv);
 856   fprintf(fmd,"# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
 857
 858         fmjdsp = xvgropen(opt2fn("-dsp",NFILE,fnm),
 859                  "MSD of the squared translational dipole moment M",
 860                  output_env_get_xvgr_tlabel(oenv),
 861                  "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
 862                  oenv);
 863
 864
 865   /* System information is read and prepared, dielectric() processes the frames
 866    * and calculates the requested quantities */
 867
 868   dielectric(fmj,fmd,outf,fcur,mcor,fmjdsp,bNoJump,bACF,bINT,ePBC,top,fr,
 869              temp,trust, bfit,efit,bvit,evit,status,isize,nmols,nshift,
 870              index0,indexm,mass2,qmol,eps_rf,oenv);
 871
 872   ffclose(fmj);
 873   ffclose(fmd);
 874   ffclose(fmjdsp);
 875   if (bACF)
 876     ffclose(outf);
 877   if(bINT)
 878       ffclose(mcor);
 879
 880   thanx(stderr);
 881
 882   return 0;
 883 }
 884