/* $Id: template.c,v 1.4 2001/07/23 15:28:29 lindahl Exp $ * * * This source code is part of * * G R O M A C S * * GROningen MAchine for Chemical Simulations * VERSION 3.0 * * Copyright (c) 1991-2001 * BIOSON Research Institute, Dept. of Biophysical Chemistry * University of Groningen, The Netherlands * * This program is free software; you can redistribute it and/or * * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * If you want to redistribute modifications, please consider that * scientific software is very special. Version control is crucial - * bugs must be traceable. We will be happy to consider code for * inclusion in the official distribution, but derived work must not * be called official GROMACS. Details are found in the README & COPYING * files - if they are missing, get the official version at www.gromacs.org. * * To help us fund GROMACS development, we humbly ask that you cite * the papers on the package - you can find them in the top README file. * * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org . * * And Hey: * Gyas ROwers Mature At Cryogenic Speed * * finished FD 09/07/08 * */ #include "statutil.h" #include "typedefs.h" #include "smalloc.h" #include "vec.h" #include "copyrite.h" #include "statutil.h" #include "tpxio.h" #include "xvgr.h" #include "rmpbc.h" #include "pbc.h" #include "physics.h" #include "index.h" #define SQR(x) (pow(x,2.0)) #define EPSI0 (EPSILON0*E_CHARGE*E_CHARGE*AVOGADRO/(KILO*NANO)) /* EPSILON0 in SI units */ static bool precalc(t_topology top, t_trxframe fr,real mass2[],real qmol[],int cmol[]){ real mtot; real qtot; real qall; int i,j,k,l; int ai,ci; bool bNEU; ai=0; ci=0; qall=0.0; for(i=0;i0.0) ci++; } cmol[0]=ci; cmol[1]=ai; if(abs(qall)>0.01){ printf("\n\nSystem not neutral (q=%f) will not calculate translational part of the dipole moment.\n",qall); bNEU=FALSE; }else bNEU=TRUE; return bNEU; } static void calc_molcharge(t_topology top,real qmol[],int catin[],int anin[],int cmol[]){ int i; int ai,ci; ai=0; ci=0; for(i=0;i0.0) catin[ci]=i; ci++; } } static void calc_mj(t_topology top,int ePBC,matrix box,bool bNoJump,int isize,int index0[],rvec fr[], rvec mj,real mass2[],real qmol[],int ani[],int cati[],int cmol[]){ int i,j,k,l,d,m; rvec tmp,tmpp,tmpm; rvec center; t_pbc pbc; clear_rvec(tmp); calc_box_center(ecenterRECT,box,center); for(i=0;i=0; m--) { while (x[i][m]-xp[i][m] <= -hbox[m]) for(d=0; d<=m; d++) x[i][d] += box[m][d]; while (x[i][m]-xp[i][m] > hbox[m]) for(d=0; d<=m; d++) x[i][d] -= box[m][d]; } } static real calc_cacf(FILE *fcacf,real prefactor,real cacf[],real time[],int nfr,int vfr[],int ei,int nshift){ int i; real corint; real deltat=0.0; real rfr; real sigma=0.0; real sigma_ret=0.0; corint=0.0; if(nfr>1){ i=0; //for(i=0;i= nalloc){ nalloc+=100; srenew(time,nalloc); srenew(mu,nalloc); if (bDSP) srenew(mjdsp,nalloc); srenew(mtrans,nalloc); for(i=nfr;i=0;j--){ rvec_sub(mtrans[nfr],mtrans[j],tmp); mjdsp[nfr-j]+=norm2(tmp); } xshfr++; } } if (fr.bV){ if(nvfr >= valloc){ valloc+=100; srenew(vfr,valloc); if(bINT) srenew(djc,valloc); srenew(v0,valloc); if(bACF) srenew(cacf,valloc); } if(time[nfr]<=bvit) ii=nvfr; if(time[nfr]<=evit) ie=nvfr; vfr[nvfr]=nfr; clear_rvec(v0[nvfr]); if(bACF) cacf[nvfr]=0.0; if(bINT) djc[nvfr]=0.0; for(i=0;i=0;j--){ if(bACF) cacf[nvfr-j]+=iprod(v0[nvfr],v0[j]); if(bINT) djc[nvfr-j]+=iprod(mu[vfr[j]],v0[nvfr]); } vshfr++; } nvfr++; } volume = det(fr.box); volume_av += volume; rvec_inc(mja_tmp,mtrans[nfr]); mjd+=iprod(mu[nfr],mtrans[nfr]); rvec_inc(mdvec,mu[nfr]); mj2+=iprod(mtrans[nfr],mtrans[nfr]); md2+=iprod(mu[nfr],mu[nfr]); 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); fprintf(fmd,"%.3f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\t%8.5f\n", \ time[nfr],mu[nfr][XX],mu[nfr][YY],mu[nfr][ZZ],md2/refr,norm(mdvec)/refr); nfr++; }while(read_next_frame(status,&fr)); /* * Now we can average and calculate the correlation functions */ mj2/=refr; mjd/=refr; md2/=refr; svmul(1.0/refr,mdvec,mdvec); svmul(1.0/refr,mja_tmp,mja_tmp); mdav2=norm2(mdvec); mj=norm2(mja_tmp); mjdav=iprod(mdvec,mja_tmp); 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); 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); volume_av/=refr; //prefactor=DEBYE2ENM*DEBYE2ENM; prefactor=1.0/(3.0*EPSILON0*volume_av*BOLTZ*temp); //prefactorav=DEBYE2ENM*DEBYE2ENM*E_CHARGE*E_CHARGE; prefactorav=E_CHARGE*E_CHARGE; prefactorav/=volume_av*BOLTZMANN*temp*NANO*6.0; fprintf(stderr,"Prefactor fit E-H: 1 / 6.0*V*k_B*T: %g\n",prefactorav); if(bDSP) calc_mjdsp(fmjdsp,volume_av,temp,prefactorav,mjdsp,time,nfr,nshift); if (v0!=NULL){ trust*=(real) nvfr; itrust=(int) trust; if(bINT){ printf("\nCalculating M_D - current correlation integral ... \n"); corint=calc_cacf(mcor,prefactorav/EPSI0,djc,time,nvfr,vfr,ie,nshift); } if(bACF){ printf("\nCalculating current autocorrelation ... \n"); sgk=calc_cacf(outf,prefactorav/PICO,cacf,time,nvfr,vfr,ie,nshift); if (ie>ii){ snew(xfit,ie-ii+1); snew(yfit,ie-ii+1); for(i=ii;i<=ie;i++){ xfit[i-ii]=log(time[vfr[i]]); rtmp=fabs(cacf[i]); yfit[i-ii]=log(rtmp); } rest=lsq_y_ax_b(ie-ii,xfit,yfit,&sigma,&malt,&err); malt=exp(malt); sigma+=1.0; malt*=prefactorav*2.0e12/sigma; sfree(xfit); sfree(yfit); } } } /* Calculation of the dielectric constant */ fprintf(stderr,"\n********************************************\n"); dk_s=calceps(prefactor,md2,mj2,mjd,eps_rf,FALSE); fprintf(stderr,"\nAbsolute values:\n epsilon=%f\n",dk_s); fprintf(stderr," , , <(M_J*M_D)^2>: (%f, %f, %f)\n\n",md2,mj2,mjd); fprintf(stderr,"********************************************\n"); dk_f=calceps(prefactor,md2-mdav2,mj2-mj,mjd-mjdav,eps_rf,FALSE); fprintf(stderr,"\n\nFluctuations:\n epsilon=%f\n\n",dk_f); fprintf(stderr,"Rotational contribution to epsilon; %f\n" , prefactor*(md2-mdav2)); fprintf(stderr,"Translational contribution to epsilon; %f\n" , prefactor*(mj2-mj)); fprintf(stderr,"Correlation contribution to epsilon; %f\n" , 2.0*prefactor*(mjd-mjdav)); fprintf(stderr,"\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n",md2-mdav2,mj2-mj,mjd-mjdav); fprintf(stderr,"\n********************************************\n"); if (bINT){ dk_d=calceps(prefactor,md2-mdav2,mj2-mj,corint,eps_rf,TRUE); fprintf(stderr,"\nStatic dielectric constant using integral and fluctuations: %f\n",dk_d); fprintf(stderr,"\n deltaM_D , deltaM_J, deltaM_JD: (%f, %f, %f)\n",md2-mdav2,mj2-mj,corint); } fprintf(stderr,"\n***************************************************"); fprintf(stderr,"\n\nAverage volume V=%f nm^3 at T=%f K\n",volume_av,temp); fprintf(stderr,"and corresponding refactor 1.0 / 3.0*V*k_B*T*EPSILON_0: %f \n",prefactor); if(bACF){ fprintf(stderr,"Integral and integrated fit to the current acf yields at t=%f:\n",time[vfr[ii]]); fprintf(stderr,"sigma=%8.3f (pure integral: %.3f)\n",sgk-malt*pow(time[vfr[ii]],sigma),sgk); } if(bDSP){ if (ei>bi){ fprintf(stderr,"\nStart fit at %f ps (%f).\n",time[bi],bfit); fprintf(stderr,"End fit at %f ps (%f).\n\n",time[ei],efit); snew(xfit,ei-bi+1); snew(yfit,ei-bi+1); for(i=bi;i<=ei;i++){ xfit[i-bi]=time[i]; yfit[i-bi]=mjdsp[i]; } rest=lsq_y_ax_b(ei-bi,xfit,yfit,&sigma,&malt,&err); sigma*=1e12; fprintf(stderr,"Einstein-Helfand fit to the MSD of the translational dipole moment yields:\n\n"); fprintf(stderr,"sigma=%.4f\n",sigma); fprintf(stderr,"translational fraction of M^2: %.4f\n",0.5*malt/prefactorav); sfree(xfit); sfree(yfit); } else{ fprintf(stderr,"Too less points for a fit.\n"); } } if (v0!=NULL) sfree(v0); if(bACF) sfree(cacf); if(bINT) sfree(djc); sfree(time); if(bDSP) sfree(mjdsp); sfree(mu); } int gmx_current(int argc,char *argv[]) { static int nshift=1000; static real temp=300.0; static real trust=0.25; static real eps_rf=0.0; static bool bNoJump=FALSE; static real bfit=100.0; static real bvit=0.5; static real efit=400.0; static real evit=5.0; t_pargs pa[] = { { "-sh", FALSE, etINT, {&nshift}, "Shift of the frames for averaging the correlation functions and the mean-square displacement."}, { "-nojump", FALSE, etBOOL, {&bNoJump}, "Removes jumps of atoms across the box."}, { "-eps", FALSE, etREAL, {&eps_rf}, "Dielectric constant of the surrounding medium. eps=0.0 corresponds to eps=infinity (thinfoil boundary conditions)."}, { "-bfit", FALSE, etREAL, {&bfit}, "Begin of the fit of the straight line to the MSD of the translational fraction of the dipole moment."}, { "-efit", FALSE, etREAL, {&efit}, "End of the fit of the straight line to the MSD of the translational fraction of the dipole moment."}, { "-bvit", FALSE, etREAL, {&bvit}, "Begin of the fit of the current autocorrelation function to a*t^b."}, { "-evit", FALSE, etREAL, {&evit}, "End of the fit of the current autocorrelation function to a*t^b."}, { "-tr", FALSE, etREAL, {&trust}, "Fraction of the trajectory taken into account for the integral."}, { "-temp", FALSE, etREAL, {&temp}, "Temperature for calculating epsilon." } }; t_topology top; char title[STRLEN]; char **grpname=NULL; char *indexfn; t_trxframe fr; real *mass2=NULL; rvec *xtop,*vtop; matrix box; atom_id *index0=NULL; int isize; int status; int flags = 0; int *cmol; bool bTop; bool bNEU; bool bDSP; bool bACF; bool bINT; int ePBC=-1; int natoms; int i,j,k=0,l; int step; int *ani,*cati; real t; real lambda; real *qmol; FILE *outf=NULL; FILE *outi=NULL; FILE *tfile=NULL; FILE *mcor=NULL; FILE *fmj=NULL; FILE *fmd=NULL; FILE *fmjdsp=NULL; FILE *fcur=NULL; t_filenm fnm[] = { { efTPS, NULL, NULL, ffREAD }, /* this is for the topology */ { efNDX, NULL, NULL, ffOPTRD }, { efTRX, "-f", NULL, ffREAD }, /* and this for the trajectory */ { efXVG, "-o", "current.xvg", ffWRITE }, { efXVG, "-caf", "caf.xvg", ffOPTWR }, { efXVG, "-dsp", "dsp.xvg", ffOPTWR }, { efXVG, "-md", "md.xvg", ffWRITE }, { efXVG, "-mj", "mj.xvg", ffWRITE}, { efXVG, "-mc", "mc.xvg", ffOPTWR } }; #define NFILE asize(fnm) static char *desc[] = { "This is a small tool for calculating the current autocorrelation function, the correlation", "of the rotational and translational dipole moment of the system, and the resulting static", "dielectric constant. To obtain a reasonable result the index group has to be neutral.", "[PAR]", "Options -sh and -tr are responsible for the averaging and integration of the", "autocorrelation functions. Since averaging proceeds by shifting the starting point", "through the trajectory, the shift can be modified with -sh to enable the choice of uncorrelated", "starting points. Towards the end, statistical inaccuracy grows and integrating the", "correlation function only yields reliable values until a certain point, depending on", "the number of frames. The option -tr controls the region of the integral taken into account", "for calculating the static dielectric constant.", "[PAR]", "Option -temp sets the temperature required for the computation of the static dielectric constant." "[PAR]", "Option -eps controls the dielectric constant of the surrounding medium for simulations using", "a Reaction Field or dipole corrections of the Ewald summation (eps=0 corresponds to", "thinfoil boundary conditions).", "[PAR]", "-[no]rmjump allows the molecules to jump over the box edges (this is implied if -dsp is set)" }; /* At first the arguments will be parsed and the system information processed */ CopyRight(stderr,argv[0]); parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW, NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL); bDSP = opt2bSet("-dsp",NFILE,fnm); bACF = opt2bSet("-caf",NFILE,fnm); bINT = opt2bSet("-mc",NFILE,fnm); bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,&vtop,box,TRUE); sfree(xtop); sfree(vtop); indexfn = ftp2fn_null(efNDX,NFILE,fnm); snew(grpname,1); get_index(&(top.atoms),indexfn,1,&isize,&index0,grpname); flags = flags | TRX_READ_X | TRX_READ_V; read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags); snew(mass2,top.atoms.nr); snew(qmol,top.mols.nr); snew(cmol,2); bNEU=precalc(top,fr,mass2,qmol,cmol); printf("\n\n %d cations and %d anions out of %d molecules\n",cmol[0],cmol[1],top.mols.nr); snew(cati,cmol[0]); snew(ani,cmol[1]); calc_molcharge(top,qmol,cati,ani,cmol); if (fr.bV){ if(bACF){ outf =xvgropen(opt2fn("-caf",NFILE,fnm), "Current autocorrelation function",xvgr_tlabel(),"ACF (e nm/ps)\\S2"); fprintf(outf,"# time\t acf\t average \t std.dev\n"); } fcur =xvgropen(opt2fn("-o",NFILE,fnm), "Current",xvgr_tlabel(),"J(t) (e nm/ps)"); fprintf(fcur,"# time\t Jx\t Jy \t J_z \n"); if(bINT){ mcor = xvgropen(opt2fn("-mc",NFILE,fnm), "M\\sD\\N - current autocorrelation function",xvgr_tlabel(), "< M\\sD\\N (0)\\c7\\CJ(t) > (e nm/ps)\\S2"); fprintf(mcor,"# time\t M_D(0) J(t) acf \t Integral acf\n"); } } fmj = xvgropen(opt2fn("-mj",NFILE,fnm), "Averaged translational part of M",xvgr_tlabel(),"< M\\sJ\\N > (enm)"); fprintf(fmj,"# time\t x\t y \t z \t average of M_J^2 \t std.dev\n"); fmd = xvgropen(opt2fn("-md",NFILE,fnm), "Averaged rotational part of M",xvgr_tlabel(),"< M\\sD\\N > (enm)"); fprintf(fmd,"# time\t x\t y \t z \t average of M_D^2 \t std.dev\n"); if(bDSP) fmjdsp = xvgropen(opt2fn("-dsp",NFILE,fnm), "MSD of the squared translational dipole moment M",xvgr_tlabel(),"<|M\\sJ\\N(t)-M\\sJ\\N(0)|\\S2\\N > / 6.0*V*k\\sB\\N*T / Sm\\S-1\\Nps\\S-1\\N"); /* System information is read and prepared, dielectric() processes the frames and calculates the dielectric constant */ dielectric(fmj,fmd,outf,fcur,mcor,fmjdsp,bNoJump,bDSP,bACF,bINT,ePBC,top,fr,temp,trust,bfit,efit,bvit,evit,status,isize,nshift,ani,cati,cmol,index0,mass2,qmol,eps_rf); fclose(fmj); fclose(fmd); if(bDSP) fclose(fmjdsp); if (bACF) fclose(outf); if(bINT) fclose(mcor); thanx(stderr); return 0; }