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

static void make_dist_leg(FILE *fp,int gnx,atom_id index[],t_atoms *atoms){  char **leg;  int  i;  snew(leg,gnx/2);  for(i=0; i<gnx; i+=2) {    snew(leg[i/2],256);    sprintf(leg[i/2],"%s %s%d - %s %s%d",	    *(atoms->atomname[index[i]]),	    *(atoms->resname[atoms->atom[index[i]].resnr]),            atoms->atom[index[i]].resnr+1,	    *(atoms->atomname[index[i+1]]),	    *(atoms->resname[atoms->atom[index[i+1]].resnr]),            atoms->atom[index[i+1]].resnr+1);  }  xvgr_legend(fp,gnx/2,leg);  for(i=0; i<gnx/2; i++)    sfree(leg[i]);  sfree(leg);}

void dump_ana_ener(t_ana_ener *ae,int nsim,real dt,char *edump,		   t_ana_struct *total){  FILE *fp;  int  i,j;  real fac;    fac = 1.0/(nsim*ELECTRONVOLT);  fp=xvgropen(edump,"Energies","Time (fs)","E (eV)");  xvgr_legend(fp,eNR,enms);  fprintf(fp,"@ s%d legend /"Ek/Nelec/"/n",eNR);  fprintf(fp,"@ type nxy/n");  for(i=0; (i<ae->nx); i++) {    fprintf(fp,"%10f",1000.0*dt*i);    for(j=0; (j<eNR); j++)      fprintf(fp,"  %8.3f",ae->e[i][j]*fac);    fprintf(fp,"  %8.3f/n",ae->e[i][eELECTRON]/(ELECTRONVOLT*total->nion[i]));  }      fprintf(fp,"&/n");  ffclose(fp);}

void write_xvg(const char *fn, const char *title, int nx, int ny, real **y,               const char **leg, const 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);}

static void pr_dev(t_coupl_rec *tcr,		   real t,real dev[eoObsNR],t_commrec *cr,int nfile,t_filenm fnm[]){  static FILE *fp=NULL;  char   **ptr;  int    i,j;    if (!fp) {    fp=xvgropen(opt2fn("-devout",nfile,fnm),		"Deviations from target value","Time (ps)","");    snew(ptr,eoObsNR);    for(i=j=0; (i<eoObsNR); i++)      if (tcr->bObsUsed[i])	ptr[j++] = eoNames[i];    xvgr_legend(fp,j,ptr);    sfree(ptr);  }  fprintf(fp,"%10.3f",t);  for(i=0; (i<eoObsNR); i++)    if (tcr->bObsUsed[i])      fprintf(fp,"  %10.3e",dev[i]);  fprintf(fp,"/n");  fflush(fp);}

//.........这里部分代码省略.........                        "Min. angle: %f, Max. angle: %f/n", bond_min, bond_max);            }            rlo.r = 1; rlo.g = 1; rlo.b = 1;            rhi.r = 0; rhi.g = 0; rhi.b = 0;            sprintf(buf, "%s av. bond angle deviation", gn_rms[0]);            write_xpm(opt2FILE("-bm", NFILE, fnm, "w"), 0, buf, "degrees",                      output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,                      axis, axis2, bond_mat, bond_min, bond_max, rlo, rhi, &nlevels);        }    }    bAv = opt2bSet("-a", NFILE, fnm);    /* Write the RMSD's to file */    if (!bPrev)    {        sprintf(buf, "%s", whatxvgname[ewhat]);    }    else    {        sprintf(buf, "%s with frame %g %s ago", whatxvgname[ewhat],                time[prev*freq]-time[0], output_env_get_time_label(oenv));    }    fp = xvgropen(opt2fn("-o", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),                  whatxvglabel[ewhat], oenv);    if (output_env_get_print_xvgr_codes(oenv))    {        fprintf(fp, "@ subtitle /"%s%s after %s%s%s/"/n",                (nrms == 1) ? "" : "of ", gn_rms[0], fitgraphlabel[efit],                bFit     ? " to " : "", bFit ? gn_fit : "");    }    if (nrms != 1)    {        xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);    }    for (i = 0; (i < teller); i++)    {        if (bSplit && i > 0 &&            abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)        {            fprintf(fp, "&/n");        }        fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]);        for (j = 0; (j < nrms); j++)        {            fprintf(fp, " %12.7f", rls[j][i]);            if (bAv)            {                rlstot += rls[j][i];            }        }        fprintf(fp, "/n");    }    gmx_ffclose(fp);    if (bMirror)    {        /* Write the mirror RMSD's to file */        sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);        sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);        fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),                      buf2, oenv);        if (nrms == 1)        {            if (output_env_get_print_xvgr_codes(oenv))            {

//.........这里部分代码省略.........    {        fprintf(stderr, "/nOne of the lowest 6 eigenvalues has a non-zero value./n");        fprintf(stderr, "This could mean that the reference structure was not/n");        fprintf(stderr, "properly energy minimized./n");    }    /* now write the output */    fprintf (stderr, "Writing eigenvalues.../n");    out = xvgropen(opt2fn("-ol", NFILE, fnm),                   "Eigenvalues", "Eigenvalue index", "Eigenvalue [Gromacs units]",                   oenv);    if (output_env_get_print_xvgr_codes(oenv))    {        if (bM)        {            fprintf(out, "@ subtitle /"mass weighted/"/n");        }        else        {            fprintf(out, "@ subtitle /"not mass weighted/"/n");        }    }    for (i = 0; i <= (end-begin); i++)    {        fprintf (out, "%6d %15g/n", begin+i, eigenvalues[i]);    }    ffclose(out);    if (opt2bSet("-qc", NFILE, fnm))    {        qc = xvgropen(opt2fn("-qc", NFILE, fnm), "Quantum Corrections", "Eigenvector index", "", oenv);        xvgr_legend(qc, asize(qcleg), qcleg, oenv);        qcvtot = qutot = 0;    }    else    {        qc = NULL;    }    printf("Writing eigenfrequencies - negative eigenvalues will be set to zero./n");    out = xvgropen(opt2fn("-of", NFILE, fnm),                   "Eigenfrequencies", "Eigenvector index", "Wavenumber [cm//S-1//N]",                   oenv);    if (output_env_get_print_xvgr_codes(oenv))    {        if (bM)        {            fprintf(out, "@ subtitle /"mass weighted/"/n");        }        else        {            fprintf(out, "@ subtitle /"not mass weighted/"/n");        }    }    /* Spectrum ? */    spec = NULL;    if (opt2bSet("-os", NFILE, fnm) && (maxspec > 0))    {        snew(spectrum, maxspec);        spec = xvgropen(opt2fn("-os", NFILE, fnm),                        "Vibrational spectrum based on harmonic approximation",                        "//f{12}w//f{4} (cm//S-1//N)",                        "Intensity [Gromacs units]",                        oenv);

//.........这里部分代码省略.........    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;}

//.........这里部分代码省略.........            dos[DOS][j] = bfac*dos[DOS][2*j];        }    }    /* Normalize it */    dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);    if (bNormalizeDos)    {        for (j = 0; (j < nframes/4); j++)        {            dos[DOS][j] *= 3*Natom/dostot;        }    }    /* Now analyze it */    DoS0 = dos[DOS][0];    /* Note this eqn. is incorrect in Pascal2011a! */    Delta = ((2*DoS0/(9*Natom))*std::sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*             std::pow((Natom/V), 1.0/3.0)*std::pow(6.0/M_PI, 2.0/3.0));    f     = calc_fluidicity(Delta, toler);    y     = calc_y(f, Delta, toler);    z     = calc_compress(y);    Sig   = BOLTZ*(5.0/2.0+std::log(2*M_PI*BOLTZ*Temp/(gmx::square(PLANCK))*V/(f*Natom)));    Shs   = Sig+calc_Shs(f, y);    rho   = (tmass*AMU)/(V*NANO*NANO*NANO);    sigHS = std::cbrt(6*y*V/(M_PI*Natom));    fprintf(fplog, "System = /"%s/"/n", *top.name);    fprintf(fplog, "Nmol = %d/n", Nmol);    fprintf(fplog, "Natom = %d/n", Natom);    fprintf(fplog, "dt = %g ps/n", dt);    fprintf(fplog, "tmass = %g amu/n", tmass);    fprintf(fplog, "V = %g nm^3/n", V);    fprintf(fplog, "rho = %g g/l/n", rho);    fprintf(fplog, "T = %g K/n", Temp);    fprintf(fplog, "beta = %g mol/kJ/n", beta);    fprintf(fplog, "/nDoS parameters/n");    fprintf(fplog, "Delta = %g/n", Delta);    fprintf(fplog, "fluidicity = %g/n", f);    fprintf(fplog, "hard sphere packing fraction = %g/n", y);    fprintf(fplog, "hard sphere compressibility = %g/n", z);    fprintf(fplog, "ideal gas entropy = %g/n", Sig);    fprintf(fplog, "hard sphere entropy = %g/n", Shs);    fprintf(fplog, "sigma_HS = %g nm/n", sigHS);    fprintf(fplog, "DoS0 = %g/n", DoS0);    fprintf(fplog, "Dos2 = %g/n", dos2);    fprintf(fplog, "DoSTot = %g/n", dostot);    /* Now compute solid (2) and diffusive (3) components */    fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",                  bRecip ? "E (cm//S-1//N)" : "//f{12}n//f{4} (1/ps)",                  "//f{4}S(//f{12}n//f{4})", oenv);    xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);    recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;    for (j = 0; (j < nframes/4); j++)    {        dos[DOS_DIFF][j]  = DoS0/(1+gmx::square(DoS0*M_PI*nu[j]/(6*f*Natom)));        dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];        fprintf(fp, "%10g  %10g  %10g  %10g/n",                recip_fac*nu[j],                dos[DOS][j]/recip_fac,                dos[DOS_SOLID][j]/recip_fac,                dos[DOS_DIFF][j]/recip_fac);    }    xvgrclose(fp);    /* Finally analyze the results! */    wCdiff = 0.5;    wSdiff = Shs/(3*BOLTZ); /* Is this correct? */    wEdiff = 0.5;    wAdiff = wEdiff-wSdiff;    for (j = 0; (j < nframes/4); j++)    {        dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +                          dos[DOS_SOLID][j]*wCsolid(nu[j], beta));        dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +                          dos[DOS_SOLID][j]*wSsolid(nu[j], beta));        dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +                          dos[DOS_SOLID][j]*wAsolid(nu[j], beta));        dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +                          dos[DOS_SOLID][j]*wEsolid(nu[j], beta));    }    DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);    DiffCoeff = 1000*DiffCoeff/3.0;    fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s/n",            DiffCoeff);    fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s/n",            1000*DoS0/(12*tmass*beta));    cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,                                   nframes/4, &stddev);    fprintf(fplog, "Heat capacity %g J/mol K/n", 1000*cP/Nmol);    fprintf(fplog, "/nArrivederci!/n");    gmx_fio_fclose(fplog);    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");    return 0;}

//.........这里部分代码省略.........    }    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, NULL, box, TRUE);    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);    if (nmol > gnx || gnx % nmol != 0)    {        gmx_fatal(FARGS, "The number of atoms in the group (%d) is not a multiple of nmol (%d)", gnx, nmol);    }    nam = gnx/nmol;    natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);    snew(x_s, natoms);    j  = 0;    t0 = t;    if (bQ)    {        out = xvgropen(ftp2fn(efXVG, NFILE, fnm),                       "Radius of Charge (total and around axes)", "Time (ps)", "Rg (nm)", oenv);    }    else if (bMOI)    {        out = xvgropen(ftp2fn(efXVG, NFILE, fnm),                       "Moments of inertia (total and around axes)", "Time (ps)", "I (a.m.u. nm//S2//N)", oenv);    }    else    {        out = xvgropen(ftp2fn(efXVG, NFILE, fnm),                       "Radius of gyration (total and around axes)", "Time (ps)", "Rg (nm)", oenv);    }    if (bMOI)    {        xvgr_legend(out, NLEG, legI, oenv);    }    else    {        if (bRot)        {            if (output_env_get_print_xvgr_codes(oenv))            {                fprintf(out, "@ subtitle /"Axes are principal component axes/"/n");            }        }        xvgr_legend(out, NLEG, leg, oenv);    }    if (nz == 0)    {        gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);    }    do    {        if (nz == 0)        {            gmx_rmpbc_copy(gpbc, natoms, box, x, x_s);        }        gyro = 0;        clear_rvec(gvec);        clear_rvec(gvec1);        clear_rvec(d);        clear_rvec(d1);        for (mol = 0; mol < nmol; mol++)        {            tm    = sub_xcm(nz == 0 ? x_s : x, nam, index+mol*nam, top.atoms.atom, xcm, bQ);            if (nz == 0)            {

//.........这里部分代码省略.........        }        else        {            for (i = 0; i < ndon; i++)            {                if (donindex[i] >= natoms)                {                    indexOK = FALSE;                }            }            for (i = 0; i < nacc; i++)            {                if (accindex[i] >= natoms)                {                    indexOK = FALSE;                }            }        }        if (indexOK)        {            if (bDatout)            {                datfp = gmx_ffopen(out_datfile, "w");            }            if (bRKout)            {                rkfp = xvgropen(out_xvgrkfile,                                "Distance and //f{Symbol}k//f{}//S2//N trajectory",                                "Time (ps)", "Distance (nm) / //f{Symbol}k//f{}//S2//N",                                oenv);                xvgr_legend(rkfp, 2, rkleg, oenv);            }            if (bInstEffout)            {                iefp = xvgropen(out_xvginstefffile,                                "Instantaneous RET Efficiency",                                "Time (ps)", "RET Efficiency",                                oenv);                xvgr_legend(iefp, 1, ieleg, oenv);            }            if (bRhistout)            {                snew(rvalues, allocblock);                rblocksallocated += 1;                snew(rhist, histbins);            }            if (bKhistout)            {                snew(kappa2values, allocblock);                kblocksallocated += 1;                snew(khist, histbins);            }            do            {                clear_rvec(donvec);                clear_rvec(accvec);                clear_rvec(donpos);                clear_rvec(accpos);

//.........这里部分代码省略.........    DoS0 = dos[DOS][0];    /* Note this eqn. is incorrect in Pascal2011a! */    Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*             pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));    f     = calc_fluidicity(Delta, toler);    y     = calc_y(f, Delta, toler);    z     = calc_compress(y);    Sig   = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));    Shs   = Sig+calc_Shs(f, y);    rho   = (tmass*AMU)/(V*NANO*NANO*NANO);    sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);    fprintf(fplog, "System = /"%s/"/n", title);    fprintf(fplog, "Nmol = %d/n", Nmol);    fprintf(fplog, "Natom = %d/n", Natom);    fprintf(fplog, "dt = %g ps/n", dt);    fprintf(fplog, "tmass = %g amu/n", tmass);    fprintf(fplog, "V = %g nm^3/n", V);    fprintf(fplog, "rho = %g g/l/n", rho);    fprintf(fplog, "T = %g K/n", Temp);    fprintf(fplog, "beta = %g mol/kJ/n", beta);    fprintf(fplog, "/nDoS parameters/n");    fprintf(fplog, "Delta = %g/n", Delta);    fprintf(fplog, "fluidicity = %g/n", f);    fprintf(fplog, "hard sphere packing fraction = %g/n", y);    fprintf(fplog, "hard sphere compressibility = %g/n", z);    fprintf(fplog, "ideal gas entropy = %g/n", Sig);    fprintf(fplog, "hard sphere entropy = %g/n", Shs);    fprintf(fplog, "sigma_HS = %g nm/n", sigHS);    fprintf(fplog, "DoS0 = %g/n", DoS0);    fprintf(fplog, "Dos2 = %g/n", dos2);    fprintf(fplog, "DoSTot = %g/n", dostot);    /* Now compute solid (2) and diffusive (3) components */    fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",                  bRecip ? "E (cm//S-1//N)" : "//f{12}n//f{4} (1/ps)",                  "//f{4}S(//f{12}n//f{4})", oenv);    xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);    recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;    for (j = 0; (j < nframes/4); j++)    {        dos[DOS_DIFF][j]  = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));        dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];        fprintf(fp, "%10g  %10g  %10g  %10g/n",                recip_fac*nu[j],                dos[DOS][j]/recip_fac,                dos[DOS_SOLID][j]/recip_fac,                dos[DOS_DIFF][j]/recip_fac);    }    xvgrclose(fp);    /* Finally analyze the results! */    wCdiff = 0.5;    wSdiff = Shs/(3*BOLTZ); /* Is this correct? */    wEdiff = 0.5;    wAdiff = wEdiff-wSdiff;    for (j = 0; (j < nframes/4); j++)    {        dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +                          dos[DOS_SOLID][j]*wCsolid(nu[j], beta));        dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +                          dos[DOS_SOLID][j]*wSsolid(nu[j], beta));        dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +                          dos[DOS_SOLID][j]*wAsolid(nu[j], beta));        dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +                          dos[DOS_SOLID][j]*wEsolid(nu[j], beta));    }    DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);    DiffCoeff = 1000*DiffCoeff/3.0;    fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s/n",            DiffCoeff);    fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s/n",            1000*DoS0/(12*tmass*beta));    cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,                                   nframes/4, &stddev);    fprintf(fplog, "Heat capacity %g J/mol K/n", 1000*cP/Nmol);    /*       S  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,                                   nframes/4,&stddev);       fprintf(fplog,"Entropy %g J/mol K/n",1000*S/Nmol);       A  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,                                   nframes/4,&stddev);       fprintf(fplog,"Helmholtz energy %g kJ/mol/n",A/Nmol);       E  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,                                   nframes/4,&stddev);       fprintf(fplog,"Internal energy %g kJ/mol/n",E/Nmol);     */    fprintf(fplog, "/nArrivederci!/n");    gmx_fio_fclose(fplog);    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");    thanx(stderr);    return 0;}

static void pr_ff(t_coupl_rec *tcr,real time,t_idef *idef,		  t_commrec *cr,int nfile,t_filenm fnm[]){  static FILE *prop;  static FILE **out=NULL;  static FILE **qq=NULL;  static FILE **ip=NULL;  t_coupl_LJ  *tclj;  t_coupl_BU  *tcbu;  char        buf[256];  char        *leg[]  =  { "C12", "C6" };  char        *eleg[] =  { "Epsilon", "Sigma" };  char        *bleg[] = { "A", "B", "C" };  char        **raleg;  int         i,j,index;    if ((prop == NULL) && (out == NULL) && (qq == NULL) && (ip == NULL)) {    prop=xvgropen(opt2fn("-runav",nfile,fnm),		  "Properties and Running Averages","Time (ps)","");    snew(raleg,2*eoObsNR);    for(i=j=0; (i<eoObsNR); i++) {      if (tcr->bObsUsed[i]) {	raleg[j++] = strdup(eoNames[i]);	sprintf(buf,"RA-%s",eoNames[i]);	raleg[j++] = strdup(buf);      }    }    xvgr_legend(prop,j,raleg);    for(i=0; (i<j); i++)       sfree(raleg[i]);    sfree(raleg);          if (tcr->nLJ) {      snew(out,tcr->nLJ);      for(i=0; (i<tcr->nLJ); i++) {	if (tcr->tcLJ[i].bPrint) {	  tclj   = &(tcr->tcLJ[i]);	  out[i] = 	    xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),			       efXVG,tclj->at_i,tclj->at_j),		     "General Coupling Lennard Jones","Time (ps)",		     "Force constant (units)");	  fprintf(out[i],"@ subtitle /"Interaction between types %d and %d/"/n",		  tclj->at_i,tclj->at_j);	  if (tcr->combrule == 1)	    xvgr_legend(out[i],asize(leg),leg);	  else	    xvgr_legend(out[i],asize(eleg),eleg);	  fflush(out[i]);	}      }    }    else if (tcr->nBU) {      snew(out,tcr->nBU);      for(i=0; (i<tcr->nBU); i++) {	if (tcr->tcBU[i].bPrint) {	  tcbu=&(tcr->tcBU[i]);	  out[i] = 	    xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,			       tcbu->at_i,tcbu->at_j),		     "General Coupling Buckingham","Time (ps)",		     "Force constant (units)");	  fprintf(out[i],"@ subtitle /"Interaction between types %d and %d/"/n",		  tcbu->at_i,tcbu->at_j);	  xvgr_legend(out[i],asize(bleg),bleg);	  fflush(out[i]);	}      }    }    snew(qq,tcr->nQ);    for(i=0; (i<tcr->nQ); i++) {      if (tcr->tcQ[i].bPrint) {	qq[i] = xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,				   tcr->tcQ[i].at_i,-1),			 "General Coupling Charge","Time (ps)","Charge (e)");	fprintf(qq[i],"@ subtitle /"Type %d/"/n",tcr->tcQ[i].at_i);	fflush(qq[i]);      }    }    snew(ip,tcr->nIP);    for(i=0; (i<tcr->nIP); i++) {      sprintf(buf,"gctIP%d",tcr->tIP[i].type);      ip[i]=xvgropen(mk_gct_nm(opt2fn("-ffout",nfile,fnm),efXVG,0,-1),		     "General Coupling iparams","Time (ps)","ip ()");      index=tcr->tIP[i].type;      fprintf(ip[i],"@ subtitle /"Coupling to %s/"/n",	      interaction_function[idef->functype[index]].longname);      fflush(ip[i]);    }  }  /* Write properties to file */  fprintf(prop,"%10.3f",time);  for(i=0; (i<eoObsNR); i++)     if (tcr->bObsUsed[i])      fprintf(prop,"  %10.3e  %10.3e",tcr->act_value[i],tcr->av_value[i]);  fprintf(prop,"/n");  fflush(prop);    for(i=0; (i<tcr->nLJ); i++) {    tclj=&(tcr->tcLJ[i]);//.........这里部分代码省略.........

//.........这里部分代码省略.........                {                    sprintf(buf, "sb-%s:%s.xvg", cg[i].label, cg[j].label);                    fp = xvgropen(buf, buf, "Time (ps)", "Distance (nm)", oenv);                    for (k = 0; (k < teller); k++)                    {                        fprintf(fp, "%10g  %10g/n", time[k], cgdist[i][j][k]);                    }                    ffclose(fp);                }            }        }        sfree(buf);    }    else    {        for (m = 0; (m < 3); m++)        {            out[m] = xvgropen(fn[m], title[m], "Time (ps)", "Distance (nm)", oenv);        }        snew(buf, 256);        for (i = 0; (i < ncg); i++)        {            qi = cg[i].q;            for (j = i+1; (j < ncg); j++)            {                qj = cg[j].q;                if (nWithin[i][j])                {                    sprintf(buf, "%s:%s", cg[i].label, cg[j].label);                    if (qi*qj < 0)                    {                        nnn = 2;                    }                    else if (qi+qj > 0)                    {                        nnn = 0;                    }                    else                    {                        nnn = 1;                    }                    if (nset[nnn] == 0)                    {                        xvgr_legend(out[nnn], 1, (const char**)&buf, oenv);                    }                    else                    {                        if (output_env_get_xvg_format(oenv) == exvgXMGR)                        {                            fprintf(out[nnn], "@ legend string %d /"%s/"/n", nset[nnn], buf);                        }                        else if (output_env_get_xvg_format(oenv) == exvgXMGRACE)                        {                            fprintf(out[nnn], "@ s%d legend /"%s/"/n", nset[nnn], buf);                        }                    }                    nset[nnn]++;                    nWithin[i][j] = nnn+1;                }            }        }        for (k = 0; (k < teller); k++)        {            for (m = 0; (m < 3); m++)            {                fprintf(out[m], "%10g", time[k]);            }            for (i = 0; (i < ncg); i++)            {                for (j = i+1; (j < ncg); j++)                {                    nnn = nWithin[i][j];                    if (nnn > 0)                    {                        fprintf(out[nnn-1], "  %10g", cgdist[i][j][k]);                    }                }            }            for (m = 0; (m < 3); m++)            {                fprintf(out[m], "/n");            }        }        for (m = 0; (m < 3); m++)        {            ffclose(out[m]);            if (nset[m] == 0)            {                remove(fn[m]);            }        }    }    thanx(stderr);    return 0;}

//.........这里部分代码省略.........        { "-mw", FALSE, etBOOL, {&bMW},          "Use mass weighted fitting" }    };    FILE           *out;    t_trxstatus    *status;    t_topology      top;    int             ePBC;    rvec           *x_ref, *x;    matrix          box, R;    real            t;    int             natoms, i;    char           *grpname, title[256];    int             gnx;    gmx_rmpbc_t     gpbc = NULL;    atom_id        *index;    output_env_t    oenv;    real           *w_rls;    const char     *leg[]  = { "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" };#define NLEG asize(leg)    t_filenm        fnm[] = {        { efTRX, "-f",   NULL,       ffREAD },        { efTPS, NULL,   NULL,       ffREAD },        { efNDX, NULL,   NULL,       ffOPTRD },        { efXVG, NULL,   "rotmat",   ffWRITE }    };#define NFILE asize(fnm)    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))    {        return 0;    }    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x_ref, NULL, box, bMW);    gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);    gmx_rmpbc(gpbc, top.atoms.nr, box, x_ref);    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);    if (reffit[0][0] != 'n')    {        get_refx(oenv, ftp2fn(efTRX, NFILE, fnm), reffit[0][2] == 'z' ? 3 : 2, skip,                 gnx, index, bMW, &top, ePBC, x_ref);    }    natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);    snew(w_rls, natoms);    for (i = 0; i < gnx; i++)    {        if (index[i] >= natoms)        {            gmx_fatal(FARGS, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[i]+1, natoms);        }        w_rls[index[i]] = (bMW ? top.atoms.atom[index[i]].m : 1.0);    }    if (reffit[0][0] == 'n')    {        reset_x(gnx, index, natoms, NULL, x_ref, w_rls);    }    out = xvgropen(ftp2fn(efXVG, NFILE, fnm),                   "Fit matrix", "Time (ps)", "", oenv);    xvgr_legend(out, NLEG, leg, oenv);    do    {        gmx_rmpbc(gpbc, natoms, box, x);        reset_x(gnx, index, natoms, NULL, x, w_rls);        if (bFitXY)        {            do_fit_ndim(2, natoms, w_rls, x_ref, x);        }        calc_fit_R(DIM, natoms, w_rls, x_ref, x, R);        fprintf(out,                "%7g %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f/n",                t,                R[XX][XX], R[XX][YY], R[XX][ZZ],                R[YY][XX], R[YY][YY], R[YY][ZZ],                R[ZZ][XX], R[ZZ][YY], R[ZZ][ZZ]);    }    while (read_next_x(oenv, status, &t, x, box));    gmx_rmpbc_done(gpbc);    close_trj(status);    gmx_ffclose(out);    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");    return 0;}

//.........这里部分代码省略.........  if ((tmpf = fopen(pdbfile,"w")) == NULL) {    sprintf(pdbfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);    gmx_tmpnam(pdbfile);    if ((tmpf = fopen(pdbfile,"w")) == NULL)       gmx_fatal(FARGS,"Can not open tmp file %s",pdbfile);  }  else    fclose(tmpf);      strcpy(tmpfile,"ddXXXXXX");  gmx_tmpnam(tmpfile);  if ((tmpf = fopen(tmpfile,"w")) == NULL) {    sprintf(tmpfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);    gmx_tmpnam(tmpfile);    if ((tmpf = fopen(tmpfile,"w")) == NULL)       gmx_fatal(FARGS,"Can not open tmp file %s",tmpfile);  }  else    fclose(tmpf);    if ((dptr=getenv("DSSP")) == NULL)    dptr="/usr/local/bin/dssp";  if (!fexist(dptr))    gmx_fatal(FARGS,"DSSP executable (%s) does not exist (use setenv DSSP)",		dptr);  sprintf(dssp,"%s %s %s %s > /dev/null %s",	  dptr,bDoAccSurf?"":"-na",pdbfile,tmpfile,bVerbose?"":"2> /dev/null");  if (bVerbose)    fprintf(stderr,"dssp cmd='%s'/n",dssp);    if (fnTArea) {    fTArea=xvgropen(fnTArea,"Solvent Accessible Surface Area",		    xvgr_tlabel(),"Area (nm//S2//N)");    xvgr_legend(fTArea,2,leg);  } else    fTArea=NULL;    mat.map=NULL;  mat.nmap=getcmap(libopen(opt2fn("-map",NFILE,fnm)),		   opt2fn("-map",NFILE,fnm),&(mat.map));    natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);  if (natoms > atoms->nr)     gmx_fatal(FARGS,"/nTrajectory does not match topology!");  if (gnx > natoms)    gmx_fatal(FARGS,"/nTrajectory does not match selected group!");    snew(average_area,atoms->nres+10);  snew(av_area,atoms->nres+10);  snew(norm_av_area,atoms->nres+10);  accr=NULL;  naccr=0;  do {    t = convert_time(t);    if (nframe>=naccr) {      naccr+=10;      srenew(accr,naccr);      for(i=naccr-10; i<naccr; i++)	snew(accr[i],atoms->nres+10);    }    rm_pbc(&(top.idef),ePBC,natoms,box,x,x);    tapein=ffopen(pdbfile,"w");    write_pdbfile_indexed(tapein,NULL,atoms,x,ePBC,box,0,-1,gnx,index);    fclose(tapein);#ifdef GMX_NO_SYSTEM

//.........这里部分代码省略.........                    else if (emax <= emid)                    {                        write_xpm(out, 0, label, "Energy (kJ/mol)",                                  "Residue Index", "Residue Index",                                  ngroups, ngroups, groupnr, groupnr, emat[m],                                  emin, emid, rlo, rmid, &nlevels);                    }                    else                    {                        write_xpm3(out, 0, label, "Energy (kJ/mol)",                                   "Residue Index", "Residue Index",                                   ngroups, ngroups, groupnr, groupnr, emat[m],                                   emin, emid, emax, rlo, rmid, rhi, &nlevels);                    }                    gmx_ffclose(out);                }            }        }        snew(etot, egNR+egSP);        for (m = 0; (m < egNR+egSP); m++)        {            snew(etot[m], ngroups);            for (i = 0; (i < ngroups); i++)            {                for (j = 0; (j < ngroups); j++)                {                    etot[m][i] += emat[m][i][j];                }            }        }        out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",                       oenv);        xvgr_legend(out, 0, NULL, oenv);        j = 0;        if (output_env_get_print_xvgr_codes(oenv))        {            char str1[STRLEN], str2[STRLEN];            if (output_env_get_xvg_format(oenv) == exvgXMGR)            {                sprintf(str1, "@ legend string ");                sprintf(str2, " ");            }            else            {                sprintf(str1, "@ s");                sprintf(str2, " legend ");            }            for (m = 0; (m < egNR+egSP); m++)            {                if (egrp_use[m])                {                    fprintf(out, "%s%d%s /"%s/"/n", str1, j++, str2, egrp_nm[m]);                }            }            if (bFree)            {                fprintf(out, "%s%d%s /"%s/"/n", str1, j++, str2, "Free");            }            if (bFree)            {                fprintf(out, "%s%d%s /"%s/"/n", str1, j++, str2, "Diff");            }            fprintf(out, "@TYPE xy/n");            fprintf(out, "#%3s", "grp");

//.........这里部分代码省略.........    {        flags = flags | TRX_READ_X;        outx  = xvgropen(opt2fn("-ox", NFILE, fnm),                         bCom ? "Center of mass" : "Coordinate",                         output_env_get_xvgr_tlabel(oenv), "Coordinate (nm)", oenv);        make_legend(outx, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);    }    if (bOXT)    {        flags      = flags | TRX_READ_X;        status_out = open_trx(opt2fn("-oxt", NFILE, fnm), "w");    }    if (bOV)    {        flags = flags | TRX_READ_V;        outv  = xvgropen(opt2fn("-ov", NFILE, fnm),                         bCom ? "Center of mass velocity" : "Velocity",                         output_env_get_xvgr_tlabel(oenv), "Velocity (nm/ps)", oenv);        make_legend(outv, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);    }    if (bOF)    {        flags = flags | TRX_READ_F;        outf  = xvgropen(opt2fn("-of", NFILE, fnm), "Force",                         output_env_get_xvgr_tlabel(oenv), "Force (kJ mol//S-1//N nm//S-1//N)",                         oenv);        make_legend(outf, ngroups, isize0[0], index0[0], grpname, bCom, bMol, bDim, oenv);    }    if (bOB)    {        outb = xvgropen(opt2fn("-ob", NFILE, fnm), "Box vector elements",                        output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);        xvgr_legend(outb, 6, box_leg, oenv);    }    if (bOT)    {        bDum[XX]  = FALSE;        bDum[YY]  = FALSE;        bDum[ZZ]  = FALSE;        bDum[DIM] = TRUE;        flags     = flags | TRX_READ_V;        outt      = xvgropen(opt2fn("-ot", NFILE, fnm), "Temperature",                             output_env_get_xvgr_tlabel(oenv), "(K)", oenv);        make_legend(outt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);    }    if (bEKT)    {        bDum[XX]  = FALSE;        bDum[YY]  = FALSE;        bDum[ZZ]  = FALSE;        bDum[DIM] = TRUE;        flags     = flags | TRX_READ_V;        outekt    = xvgropen(opt2fn("-ekt", NFILE, fnm), "Center of mass translation",                             output_env_get_xvgr_tlabel(oenv), "Energy (kJ mol//S-1//N)", oenv);        make_legend(outekt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);    }    if (bEKR)    {        bDum[XX]  = FALSE;        bDum[YY]  = FALSE;        bDum[ZZ]  = FALSE;        bDum[DIM] = TRUE;        flags     = flags | TRX_READ_X | TRX_READ_V;        outekr    = xvgropen(opt2fn("-ekr", NFILE, fnm), "Center of mass rotation",                             output_env_get_xvgr_tlabel(oenv), "Energy (kJ mol//S-1//N)", oenv);

//.........这里部分代码省略.........      }    }    if (npcheck != nphobic)      gmx_fatal(FARGS,"Consistency check failed: not all %d atoms in the hydrophobic index/n"		  "found in the normal index selection (%d atoms)",nphobic,npcheck);  }  else    nphobic = 0;      for(i=0; (i<nx[0]); i++) {    ii = index[0][i];    if (!bFindex) {      bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;      if (bPhobic[i] && bOut[ii])	nphobic++;    }    if (bDGsol)      if (!gmx_atomprop_query(aps,epropDGsol,			      *(atoms->resinfo[atoms->atom[ii].resind].name),			      *(atoms->atomtype[ii]),&(dgs_factor[i])))	dgs_factor[i] = dgs_default;    if (debug)      fprintf(debug,"Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s/n",	      ii+1,*(atoms->resinfo[atoms->atom[ii].resind].name),	      *(atoms->atomname[ii]),	      atoms->atom[ii].q,radius[ii]-solsize,dgs_factor[i],	      BOOL(bPhobic[i]));  }  fprintf(stderr,"%d out of %d atoms were classified as hydrophobic/n",	  nphobic,nx[1]);    fp=xvgropen(opt2fn("-o",nfile,fnm),"Solvent Accessible Surface","Time (ps)",	      "Area (nm//S2//N)",oenv);  xvgr_legend(fp,asize(flegend) - (bDGsol ? 0 : 1),flegend,oenv);  vfile = opt2fn_null("-tv",nfile,fnm);  if (vfile) {    if (!bTop) {      gmx_fatal(FARGS,"Need a tpr file for option -tv");    }    vp=xvgropen(vfile,"Volume and Density","Time (ps)","",oenv);    xvgr_legend(vp,asize(vlegend),vlegend,oenv);    totmass  = 0;    ndefault = 0;    for(i=0; (i<nx[0]); i++) {      real mm;      ii = index[0][i];      /*      if (!query_atomprop(atomprop,epropMass,			  *(top->atoms.resname[top->atoms.atom[ii].resnr]),			  *(top->atoms.atomname[ii]),&mm))	ndefault++;      totmass += mm;      */      totmass += atoms->atom[ii].m;    }    if (ndefault)      fprintf(stderr,"WARNING: Using %d default masses for density calculation, which most likely are inaccurate/n",ndefault);  }  else    vp = NULL;      gmx_atomprop_destroy(aps);  if (bPBC)    gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);  

//.........这里部分代码省略.........    /* read index files */  ngrps = 2;  snew(com,ngrps);  snew(grpname,ngrps);  snew(index,ngrps);  snew(isize,ngrps);  get_index(&top->atoms,ftp2fn(efNDX,NFILE,fnm),ngrps,isize,index,grpname);    /* calculate mass */  max=0;  snew(mass,ngrps);  for(g=0;(g<ngrps);g++) {    mass[g]=0;    for(i=0;(i<isize[g]);i++) {      if (index[g][i]>max)	max=index[g][i];      if (index[g][i] >= top->atoms.nr)	gmx_fatal(FARGS,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)/n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);      mass[g]+=top->atoms.atom[index[g][i]].m;    }  }  natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);  t0 = t;  if (max>=natoms)    gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)/n",(int)max+1,natoms);  if (!bCutoff) {    /* open output file */    fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),		  "Distance","Time (ps)","Distance (nm)",oenv);    xvgr_legend(fp,4,leg,oenv);  } else {    ngrps = 1;    if (bLifeTime)      snew(contact_time,isize[1]);  }  if (ePBC != epbcNONE)    snew(pbc,1);  else    pbc = NULL;      gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);  do {    /* initialisation for correct distance calculations */    if (pbc) {      set_pbc(pbc,ePBC,box);      /* make molecules whole again */      gmx_rmpbc(gpbc,natoms,box,x);    }    /* calculate center of masses */    for(g=0;(g<ngrps);g++) {      if (isize[g] == 1) {	copy_rvec(x[index[g][0]],com[g]);      } else {	for(d=0;(d<DIM);d++) {	  com[g][d]=0;	  for(i=0;(i<isize[g]);i++) {	    com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;	  }	  com[g][d] /= mass[g];	}      }    }

//.........这里部分代码省略.........        leg[e++] = strdup(str);        sprintf(str, "f. <Ue//S-//betaU//N>");        leg[e++] = strdup(str);        for (i = 0; i < ngid; i++)        {            sprintf(str, "f. <U//sVdW %s//Ne//S-//betaU//N>",                    *(groups->grpname[groups->grps[egcENER].nm_ind[i]]));            leg[e++] = strdup(str);        }        if (bDispCorr)        {            sprintf(str, "f. <U//sdisp c//Ne//S-//betaU//N>");            leg[e++] = strdup(str);        }        if (bCharge)        {            for (i = 0; i < ngid; i++)            {                sprintf(str, "f. <U//sCoul %s//Ne//S-//betaU//N>",                        *(groups->grpname[groups->grps[egcENER].nm_ind[i]]));                leg[e++] = strdup(str);            }            if (bRFExcl)            {                sprintf(str, "f. <U//sRF excl//Ne//S-//betaU//N>");                leg[e++] = strdup(str);            }            if (EEL_FULL(fr->eeltype))            {                sprintf(str, "f. <U//sCoul recip//Ne//S-//betaU//N>");                leg[e++] = strdup(str);            }        }        xvgr_legend(fp_tpi, 4+nener, (const char**)leg, oenv);        for (i = 0; i < 4+nener; i++)        {            sfree(leg[i]);        }        sfree(leg);    }    clear_rvec(x_init);    V_all     = 0;    VembU_all = 0;    invbinw = 10;    nbin    = 10;    snew(bin, nbin);    /* Avoid frame step numbers <= -1 */    frame_step_prev = -1;    bNotLastFrame = read_first_frame(oenv, &status, opt2fn("-rerun", nfile, fnm),                                     &rerun_fr, TRX_NEED_X);    frame = 0;    if (rerun_fr.natoms - (bCavity ? nat_cavity : 0) !=        mdatoms->nr - (a_tp1 - a_tp0))    {        gmx_fatal(FARGS, "Number of atoms in trajectory (%d)%s "                  "is not equal the number in the run input file (%d) "                  "minus the number of atoms to insert (%d)/n",                  rerun_fr.natoms, bCavity ? " minus one" : "",                  mdatoms->nr, a_tp1-a_tp0);    }    refvolshift = log(det(rerun_fr.box));

//.........这里部分代码省略.........                break;            case edHdLPrintEnergyTOTAL:            case edHdLPrintEnergyYES:            default:                sprintf(buf, "%s (%s)", "Total Energy", unit_energy);        }        setname[s] = gmx_strdup(buf);        s         += 1;    }    if (fep->dhdl_derivatives == edhdlderivativesYES)    {        for (i = 0; i < efptNR; i++)        {            if (fep->separate_dvdl[i])            {                if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))                {                    /* compatibility output */                    sprintf(buf, "%s %s %.4f", dhdl, lambda, fep->init_lambda);                }                else                {                    double lam = fep->init_lambda;                    if (fep->init_lambda < 0)                    {                        lam = fep->all_lambda[i][fep->init_fep_state];                    }                    sprintf(buf, "%s %s = %.4f", dhdl, efpt_singular_names[i],                            lam);                }                setname[s] = gmx_strdup(buf);                s         += 1;            }        }    }    if (fep->n_lambda > 0)    {        /* g_bar has to determine the lambda values used in this simulation         * from this xvg legend.         */        if (expand->elmcmove > elmcmoveNO)        {            nsetsbegin = 1;  /* for including the expanded ensemble */        }        else        {            nsetsbegin = 0;        }        if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)        {            nsetsbegin += 1;        }        nsetsbegin += nsets_dhdl;        for (i = fep->lambda_start_n; i < fep->lambda_stop_n; i++)        {            print_lambda_vector(fep, i, FALSE, FALSE, lambda_vec_str);            if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))            {                /* for compatible dhdl.xvg files */                nps = sprintf(buf, "%s %s %s", deltag, lambda, lambda_vec_str);            }            else            {                nps = sprintf(buf, "%s %s to %s", deltag, lambda, lambda_vec_str);            }            if (ir->bSimTemp)            {                /* print the temperature for this state if doing simulated annealing */                sprintf(&buf[nps], "T = %g (%s)",                        ir->simtempvals->temperatures[s-(nsetsbegin)],                        unit_temp_K);            }            setname[s] = gmx_strdup(buf);            s++;        }        if (write_pV)        {            sprintf(buf, "pV (%s)", unit_energy);            setname[nsetsextend-1] = gmx_strdup(buf);  /* the first entry after                                                          nsets */        }        xvgr_legend(fp, nsetsextend, (const char **)setname, oenv);        for (s = 0; s < nsetsextend; s++)        {            sfree(setname[s]);        }        sfree(setname);    }    return fp;}

void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,                const output_env_t oenv){    FILE        *fp;    t_mapping   *map;    int          f, r, *count, *total, ss_count, total_count;    size_t       s;    const char** leg;    map = mat->map;    snew(count, mat->nmap);    snew(total, mat->nmap);    snew(leg, mat->nmap+1);    leg[0] = "Structure";    for (s = 0; s < (size_t)mat->nmap; s++)    {        leg[s+1] = gmx_strdup(map[s].desc);    }    fp = xvgropen(outfile, "Secondary Structure",                  output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);    if (output_env_get_print_xvgr_codes(oenv))    {        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, oenv);    total_count = 0;    for (s = 0; s < (size_t)mat->nmap; s++)    {        total[s] = 0;    }    for (f = 0; f < mat->nx; f++)    {        ss_count = 0;        for (s = 0; s < (size_t)mat->nmap; s++)        {            count[s] = 0;        }        for (r = 0; r < mat->ny; r++)        {            count[mat->matrix[f][r]]++;            total[mat->matrix[f][r]]++;        }        for (s = 0; s < (size_t)mat->nmap; s++)        {            if (strchr(ss_string, map[s].code.c1))            {                ss_count    += count[s];                total_count += count[s];            }        }        fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);        for (s = 0; s < (size_t)mat->nmap; s++)        {            fprintf(fp, " %5d", count[s]);        }        fprintf(fp, "/n");    }    /* now print column totals */    fprintf(fp, "%-8s %5d", "# Totals", total_count);    for (s = 0; s < (size_t)mat->nmap; s++)    {        fprintf(fp, " %5d", total[s]);    }    fprintf(fp, "/n");    /* now print percentages */    fprintf(fp, "%-8s %5.2f", "# SS %", total_count / (real) (mat->nx * mat->ny));    for (s = 0; s < (size_t)mat->nmap; s++)    {        fprintf(fp, " %5.2f", total[s] / (real) (mat->nx * mat->ny));    }    fprintf(fp, "/n");    xvgrclose(fp);    sfree(leg);    sfree(count);}

int gmx_sigeps(int argc, char *argv[]){    const char   *desc[] = {        "[TT]g_sigeps[tt] is a simple utility that converts C6/C12 or C6/Cn combinations",        "to [GRK]sigma[grk] and [GRK]epsilon[grk], or vice versa. It can also plot the potential",        "in  file. In addition, it makes an approximation of a Buckingham potential",        "to a Lennard-Jones potential."    };    static real   c6   = 1.0e-3, cn = 1.0e-6, qi = 0, qj = 0, sig = 0.3, eps = 1, sigfac = 0.7;    static real   Abh  = 1e5, Bbh = 32, Cbh = 1e-3;    static int    npow = 12;    t_pargs       pa[] = {        { "-c6",   FALSE,  etREAL,  {&c6},  "C6"   },        { "-cn",   FALSE,  etREAL,  {&cn},  "Constant for repulsion"   },        { "-pow",  FALSE,  etINT,   {&npow}, "Power of the repulsion term" },        { "-sig",  FALSE,  etREAL,  {&sig}, "[GRK]sigma[grk]"  },        { "-eps",  FALSE,  etREAL,  {&eps}, "[GRK]epsilon[grk]"  },        { "-A",    FALSE,  etREAL,  {&Abh}, "Buckingham A" },        { "-B",    FALSE,  etREAL,  {&Bbh}, "Buckingham B" },        { "-C",    FALSE,  etREAL,  {&Cbh}, "Buckingham C" },        { "-qi",   FALSE,  etREAL,  {&qi},  "qi"   },        { "-qj",   FALSE,  etREAL,  {&qj},  "qj"   },        { "-sigfac", FALSE, etREAL, {&sigfac}, "Factor in front of [GRK]sigma[grk] for starting the plot" }    };    t_filenm      fnm[] = {        { efXVG, "-o", "potje", ffWRITE }    };    output_env_t  oenv;#define NFILE asize(fnm)    const char   *legend[] = { "Lennard-Jones", "Buckingham" };    FILE         *fp;    int           i;    gmx_bool      bBham;    real          qq, x, oldx, minimum, mval, dp[2], pp[2];    int           cur = 0;#define next (1-cur)    parse_common_args(&argc, argv, PCA_CAN_VIEW,                      NFILE, fnm, asize(pa), pa, asize(desc),                      desc, 0, NULL, &oenv);    bBham = (opt2parg_bSet("-A", asize(pa), pa) ||             opt2parg_bSet("-B", asize(pa), pa) ||             opt2parg_bSet("-C", asize(pa), pa));    if (bBham)    {        c6  = Cbh;        sig = pow((6.0/npow)*pow(npow/Bbh, npow-6.0), 1.0/(npow-6.0));        eps = c6/(4*pow(sig, 6.0));        cn  = 4*eps*pow(sig, npow);    }    else    {        if (opt2parg_bSet("-sig", asize(pa), pa) ||            opt2parg_bSet("-eps", asize(pa), pa))        {            c6  = 4*eps*pow(sig, 6);            cn  = 4*eps*pow(sig, npow);        }        else if (opt2parg_bSet("-c6", asize(pa), pa) ||                 opt2parg_bSet("-cn", asize(pa), pa) ||                 opt2parg_bSet("-pow", asize(pa), pa))        {            sig = pow(cn/c6, 1.0/(npow-6.0));            eps = 0.25*c6*pow(sig, -6.0);        }        else        {            sig = eps = 0;        }        printf("c6    = %12.5e, c%d    = %12.5e/n", c6, npow, cn);        printf("sigma = %12.5f, epsilon = %12.5f/n", sig, eps);        minimum = pow(npow/6.0*pow(sig, npow-6.0), 1.0/(npow-6));        printf("Van der Waals minimum at %g, V = %g/n/n",               minimum, pot(minimum, 0, c6, cn, npow));        printf("Fit of Lennard Jones (%d-6) to Buckingham:/n", npow);        Bbh = npow/minimum;        Cbh = c6;        Abh = 4*eps*pow(sig/minimum, npow)*exp(npow);        printf("A = %g, B = %g, C = %g/n", Abh, Bbh, Cbh);    }    qq = qi*qj;    fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Potential", "r (nm)", "E (kJ/mol)",                  oenv);    xvgr_legend(fp, asize(legend), legend,                oenv);    if (sig == 0)    {        sig = 0.25;    }    minimum = -1;    mval    = 0;    oldx    = 0;    for (i = 0; (i < 100); i++)    {        x        = sigfac*sig+sig*i*0.02;        dp[next] = dpot(x, qq, c6, cn, npow);//.........这里部分代码省略.........

//.........这里部分代码省略.........    }    if ((dptr = getenv("DSSP")) == NULL)    {        dptr = "/usr/local/bin/dssp";    }    if (!gmx_fexist(dptr))    {        gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)",                  dptr);    }    if (dsspVersion >= 2)    {        if (dsspVersion > 2)        {            printf("/nWARNING: You use DSSP version %d, which is not explicitly/nsupported by do_dssp. Assuming version 2 syntax./n/n", dsspVersion);        }        sprintf(dssp, "%s -i %s -o %s > /dev/null %s",                dptr, pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");    }    else    {        sprintf(dssp, "%s %s %s %s > /dev/null %s",                dptr, bDoAccSurf ? "" : "-na", pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");    }    fprintf(stderr, "dssp cmd='%s'/n", dssp);    if (fnTArea)    {        fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",                          output_env_get_xvgr_tlabel(oenv), "Area (nm//S2//N)", oenv);        xvgr_legend(fTArea, 2, leg, oenv);    }    else    {        fTArea = NULL;    }    mat.map  = NULL;    mat.nmap = readcmap(opt2fn("-map", NFILE, fnm), &(mat.map));    natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);    if (natoms > atoms->nr)    {        gmx_fatal(FARGS, "/nTrajectory does not match topology!");    }    if (gnx > natoms)    {        gmx_fatal(FARGS, "/nTrajectory does not match selected group!");    }    snew(average_area, atoms->nres);    snew(av_area, atoms->nres);    snew(norm_av_area, atoms->nres);    accr  = NULL;    naccr = 0;    gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);    do    {        t = output_env_conv_time(oenv, t);        if (bDoAccSurf && nframe >= naccr)        {            naccr += 10;

//.........这里部分代码省略.........    get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm),              1, &isize, &index, &grpname);    snew(molind, top->mols.nr+1);    nmol = 0;    mol  = -1;    for (i = 0; i < isize; i++)    {        if (i == 0 || index[i] >= top->mols.index[mol+1])        {            molind[nmol++] = i;            do            {                mol++;            }            while (index[i] >= top->mols.index[mol+1]);        }    }    molind[nmol] = i;    nat_min      = top->atoms.nr;    nat_max      = 0;    for (mol = 0; mol < nmol; mol++)    {        nat_min = std::min(nat_min, molind[mol+1]-molind[mol]);        nat_max = std::max(nat_max, molind[mol+1]-molind[mol]);    }    fprintf(stderr, "Group %s consists of %d molecules/n", grpname, nmol);    fprintf(stderr, "Group size per molecule, min: %d atoms, max %d atoms/n",            nat_min, nat_max);    sprintf(title, "Size of %d polymers", nmol);    out = xvgropen(opt2fn("-o", NFILE, fnm), title, output_env_get_xvgr_tlabel(oenv), "(nm)",                   oenv);    xvgr_legend(out, bPC ? 8 : 5, leg, oenv);    if (opt2bSet("-v", NFILE, fnm))    {        outv = xvgropen(opt2fn("-v", NFILE, fnm), "Principal components",                        output_env_get_xvgr_tlabel(oenv), "(nm)", oenv);        snew(legp, DIM*DIM);        for (d = 0; d < DIM; d++)        {            for (d2 = 0; d2 < DIM; d2++)            {                sprintf(buf, "eig%d %c", d+1, 'x'+d2);                legp[d*DIM+d2] = gmx_strdup(buf);            }        }        xvgr_legend(outv, DIM*DIM, (const char**)legp, oenv);    }    else    {        outv = nullptr;    }    if (opt2bSet("-p", NFILE, fnm))    {        outp = xvgropen(opt2fn("-p", NFILE, fnm), "Persistence length",                        output_env_get_xvgr_tlabel(oenv), "bonds", oenv);        snew(bond, nat_max-1);        snew(sum_inp, nat_min/2);        snew(ninp, nat_min/2);    }    else    {        outp = nullptr;

//.........这里部分代码省略.........	    pt[fbin]++;	}	if (matfile)	  mcount[mbin]++;	if (otfile)	  tcount[fbin]++;      }    }    if (orfile) {      for(fbin=0; fbin<nr; fbin++) {	ff = f - (fbin + 1)*fshift;	if (ff >= 0) {	  for(i=0; i<isize; i++) {	    d2 = distance2(sx[f][i],sx[ff][i]);	    bin = (int)(sqrt(d2)*invbin);	    if (bin >= nalloc) {	      nallocn = 10*(bin/10) + 11;	      for(m=0; m<nr; m++) {		srenew(pr[m],nallocn);		for(i=nalloc; i<nallocn; i++)		  pr[m][i] = 0;	      }	      nalloc = nallocn;	    }	    pr[fbin][bin]++;	  }	  rcount[fbin]++;	}      }    }  }  fprintf(stderr,"/n");    if (matfile) {    matmax = 0;    for(f=0; f<mat_nx; f++) {      normfac = 1.0/(mcount[f]*isize*rbin);      for(i=0; i<nbin; i++) {	mat[f][i] *= normfac;	if (mat[f][i] > matmax && (f!=0 || i!=0))	  matmax = mat[f][i];      }    }    fprintf(stdout,"Value at (0,0): %.3f, maximum of the rest %.3f/n",	    mat[0][0],matmax);    if (mmax > 0)      matmax = mmax;    snew(tickx,mat_nx);    for(f=0; f<mat_nx; f++) {      if (sbin == 0)	tickx[f] = f*dt;      else	tickx[f] = f*sbin;    }    snew(ticky,nbin+1);    for(i=0; i<=nbin; i++)      ticky[i] = i*rbin;    fp = ffopen(matfile,"w");    write_xpm(fp,MAT_SPATIAL_Y,"Van Hove function","G (1/nm)",	      sbin==0 ? "time (ps)" : "sqrt(time) (ps^1/2)","r (nm)",	      mat_nx,nbin,tickx,ticky,mat,0,matmax,rlo,rhi,&nlev);         ffclose(fp);  }    if (orfile) {    fp = xvgropen(orfile,"Van Hove function","r (nm)","G (nm//S-1//N)",oenv);    fprintf(fp,"@ subtitle /"for particles in group %s/"/n",grpname);    snew(legend,nr);    for(fbin=0; fbin<nr; fbin++) {      sprintf(buf,"%g ps",(fbin + 1)*fshift*dt);      legend[fbin] = strdup(buf);    }    xvgr_legend(fp,nr,(const char**)legend,oenv);    for(i=0; i<nalloc; i++) {      fprintf(fp,"%g",i*rbin);      for(fbin=0; fbin<nr; fbin++)	fprintf(fp," %g",		(real)pr[fbin][i]/(rcount[fbin]*isize*rbin*(i==0 ? 0.5 : 1)));      fprintf(fp,"/n");    }    ffclose(fp);  }    if (otfile) {    sprintf(buf,"Probability of moving less than %g nm",rint);    fp = xvgropen(otfile,buf,"t (ps)","",oenv);    fprintf(fp,"@ subtitle /"for particles in group %s/"/n",grpname);    for(f=0; f<=ftmax; f++)      fprintf(fp,"%g %g/n",f*dt,(real)pt[f]/(tcount[f]*isize));    ffclose(fp);  }  do_view(oenv, matfile,NULL);  do_view(oenv, orfile,NULL);  do_view(oenv, otfile,NULL);  thanx(stderr);    return 0;}

//.........这里部分代码省略.........        label = "Residue";    }    else    {        label = "Atom";    }    for (i = 0; i < isize; i++)    {        rmsf[i] = U[i][XX*DIM + XX] + U[i][YY*DIM + YY] + U[i][ZZ*DIM + ZZ];    }    if (dirfn)    {        fprintf(stdout, "/n");        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++)

//.........这里部分代码省略.........      ind_fit[kkk] = kkk;    }        snew(atoms,1);    *atoms = gmx_mtop_global_atoms(&mtop);        if (atoms->pdbinfo == NULL) {      snew(atoms->pdbinfo,atoms->nr);    }  }   top = gmx_mtop_generate_local_top(&mtop,&ir);  g = NULL;  pbc_null = NULL;  if (ir.ePBC != epbcNONE) {    if (ir.bPeriodicMols)      pbc_null = &pbc;    else      g = mk_graph(fplog,&top->idef,0,mtop.natoms,FALSE,FALSE);  }    if (ftp2bSet(efNDX,NFILE,fnm)) {    rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);    xvg=xvgropen(opt2fn("-dr",NFILE,fnm),"Inidividual Restraints","Time (ps)",		 "nm",oenv);    snew(vvindex,isize);    snew(leg,isize);    for(i=0; (i<isize); i++) {      index[i]++;      snew(leg[i],12);      sprintf(leg[i],"index %d",index[i]);    }    xvgr_legend(xvg,isize,(const char**)leg,oenv);  }  else     isize=0;  ir.dr_tau=0.0;  init_disres(fplog,&mtop,&ir,NULL,FALSE,&fcd,NULL);  natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);  snew(f,5*natoms);    init_dr_res(&dr,fcd.disres.nres);  if (opt2bSet("-c",NFILE,fnm)) {    clust = cluster_index(fplog,opt2fn("-c",NFILE,fnm));    snew(dr_clust,clust->clust->nr+1);    for(i=0; (i<=clust->clust->nr); i++)      init_dr_res(&dr_clust[i],fcd.disres.nres);  }  else {	    out =xvgropen(opt2fn("-ds",NFILE,fnm),		  "Sum of Violations","Time (ps)","nm",oenv);    aver=xvgropen(opt2fn("-da",NFILE,fnm),		  "Average Violation","Time (ps)","nm",oenv);    numv=xvgropen(opt2fn("-dn",NFILE,fnm),		  "# Violations","Time (ps)","#",oenv);    maxxv=xvgropen(opt2fn("-dm",NFILE,fnm),		   "Largest Violation","Time (ps)","nm",oenv);  }  mdatoms = init_mdatoms(fplog,&mtop,ir.efep!=efepNO);  atoms2md(&mtop,&ir,0,NULL,0,mtop.natoms,mdatoms);  update_mdatoms(mdatoms,ir.init_lambda);  fr      = mk_forcerec();