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                        t_trxstatus *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   gmx_rmpbc_t  gpbc=NULL;
 348
 349   /*
 350    * indices for EH
 351    */
 352
 353   ei=0;
 354   bi=0;
 355
 356   /*
 357    * indices for GK
 358    */
 359
 360   ii=0;
 361   ie=0;
 362   t0=0;
 363   sgk=0.0;
 364
 365
 366   /* This is the main loop over frames */
 367
 368
 369   nfr = 0;
 370
 371
 372   nvfr = 0;
 373   vshfr = 0;
 374   nalloc = 0;
 375   valloc = 0;
 376
 377   clear_rvec(mja_tmp);
 378   clear_rvec(mjd_tmp);
 379   clear_rvec(mdvec);
 380   clear_rvec(tmp);
 381   gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box);
 382
 383   do{
 384
 385     refr=(real)(nfr+1);
 386
 387     if(nfr >= nalloc){
 388       nalloc+=100;
 389       srenew(time,nalloc);
 390       srenew(mu,nalloc);
 391       srenew(mjdsp,nalloc);
 392       srenew(dsp2,nalloc);
 393       srenew(mtrans,nalloc);
 394       srenew(xshfr,nalloc);
 395
 396       for(i=nfr;i<nalloc;i++){
 397         clear_rvec(mjdsp[i]);
 398         clear_rvec(mu[i]);
 399         clear_rvec(mtrans[i]);
 400         dsp2[i]=0.0;
 401         xshfr[i]=0.0;
 402       }
 403     }
 404
 405
 406
 407     if (nfr==0){
 408       t0=fr.time;
 409
 410     }
 411
 412     time[nfr]=fr.time-t0;
 413
 414     if(time[nfr]<=bfit)
 415         bi=nfr;
 416                 if(time[nfr]<=efit)
 417                         ei=nfr;
 418
 419                 if(bNoJump){
 420
 421     if (xp)
 422       remove_jump(fr.box,fr.natoms,xp,fr.x);
 423                         else
 424       snew(xp,fr.natoms);
 425
 426                         for(i=0; i<fr.natoms;i++)
 427       copy_rvec(fr.x[i],xp[i]);
 428
 429     }
 430
 431                 gmx_rmpbc_trxfr(gpbc,&fr);
 432
 433                 calc_mj(top,ePBC,fr.box,bNoJump,nmols,indexm,fr.x,mtrans[nfr],mass2,qmol);
 434
 435     for(i=0;i<isize;i++){
 436       j=index0[i];
 437       svmul(top.atoms.atom[j].q,fr.x[j],fr.x[j]);
 438       rvec_inc(mu[nfr],fr.x[j]);
 439     }
 440
 441     /*if(mod(nfr,nshift)==0){*/
 442                 if(nfr%nshift==0){
 443                         for(j=nfr;j>=0;j--){
 444                                         rvec_sub(mtrans[nfr],mtrans[j],tmp);
 445                                         dsp2[nfr-j]+=norm2(tmp);
 446                                         xshfr[nfr-j]+=1.0;
 447                         }
 448                 }
 449
 450                 if (fr.bV){
 451         if(nvfr >= valloc){
 452                 valloc+=100;
 453                 srenew(vfr,valloc);
 454                 if(bINT)
 455                         srenew(djc,valloc);
 456                 srenew(v0,valloc);
 457                 if(bACF)
 458                         srenew(cacf,valloc);
 459         }
 460         if(time[nfr]<=bvit)
 461                 ii=nvfr;
 462         if(time[nfr]<=evit)
 463                                 ie=nvfr;
 464         vfr[nvfr]=nfr;
 465         clear_rvec(v0[nvfr]);
 466         if(bACF)
 467                 cacf[nvfr]=0.0;
 468         if(bINT)
 469                 djc[nvfr]=0.0;
 470         for(i=0;i<isize;i++){
 471                 j=index0[i];
 472                 svmul(mass2[j],fr.v[j],fr.v[j]);
 473                 svmul(qmol[j],fr.v[j],fr.v[j]);
 474                 rvec_inc(v0[nvfr],fr.v[j]);
 475       }
 476
 477         fprintf(fcur,"%.3f\t%.6f\t%.6f\t%.6f\n",time[nfr],v0[nfr][XX],v0[nfr][YY],v0[nfr][ZZ]);
 478         if(bACF || bINT)
 479                 /*if(mod(nvfr,nshift)==0){*/
 480                                 if(nvfr%nshift==0){
 481                         for(j=nvfr;j>=0;j--){
 482                                 if(bACF)
 483                                         cacf[nvfr-j]+=iprod(v0[nvfr],v0[j]);
 484                                 if(bINT)
 485                                         djc[nvfr-j]+=iprod(mu[vfr[j]],v0[nvfr]);
 486                         }
 487                         vshfr++;
 488                 }
 489                         nvfr++;
 490     }
 491
 492     volume = det(fr.box);
 493     volume_av += volume;
 494
 495     rvec_inc(mja_tmp,mtrans[nfr]);
 496     mjd+=iprod(mu[nfr],mtrans[nfr]);
 497     rvec_inc(mdvec,mu[nfr]);
 498
 499     mj2+=iprod(mtrans[nfr],mtrans[nfr]);
 500     md2+=iprod(mu[nfr],mu[nfr]);
 501
 502     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);
 503     fprintf(fmd,"%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n",    \
 504             time[nfr],mu[nfr][XX],mu[nfr][YY],mu[nfr][ZZ],md2/refr,norm(mdvec)/refr);
 505
 506     nfr++;
 507
 508   }while(read_next_frame(oenv,status,&fr));
 509
 510   gmx_rmpbc_done(gpbc);
 511
 512   volume_av/=refr;
 513
 514         prefactor=1.0;
 515         prefactor/=3.0*EPSILON0*volume_av*BOLTZ*temp;
 516
 517
 518         prefactorav=E_CHARGE*E_CHARGE;
 519         prefactorav/=volume_av*BOLTZMANN*temp*NANO*6.0;
 520
 521         fprintf(stderr,"Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n",prefactorav);
 522
 523         calc_mjdsp(fmjdsp,volume_av,temp,prefactorav,mjdsp,dsp2,time,nfr,xshfr);
 524
 525   /*
 526    * Now we can average and calculate the correlation functions
 527    */
 528
 529
 530   mj2/=refr;
 531   mjd/=refr;
 532   md2/=refr;
 533
 534   svmul(1.0/refr,mdvec,mdvec);
 535   svmul(1.0/refr,mja_tmp,mja_tmp);
 536
 537   mdav2=norm2(mdvec);
 538   mj=norm2(mja_tmp);
 539   mjdav=iprod(mdvec,mja_tmp);
 540
 541
 542   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);
 543   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);
 544
 545   if (v0!=NULL){
 546         trust*=(real) nvfr;
 547                 itrust=(int) trust;
 548
 549         if(bINT){
 550
 551                 printf("\nCalculating M_D - current correlation integral ... \n");
 552                 corint=calc_cacf(mcor,prefactorav/EPSI0,djc,time,nvfr,vfr,ie,nshift);
 553
 554         }
 555
 556                 if(bACF){
 557
 558     printf("\nCalculating current autocorrelation ... \n");
 559                         sgk=calc_cacf(outf,prefactorav/PICO,cacf,time,nvfr,vfr,ie,nshift);
 560
 561                         if (ie>ii){
 562
 563                                 snew(xfit,ie-ii+1);
 564                                 snew(yfit,ie-ii+1);
 565
 566                                 for(i=ii;i<=ie;i++){
 567
 568                                         xfit[i-ii]=log(time[vfr[i]]);
 569                                         rtmp=fabs(cacf[i]);
 570                                         yfit[i-ii]=log(rtmp);
 571
 572   }
 573
 574                                 lsq_y_ax_b(ie-ii,xfit,yfit,&sigma,&malt,&err,&rest);
 575
 576                                 malt=exp(malt);
 577
 578                                 sigma+=1.0;
 579
 580                                 malt*=prefactorav*2.0e12/sigma;
 581
 582                                 sfree(xfit);
 583                                 sfree(yfit);
 584
 585                         }
 586                 }
 587   }
 588
 589
 590   /* Calculation of the dielectric constant */
 591
 592   fprintf(stderr,"\n********************************************\n");
 593         dk_s=calceps(prefactor,md2,mj2,mjd,eps_rf,FALSE);
 594         fprintf(stderr,"\nAbsolute values:\n epsilon=%f\n",dk_s);
 595         fprintf(stderr," <M_D^2> , <M_J^2>, <(M_J*M_D)^2>:  (%f, %f, %f)\n\n",md2,mj2,mjd);
 596         fprintf(stderr,"********************************************\n");
 597
 598
 599         dk_f=calceps(prefactor,md2-mdav2,mj2-mj,mjd-mjdav,eps_rf,FALSE);
 600         fprintf(stderr,"\n\nFluctuations:\n epsilon=%f\n\n",dk_f);
 601         fprintf(stderr,"\n deltaM_D , deltaM_J, deltaM_JD:  (%f, %f, %f)\n",md2-mdav2,mj2-mj,mjd-mjdav);
 602         fprintf(stderr,"\n********************************************\n");
 603         if (bINT){
 604                 dk_d=calceps(prefactor,md2-mdav2,mj2-mj,corint,eps_rf,TRUE);
 605                 fprintf(stderr,"\nStatic dielectric constant using integral and fluctuations: %f\n",dk_d);
 606                 fprintf(stderr,"\n < M_JM_D > via integral:  %.3f\n",-1.0*corint);
 607         }
 608
 609         fprintf(stderr,"\n***************************************************");
 610         fprintf(stderr,"\n\nAverage volume V=%f nm^3 at T=%f K\n",volume_av,temp);
 611         fprintf(stderr,"and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n",prefactor);
 612
 613
 614
 615         if(bACF){
 616                 fprintf(stderr,"Integral and integrated fit to the current acf yields at t=%f:\n",time[vfr[ii]]);
 617                 fprintf(stderr,"sigma=%8.3f (pure integral: %.3f)\n",sgk-malt*pow(time[vfr[ii]],sigma),sgk);
 618         }
 619
 620         if (ei>bi){
 621                         fprintf(stderr,"\nStart fit at %f ps (%f).\n",time[bi],bfit);
 622                         fprintf(stderr,"End fit at %f ps (%f).\n\n",time[ei],efit);
 623
 624                         snew(xfit,ei-bi+1);
 625                         snew(yfit,ei-bi+1);
 626
 627                         for(i=bi;i<=ei;i++){
 628                                 xfit[i-bi]=time[i];
 629                                 yfit[i-bi]=dsp2[i];
 630                         }
 631
 632                         lsq_y_ax_b(ei-bi,xfit,yfit,&sigma,&malt,&err,&rest);
 633
 634                         sigma*=1e12;
 635                         dk_d=calceps(prefactor,md2,0.5*malt/prefactorav,corint,eps_rf,TRUE);
 636
 637                         fprintf(stderr,"Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n");
 638                         fprintf(stderr,"sigma=%.4f\n",sigma);
 639                         fprintf(stderr,"translational fraction of M^2: %.4f\n",0.5*malt/prefactorav);
 640                         fprintf(stderr,"Dielectric constant using EH: %.4f\n",dk_d);
 641
 642                         sfree(xfit);
 643                         sfree(yfit);
 644   }
 645   else{
 646                         fprintf(stderr,"Too less points for a fit.\n");
 647   }
 648
 649
 650         if (v0!=NULL)
 651     sfree(v0);
 652         if(bACF)
 653                 sfree(cacf);
 654         if(bINT)
 655                 sfree(djc);
 656
 657   sfree(time);
 658
 659
 660         sfree(mjdsp);
 661   sfree(mu);
 662 }
 663
 664
 665
 666 int gmx_current(int argc,char *argv[])
 667 {
 668
 669   static int  nshift=1000;
 670   static real temp=300.0;
 671   static real trust=0.25;
 672   static real eps_rf=0.0;
 673   static bool bNoJump=TRUE;
 674   static real bfit=100.0;
 675   static real bvit=0.5;
 676   static real efit=400.0;
 677   static real evit=5.0;
 678   t_pargs pa[] = {
 679     { "-sh", FALSE, etINT, {&nshift},
 680       "Shift of the frames for averaging the correlation functions and the mean-square displacement."},
 681     { "-nojump", FALSE, etBOOL, {&bNoJump},
 682       "Removes jumps of atoms across the box."},
 683     { "-eps", FALSE, etREAL, {&eps_rf},
 684                 "Dielectric constant of the surrounding medium. eps=0.0 corresponds to eps=infinity (thinfoil boundary conditions)."},
 685         { "-bfit", FALSE, etREAL, {&bfit},
 686                 "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
 687         { "-efit", FALSE, etREAL, {&efit},
 688                 "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."},
 689         { "-bvit", FALSE, etREAL, {&bvit},
 690                                 "Begin of the fit of the current autocorrelation function to a*t^b."},
 691         { "-evit", FALSE, etREAL, {&evit},
 692                                 "End of the fit of the current autocorrelation function to a*t^b."},
 693     { "-tr", FALSE, etREAL, {&trust},
 694                 "Fraction of the trajectory taken into account for the integral."},
 695     { "-temp", FALSE, etREAL, {&temp},
 696                 "Temperature for calculating epsilon."
 697     }
 698   };
 699
 700   output_env_t oenv;
 701   t_topology top;
 702   char       title[STRLEN];
 703   char       **grpname=NULL;
 704   const char *indexfn;
 705   t_trxframe fr;
 706   real       *mass2=NULL;
 707   rvec       *xtop,*vtop;
 708   matrix     box;
 709   atom_id    *index0=NULL;
 710   int                                   *indexm=NULL;
 711   int        isize;
 712   t_trxstatus *status;
 713   int        flags = 0;
 714   bool       bTop;
 715   bool           bNEU;
 716   bool           bACF;
 717   bool           bINT;
 718   int        ePBC=-1;
 719   int        natoms;
 720   int                   nmols;
 721   int        i,j,k=0,l;
 722   int        step;
 723   real       t;
 724   real       lambda;
 725   real                   *qmol;
 726   FILE       *outf=NULL;
 727   FILE       *outi=NULL;
 728   FILE       *tfile=NULL;
 729   FILE       *mcor=NULL;
 730   FILE       *fmj=NULL;
 731   FILE       *fmd=NULL;
 732   FILE           *fmjdsp=NULL;
 733   FILE          *fcur=NULL;
 734   t_filenm fnm[] = {
 735     { efTPS,  NULL,  NULL, ffREAD },   /* this is for the topology */
 736     { efNDX, NULL, NULL, ffOPTRD },
 737     { efTRX, "-f", NULL, ffREAD },     /* and this for the trajectory */
 738     { efXVG, "-o", "current.xvg", ffWRITE },
 739     { efXVG, "-caf", "caf.xvg", ffOPTWR },
 740     { efXVG, "-dsp", "dsp.xvg", ffWRITE },
 741     { efXVG, "-md", "md.xvg", ffWRITE },
 742     { efXVG, "-mj", "mj.xvg", ffWRITE},
 743     { efXVG, "-mc", "mc.xvg", ffOPTWR }
 744   };
 745
 746 #define NFILE asize(fnm)
 747
 748
 749     const char *desc[] = {
 750     "This is a tool for calculating the current autocorrelation function, the correlation",
 751     "of the rotational and translational dipole moment of the system, and the resulting static",
 752     "dielectric constant. To obtain a reasonable result the index group has to be neutral.",
 753     "Furthermore the routine is capable of extracting the static conductivity from the current ",
 754     "autocorrelation function, if velocities are given. Additionally an Einstein-Helfand fit also",
 755     "allows to get the static conductivity."
 756     "[PAR]",
 757     "The flag [TT]-caf[tt] is for the output of the current autocorrelation function and [TT]-mc[tt] writes the",
 758     "correlation of the rotational and translational part of the dipole moment in the corresponding",
 759     "file. However this option is only available for trajectories containing velocities.",
 760     "Options [TT]-sh[tt] and [TT]-tr[tt] are responsible for the averaging and integration of the",
 761     "autocorrelation functions. Since averaging proceeds by shifting the starting point",
 762     "through the trajectory, the shift can be modified with [TT]-sh[tt] to enable the choice of uncorrelated",
 763     "starting points. Towards the end, statistical inaccuracy grows and integrating the",
 764     "correlation function only yields reliable values until a certain point, depending on",
 765     "the number of frames. The option [TT]-tr[tt] controls the region of the integral taken into account",
 766     "for calculating the static dielectric constant.",
 767     "[PAR]",
 768     "Option [TT]-temp[tt] sets the temperature required for the computation of the static dielectric constant.",
 769     "[PAR]",
 770     "Option [TT]-eps[tt] controls the dielectric constant of the surrounding medium for simulations using",
 771     "a Reaction Field or dipole corrections of the Ewald summation (eps=0 corresponds to",
 772     "tin-foil boundary conditions).",
 773     "[PAR]",
 774     "[TT]-[no]nojump[tt] unfolds the coordinates to allow free diffusion. This is required to get a continuous",
 775     "translational dipole moment, required for the Einstein-Helfand fit. The resuls from the fit allow to",
 776     "determine the dielectric constant for system of charged molecules. However it is also possible to extract",
 777     "the dielectric constant from the fluctuations of the total dipole moment in folded coordinates. But this",
 778     "options has to be used with care, since only very short time spans fulfill the approximation, that the density",
 779     "of the molecules is approximately constant and the averages are already converged. To be on the safe side,",
 780     "the dielectric constant should be calculated with the help of the Einstein-Helfand method for",
 781     "the translational part of the dielectric constant."
 782   };
 783
 784
 785   /* At first the arguments will be parsed and the system information processed */
 786
 787   CopyRight(stderr,argv[0]);
 788
 789   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,
 790                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 791
 792   bACF = opt2bSet("-caf",NFILE,fnm);
 793   bINT = opt2bSet("-mc",NFILE,fnm);
 794
 795   bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,&vtop,box,TRUE);
 796
 797
 798
 799   sfree(xtop);
 800   sfree(vtop);
 801   indexfn = ftp2fn_null(efNDX,NFILE,fnm);
 802   snew(grpname,1);
 803
 804   get_index(&(top.atoms),indexfn,1,&isize,&index0,grpname);
 805
 806   flags = flags | TRX_READ_X | TRX_READ_V;
 807
 808   read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
 809
 810   snew(mass2,top.atoms.nr);
 811   snew(qmol,top.atoms.nr);
 812
 813   bNEU=precalc(top,mass2,qmol);
 814
 815
 816   snew(indexm,isize);
 817
 818   for(i=0;i<isize;i++)
 819         indexm[i]=index0[i];
 820
 821   nmols=isize;
 822
 823
 824   index_atom2mol(&nmols,indexm,&top.mols);
 825
 826   if (fr.bV){
 827       if(bACF){
 828           outf =xvgropen(opt2fn("-caf",NFILE,fnm),
 829                   "Current autocorrelation function",output_env_get_xvgr_tlabel(oenv),
 830                   "ACF (e nm/ps)\\S2",oenv);
 831           fprintf(outf,"# time\t acf\t average \t std.dev\n");
 832       }
 833       fcur =xvgropen(opt2fn("-o",NFILE,fnm),
 834               "Current",output_env_get_xvgr_tlabel(oenv),"J(t) (e nm/ps)",oenv);
 835       fprintf(fcur,"# time\t Jx\t Jy \t J_z \n");
 836       if(bINT){
 837           mcor = xvgropen(opt2fn("-mc",NFILE,fnm),
 838                   "M\\sD\\N - current  autocorrelation function",
 839                   output_env_get_xvgr_tlabel(oenv),
 840                   "< M\\sD\\N (0)\\c7\\CJ(t) >  (e nm/ps)\\S2",oenv);
 841           fprintf(mcor,"# time\t M_D(0) J(t) acf \t Integral acf\n");
 842     }
 843   }
 844
 845   fmj = xvgropen(opt2fn("-mj",NFILE,fnm),
 846                  "Averaged translational part of M",output_env_get_xvgr_tlabel(oenv),
 847                  "< M\\sJ\\N > (enm)",oenv);
 848   fprintf(fmj,"# time\t x\t y \t z \t average of M_J^2 \t std.dev\n");
 849   fmd = xvgropen(opt2fn("-md",NFILE,fnm),
 850                  "Averaged rotational part of M",output_env_get_xvgr_tlabel(oenv),
 851                  "< M\\sD\\N > (enm)",oenv);
 852   fprintf(fmd,"# time\t x\t y \t z \t average of M_D^2 \t std.dev\n");
 853
 854         fmjdsp = xvgropen(opt2fn("-dsp",NFILE,fnm),
 855                  "MSD of the squared translational dipole moment M",
 856                  output_env_get_xvgr_tlabel(oenv),
 857                  "<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N",
 858                  oenv);
 859
 860
 861   /* System information is read and prepared, dielectric() processes the frames
 862    * and calculates the requested quantities */
 863
 864   dielectric(fmj,fmd,outf,fcur,mcor,fmjdsp,bNoJump,bACF,bINT,ePBC,top,fr,
 865              temp,trust, bfit,efit,bvit,evit,status,isize,nmols,nshift,
 866              index0,indexm,mass2,qmol,eps_rf,oenv);
 867
 868   ffclose(fmj);
 869   ffclose(fmd);
 870   ffclose(fmjdsp);
 871   if (bACF)
 872     ffclose(outf);
 873   if(bINT)
 874       ffclose(mcor);
 875
 876   thanx(stderr);
 877
 878   return 0;
 879 }
 880