自学教程：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);}

voidAbstractPlotModule::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);    }}

static void corr_print(t_corr *curr, gmx_bool bTen, const char *fn, const char *title,                       const char *yaxis,                       real msdtime, real beginfit, real endfit,                       real *DD, real *SigmaD, char *grpname[],                       const output_env_t oenv){    FILE *out;    int   i, j;    out = xvgropen(fn, title, output_env_get_xvgr_tlabel(oenv), yaxis, oenv);    if (DD)    {        fprintf(out, "# MSD gathered over %g %s with %d restarts/n",                msdtime, output_env_get_time_unit(oenv), curr->nrestart);        fprintf(out, "# Diffusion constants fitted from time %g to %g %s/n",                beginfit, endfit, output_env_get_time_unit(oenv));        for (i = 0; i < curr->ngrp; i++)        {            fprintf(out, "# D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)/n",                    grpname[i], DD[i], SigmaD[i]);        }    }    for (i = 0; i < curr->nframes; i++)    {        fprintf(out, "%10g", output_env_conv_time(oenv, curr->time[i]));        for (j = 0; j < curr->ngrp; j++)        {            fprintf(out, "  %10g", curr->data[j][i]);            if (bTen)            {                fprintf(out, " %10g %10g %10g %10g %10g %10g",                        curr->datam[j][i][XX][XX],                        curr->datam[j][i][YY][YY],                        curr->datam[j][i][ZZ][ZZ],                        curr->datam[j][i][YY][XX],                        curr->datam[j][i][ZZ][XX],                        curr->datam[j][i][ZZ][YY]);            }        }        fprintf(out, "/n");    }    xvgrclose(out);}

static void dump_fy(output_env_t oenv, real toler){    FILE       *fp;    double      Delta, f, y, DD;    const char *leg[] = { "f", "fy", "y" };    DD = pow(10.0, 0.125);    fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);    xvgr_legend(fp, asize(leg), leg, oenv);    fprintf(fp, "@    world 1e-05, 0, 1000, 1/n");    fprintf(fp, "@    xaxes scale Logarithmic/n");    for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)    {        f = calc_fluidicity(Delta, toler);        y = calc_y(f, Delta, toler);        fprintf(fp, "%10g  %10g  %10g  %10g/n", Delta, f, f*y, y);    }    xvgrclose(fp);}

static void correlate_aniso(const char *fn, t_atoms *ref, t_atoms *calc,                            const gmx_output_env_t *oenv){    FILE *fp;    int   i, j;    fp = xvgropen(fn, "Correlation between X-Ray and Computed Uij", "X-Ray",                  "Computed", oenv);    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]);            }        }    }    xvgrclose(fp);}

static void plot_coscont(const char *ccfile, int n, int nset, real **val,                         const gmx_output_env_t *oenv){    FILE *fp;    int   s;    real  cc;    fp = xvgropen(ccfile, "Cosine content", "set / half periods", "cosine content",                  oenv);    for (s = 0; s < nset; s++)    {        cc = cosine_content(s+1, n, val[s]);        fprintf(fp, " %d %g/n", s+1, cc);        fprintf(stdout, "Cosine content of set %d with %.1f periods: %g/n",                s+1, 0.5*(s+1), cc);    }    fprintf(stdout, "/n");    xvgrclose(fp);}

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

void h2order_plot(rvec dipole[], real order[], const char *afile,                  int nslices, real slWidth, const gmx_output_env_t *oenv){    FILE       *ord;              /* xvgr files with order parameters  */    int         slice;            /* loop index     */    char        buf[256];         /* for xvgr title */    real        factor;           /* conversion to Debye from electron*nm */    /*  factor = 1e-9*1.60217733e-19/3.336e-30 */    factor = 1.60217733/3.336e-2;    fprintf(stderr, "%d slices/n", nslices);    sprintf(buf, "Water orientation with respect to normal");    ord = xvgropen(afile, buf,                   "box (nm)", "mu_x, mu_y, mu_z (D), cosine with normal", oenv);    for (slice = 0; slice < nslices; slice++)    {        fprintf(ord, "%8.3f %8.3f %8.3f %8.3f %e/n", slWidth*slice,                factor*dipole[slice][XX], factor*dipole[slice][YY],                factor*dipole[slice][ZZ], order[slice]);    }    xvgrclose(ord);}

//.........这里部分代码省略.........        print_dir(stdout, Uaver);        fp = gmx_ffopen(dirfn, "w");        print_dir(fp, Uaver);        gmx_ffclose(fp);    }    for (i = 0; i < isize; i++)    {        sfree(U[i]);    }    sfree(U);    /* Write RMSF output */    if (bReadPDB)    {        bfac = 8.0*M_PI*M_PI/3.0*100;        fp   = xvgropen(ftp2fn(efXVG, NFILE, fnm), "B-Factors",                        label, "(A//b//S//So//N//S2//N)", oenv);        xvgr_legend(fp, 2, leg, oenv);        for (i = 0; (i < isize); i++)        {            if (!bRes || i+1 == isize ||                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)            {                resind    = top.atoms.atom[index[i]].resind;                pdb_bfac  = find_pdb_bfac(pdbatoms, &top.atoms.resinfo[resind],                                          *(top.atoms.atomname[index[i]]));                fprintf(fp, "%5d  %10.5f  %10.5f/n",                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, rmsf[i]*bfac,                        pdb_bfac);            }        }        xvgrclose(fp);    }    else    {        fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS fluctuation", label, "(nm)", oenv);        for (i = 0; i < isize; i++)        {            if (!bRes || i+1 == isize ||                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)            {                fprintf(fp, "%5d %8.4f/n",                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, std::sqrt(rmsf[i]));            }        }        xvgrclose(fp);    }    for (i = 0; i < isize; i++)    {        pdbatoms->pdbinfo[index[i]].bfac = 800*M_PI*M_PI/3.0*rmsf[i];    }    if (devfn)    {        for (i = 0; i < isize; i++)        {            rmsf[i] = (rmsd_x[i][XX]+rmsd_x[i][YY]+rmsd_x[i][ZZ])/count;        }        if (bRes)        {            average_residues(rmsf, nullptr, 0, isize, index, w_rls, &top.atoms);        }        /* Write RMSD output */

static void process_tcaf(int nframes, real dt, int nkc, real **tc, rvec *kfac,                         real rho, real wt, const char *fn_trans,                         const char *fn_tca, const char *fn_tc,                         const char *fn_tcf, const char *fn_cub,                         const char *fn_vk, const gmx_output_env_t *oenv){    FILE  *fp, *fp_vk, *fp_cub = NULL;    int    nk, ntc;    real **tcaf, **tcafc = NULL, eta, *sig;    int    i, j, k, kc;    int    ncorr;    double fitparms[3];    nk  = kset_c[nkc];    ntc = nk*NPK;    if (fn_trans)    {        fp = xvgropen(fn_trans, "Transverse Current", "Time (ps)", "TC (nm/ps)",                      oenv);        for (i = 0; i < nframes; i++)        {            fprintf(fp, "%g", i*dt);            for (j = 0; j < ntc; j++)            {                fprintf(fp, " %g", tc[j][i]);            }            fprintf(fp, "/n");        }        xvgrclose(fp);        do_view(oenv, fn_trans, "-nxy");    }    ncorr = (nframes+1)/2;    if (ncorr > static_cast<int>(5*wt/dt+0.5))    {        ncorr = static_cast<int>(5*wt/dt+0.5)+1;    }    snew(tcaf, nk);    for (k = 0; k < nk; k++)    {        snew(tcaf[k], ncorr);    }    if (fn_cub)    {        snew(tcafc, nkc);        for (k = 0; k < nkc; k++)        {            snew(tcafc[k], ncorr);        }    }    snew(sig, ncorr);    for (i = 0; i < ncorr; i++)    {        sig[i] = std::exp(0.5*i*dt/wt);    }    low_do_autocorr(fn_tca, oenv, "Transverse Current Autocorrelation Functions",                    nframes, ntc, ncorr, tc, dt, eacNormal,                    1, FALSE, FALSE, FALSE, 0, 0, 0);    do_view(oenv, fn_tca, "-nxy");    fp = xvgropen(fn_tc, "Transverse Current Autocorrelation Functions",                  "Time (ps)", "TCAF", oenv);    for (i = 0; i < ncorr; i++)    {        kc = 0;        fprintf(fp, "%g", i*dt);        for (k = 0; k < nk; k++)        {            for (j = 0; j < NPK; j++)            {                tcaf[k][i] += tc[NPK*k+j][i];            }            if (fn_cub)            {                for (j = 0; j < NPK; j++)                {                    tcafc[kc][i] += tc[NPK*k+j][i];                }            }            if (i == 0)            {                fprintf(fp, " %g", 1.0);            }            else            {                tcaf[k][i] /= tcaf[k][0];                fprintf(fp, " %g", tcaf[k][i]);            }            if (k+1 == kset_c[kc+1])            {                kc++;            }        }        fprintf(fp, "/n");    }    xvgrclose(fp);    do_view(oenv, fn_tc, "-nxy");//.........这里部分代码省略.........

//.........这里部分代码省略.........        sum_gyro_tot += sum_gyro;        if (outp)        {            i = -1;            for (j = 0; j < nat_min/2; j += 2)            {                sum_inp[j] /= ninp[j];                if (i == -1 && sum_inp[j] <= std::exp(-1.0))                {                    i = j;                }            }            if (i == -1)            {                pers = j;            }            else            {                /* Do linear interpolation on a log scale */                pers = i - 2.0                    + 2.0*(std::log(sum_inp[i-2]) + 1.0)/(std::log(sum_inp[i-2]) - std::log(sum_inp[i]));            }            fprintf(outp, "%10.3f %8.4f/n", t*output_env_get_time_factor(oenv), pers);            sum_pers_tot += pers;        }        frame++;    }    while (read_next_x(oenv, status, &t, x, box));    gmx_rmpbc_done(gpbc);    close_trx(status);    xvgrclose(out);    if (outv)    {        xvgrclose(outv);    }    if (outp)    {        xvgrclose(outp);    }    sum_eed2_tot /= frame;    sum_gyro_tot /= frame;    sum_pers_tot /= frame;    fprintf(stdout, "/nAverage end to end distance: %.3f (nm)/n",            std::sqrt(sum_eed2_tot));    fprintf(stdout, "/nAverage radius of gyration:  %.3f (nm)/n",            std::sqrt(sum_gyro_tot));    if (opt2bSet("-p", NFILE, fnm))    {        fprintf(stdout, "/nAverage persistence length:  %.2f bonds/n",                sum_pers_tot);    }    /* Handle printing of internal distances. */    if (outi)    {        if (output_env_get_print_xvgr_codes(oenv))        {            fprintf(outi, "@    xaxes scale Logarithmic/n");        }        ymax = -1;        ymin = 1e300;        j    = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */        for (i = 0; i < j; i++)        {            intd[i] /= (i + 1) * frame * nmol;            if (intd[i] > ymax)            {                ymax = intd[i];            }            if (intd[i] < ymin)            {                ymin = intd[i];            }        }        xvgr_world(outi, 1, ymin, j, ymax, oenv);        for (i = 0; i < j; i++)        {            fprintf(outi, "%d  %8.4f/n", i+1, intd[i]);        }        xvgrclose(outi);    }    do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");    if (opt2bSet("-v", NFILE, fnm))    {        do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");    }    if (opt2bSet("-p", NFILE, fnm))    {        do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");    }    return 0;}

//.........这里部分代码省略.........                }                if (bDatout)                {                    fprintf(datfp, "%12.7f %12.7f %12.7f/n", fr.time, R, kappa2);                }                if (bRhistout)                {                    rvalues[rkcount-1] = R;                    if (rkcount % allocblock == 0)                    {                        srenew(rvalues, allocblock*(rblocksallocated+1));                        rblocksallocated += 1;                    }                }                if (bKhistout)                {                    kappa2values[rkcount-1] = kappa2;                    if (rkcount % allocblock == 0)                    {                        srenew(kappa2values, allocblock*(kblocksallocated+1));                        kblocksallocated += 1;                    }                }                bHaveNextFrame = read_next_frame(oenv, status, &fr);            }            while (bHaveNextFrame);            if (bRKout)            {                xvgrclose(rkfp);            }            if (bDatout)            {                gmx_ffclose(datfp);            }            if (bInstEffout)            {                xvgrclose(iefp);            }            if (bRhistout)            {                printf("Writing R-Histogram/n");                rmin = rvalues[0];                rmax = rvalues[0];                for (i = 1; i < rkcount; i++)                {                    if (rvalues[i] < rmin)                    {                        rmin = rvalues[i];                    }                    else if (rvalues[i] > rmax)                    {                        rmax = rvalues[i];                    }                }                rmin -= histexpand;                rmax += histexpand;

//.........这里部分代码省略.........        printf("Going to do %d fourier transforms of length %d. Hang on./n",               gnx, nframes);    }    low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,                    FALSE, FALSE, -1, -1, 0, 0);    snew(dos, DOS_NR);    for (j = 0; (j < DOS_NR); j++)    {        snew(dos[j], nframes+4);    }    if (bVerbose)    {        printf("Going to merge the ACFs into the mass-weighted and plain ACF/n");    }    for (i = 0; (i < gnx); i += DIM)    {        mi = top.atoms.atom[i/DIM].m;        for (j = 0; (j < nframes/2); j++)        {            c1j            = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);            dos[VACF][j]  += c1j/Natom;            dos[MVACF][j] += mi*c1j;        }    }    fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",                  "Time (ps)", "C(t)", oenv);    snew(tt, nframes/2);    for (j = 0; (j < nframes/2); j++)    {        tt[j] = j*dt;        fprintf(fp, "%10g  %10g/n", tt[j], dos[VACF][j]);    }    xvgrclose(fp);    fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",                  "Time (ps)", "C(t)", oenv);    for (j = 0; (j < nframes/2); j++)    {        fprintf(fp, "%10g  %10g/n", tt[j], dos[MVACF][j]);    }    xvgrclose(fp);    if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,                                        GMX_FFT_FLAG_NONE)) != 0)    {        gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);    }    if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,                                   (void *)dos[MVACF], (void *)dos[DOS])) != 0)    {        gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);    }    /* First compute the DoS */    /* Magic factor of 8 included now. */    bfac = 8*dt*beta/2;    dos2 = 0;    snew(nu, nframes/4);    for (j = 0; (j < nframes/4); j++)    {        nu[j] = 2*j/(t1-t0);        dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);        if (bAbsolute)        {            dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));        }

//.........这里部分代码省略.........            j++;        }    }    nfitpnts = j;    integral = 0;    if (nfitpnts < effnNparams(eFitFn))    {        fprintf(stderr, "Not enough (%d) data points for fitting, dt = %g!/n",                nfitpnts, dt);    }    else    {        gmx_bool bSuccess;        if (bVerbose)        {            print_chi2_params(stdout, eFitFn, fitparms, "initial", nfitpnts, x, y);        }        initiate_fit_params(eFitFn, fitparms);        bSuccess = lmfit_exp(nfitpnts, x, y, dy, fitparms, bVerbose, eFitFn, fix);        extract_fit_params(eFitFn, fitparms);        if (!bSuccess)        {            fprintf(stderr, "Fit failed!/n");        }        else        {            if (bVerbose)            {                print_chi2_params(stdout, eFitFn, fitparms, "final", nfitpnts, x, y);            }            switch (eFitFn)            {                case effnEXP1:                    integral = fitparms[0]*myexp(begintimefit, 1,  fitparms[0]);                    break;                case effnEXP2:                    integral = fitparms[0]*myexp(begintimefit, fitparms[1],  fitparms[0]);                    break;                case effnEXPEXP:                    integral = (fitparms[0]*myexp(begintimefit, fitparms[1],  fitparms[0]) +                                fitparms[2]*myexp(begintimefit, 1-fitparms[1], fitparms[2]));                    break;                case effnEXP5:                case effnEXP7:                case effnEXP9:                    integral = 0;                    for (i = 0; (i < (effnNparams(eFitFn)-1)/2); i++)                    {                        integral += fitparms[2*i]*myexp(begintimefit, fitparms[2*i+1], fitparms[2*i]);                    }                    break;                default:                    /* Do numerical integration */                    integral = 0;                    for (i = 0; (i < nfitpnts-1); i++)                    {                        double y0 = lmcurves[eFitFn](x[i], fitparms);                        double y1 = lmcurves[eFitFn](x[i+1], fitparms);                        integral += (x[i+1]-x[i])*(y1+y0)*0.5;                    }                    break;            }            if (bVerbose)            {                printf("FIT: Integral of fitted function: %g/n", integral);                if ((effnEXP5 == eFitFn) || (effnEXP7 == eFitFn) || (effnEXP9 == eFitFn))                {                    printf("FIT: Note that the constant term is not taken into account when computing integral./n");                }            }            /* Generate debug output */            if (NULL != fn_fitted)            {                fp = xvgropen(fn_fitted, "Data + Fit", "Time (ps)",                              "Data (t)", oenv);                for (i = 0; (i < effnNparams(eFitFn)); i++)                {                    fprintf(fp, "# fitparms[%d] = %g/n", i, fitparms[i]);                }                for (j = 0; (j < nfitpnts); j++)                {                    real ttt = x0 ? x0[j] : dt*j;                    fprintf(fp, "%10.5e  %10.5e  %10.5e/n",                            x[j], y[j], (lmcurves[eFitFn])(ttt, fitparms));                }                xvgrclose(fp);            }        }    }    sfree(x);    sfree(y);    sfree(dy);    return integral;}

//.........这里部分代码省略.........                dy[j] = 1.0e-7;            }            else            {                dy[j] = sig[i];            }            if (debug)            {                fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g/n",                        j, i, x[j], y[j], dy[j]);            }            j++;        }    }    nfitpnts = j;    integral = 0;    if (nfitpnts < nparm)    {        fprintf(stderr, "Not enough data points for fitting!/n");    }    else    {        snew(parm, nparm);        snew(dparm, nparm);        if (fitparms)        {            for (i = 0; (i < nparm); i++)            {                parm[i] = fitparms[i];            }        }        if (!lmfit_exp(nfitpnts, x, y, dy, ftol, parm, dparm, bVerbose, eFitFn, fix))        {            fprintf(stderr, "Fit failed!/n");        }        else if (nparm <= 3)        {            /* Compute the integral from begintimefit */            if (nparm == 3)            {                integral = (parm[0]*myexp(begintimefit, parm[1],  parm[0]) +                            parm[2]*myexp(begintimefit, 1-parm[1], parm[2]));            }            else if (nparm == 2)            {                integral = parm[0]*myexp(begintimefit, parm[1],  parm[0]);            }            else if (nparm == 1)            {                integral = parm[0]*myexp(begintimefit, 1,  parm[0]);            }            else            {                gmx_fatal(FARGS, "nparm = %d in file %s, line %d",                          nparm, __FILE__, __LINE__);            }            /* Generate THE output */            if (bVerbose)            {                fprintf(stderr, "FIT: # points used in fit is: %d/n", nfitpnts);                fprintf(stderr, "FIT: %21s%21s%21s/n",                        "parm0     ", "parm1 (ps)   ", "parm2 (ps)    ");                fprintf(stderr, "FIT: ------------------------------------------------------------/n");                fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g/n",                        parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]);                fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g/n",                        begintimefit, integral);                sprintf(buf, "test%d.xvg", nfitpnts);                fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv);                fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g/n",                        parm[0], parm[1], parm[2]);                for (j = 0; (j < nfitpnts); j++)                {                    ttt = x0 ? x0[j] : dt*j;                    fprintf(fp, "%10.5e  %10.5e  %10.5e/n",                            ttt, c1[j], fit_function(eFitFn, parm, ttt));                }                xvgrclose(fp);            }        }        if (fitparms)        {            for (i = 0; (i < nparm); i++)            {                fitparms[i] = parm[i];            }        }        sfree(parm);        sfree(dparm);    }    sfree(x);    sfree(y);    sfree(dy);    return integral;}

//.........这里部分代码省略.........        }        pr->grn      = prframecurrent->grn;        pr->binwidth = prframecurrent->binwidth;        /* summ up gr and fill r */        for (i = 0; i < prframecurrent->grn; i++)        {            pr->gr[i] += prframecurrent->gr[i];            pr->r[i]   = prframecurrent->r[i];        }        /* normalize histo */        normalize_probability(prframecurrent->grn, prframecurrent->gr);        /* convert p(r) to sq */        sqframecurrent = convert_histogram_to_intensity_curve(prframecurrent, start_q, end_q, q_step);        /* print frame data if needed */        if (opt2fn_null("-prframe", NFILE, fnm))        {            snew(hdr, 25);            snew(suffix, GMX_PATH_MAX);            /* prepare header */            sprintf(hdr, "g(r), t = %f", t);            /* prepare output filename */            fnmdup = dup_tfn(NFILE, fnm);            sprintf(suffix, "-t%.2f", t);            add_suffix_to_output_names(fnmdup, NFILE, suffix);            fp = xvgropen(opt2fn_null("-prframe", NFILE, fnmdup), hdr, "Distance (nm)", "Probability", oenv);            for (i = 0; i < prframecurrent->grn; i++)            {                fprintf(fp, "%10.6f%10.6f/n", prframecurrent->r[i], prframecurrent->gr[i]);            }            done_filenms(NFILE, fnmdup);            fclose(fp);            sfree(hdr);            sfree(suffix);            sfree(fnmdup);        }        if (opt2fn_null("-sqframe", NFILE, fnm))        {            snew(hdr, 25);            snew(suffix, GMX_PATH_MAX);            /* prepare header */            sprintf(hdr, "I(q), t = %f", t);            /* prepare output filename */            fnmdup = dup_tfn(NFILE, fnm);            sprintf(suffix, "-t%.2f", t);            add_suffix_to_output_names(fnmdup, NFILE, suffix);            fp = xvgropen(opt2fn_null("-sqframe", NFILE, fnmdup), hdr, "q (nm^-1)", "s(q)/s(0)", oenv);            for (i = 0; i < sqframecurrent->qn; i++)            {                fprintf(fp, "%10.6f%10.6f/n", sqframecurrent->q[i], sqframecurrent->s[i]);            }            done_filenms(NFILE, fnmdup);            fclose(fp);            sfree(hdr);            sfree(suffix);            sfree(fnmdup);        }        /* free pr structure */        sfree(prframecurrent->gr);        sfree(prframecurrent->r);        sfree(prframecurrent);        /* free sq structure */        sfree(sqframecurrent->q);        sfree(sqframecurrent->s);        sfree(sqframecurrent);    }    while (read_next_x(oenv, status, &t, x, box));    close_trj(status);    /* normalize histo */    normalize_probability(pr->grn, pr->gr);    sq = convert_histogram_to_intensity_curve(pr, start_q, end_q, q_step);    /* prepare pr.xvg */    fp = xvgropen(opt2fn_null("-pr", NFILE, fnm), "G(r)", "Distance (nm)", "Probability", oenv);    for (i = 0; i < pr->grn; i++)    {        fprintf(fp, "%10.6f%10.6f/n", pr->r[i], pr->gr[i]);    }    xvgrclose(fp);    /* prepare sq.xvg */    fp = xvgropen(opt2fn_null("-sq", NFILE, fnm), "I(q)", "q (nm^-1)", "s(q)/s(0)", oenv);    for (i = 0; i < sq->qn; i++)    {        fprintf(fp, "%10.6f%10.6f/n", sq->q[i], sq->s[i]);    }    xvgrclose(fp);    /*     * Clean up memory     */    sfree(pr->gr);    sfree(pr->r);    sfree(pr);    sfree(sq->q);    sfree(sq->s);    sfree(sq);    please_cite(stdout, "Garmay2012");    return 0;}

//.........这里部分代码省略.........    trxnat = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);    nframes = 0;    rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;    rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;    gpbc = gmx_rmpbc_init(&top.idef, ePBC, trxnat);    if (bFrames)    {        out = opt2FILE("-frames", NFILE, fnm, "w");    }    do    {        gmx_rmpbc(gpbc, trxnat, box, x);        nframes++;        calc_mat(nres, natoms, rndx, x, index, truncate, mdmat, nmat, ePBC, box);        for (i = 0; (i < nres); i++)        {            for (j = 0; (j < natoms); j++)            {                if (nmat[i][j])                {                    totnmat[i][j]++;                }            }        }        for (i = 0; (i < nres); i++)        {            for (j = 0; (j < nres); j++)            {                totmdmat[i][j] += mdmat[i][j];            }        }        if (bFrames)        {            sprintf(label, "t=%.0f ps", t);            write_xpm(out, 0, label, "Distance (nm)", "Residue Index", "Residue Index",                      nres, nres, resnr, resnr, mdmat, 0, truncate, rlo, rhi, &nlevels);        }    }    while (read_next_x(oenv, status, &t, x, box));    fprintf(stderr, "/n");    close_trj(status);    gmx_rmpbc_done(gpbc);    if (bFrames)    {        gmx_ffclose(out);    }    fprintf(stderr, "Processed %d frames/n", nframes);    for (i = 0; (i < nres); i++)    {        for (j = 0; (j < nres); j++)        {            totmdmat[i][j] /= nframes;        }    }    write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean smallest distance",              "Distance (nm)", "Residue Index", "Residue Index",              nres, nres, resnr, resnr, totmdmat, 0, truncate, rlo, rhi, &nlevels);    if (bCalcN)    {        char **legend;        snew(legend, 5);        for (i = 0; i < 5; i++)        {            snew(legend[i], STRLEN);        }        tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);        fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),                      "Increase in number of contacts", "Residue", "Ratio", oenv);        sprintf(legend[0], "Total/mean");        sprintf(legend[1], "Total");        sprintf(legend[2], "Mean");        sprintf(legend[3], "# atoms");        sprintf(legend[4], "Mean/# atoms");        xvgr_legend(fp, 5, (const char**)legend, oenv);        for (i = 0; (i < nres); i++)        {            if (mean_n[i] == 0)            {                ratio = 1;            }            else            {                ratio = tot_n[i]/mean_n[i];            }            fprintf(fp, "%3d  %8.3f  %3d  %8.3f  %3d  %8.3f/n",                    i+1, ratio, tot_n[i], mean_n[i], natm[i], mean_n[i]/natm[i]);        }        xvgrclose(fp);    }    return 0;}