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