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