C++ xvgrclose函数代码示例

void low_rmsd_dist(const char *fn, real maxrms, int nn, real **mat,
                   const gmx_output_env_t *oenv)
{
    FILE   *fp;
    int     i, j, *histo, x;
    real    fac;

    fac = 100/maxrms;
    snew(histo, 101);
    for (i = 0; i < nn; i++)
    {
        for (j = i+1; j < nn; j++)
        {
            x = static_cast<int>(fac*mat[i][j]+0.5);
            if (x <= 100)
            {
                histo[x]++;
            }
        }
    }

    fp = xvgropen(fn, "RMS Distribution", "RMS (nm)", "a.u.", oenv);
    for (i = 0; (i < 101); i++)
    {
        fprintf(fp, "%10g  %10d\n", i/fac, histo[i]);
    }
    xvgrclose(fp);
    sfree(histo);
}

void histogram(const char *distfile, real binwidth, int n, int nset, real **val,
               const gmx_output_env_t *oenv)
{
    FILE          *fp;
    int            i, s;
    double         minval, maxval;
    int            nbin;
    gmx_int64_t   *histo;

    minval = val[0][0];
    maxval = val[0][0];
    for (s = 0; s < nset; s++)
    {
        for (i = 0; i < n; i++)
        {
            if (val[s][i] < minval)
            {
                minval = val[s][i];
            }
            else if (val[s][i] > maxval)
            {
                maxval = val[s][i];
            }
        }
    }

    minval = binwidth*std::floor(minval/binwidth);
    maxval = binwidth*std::ceil(maxval/binwidth);
    if (minval != 0)
    {
        minval -= binwidth;
    }
    maxval += binwidth;

    nbin = static_cast<int>(((maxval - minval)/binwidth + 0.5) + 1);
    fprintf(stderr, "Making distributions with %d bins\n", nbin);
    snew(histo, nbin);
    fp = xvgropen(distfile, "Distribution", "", "", oenv);
    for (s = 0; s < nset; s++)
    {
        for (i = 0; i < nbin; i++)
        {
            histo[i] = 0;
        }
        for (i = 0; i < n; i++)
        {
            histo[static_cast<int>((val[s][i] - minval)/binwidth + 0.5)]++;
        }
        for (i = 0; i < nbin; i++)
        {
            fprintf(fp, " %g  %g\n", minval+i*binwidth, static_cast<double>(histo[i])/(n*binwidth));
        }
        if (s < nset-1)
        {
            fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
        }
    }
    xvgrclose(fp);
}

void histogram(const char *distfile, real binwidth, int n, int nset, real **val,
               const output_env_t oenv)
{
    FILE          *fp;
    int            i, s;
    double         min, max;
    int            nbin;
    gmx_int64_t   *histo;

    min = val[0][0];
    max = val[0][0];
    for (s = 0; s < nset; s++)
    {
        for (i = 0; i < n; i++)
        {
            if (val[s][i] < min)
            {
                min = val[s][i];
            }
            else if (val[s][i] > max)
            {
                max = val[s][i];
            }
        }
    }

    min = binwidth*floor(min/binwidth);
    max = binwidth*ceil(max/binwidth);
    if (min != 0)
    {
        min -= binwidth;
    }
    max += binwidth;

    nbin = (int)((max - min)/binwidth + 0.5) + 1;
    fprintf(stderr, "Making distributions with %d bins\n", nbin);
    snew(histo, nbin);
    fp = xvgropen(distfile, "Distribution", "", "", oenv);
    for (s = 0; s < nset; s++)
    {
        for (i = 0; i < nbin; i++)
        {
            histo[i] = 0;
        }
        for (i = 0; i < n; i++)
        {
            histo[(int)((val[s][i] - min)/binwidth + 0.5)]++;
        }
        for (i = 0; i < nbin; i++)
        {
            fprintf(fp, " %g  %g\n", min+i*binwidth, (double)histo[i]/(n*binwidth));
        }
        if (s < nset-1)
        {
            fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
        }
    }
    xvgrclose(fp);
}

static void do_ballistic(const char *balFile, int nData,
                         real *t, real **val, int nSet,
                         real balTime, int nBalExp,
                         const output_env_t oenv)
{
    double     **ctd   = NULL, *td = NULL;
    t_gemParams *GP    = init_gemParams(0, 0, t, nData, 0, 0, 0, balTime, nBalExp);
    static char *leg[] = {"Ac'(t)"};
    FILE        *fp;
    int          i, set;

    if (GP->ballistic/GP->tDelta >= GP->nExpFit*2+1)
    {
        snew(ctd, nSet);
        snew(td,  nData);

        fp = xvgropen(balFile, "Hydrogen Bond Autocorrelation", "Time (ps)", "C'(t)", oenv);
        xvgr_legend(fp, asize(leg), (const char**)leg, oenv);

        for (set = 0; set < nSet; set++)
        {
            snew(ctd[set], nData);
            for (i = 0; i < nData; i++)
            {
                ctd[set][i] = (double)val[set][i];
                if (set == 0)
                {
                    td[i] = (double)t[i];
                }
            }

            takeAwayBallistic(ctd[set], td, nData, GP->ballistic, GP->nExpFit, GP->bDt);
        }

        for (i = 0; i < nData; i++)
        {
            fprintf(fp, "  %g", t[i]);
            for (set = 0; set < nSet; set++)
            {
                fprintf(fp, "  %g", ctd[set][i]);
            }
            fprintf(fp, "\n");
        }


        for (set = 0; set < nSet; set++)
        {
            sfree(ctd[set]);
        }
        sfree(ctd);
        sfree(td);
        xvgrclose(fp);
    }
    else
    {
        printf("Number of data points is less than the number of parameters to fit\n."
               "The system is underdetermined, hence no ballistic term can be found.\n\n");
    }
}

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 = std::min(enerT[0][j], bmin);
        bmax = std::max(enerT[0][j], bmax);
    }
    bwidth  = 1.0;
    blength = static_cast<int>((bmax - bmin)/bwidth + 2);
    snew(histo, nbin);
    for (i = 0; (i < nbin); i++)
    {
        snew(histo[i], blength);
    }
    for (j = 0; (j < n); j++)
    {
        k = static_cast<int>((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");
    }
    xvgrclose(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 gmx_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);
}

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                   = static_cast<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");
        }
    }

    xvgrclose (fp);
}

void
AbstractPlotModule::Impl::closeFile()
{
    if (fp_ != NULL)
    {
        if (bPlain_)
        {
            gmx_fio_fclose(fp_);
        }
        else
        {
            xvgrclose(fp_);
        }
        fp_ = NULL;
    }
}

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");
  }
  xvgrclose(fp);
}

static void dump_w(output_env_t oenv, real beta)
{
    FILE       *fp;
    double      nu;
    const char *leg[] = { "wCv", "wS", "wA", "wE" };

    fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
                  "w", oenv);
    xvgr_legend(fp, asize(leg), leg, oenv);
    for (nu = 1; (nu < 100); nu += 0.05)
    {
        fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n", beta*PLANCK*nu,
                wCsolid(nu, beta), wSsolid(nu, beta),
                wAsolid(nu, beta), wEsolid(nu, beta));
    }
    xvgrclose(fp);
}

void print_one(const output_env_t oenv, const char *base, const char *name,
               const char *title, const char *ylabel, int nf, real time[],
               real data[])
{
    FILE *fp;
    char  buf[256], t2[256];
    int   k;

    sprintf(buf, "%s%s.xvg", base, name);
    fprintf(stderr, "\rPrinting %s  ", buf);
    sprintf(t2, "%s %s", title, name);
    fp = xvgropen(buf, t2, "Time (ps)", ylabel, oenv);
    for (k = 0; (k < nf); k++)
    {
        fprintf(fp, "%10g  %10g\n", time[k], data[k]);
    }
    xvgrclose(fp);
}

static void analyse_em_all(int npdb, t_pdbfile *pdbf[], const char *edocked,
                           const char *efree, const gmx_output_env_t *oenv)
{
    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)", oenv);
        for (i = 0; (i < npdb); i++)
        {
            fprintf(fp, "%12lf\n", bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
        }
        xvgrclose(fp);
    }
}