Added warning in g_dipoles
authorDavid van der Spoel <spoel@anfinsen.bmc.uu.se>
Sat, 25 Dec 2010 08:55:12 +0000 (09:55 +0100)
committerDavid van der Spoel <spoel@anfinsen.bmc.uu.se>
Sat, 25 Dec 2010 08:55:12 +0000 (09:55 +0100)
Removed heat capacity calculation from g_energy.

src/tools/gmx_dipoles.c
src/tools/gmx_energy.c

index 1d1e3c48ee5fbdc8529b7ebdfa8598d665aa5ae9..e159f736d9c8269f96ee5a1ee1410c7102acfecf 100644 (file)
@@ -683,7 +683,7 @@ static void do_dip(t_topology *top,int ePBC,real volume,
     real       rcut=0,t,t0,t1,dt,lambda,dd,rms_cos;
     rvec       dipaxis;
     matrix     box;
-    gmx_bool       bCorr,bTotal,bCont;
+    gmx_bool   bCorr,bTotal,bCont;
     double     M_diff=0,epsilon,invtel,vol_aver;
     double     mu_ave,mu_mol,M2_ave=0,M_ave2=0,M_av[DIM],M_av2[DIM];
     double     M[3],M2[3],M4[3],Gk=0,g_k=0;
@@ -726,7 +726,7 @@ static void do_dip(t_topology *top,int ePBC,real volume,
         mols = &(top->mols);
     }
   
-    if (iVol == -1)
+    if ((iVol == -1) && bMU)
         printf("Using Volume from topology: %g nm^3\n",volume);
 
     /* Correlation stuff */ 
@@ -1349,7 +1349,7 @@ int gmx_dipoles(int argc,char *argv[])
     int          nFF[2];
     atom_id      **grpindex;
     char         **grpname=NULL;
-    gmx_bool         bCorr,bQuad,bGkr,bMU,bSlab;  
+    gmx_bool     bCorr,bQuad,bGkr,bMU,bSlab;  
     t_filenm fnm[] = {
         { efEDR, "-en", NULL,         ffOPTRD },
         { efTRX, "-f", NULL,           ffREAD },
index cd299645b164f26c58635b1cae90ebae8cfa1193..1c72609710832b452c1ac2c0d777300a7d040745 100644 (file)
@@ -841,9 +841,8 @@ static char *ee_pr(double ee,char *buf)
 }
 
 static void analyse_ener(gmx_bool bCorr,const char *corrfn,
-                         gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
+                         gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,
                          gmx_bool bVisco,const char *visfn,int nmol,
-                         int nconstr,
                          gmx_large_int_t start_step,double start_t,
                          gmx_large_int_t step,double t,
                          double time[], real reftemp,
@@ -864,7 +863,7 @@ static void analyse_ener(gmx_bool bCorr,const char *corrfn,
   gmx_large_int_t nsteps;
   int  nexact,nnotexact;
   double x1m,x1mk;
-  real Temp=-1,Pres=-1,VarV=-1,VarT=-1,VarEtot=-1,AvEtot=0,VarEnthalpy=-1;
+  real Temp=-1,Pres=-1,VarV=-1,VarT=-1,VarEtot=-1,AvEtot=0;
   int  i,j,nout;
   real chi2;
   char buf[256],eebuf[100];
@@ -962,9 +961,7 @@ static void analyse_ener(gmx_bool bCorr,const char *corrfn,
       } else if (strstr(leg[i],"otal") != NULL) {
        VarEtot = sqr(stddev);
        AvEtot = aver;
-      } else if (strstr(leg[i],"nthalpy") != NULL) {
-       VarEnthalpy = sqr(stddev);
-      }
+      } 
       if (bIsEner[i]) {
        pr_aver   = aver/nmol-ezero;
        pr_stddev = stddev/nmol;
@@ -1007,35 +1004,6 @@ static void analyse_ener(gmx_bool bCorr,const char *corrfn,
       else
        fprintf(stdout,"\n");
     }
-    if (bTempFluct && Temp != -1) {
-      printf("\nTemperature dependent fluctuation properties at T = %g. #constr/mol = %d\n",Temp,nconstr);
-      if (nmol < 2)
-       printf("Warning: nmol = %d, this may not be what you want.\n",
-              nmol);
-      if (VarV != -1) {
-       real tmp = VarV/(Vaver*BOLTZ*Temp*PRESFAC);
-       
-       printf("Isothermal Compressibility: %10g /%s\n",
-              tmp,unit_pres_bar);
-       printf("Adiabatic bulk modulus:     %10g  %s\n",
-              1.0/tmp,unit_pres_bar);
-      }
-      if (VarEnthalpy != -1) {
-       real Cp = 1000*((VarEnthalpy/nmol)/(BOLTZ*Temp*Temp) - 
-                       0.5*BOLTZ*nconstr);
-       printf("Heat capacity at constant pressure Cp: %10g J/mol K\n",Cp);
-      }
-      if ((VarV != -1) && (VarEnthalpy != -1)) {
-       real aP = (sqrt(VarEnthalpy*VarV/nmol))/(BOLTZ*Vaver*Temp*Temp);
-       printf("Thermal expansion coefficient alphaP: %10g 1/K\n",aP);
-      }
-      if ((VarV == -1) && (VarEtot != -1)) {
-       real Cv = 1000*((VarEtot/nmol)/(BOLTZ*Temp*Temp) - 
-                       0.5*BOLTZ*nconstr);
-       printf("Heat capacity at constant volume Cv: %10g J/mol K\n",Cv);
-      }
-      please_cite(stdout,"Allen1987a");
-    }
     /* Do correlation function */
     if (edat->nframes > 1)
     {
@@ -1525,24 +1493,6 @@ int gmx_energy(int argc,char *argv[])
 
     "The term fluctuation gives the RMSD around the LSQ fit.[PAR]",
     
-    "Some fluctuation-dependent properties can be calculated provided",
-    "the correct energy terms are selected. The following properties",
-    "will be computed:[BR]",
-    "Property                        Energy terms needed[BR]",
-    "---------------------------------------------------[BR]",
-    "Heat capacity Cp (NPT sims):    Enthalpy, Temp     [BR]",
-    "Heat capacity Cv (NVT sims):    Etot, Temp         [BR]",
-    "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
-    "Isothermal compressibility:     Vol, Temp          [BR]",
-    "Adiabatic bulk modulus:         Vol, Temp          [BR]",
-    "---------------------------------------------------[BR]",
-    "You always need to set the number of molecules [TT]-nmol[tt], and,",
-    "if you used constraints in your simulations you will need to give",
-    "the number of constraints per molecule [TT]-nconstr[tt] in order to",
-    "correct for this: (nconstr/2) kB is subtracted from the heat",
-    "capacity in this case. For instance in the case of rigid water",
-    "you need to give the value 3 to this option.[PAR]",
-    
     "When the [TT]-viol[tt] option is set, the time averaged",
     "violations are plotted and the running time-averaged and",
     "instantaneous sum of violations are recalculated. Additionally",
@@ -1594,7 +1544,7 @@ int gmx_energy(int argc,char *argv[])
   };
   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE;
   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
-  static int  skip=0,nmol=1,nconstr=0,nbmin=5,nbmax=5;
+  static int  skip=0,nmol=1,nbmin=5,nbmax=5;
   static real reftemp=300.0,ezero=0;
   t_pargs pa[] = {
     { "-fee",   FALSE, etBOOL,  {&bFee},
@@ -1619,8 +1569,6 @@ int gmx_energy(int argc,char *argv[])
       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
     { "-nmol", FALSE, etINT,  {&nmol},
       "Number of molecules in your sample: the energies are divided by this number" },
-    { "-nconstr",  FALSE, etINT,  {&nconstr},
-      "Number of constraints per molecule. Necessary for calculating the heat capacity" },
     { "-fluc", FALSE, etBOOL, {&bFluct},
       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
     { "-orinst", FALSE, etBOOL, {&bOrinst},
@@ -2336,9 +2284,9 @@ int gmx_energy(int argc,char *argv[])
   else
   {
       analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
-                   bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
+                   bFee,bSum,bFluct,
                    bVisco,opt2fn("-vis",NFILE,fnm),
-                   nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
+                   nmol,start_step,start_t,frame[cur].step,frame[cur].t,
                    time,reftemp,&edat,
                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
                    oenv);