/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
*
- *
+ *
* This source code is part of
- *
+ *
* G R O M A C S
- *
+ *
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.2.0
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* 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.
- *
+ *
* For more info, check our website at http://www.gromacs.org
- *
+ *
* And Hey:
* Green Red Orange Magenta Azure Cyan Skyblue
*/
#include "viewit.h"
#include "mtop_util.h"
#include "gmx_ana.h"
-
+#include "mdebin.h"
static real minthird=-1.0/3.0,minsixth=-1.0/6.0;
{
if (x > 0)
return pow(x,y);
- else
+ else
return 0.0;
}
int n,k,j,i;
int *set;
gmx_bool bVerbose = TRUE;
-
+
if ((getenv("VERBOSE")) != NULL)
bVerbose = FALSE;
-
+
fprintf(stderr,"Select the terms you want from the following list\n");
fprintf(stderr,"End your selection with 0\n");
if ( bVerbose ) {
for(k=0; (k<nre); ) {
- for(j=0; (j<4) && (k<nre); j++,k++)
+ for(j=0; (j<4) && (k<nre); j++,k++)
fprintf(stderr," %3d=%14s",k+1,nm[k]);
fprintf(stderr,"\n");
}
for(i=(*nset)=0; (i<nre); i++)
if (bE[i])
set[(*nset)++]=i;
-
+
sfree(bE);
-
- return set;
-}
-static int strcount(const char *s1,const char *s2)
-{
- int n=0;
- while (s1 && s2 && (toupper(s1[n]) == toupper(s2[n])))
- n++;
- return n;
+ return set;
}
static void chomp(char *buf)
{
int len = strlen(buf);
-
+
while ((len > 0) && (buf[len-1] == '\n')) {
buf[len-1] = '\0';
len--;
const char *fm4="%3d %-14s";
const char *fm2="%3d %-34s";
char **newnm=NULL;
-
+
if ((getenv("VERBOSE")) != NULL)
bVerbose = FALSE;
-
+
fprintf(stderr,"\n");
fprintf(stderr,"Select the terms you want from the following list by\n");
fprintf(stderr,"selecting either (part of) the name or the number or a combination.\n");
fprintf(stderr,"End your selection with an empty line or a zero.\n");
fprintf(stderr,"-------------------------------------------------------------------\n");
-
+
snew(newnm,nre);
j = 0;
for(k=0; k<nre; k++) {
if ( bVerbose ) {
fprintf(stderr,"\n\n");
}
-
+
snew(bE,nre);
-
+
bEOF = FALSE;
while (!bEOF && (fgets2(buf,STRLEN-1,stdin))) {
/* Remove newlines */
chomp(buf);
-
+
/* Remove spaces */
trim(buf);
-
+
/* Empty line means end of input */
bEOF = (strlen(buf) == 0);
if (!bEOF) {
} while (!bEOF && (ptr && (strlen(ptr) > 0)));
}
}
-
+
snew(set,nre);
for(i=(*nset)=0; (i<nre); i++)
if (bE[i])
set[(*nset)++]=i;
-
+
sfree(bE);
-
+
if (*nset == 0)
gmx_fatal(FARGS,"No energy terms selected");
- for(i=0; (i<nre); i++)
+ for(i=0; (i<nre); i++)
sfree(newnm[i]);
sfree(newnm);
-
+
return set;
}
+static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
+{
+ gmx_mtop_t mtop;
+ int natoms;
+ t_iatom *iatom;
+ matrix box;
+
+ /* all we need is the ir to be able to write the label */
+ read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,&mtop);
+}
+
static void get_orires_parms(const char *topnm,
- int *nor,int *nex,int **label,real **obs)
+ int *nor,int *nex,int **label,real **obs)
{
gmx_mtop_t mtop;
gmx_localtop_t *top;
ip = top->idef.iparams;
iatom = top->idef.il[F_ORIRES].iatoms;
-
+
/* Count how many distance restraint there are... */
nb = top->idef.il[F_ORIRES].nr;
if (nb == 0)
gmx_fatal(FARGS,"No orientation restraints in topology!\n");
-
+
*nor = nb/3;
*nex = 0;
snew(*label,*nor);
functype = top->idef.functype;
ip = top->idef.iparams;
-
+
/* Count how many distance restraint there are... */
nb=top->idef.il[F_DISRES].nr;
if (nb == 0)
gmx_fatal(FARGS,"No distance restraints in topology!\n");
-
+
/* Allocate memory */
snew(b,nb);
snew(ind,nb);
snew(pair,nb+1);
-
+
/* Fill the bound array */
nb=0;
for(i=0; (i<top->idef.ntypes); i++) {
}
}
*bounds = b;
-
+
/* Fill the index array */
label1 = -1;
disres = &(top->idef.il[F_DISRES]);
natom = interaction_function[ftype].nratoms+1;
if (label1 != top->idef.iparams[type].disres.label) {
pair[j] = k;
- label1 = top->idef.iparams[type].disres.label;
+ label1 = top->idef.iparams[type].disres.label;
j ++;
}
k++;
*index = ind;
*dr_pair = pair;
-
+
return nb;
}
const real sixth=1.0/6.0;
int i,j;
double rsum,rav,sumaver,sumt;
-
+
sumaver = 0;
sumt = 0;
for(i=0; (i<nb); i++) {
}
rsum = max(0.0,mypow(rsum,-sixth)-bounds[i]);
rav = max(0.0,mypow(rav, -sixth)-bounds[i]);
-
+
sumt += rsum;
sumaver += rav;
}
sumaver);
fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",
sumt);
-#endif
+#endif
vout=xvgropen(voutfn,"r\\S-3\\N average violations","DR Index","nm",
oenv);
sum = 0.0;
for(i=0; (i<nbounds); i++) {
/* Do ensemble averaging */
sumaver = 0;
- for(j=pair[i]; (j<pair[i+1]); j++)
- sumaver += sqr(violaver[j]/nframes);
+ for(j=pair[i]; (j<pair[i+1]); j++)
+ sumaver += sqr(violaver[j]/nframes);
sumaver = max(0.0,mypow(sumaver,minsixth)-bounds[i]);
sumt += sumaver;
static void add_ee_av(ee_sum_t *ees)
{
double av;
-
+
av = ees->sum/ees->np;
ees->sav += av;
ees->sav2 += av*av;
for(i=0; i<nset; i++)
{
ed = &edat->s[i];
-
+
sum = 0;
sum2 = 0;
np = 0;
sump = ed->ener[f];
sum2 += dsqr(sump);
}
-
+
/* sum has to be increased after sum2 */
np += p;
sum += sump;
enerdat_t *s;
int f,i;
double sum;
-
+
snew(esum,1);
*esum = *edat;
snew(esum->s,1);
s = &esum->s[0];
snew(s->ener,esum->nframes);
snew(s->es ,esum->nframes);
-
+
s->bExactStat = TRUE;
s->slope = 0;
for(i=0; i<nset; i++)
}
s->slope += edat->s[i].slope;
}
-
+
for(f=0; f<edat->nframes; f++)
{
sum = 0;
s->es[f].sum = sum;
s->es[f].sum2 = 0;
}
-
+
calc_averages(1,esum,nbmin,nbmax);
return esum;
{
char tmp[100];
double rnd;
-
+
if (ee < 0)
{
sprintf(buf,"%s","--");
edat->npoints = edat->nframes;
edat->nsteps = edat->nframes;
-
+
for(k=0; (k<5); k++)
{
- for(i=0; (i<nset); i++)
+ for(i=0; (i<nset); i++)
{
delta = edat->s[i].slope*dt;
-
+
if (NULL != debug)
fprintf(debug,"slope for set %d is %g\n",i,edat->s[i].slope);
-
+
for(j=0; (j<edat->nframes); j++)
{
edat->s[i].ener[j] -= j*delta;
enum { eVol, eEnth, eTemp, eEtot, eNR };
const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
int ii[eNR];
-
+
NANO3 = NANO*NANO*NANO;
if (!bDriftCorr)
{
{
remove_drift(nset,nbmin,nbmax,dt,edat);
}
- for(i=0; (i<eNR); i++)
+ for(i=0; (i<eNR); i++)
{
- for(ii[i]=0; (ii[i]<nset &&
+ for(ii[i]=0; (ii[i]<nset &&
(gmx_strcasecmp(leg[ii[i]],my_ener[i]) != 0)); ii[i]++)
;
/* if (ii[i] < nset)
*/ }
/* Compute it all! */
vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
-
+
/* Temperature */
tt = NOTSET;
if (ii[eTemp] < nset)
if (debug != NULL)
{
- if (varv != NOTSET)
+ if (varv != NOTSET)
fprintf(fp,"varv = %10g (m^6)\n",varv*AVOGADRO/nmol);
if (vvhh != NOTSET)
fprintf(fp,"vvhh = %10g (m^3 J)\n",vvhh);
dcp);
please_cite(fp,"Allen1987a");
}
- else
+ else
{
fprintf(fp,"You should select the temperature in order to obtain fluctuation properties.\n");
}
else {
/* Calculate the time difference */
delta_t = t - start_t;
-
+
fprintf(stdout,"\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
gmx_step_str(nsteps,buf),start_t,t,nset);
calc_averages(nset,edat,nbmin,nbmax);
-
+
if (bSum) {
esum = calc_sum(nset,edat,nbmin,nbmax);
}
}
}
}
-
+
if (nnotexact == 0) {
fprintf(stdout,"All statistics are over %s points\n",
gmx_step_str(edat->npoints,buf));
else
fprintf(stdout,"\n");
fprintf(stdout,"-------------------------------------------------------------------------------\n");
-
+
/* Initiate locals, only used with -sum */
expEtot=0;
if (bFee) {
for(j=0; (j<edat->nframes); j++) {
expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
}
- if (bSum)
+ if (bSum)
expEtot+=expE/edat->nframes;
-
+
fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
}
if (strstr(leg[i],"empera") != NULL) {
Vaver= aver;
} else if (strstr(leg[i],"essure") != NULL) {
Pres = aver;
- }
+ }
if (bIsEner[i]) {
pr_aver = aver/nmol-ezero;
pr_stddev = stddev/nmol;
fprintf(stdout,"%-24s %10g %10s %10g %10g",
leg[i],pr_aver,ee_pr(pr_errest,eebuf),pr_stddev,totaldrift);
- if (bFee)
+ if (bFee)
fprintf(stdout," %10g",fee[i]);
-
+
fprintf(stdout," (%s)\n",enm[set[i]].unit);
if (bFluct) {
"Total",esum->s[0].av/nmol,ee_pr(esum->s[0].ee/nmol,eebuf),
"--",totaldrift/nmol,enm[set[0]].unit);
/* pr_aver,pr_stddev,a,totaldrift */
- if (bFee)
+ if (bFee)
fprintf(stdout," %10g %10g\n",
log(expEtot)/beta + esum->s[0].av/nmol,log(expEtot)/beta);
else
fprintf(stdout,"\n");
}
-
+
/* Do correlation function */
if (edat->nframes > 1)
{
real factor;
real **eneset;
real **enesum;
-
+
/* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
-
- /* Symmetrise tensor! (and store in first three elements)
+
+ /* Symmetrise tensor! (and store in first three elements)
* And subtract average pressure!
*/
snew(eneset,12);
enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
}
-
+
einstein_visco("evisco.xvg","eviscoi.xvg",
3,edat->nframes,enesum,Vaver,Temp,nsteps,time,oenv);
-
+
/*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
/* Do it for shear viscosity */
strcpy(buf,"Shear Viscosity");
low_do_autocorr(corrfn,oenv,buf,edat->nframes,3,
(edat->nframes+1)/2,eneset,Dt,
eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
-
+
/* Now for bulk viscosity */
strcpy(buf,"Bulk Viscosity");
low_do_autocorr(corrfn,oenv,buf,edat->nframes,1,
(edat->nframes+1)/2,&(eneset[11]),Dt,
eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
-
+
factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)",oenv);
xvgr_legend(fp,asize(leg),leg,oenv);
-
+
/* Use trapezium rule for integration */
integral = 0;
intBulk = 0;
nout = get_acfnout();
if ((nout < 2) || (nout >= edat->nframes/2))
nout = edat->nframes/2;
- for(i=1; (i<nout); i++)
+ for(i=1; (i<nout); i++)
{
integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
fprintf(fp," %10.6f",e);
}
-static void fec(const char *ene2fn, const char *runavgfn,
- real reftemp, int nset, int set[], char *leg[],
+static void fec(const char *ene2fn, const char *runavgfn,
+ real reftemp, int nset, int set[], char *leg[],
enerdata_t *edat, double time[],
const output_env_t oenv)
{
- const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
+ const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
"<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
FILE *fp;
ener_file_t enx;
gmx_enxnm_t *enm=NULL;
t_enxframe *fr;
char buf[22];
-
+
/* read second energy file */
snew(fr,1);
enm = NULL;
enx = open_enx(ene2fn,"r");
do_enxnms(enx,&(fr->nre),&enm);
-
+
snew(eneset2,nset+1);
nenergy2=0;
maxenergy=0;
timecheck=0;
do {
- /* This loop searches for the first frame (when -b option is given),
+ /* This loop searches for the first frame (when -b option is given),
* or when this has been found it reads just one energy frame
*/
do {
bCont = do_enx(enx,fr);
-
+
if (bCont)
timecheck = check_times(fr->t);
-
+
} while (bCont && (timecheck < 0));
-
+
/* Store energies for analysis afterwards... */
if ((timecheck == 0) && bCont) {
if (fr->nre > 0) {
}
}
} while (bCont && (timecheck == 0));
-
+
/* check */
if (edat->nframes != nenergy2) {
fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
edat->nframes,nenergy2);
}
nenergy = min(edat->nframes,nenergy2);
-
+
/* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
fp=NULL;
if (runavgfn) {
dE = eneset2[i][j] - edat->s[i].ener[j];
sum += exp(-dE*beta);
if (fp)
- fprintf(fp,"%10g %10g %10g\n",
+ fprintf(fp,"%10g %10g %10g\n",
time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
}
aver = -BOLTZ*reftemp*log(sum/nenergy);
}
-static void do_dhdl(t_enxframe *fr, FILE **fp_dhdl, const char *filename,
- int *blocks, int *hists, int *samples, int *nlambdas,
- const output_env_t oenv)
+static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, const char *filename, gmx_bool bDp,
+ int *blocks, int *hists, int *samples, int *nlambdas, const output_env_t oenv)
{
const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
{
if (nblock_dh>0)
{
- sprintf(title,"%s, %s",dhdl,deltag);
- sprintf(label_x,"%s (%s)","Time",unit_time);
- sprintf(label_y,"(%s)",unit_energy);
+ /* we have standard, non-histogram data -- call open_dhdl to open the file */
+ *fp_dhdl=open_dhdl(filename,ir,oenv);
}
else
{
sprintf(title,"N(%s)",deltag);
sprintf(label_x,"%s (%s)",deltag,unit_energy);
sprintf(label_y,"Samples");
- }
- *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,
- oenv);
- if (! changing_lambda)
- {
+ *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,oenv);
sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
+ xvgr_subtitle(*fp_dhdl,buf,oenv);
}
- else
- {
- sprintf(buf,"T = %g (K)", temp);
- }
- xvgr_subtitle(*fp_dhdl,buf,oenv);
- first=TRUE;
}
-
-
(*hists)+=nblock_hist;
(*blocks)+=nblock_dh;
(*nlambdas) = nblock_hist+nblock_dh;
-
/* write the data */
if (nblock_hist > 0)
{
if (!derivative)
{
sprintf(legend, "N(%s(%s=%g) | %s=%g)",
- deltag, lambda, foreign_lambda,
+ deltag, lambda, foreign_lambda,
lambda, start_lambda);
}
else
{
- sprintf(legend, "N(%s | %s=%g)",
+ sprintf(legend, "N(%s | %s=%g)",
dhdl, lambda, start_lambda);
}
-
+
lg[0]=legend;
- xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
+ xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
setnr++;
for(k=0;k<blk->sub[j+2].nr;k++)
{
int hist;
double xmin, xmax;
-
+
hist=blk->sub[j+2].ival[k];
xmin=(x0+k)*dx;
xmax=(x0+k+1)*dx;
- fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
+ fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
xmax, hist);
sum+=hist;
}
}
}
}
-
(*samples) += (int)(sum/nblock_hist);
}
else
char **setnames=NULL;
int nnames=nblock_dh;
- if (changing_lambda)
- {
- nnames++;
- }
- if (first)
- {
- snew(setnames, nnames);
- }
- j=0;
-
- if ( changing_lambda && first)
- {
- /* lambda is a plotted value */
- setnames[j]=gmx_strdup(lambda);
- j++;
- }
-
-
for(i=0;i<fr->nblock;i++)
{
t_enxblock *blk=&(fr->block[i]);
if (blk->id == enxDH)
{
- if (first)
- {
- /* do the legends */
- int derivative;
- double foreign_lambda;
-
- derivative=blk->sub[0].ival[0];
- foreign_lambda=blk->sub[1].dval[0];
-
- if (derivative)
- {
- sprintf(buf, "%s %s %g",dhdl,lambda,start_lambda);
- }
- else
- {
- sprintf(buf, "%s %s %g",deltag,lambda, foreign_lambda);
- }
- setnames[j] = gmx_strdup(buf);
- j++;
- }
-
if (len == 0)
- {
+ {
len=blk->sub[2].nr;
}
else
}
}
}
-
-
- if (first)
- {
- xvgr_legend(*fp_dhdl, nblock_dh, (const char**)setnames, oenv);
- setnr += nblock_dh;
- for(i=0;i<nblock_dh;i++)
- {
- sfree(setnames[i]);
- }
- sfree(setnames);
- }
-
(*samples) += len;
+
for(i=0;i<len;i++)
{
double time=start_time + delta_time*i;
- fprintf(*fp_dhdl,"%.4f", time);
- if (fabs(delta_lambda) > 1e-9)
- {
- double lambda_now=i*delta_lambda + start_lambda;
- fprintf(*fp_dhdl," %.4f", lambda_now);
- }
+ fprintf(*fp_dhdl,"%.4f ", time);
+
for(j=0;j<fr->nblock;j++)
{
t_enxblock *blk=&(fr->block[j]);
{
value=blk->sub[2].dval[i];
}
- fprintf(*fp_dhdl," %g", value);
+ /* we need to decide which data type it is based on the count*/
+
+ if (j==1 && ir->bExpanded)
+ {
+ fprintf(*fp_dhdl,"%4d", (int)value); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
+ } else {
+ if (bDp) {
+ fprintf(*fp_dhdl," %#.12g", value); /* print normal precision */
+ }
+ else
+ {
+ fprintf(*fp_dhdl," %#.8g", value); /* print normal precision */
+ }
+ }
}
}
fprintf(*fp_dhdl, "\n");
int gmx_energy(int argc,char *argv[])
{
const char *desc[] = {
-
+
"[TT]g_energy[tt] extracts energy components or distance restraint",
"data from an energy file. The user is prompted to interactively",
"select the desired energy terms.[PAR]",
-
+
"Average, RMSD, and drift are calculated with full precision from the",
"simulation (see printed manual). Drift is calculated by performing",
"a least-squares fit of the data to a straight line. The reported total drift",
"energy values.[PAR]",
"The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
-
+
"Some fluctuation-dependent properties can be calculated provided",
"the correct energy terms are selected, and that the command line option",
"[TT]-fluct_props[tt] is given. The following properties",
" [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]",
" [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
"[PAR]",
-
+
"When a second energy file is specified ([TT]-f2[tt]), a free energy",
"difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
"where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
"files, and the average is over the ensemble A. The running average",
"of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
"[BB]Note[bb] that the energies must both be calculated from the same trajectory."
-
+
};
static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE,bFluctProps=FALSE;
"Print energies in high precision" },
{ "-nbmin", FALSE, etINT, {&nbmin},
"Minimum number of blocks for error estimate" },
- { "-nbmax", FALSE, etINT, {&nbmax},
+ { "-nbmax", FALSE, etINT, {&nbmax},
"Maximum number of blocks for error estimate" },
{ "-mutot",FALSE, etBOOL, {&bMutot},
"Compute the total dipole moment from the components" },
"Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
"Volume", "Pressure"
};
-
+
FILE *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
FILE *fp_dhdl=NULL;
FILE **drout;
#define NFILE asize(fnm)
int npargs;
t_pargs *ppa;
-
+
CopyRight(stderr,argv[0]);
npargs = asize(pa);
ppa = add_acf_pargs(&npargs,pa);
parse_common_args(&argc,argv,
PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
-
+
bDRAll = opt2bSet("-pairs",NFILE,fnm);
bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
bORA = opt2bSet("-ora",NFILE,fnm);
do_enxnms(fp,&nre,&enm);
Vaver = -1;
-
+
bVisco = opt2bSet("-vis",NFILE,fnm);
-
- if (!bDisRe && !bDHDL)
+
+ if ((!bDisRe) && (!bDHDL))
{
if (bVisco) {
nset=asize(setnm);
}
}
}
- else
+ else
{
set=select_by_name(nre,enm,&nset);
}
snew(violaver,npairs);
out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
"Time (ps)","nm",oenv);
- xvgr_legend(out,2,drleg,oenv);
- if (bDRAll) {
+ xvgr_legend(out,2,drleg,oenv);
+ if (bDRAll) {
fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
"Time (ps)","Distance (nm)",oenv);
if (output_env_get_print_xvgr_codes(oenv))
fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
ir.dr_tau);
}
+ } else if (bDHDL) {
+ get_dhdl_parms(ftp2fn(efTPX,NFILE,fnm),&ir);
}
-
- /* Initiate energies and set them to zero */
- edat.nsteps = 0;
- edat.npoints = 0;
- edat.nframes = 0;
- edat.step = NULL;
- edat.steps = NULL;
- edat.points = NULL;
- snew(edat.s,nset);
-
- /* Initiate counters */
- teller = 0;
- teller_disre = 0;
- bFoundStart = FALSE;
- start_step = 0;
- start_t = 0;
- do {
- /* This loop searches for the first frame (when -b option is given),
- * or when this has been found it reads just one energy frame
- */
- do {
- bCont = do_enx(fp,&(frame[NEXT]));
-
- if (bCont) {
- timecheck = check_times(frame[NEXT].t);
- }
- } while (bCont && (timecheck < 0));
-
- if ((timecheck == 0) && bCont) {
- /* We read a valid frame, so we can use it */
- fr = &(frame[NEXT]);
-
- if (fr->nre > 0) {
- /* The frame contains energies, so update cur */
- cur = NEXT;
-
- if (edat.nframes % 1000 == 0)
- {
- srenew(edat.step,edat.nframes+1000);
- memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
- srenew(edat.steps,edat.nframes+1000);
- memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
- srenew(edat.points,edat.nframes+1000);
- memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
- for(i=0; i<nset; i++)
- {
- srenew(edat.s[i].ener,edat.nframes+1000);
- memset(&(edat.s[i].ener[edat.nframes]),0,
- 1000*sizeof(edat.s[i].ener[0]));
-
- srenew(edat.s[i].es ,edat.nframes+1000);
- memset(&(edat.s[i].es[edat.nframes]),0,
- 1000*sizeof(edat.s[i].es[0]));
- }
- }
-
- nfr = edat.nframes;
- edat.step[nfr] = fr->step;
-
- if (!bFoundStart)
- {
- bFoundStart = TRUE;
- /* Initiate the previous step data */
- start_step = fr->step;
- start_t = fr->t;
- /* Initiate the energy sums */
- edat.steps[nfr] = 1;
- edat.points[nfr] = 1;
- for(i=0; i<nset; i++)
- {
- sss = set[i];
- edat.s[i].es[nfr].sum = fr->ener[sss].e;
- edat.s[i].es[nfr].sum2 = 0;
- }
- edat.nsteps = 1;
- edat.npoints = 1;
- }
- else
- {
- edat.steps[nfr] = fr->nsteps;
- {
- if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
- {
- if (fr->nsum <= 1)
- {
- edat.points[nfr] = 1;
- for(i=0; i<nset; i++)
- {
- sss = set[i];
- edat.s[i].es[nfr].sum = fr->ener[sss].e;
- edat.s[i].es[nfr].sum2 = 0;
- }
- edat.npoints += 1;
- }
- else
- {
- edat.points[nfr] = fr->nsum;
- for(i=0; i<nset; i++)
- {
- sss = set[i];
- edat.s[i].es[nfr].sum = fr->ener[sss].esum;
- edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
- }
- edat.npoints += fr->nsum;
- }
- }
- else
- {
- /* The interval does not match fr->nsteps:
- * can not do exact averages.
- */
- edat.npoints = 0;
- }
- edat.nsteps = fr->step - start_step + 1;
- }
- }
- for(i=0; i<nset; i++)
- {
- edat.s[i].ener[nfr] = fr->ener[set[i]].e;
- }
- }
- /*
- * Define distance restraint legends. Can only be done after
- * the first frame has been read... (Then we know how many there are)
- */
- blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
- if (bDisRe && bDRAll && !leg && blk_disre)
- {
- t_iatom *fa;
- t_iparams *ip;
-
- fa = top->idef.il[F_DISRES].iatoms;
- ip = top->idef.iparams;
- if (blk_disre->nsub != 2 ||
- (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
- {
- gmx_incons("Number of disre sub-blocks not equal to 2");
- }
-
- ndisre=blk_disre->sub[0].nr ;
- if (ndisre != top->idef.il[F_DISRES].nr/3)
- {
- gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
- ndisre,top->idef.il[F_DISRES].nr/3);
- }
- snew(pairleg,ndisre);
- for(i=0; i<ndisre; i++)
- {
- snew(pairleg[i],30);
- j=fa[3*i+1];
- k=fa[3*i+2];
- gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
- gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
- sprintf(pairleg[i],"%d %s %d %s (%d)",
- resnr_j,anm_j,resnr_k,anm_k,
- ip[fa[3*i]].disres.label);
- }
- set=select_it(ndisre,pairleg,&nset);
- snew(leg,2*nset);
- for(i=0; (i<nset); i++)
- {
- snew(leg[2*i],32);
- sprintf(leg[2*i], "a %s",pairleg[set[i]]);
- snew(leg[2*i+1],32);
- sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
- }
- xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
- }
-
- /*
- * Store energies for analysis afterwards...
- */
- if (!bDisRe && !bDHDL && (fr->nre > 0)) {
- if (edat.nframes % 1000 == 0) {
- srenew(time,edat.nframes+1000);
- }
- time[edat.nframes] = fr->t;
- edat.nframes++;
- }
- /*
- * Printing time, only when we do not want to skip frames
- */
- if (!skip || teller % skip == 0) {
- if (bDisRe) {
- /*******************************************
- * D I S T A N C E R E S T R A I N T S
- *******************************************/
- if (ndisre > 0)
- {
-#ifndef GMX_DOUBLE
- float *disre_rt = blk_disre->sub[0].fval;
- float *disre_rm3tav = blk_disre->sub[1].fval;
-#else
- double *disre_rt = blk_disre->sub[0].dval;
- double *disre_rm3tav = blk_disre->sub[1].dval;
-#endif
-
- print_time(out,fr->t);
- if (violaver == NULL)
- snew(violaver,ndisre);
-
- /* Subtract bounds from distances, to calculate violations */
- calc_violations(disre_rt, disre_rm3tav,
- nbounds,pair,bounds,violaver,&sumt,&sumaver);
-
- fprintf(out," %8.4f %8.4f\n",sumaver,sumt);
- if (bDRAll) {
- print_time(fp_pairs,fr->t);
- for(i=0; (i<nset); i++) {
- sss=set[i];
- fprintf(fp_pairs," %8.4f", mypow(disre_rm3tav[sss],minthird));
- fprintf(fp_pairs," %8.4f", disre_rt[sss]);
- }
- fprintf(fp_pairs,"\n");
- }
- teller_disre++;
- }
- }
- else if (bDHDL)
- {
- do_dhdl(fr, &fp_dhdl, opt2fn("-odh",NFILE,fnm),
- &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas,
- oenv);
- }
- /*******************************************
- * E N E R G I E S
- *******************************************/
- else {
- if (fr->nre > 0) {
- if (bPrAll)
- {
- /* We skip frames with single points (usually only the first frame),
- * since they would result in an average plot with outliers.
- */
- if (fr->nsum > 1) {
- print_time(out,fr->t);
- print1(out,bDp,fr->ener[set[0]].e);
- print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
- print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
- fprintf(out,"\n");
- }
- }
- else
- {
- print_time(out,fr->t);
- if (bSum)
- {
- sum = 0;
- for(i=0; i<nset; i++)
- {
- sum += fr->ener[set[i]].e;
- }
- print1(out,bDp,sum/nmol-ezero);
- }
- else
- {
- for(i=0; (i<nset); i++)
- {
- if (bIsEner[i])
- {
- print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
- }
- else
- {
- print1(out,bDp,fr->ener[set[i]].e);
- }
- }
- }
- fprintf(out,"\n");
- }
- }
-#if 0
- /* we first count the blocks that have id 0: the orire blocks */
- block_orire=0;
- for(b=0;b<fr->nblock;b++)
- {
- if (fr->block[b].id == mde_block_type_orire)
- nblock_orire++;
- }
-#endif
+ /* Initiate energies and set them to zero */
+ edat.nsteps = 0;
+ edat.npoints = 0;
+ edat.nframes = 0;
+ edat.step = NULL;
+ edat.steps = NULL;
+ edat.points = NULL;
+ snew(edat.s,nset);
+
+ /* Initiate counters */
+ teller = 0;
+ teller_disre = 0;
+ bFoundStart = FALSE;
+ start_step = 0;
+ start_t = 0;
+ do {
+ /* This loop searches for the first frame (when -b option is given),
+ * or when this has been found it reads just one energy frame
+ */
+ do {
+ bCont = do_enx(fp,&(frame[NEXT]));
+ if (bCont) {
+ timecheck = check_times(frame[NEXT].t);
+ }
+ } while (bCont && (timecheck < 0));
+
+ if ((timecheck == 0) && bCont) {
+ /* We read a valid frame, so we can use it */
+ fr = &(frame[NEXT]);
+
+ if (fr->nre > 0) {
+ /* The frame contains energies, so update cur */
+ cur = NEXT;
+
+ if (edat.nframes % 1000 == 0)
+ {
+ srenew(edat.step,edat.nframes+1000);
+ memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
+ srenew(edat.steps,edat.nframes+1000);
+ memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
+ srenew(edat.points,edat.nframes+1000);
+ memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
+
+ for(i=0; i<nset; i++)
+ {
+ srenew(edat.s[i].ener,edat.nframes+1000);
+ memset(&(edat.s[i].ener[edat.nframes]),0,
+ 1000*sizeof(edat.s[i].ener[0]));
+ srenew(edat.s[i].es ,edat.nframes+1000);
+ memset(&(edat.s[i].es[edat.nframes]),0,
+ 1000*sizeof(edat.s[i].es[0]));
+ }
+ }
+
+ nfr = edat.nframes;
+ edat.step[nfr] = fr->step;
+
+ if (!bFoundStart)
+ {
+ bFoundStart = TRUE;
+ /* Initiate the previous step data */
+ start_step = fr->step;
+ start_t = fr->t;
+ /* Initiate the energy sums */
+ edat.steps[nfr] = 1;
+ edat.points[nfr] = 1;
+ for(i=0; i<nset; i++)
+ {
+ sss = set[i];
+ edat.s[i].es[nfr].sum = fr->ener[sss].e;
+ edat.s[i].es[nfr].sum2 = 0;
+ }
+ edat.nsteps = 1;
+ edat.npoints = 1;
+ }
+ else
+ {
+ edat.steps[nfr] = fr->nsteps;
+ {
+ if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
+ {
+ if (fr->nsum <= 1)
+ {
+ edat.points[nfr] = 1;
+ for(i=0; i<nset; i++)
+ {
+ sss = set[i];
+ edat.s[i].es[nfr].sum = fr->ener[sss].e;
+ edat.s[i].es[nfr].sum2 = 0;
+ }
+ edat.npoints += 1;
+ }
+ else
+ {
+ edat.points[nfr] = fr->nsum;
+ for(i=0; i<nset; i++)
+ {
+ sss = set[i];
+ edat.s[i].es[nfr].sum = fr->ener[sss].esum;
+ edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
+ }
+ edat.npoints += fr->nsum;
+ }
+ }
+ else
+ {
+ /* The interval does not match fr->nsteps:
+ * can not do exact averages.
+ */
+ edat.npoints = 0;
+ }
+ edat.nsteps = fr->step - start_step + 1;
+ }
+ }
+ for(i=0; i<nset; i++)
+ {
+ edat.s[i].ener[nfr] = fr->ener[set[i]].e;
+ }
+ }
+ /*
+ * Define distance restraint legends. Can only be done after
+ * the first frame has been read... (Then we know how many there are)
+ */
+ blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
+ if (bDisRe && bDRAll && !leg && blk_disre)
+ {
+ t_iatom *fa;
+ t_iparams *ip;
+
+ fa = top->idef.il[F_DISRES].iatoms;
+ ip = top->idef.iparams;
+ if (blk_disre->nsub != 2 ||
+ (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
+ {
+ gmx_incons("Number of disre sub-blocks not equal to 2");
+ }
+
+ ndisre=blk_disre->sub[0].nr ;
+ if (ndisre != top->idef.il[F_DISRES].nr/3)
+ {
+ gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
+ ndisre,top->idef.il[F_DISRES].nr/3);
+ }
+ snew(pairleg,ndisre);
+ for(i=0; i<ndisre; i++)
+ {
+ snew(pairleg[i],30);
+ j=fa[3*i+1];
+ k=fa[3*i+2];
+ gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
+ gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
+ sprintf(pairleg[i],"%d %s %d %s (%d)",
+ resnr_j,anm_j,resnr_k,anm_k,
+ ip[fa[3*i]].disres.label);
+ }
+ set=select_it(ndisre,pairleg,&nset);
+ snew(leg,2*nset);
+ for(i=0; (i<nset); i++)
+ {
+ snew(leg[2*i],32);
+ sprintf(leg[2*i], "a %s",pairleg[set[i]]);
+ snew(leg[2*i+1],32);
+ sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
+ }
+ xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
+ }
+
+ /*
+ * Store energies for analysis afterwards...
+ */
+ if (!bDisRe && !bDHDL && (fr->nre > 0)) {
+ if (edat.nframes % 1000 == 0) {
+ srenew(time,edat.nframes+1000);
+ }
+ time[edat.nframes] = fr->t;
+ edat.nframes++;
+ }
+ /*
+ * Printing time, only when we do not want to skip frames
+ */
+ if (!skip || teller % skip == 0) {
+ if (bDisRe) {
+ /*******************************************
+ * D I S T A N C E R E S T R A I N T S
+ *******************************************/
+ if (ndisre > 0)
+ {
+ #ifndef GMX_DOUBLE
+ float *disre_rt = blk_disre->sub[0].fval;
+ float *disre_rm3tav = blk_disre->sub[1].fval;
+ #else
+ double *disre_rt = blk_disre->sub[0].dval;
+ double *disre_rm3tav = blk_disre->sub[1].dval;
+ #endif
+
+ print_time(out,fr->t);
+ if (violaver == NULL)
+ snew(violaver,ndisre);
+
+ /* Subtract bounds from distances, to calculate violations */
+ calc_violations(disre_rt, disre_rm3tav,
+ nbounds,pair,bounds,violaver,&sumt,&sumaver);
+
+ fprintf(out," %8.4f %8.4f\n",sumaver,sumt);
+ if (bDRAll) {
+ print_time(fp_pairs,fr->t);
+ for(i=0; (i<nset); i++) {
+ sss=set[i];
+ fprintf(fp_pairs," %8.4f", mypow(disre_rm3tav[sss],minthird));
+ fprintf(fp_pairs," %8.4f", disre_rt[sss]);
+ }
+ fprintf(fp_pairs,"\n");
+ }
+ teller_disre++;
+ }
+ }
+ else if (bDHDL)
+ {
+ do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh",NFILE,fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
+ }
+
+ /*******************************************
+ * E N E R G I E S
+ *******************************************/
+ else {
+ if (fr->nre > 0) {
+ if (bPrAll)
+ {
+ /* We skip frames with single points (usually only the first frame),
+ * since they would result in an average plot with outliers.
+ */
+ if (fr->nsum > 1) {
+ print_time(out,fr->t);
+ print1(out,bDp,fr->ener[set[0]].e);
+ print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
+ print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
+ fprintf(out,"\n");
+ }
+ }
+ else
+ {
+ print_time(out,fr->t);
+ if (bSum)
+ {
+ sum = 0;
+ for(i=0; i<nset; i++)
+ {
+ sum += fr->ener[set[i]].e;
+ }
+ print1(out,bDp,sum/nmol-ezero);
+ }
+ else
+ {
+ for(i=0; (i<nset); i++)
+ {
+ if (bIsEner[i])
+ {
+ print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
+ }
+ else
+ {
+ print1(out,bDp,fr->ener[set[i]].e);
+ }
+ }
+ }
+ fprintf(out,"\n");
+ }
+ }
blk = find_block_id_enxframe(fr, enx_i, NULL);
if (bORIRE && blk)
{
vals=blk->sub[0].dval;
#endif
- if (blk->sub[0].nr != (size_t)nor)
+ if (blk->sub[0].nr != (size_t)nor)
gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
if (bORA || bODA)
{
for(i=0; i<nor; i++)
odrms[i] += sqr(vals[i]-oobs[i]);
}
- if (bORT)
+ if (bORT)
{
fprintf(fort," %10f",fr->t);
for(i=0; i<norsel; i++)
fprintf(fort," %g",vals[orsel[i]]);
fprintf(fort,"\n");
}
- if (bODT)
+ if (bODT)
{
fprintf(fodt," %10f",fr->t);
for(i=0; i<norsel; i++)
}
norfr++;
}
- blk = find_block_id_enxframe(fr, enxORT, NULL);
- if (bOTEN && blk)
- {
+ blk = find_block_id_enxframe(fr, enxORT, NULL);
+ if (bOTEN && blk)
+ {
#ifndef GMX_DOUBLE
- xdr_datatype dt=xdr_datatype_float;
+ xdr_datatype dt=xdr_datatype_float;
#else
- xdr_datatype dt=xdr_datatype_double;
+ xdr_datatype dt=xdr_datatype_double;
#endif
- real *vals;
-
- if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
- gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
+ real *vals;
+
+ if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
+ gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
#ifndef GMX_DOUBLE
- vals=blk->sub[0].fval;
+ vals=blk->sub[0].fval;
#else
- vals=blk->sub[0].dval;
+ vals=blk->sub[0].dval;
#endif
if (blk->sub[0].nr != (size_t)(nex*12))
teller++;
}
} while (bCont && (timecheck == 0));
-
+
fprintf(stderr,"\n");
close_enx(fp);
- if (out)
+ if (out)
ffclose(out);
if (bDRAll)
ffclose(fort);
if (bODT)
ffclose(fodt);
- if (bORA)
+ if (bORA)
{
out = xvgropen(opt2fn("-ora",NFILE,fnm),
"Average calculated orientations",
if (bOTEN)
ffclose(foten);
- if (bDisRe)
+ if (bDisRe)
{
analyse_disre(opt2fn("-viol",NFILE,fnm),
teller_disre,violaver,bounds,index,pair,nbounds,oenv);
- }
+ }
else if (bDHDL)
{
if (fp_dhdl)
{
ffclose(fp_dhdl);
- printf("\n\nWrote %d lambda values with %d samples as ",
+ printf("\n\nWrote %d lambda values with %d samples as ",
dh_lambdas, dh_samples);
if (dh_hists > 0)
{
nbmin,nbmax);
}
if (opt2bSet("-f2",NFILE,fnm)) {
- fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
+ fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
reftemp, nset, set, leg, &edat, time ,oenv);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff