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 protected]    xaxis  tick major 60\[email protected]    xaxis  tick minor 30\n");
  fprintf(out,"@    yaxis  tick on\[email protected]    yaxis  tick major 60\[email protected]    yaxis  tick minor 30\n");
  fprintf(out,"@ s0 symbol 2\[email protected] s0 symbol size 0.4\[email protected] s0 symbol fill 1\n");
  
  j=0;
  do {
    plot_rama(out,xr);
    j++;
  } while (new_data(xr));
  fprintf(stderr,"\n");
  ffclose(out);
  
  do_view(oenv,ftp2fn(efXVG,NFILE,fnm),NULL);
  
  thanx(stderr);
  
  return 0;
}

static void analyse_em_all(int npdb,t_pdbfile *pdbf[],
			   char *edocked,char *efree)
{
  FILE *fp;
  int i;
  
  for(bFreeSort = FALSE; (bFreeSort <= TRUE); bFreeSort++) {
    qsort(pdbf,npdb,sizeof(pdbf[0]),pdbf_comp);
    fp = xvgropen(bFreeSort ? efree : edocked,
		  etitles[bFreeSort],"()","E (kJ/mol)");
    for(i=0; (i<npdb); i++)
      fprintf(fp,"%12lf\n",bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
    fclose(fp);
  }
}

static void ehisto(const char *fh,int n,real **enerT, const output_env_t oenv)
{
  FILE *fp;
  int  i,j,k,nbin,blength;
  int  *bindex;
  real *T,bmin,bmax,bwidth;
  int  **histo;
  
  bmin =  1e8;
  bmax = -1e8;
  snew(bindex,n);
  snew(T,n);
  nbin = 0;
  for(j=1; (j<n); j++) {
    for(k=0; (k<nbin); k++) {
      if (T[k] == enerT[1][j]) {
	bindex[j] = k;
	break;
      }
    }
    if (k == nbin) {
      bindex[j] = nbin;
      T[nbin]   = enerT[1][j];
      nbin++;
    }
    bmin = min(enerT[0][j],bmin);
    bmax = max(enerT[0][j],bmax);
}
  bwidth  = 1.0;
  blength = (bmax - bmin)/bwidth + 2;
  snew(histo,nbin);
  for(i=0; (i<nbin); i++) {
    snew(histo[i],blength);
  }
  for(j=0; (j<n); j++) {
    k = (enerT[0][j]-bmin)/bwidth;
    histo[bindex[j]][k]++;
  }
  fp = xvgropen(fh,"Energy distribution","E (kJ/mol)","",oenv);
  for(j=0; (j<blength); j++) {
    fprintf(fp,"%8.3f",bmin+j*bwidth);
    for(k=0; (k<nbin); k++) {
      fprintf(fp,"  %6d",histo[k][j]);
    }
    fprintf(fp,"\n");
  }  
  ffclose(fp);
}

void analyse_ss(char *outfile, t_matrix *mat, char *ss_string)
{
  FILE *fp;
  t_mapping *map;
  int s,f,r,*count,ss_count;
  char **leg;
  
  map=mat->map;
  snew(count,mat->nmap);
  snew(leg,mat->nmap+1);
  leg[0]="Structure";
  for(s=0; s<mat->nmap; s++)
    leg[s+1]=map[s].desc;
  
  fp=xvgropen(outfile,"Secondary Structure",
	      xvgr_tlabel(),"Number of Residues");
  if (bPrintXvgrCodes())
    fprintf(fp,"@ subtitle \"Structure = ");
  for(s=0; s<strlen(ss_string); s++) {
    if (s>0)
      fprintf(fp," + ");
    for(f=0; f<mat->nmap; f++)
      if (ss_string[s]==map[f].code.c1)
	fprintf(fp,"%s",map[f].desc);
  }
  fprintf(fp,"\"\n");
  xvgr_legend(fp,mat->nmap+1,leg);
  
  for(f=0; f<mat->nx; f++) {
    ss_count=0;
    for(s=0; s<mat->nmap; s++)
      count[s]=0;
    for(r=0; r<mat->ny; r++)
      count[mat->matrix[f][r]]++;
    for(s=0; s<mat->nmap; s++) {
      if (strchr(ss_string,map[s].code.c1))
	ss_count+=count[s];
    }
    fprintf(fp,"%8g %5d",mat->axis_x[f],ss_count);
    for(s=0; s<mat->nmap; s++)
      fprintf(fp," %5d",count[s]);
    fprintf(fp,"\n");
  }
  
  fclose(fp);
  sfree(leg);
  sfree(count);
}

void write_xvg(char *fn,char *title,int nx,int ny,real **y,char **leg)
{
  FILE *fp;
  int  i,j;
  
  fp=xvgropen(fn,title,"X","Y");
  if (leg)
    xvgr_legend(fp,ny-1,leg);
  for(i=0; (i<nx); i++) {
    for(j=0; (j<ny); j++) {
      fprintf(fp,"  %12.5e",y[j][i]);
    }
    fprintf(fp,"\n");
  }
  fclose(fp);
}

void
AbstractPlotModule::dataStarted(AbstractAnalysisData * /* data */)
{
    if (!impl_->filename_.empty())
    {
        if (impl_->bPlain_)
        {
            impl_->fp_ = gmx_fio_fopen(impl_->filename_.c_str(), "w");
        }
        else
        {
            time_unit_t  time_unit
                = static_cast<time_unit_t>(impl_->settings_.timeUnit() + 1);
            xvg_format_t xvg_format
                = (impl_->settings_.plotFormat() > 0
                   ? static_cast<xvg_format_t>(impl_->settings_.plotFormat())
                   : exvgNONE);
            output_env_t                  oenv;
            output_env_init(&oenv, getProgramContext(), time_unit, FALSE, xvg_format, 0);
            boost::shared_ptr<output_env> oenvGuard(oenv, &output_env_done);
            impl_->fp_ = xvgropen(impl_->filename_.c_str(), impl_->title_.c_str(),
                                  impl_->xlabel_.c_str(), impl_->ylabel_.c_str(),
                                  oenv);
            const SelectionCollection *selections
                = impl_->settings_.selectionCollection();
            if (selections != NULL)
            {
                selections->printXvgrInfo(impl_->fp_, oenv);
            }
            if (!impl_->subtitle_.empty())
            {
                xvgr_subtitle(impl_->fp_, impl_->subtitle_.c_str(), oenv);
            }
            if (output_env_get_print_xvgr_codes(oenv)
                && !impl_->legend_.empty())
            {
                std::vector<const char *> legend;
                legend.reserve(impl_->legend_.size());
                for (size_t i = 0; i < impl_->legend_.size(); ++i)
                {
                    legend.push_back(impl_->legend_[i].c_str());
                }
                xvgr_legend(impl_->fp_, legend.size(), &legend[0], oenv);
            }
        }
    }
}