自学教程：C++ xvgropen函数代码示例

static void print_histo(char *fn,int nhisto,int histo[],real binwidth){  FILE *fp;  int i;    fp = xvgropen(fn,"Velocity distribution","V (nm/ps)","arbitrary units");  for(i=0; (i<nhisto); i++)     fprintf(fp,"%10.3e  %10d/n",i*binwidth,histo[i]);  fclose(fp);}

void plot_potential(double *potential[], double *charge[], double *field[],                    const char *afile, const char *bfile, const char *cfile,                    int nslices, int nr_grps, const char *grpname[], double slWidth,                    const output_env_t oenv){    FILE       *pot,     /* xvgr file with potential */    *cha,                /* xvgr file with charges   */    *fie;                /* xvgr files with fields   */    char       buf[256]; /* for xvgr title */    int        slice, n;    sprintf(buf, "Electrostatic Potential");    pot = xvgropen(afile, buf, "Box (nm)", "Potential (V)", oenv);    xvgr_legend(pot, nr_grps, grpname, oenv);    sprintf(buf, "Charge Distribution");    cha = xvgropen(bfile, buf, "Box (nm)", "Charge density (q/nm//S3//N)", oenv);    xvgr_legend(cha, nr_grps, grpname, oenv);    sprintf(buf, "Electric Field");    fie = xvgropen(cfile, buf, "Box (nm)", "Field (V/nm)", oenv);    xvgr_legend(fie, nr_grps, grpname, oenv);    for (slice = cb; slice < (nslices - ce); slice++)    {        fprintf(pot, "%20.16g  ", slice * slWidth);        fprintf(cha, "%20.16g  ", slice * slWidth);        fprintf(fie, "%20.16g  ", slice * slWidth);        for (n = 0; n < nr_grps; n++)        {            fprintf(pot, "   %20.16g", potential[n][slice]);            fprintf(fie, "   %20.16g", field[n][slice]/1e9); /* convert to V/nm */            fprintf(cha, "   %20.16g", charge[n][slice]);        }        fprintf(pot, "/n");        fprintf(cha, "/n");        fprintf(fie, "/n");    }    xvgrclose(pot);    xvgrclose(cha);    xvgrclose(fie);}

static void do_geminate(const char *gemFile, int nData,                        real *t, real **val, int nSet,                        const real D, const real rcut, const real balTime,                        const int nFitPoints, const real begFit, const real endFit,                        const output_env_t oenv){    double     **ctd = NULL, **ctdGem = NULL, *td = NULL;    t_gemParams *GP  = init_gemParams(rcut, D, t, nData, nFitPoints,                                      begFit, endFit, balTime, 1);    const char  *leg[] = {"Ac//sgem//N(t)"};    FILE        *fp;    int          i, set;    snew(ctd,    nSet);    snew(ctdGem, nSet);    snew(td,  nData);    fp = xvgropen(gemFile, "Hydrogen Bond Autocorrelation", "Time (ps)", "C'(t)", oenv);    xvgr_legend(fp, asize(leg), leg, oenv);    for (set = 0; set < nSet; set++)    {        snew(ctd[set],    nData);        snew(ctdGem[set], nData);        for (i = 0; i < nData; i++)        {            ctd[set][i] = (double)val[set][i];            if (set == 0)            {                td[i] = (double)t[i];            }        }        fitGemRecomb(ctd[set], td, &(ctd[set]), nData, GP);    }    for (i = 0; i < nData; i++)    {        fprintf(fp, "  %g", t[i]);        for (set = 0; set < nSet; set++)        {            fprintf(fp, "  %g", ctdGem[set][i]);        }        fprintf(fp, "/n");    }    for (set = 0; set < nSet; set++)    {        sfree(ctd[set]);        sfree(ctdGem[set]);    }    sfree(ctd);    sfree(ctdGem);    sfree(td);    xvgrclose(fp);}

static void calc_spectrum(int n, real c[], real dt, const char *fn,                          gmx_output_env_t *oenv, gmx_bool bRecip){    FILE     *fp;    gmx_fft_t fft;    int       i, status;    real     *data;    real      nu, omega, recip_fac;    snew(data, n*2);    for (i = 0; (i < n); i++)    {        data[i] = c[i];    }    if ((status = gmx_fft_init_1d_real(&fft, n, GMX_FFT_FLAG_NONE)) != 0)    {        gmx_fatal(FARGS, "Invalid fft return status %d", status);    }    if ((status = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX, data, data)) != 0)    {        gmx_fatal(FARGS, "Invalid fft return status %d", status);    }    fp = xvgropen(fn, "Vibrational Power Spectrum",                  bRecip ? "//f{12}w//f{4} (cm//S-1//N)" :                  "//f{12}n//f{4} (ps//S-1//N)",                  "a.u.", oenv);    /* This is difficult.     * The length of the ACF is dt (as passed to this routine).     * We pass the vacf with N time steps from 0 to dt.     * That means that after FFT we have lowest frequency = 1/dt     * then 1/(2 dt) etc. (this is the X-axis of the data after FFT).     * To convert to 1/cm we need to have to realize that     * E = hbar w = h nu = h c/lambda. We want to have reciprokal cm     * on the x-axis, that is 1/lambda, so we then have     * 1/lambda = nu/c. Since nu has units of 1/ps and c has gromacs units     * of nm/ps, we need to multiply by 1e7.     * The timestep between saving the trajectory is     * 1e7 is to convert nanometer to cm     */    recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;    for (i = 0; (i < n); i += 2)    {        nu    = i/(2*dt);        omega = nu*recip_fac;        /* Computing the square magnitude of a complex number, since this is a power         * spectrum.         */        fprintf(fp, "%10g  %10g/n", omega, gmx::square(data[i])+gmx::square(data[i+1]));    }    xvgrclose(fp);    gmx_fft_destroy(fft);    sfree(data);}

void dump_ana_struct(char *rmax,char *nion,char *gyr_com,char *gyr_origin,		     t_ana_struct *anal,int nsim){  FILE *fp,*gp,*hp,*kp;  int  i,j;  real t,d2;  char *legend[] = { "Rg", "RgX", "RgY", "RgZ" };    fp = xvgropen(rmax,"rmax","Time (fs)","r (nm)");  gp = xvgropen(nion,"N ion","Time (fs)","N ions");  hp = xvgropen(gyr_com,"Radius of gyration wrt. C.O.M.",		"Time (fs)","Rg (nm)");  xvgr_legend(hp,asize(legend),legend);  kp = xvgropen(gyr_origin,"Radius of gyration wrt. Origin",		"Time (fs)","Rg (nm)");  xvgr_legend(kp,asize(legend),legend);  for(i=0; (i<anal->index); i++) {    t = 1000*anal->t[i];    fprintf(fp,"%12g  %10.3f/n",t,anal->maxdist[i]);    fprintf(gp,"%12g  %10.3f/n",t,(1.0*anal->nion[i])/nsim-1);    d2 = anal->d2_com[i][XX] + anal->d2_com[i][YY] + anal->d2_com[i][ZZ];    fprintf(hp,"%12g  %10.3f  %10.3f  %10.3f  %10.3f/n",	    t,sqrt(d2/nsim),	    sqrt(anal->d2_com[i][XX]/nsim),	    sqrt(anal->d2_com[i][YY]/nsim),	    sqrt(anal->d2_com[i][ZZ]/nsim));    d2 = anal->d2_origin[i][XX] + anal->d2_origin[i][YY] + anal->d2_origin[i][ZZ];    fprintf(kp,"%12g  %10.3f  %10.3f  %10.3f  %10.3f/n",	    t,sqrt(d2/nsim),	    sqrt(anal->d2_origin[i][XX]/nsim),	    sqrt(anal->d2_origin[i][YY]/nsim),	    sqrt(anal->d2_origin[i][ZZ]/nsim));  }  ffclose(hp);  ffclose(gp);  ffclose(fp);  ffclose(kp);}

static void ana_trans(t_clusters *clust, int nf, 		      char *transfn, char *ntransfn, FILE *log,		      t_rgb rlo,t_rgb rhi){  FILE *fp;  real **trans,*axis;  int  *ntrans;  int  i,ntranst,maxtrans;  char buf[STRLEN];    snew(ntrans,clust->ncl);  snew(trans,clust->ncl);  snew(axis,clust->ncl);  for(i=0; i<clust->ncl; i++) {    axis[i]=i+1;    snew(trans[i],clust->ncl);  }  ntranst=0;  maxtrans=0;  for(i=1; i<nf; i++)    if(clust->cl[i] != clust->cl[i-1]) {      ntranst++;      ntrans[clust->cl[i-1]-1]++;      ntrans[clust->cl[i]-1]++;      trans[clust->cl[i-1]-1][clust->cl[i]-1]++;      maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);    }  ffprintf2(stderr,log,buf,"Counted %d transitions in total, "	    "max %d between two specific clusters/n",ntranst,maxtrans);  if (transfn) {    fp=ffopen(transfn,"w");    i = min(maxtrans+1, 80);    write_xpm(fp,0,"Cluster Transitions","# transitions",	      "from cluster","to cluster", 	      clust->ncl, clust->ncl, axis, axis, trans, 	      0, maxtrans, rlo, rhi, &i);    ffclose(fp);  }  if (ntransfn) {    fp=xvgropen(ntransfn,"Cluster Transitions","Cluster #","# transitions");    for(i=0; i<clust->ncl; i++)      fprintf(fp,"%5d %5d/n",i+1,ntrans[i]);    ffclose(fp);  }  sfree(ntrans);  for(i=0; i<clust->ncl; i++)    sfree(trans[i]);  sfree(trans);  sfree(axis);}

void correlate_aniso(char *fn,t_atoms *ref,t_atoms *calc){  FILE *fp;  int  i,j;    fp = xvgropen(fn,"Correlation between X-Ray and Computed Uij","X-Ray","Computed");  for(i=0; (i<ref->nr); i++) {    if (ref->pdbinfo[i].bAnisotropic) {      for(j=U11; (j<=U23); j++)	fprintf(fp,"%10d  %10d/n",ref->pdbinfo[i].uij[j],calc->pdbinfo[i].uij[j]);    }  }  fclose(fp);}

static void print_histo(const char *fn, int nhisto, int histo[], real binwidth,                        const gmx_output_env_t *oenv){    FILE *fp;    int   i;    fp = xvgropen(fn, "Velocity distribution", "V (nm/ps)", "arbitrary units",                  oenv);    for (i = 0; (i < nhisto); i++)    {        fprintf(fp, "%10.3e  %10d/n", i*binwidth, histo[i]);    }    xvgrclose(fp);}

extern void save_data (structure_factor_t *sft, const char *file, int ngrps,                       real start_q, real end_q, const output_env_t oenv){    FILE             *fp;    int               i, g = 0;    double           *tmp, polarization_factor, A;    structure_factor *sf = (structure_factor *)sft;    fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",                   "Intensity (a.u.)", oenv);    snew (tmp, ngrps);    for (g = 0; g < ngrps; g++)    {        for (i = 0; i < sf->n_angles; i++)        {/* *          theta is half the angle between incoming and scattered vectors. * *          polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta) * *          sin(theta) = q/(2k) := A  ->  sin^2(theta) = 4*A^2 (1-A^2) -> *          -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2) */            A                   = (double) (i * sf->ref_k) / (2.0 * sf->momentum);            polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));            sf->F[g][i]        *= polarization_factor;        }    }    for (i = 0; i < sf->n_angles; i++)    {        if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q)        {            fprintf (fp, "%10.5f  ", i * sf->ref_k);            for (g = 0; g < ngrps; g++)            {                fprintf (fp, "  %10.5f ", (sf->F[g][i]) /( sf->total_n_atoms*                                                           sf->nSteps));            }            fprintf (fp, "/n");        }    }    ffclose (fp);}

int gmx_rama(int argc,char *argv[]){  const char *desc[] = {    "[TT]g_rama[tt] selects the [GRK]phi[grk]/[GRK]psi[grk] dihedral combinations from your topology file",    "and computes these as a function of time.",    "Using simple Unix tools such as [IT]grep[it] you can select out",     "specific residues."  };  FILE      *out;  t_xrama   *xr;  int       j;  output_env_t oenv;  t_filenm  fnm[] = {    { efTRX, "-f", NULL,  ffREAD },    { efTPX, NULL, NULL,  ffREAD },    { efXVG, NULL, "rama",ffWRITE }  };#define NFILE asize(fnm)  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,		    NFILE,fnm,0,NULL,asize(desc),desc,0,NULL,&oenv);		        snew(xr,1);  init_rama(oenv,ftp2fn(efTRX,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),xr,3);    out=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Ramachandran Plot","Phi","Psi",oenv);  xvgr_line_props(out,0,elNone,ecFrank,oenv);  xvgr_view(out,0.2,0.2,0.8,0.8,oenv);  xvgr_world(out,-180,-180,180,180,oenv);  fprintf(out,"@    xaxis  tick on/[email
C++ xvid_encore函数代码示例
C++ xvgrclose函数代码示例