//.........这里部分代码省略.........
    }
    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;
            srenew(accr, naccr);
            for (i = naccr-10; i < naccr; i++)
            {
                snew(accr[i], 2*atoms->nres-1);
            }
        }
        gmx_rmpbc(gpbc, natoms, box, x);
        tapein = gmx_ffopen(pdbfile, "w");
        write_pdbfile_indexed(tapein, NULL, atoms, x, ePBC, box, ' ', -1, gnx, index, NULL, TRUE);
        gmx_ffclose(tapein);

#ifdef GMX_NO_SYSTEM
        printf("Warning-- No calls to system(3) supported on this platform.");
        printf("Warning-- Skipping execution of 'system(\"%s\")'.", dssp);
        exit(1);
#else
        if (0 != system(dssp))
        {
            gmx_fatal(FARGS, "Failed to execute command: %s\n",
                      "Try specifying your dssp version with the -ver option.", dssp);

static void process_tcaf(int nframes, real dt, int nkc, real **tc, rvec *kfac,
                         real rho, real wt, const char *fn_trans,
                         const char *fn_tca, const char *fn_tc,
                         const char *fn_tcf, const char *fn_cub,
                         const char *fn_vk, const output_env_t oenv)
{
    FILE  *fp, *fp_vk, *fp_cub = NULL;
    int    nk, ntc;
    real **tcaf, **tcafc = NULL, eta;
    int    i, j, k, kc;
    int    ncorr;
    real   fitparms[3], *sig;

    nk  = kset_c[nkc];
    ntc = nk*NPK;

    if (fn_trans)
    {
        fp = xvgropen(fn_trans, "Transverse Current", "Time (ps)", "TC (nm/ps)",
                      oenv);
        for (i = 0; i < nframes; i++)
        {
            fprintf(fp, "%g", i*dt);
            for (j = 0; j < ntc; j++)
            {
                fprintf(fp, " %g", tc[j][i]);
            }
            fprintf(fp, "\n");
        }
        ffclose(fp);
        do_view(oenv, fn_trans, "-nxy");
    }

    ncorr = (nframes+1)/2;
    if (ncorr > (int)(5*wt/dt+0.5))
    {
        ncorr = (int)(5*wt/dt+0.5)+1;
    }
    snew(tcaf, nk);
    for (k = 0; k < nk; k++)
    {
        snew(tcaf[k], ncorr);
    }
    if (fn_cub)
    {
        snew(tcafc, nkc);
        for (k = 0; k < nkc; k++)
        {
            snew(tcafc[k], ncorr);
        }
    }
    snew(sig, ncorr);
    for (i = 0; i < ncorr; i++)
    {
        sig[i] = exp(0.5*i*dt/wt);
    }

    low_do_autocorr(fn_tca, oenv, "Transverse Current Autocorrelation Functions",
                    nframes, ntc, ncorr, tc, dt, eacNormal,
                    1, FALSE, FALSE, FALSE, 0, 0, 0);
    do_view(oenv, fn_tca, "-nxy");

    fp = xvgropen(fn_tc, "Transverse Current Autocorrelation Functions",
                  "Time (ps)", "TCAF", oenv);
    for (i = 0; i < ncorr; i++)
    {
        kc = 0;
        fprintf(fp, "%g", i*dt);
        for (k = 0; k < nk; k++)
        {
            for (j = 0; j < NPK; j++)
            {
                tcaf[k][i] += tc[NPK*k+j][i];
            }
            if (fn_cub)
            {
                for (j = 0; j < NPK; j++)
                {
                    tcafc[kc][i] += tc[NPK*k+j][i];
                }
            }
            if (i == 0)
            {
                fprintf(fp, " %g", 1.0);
            }
            else
            {
                tcaf[k][i] /= tcaf[k][0];
                fprintf(fp, " %g", tcaf[k][i]);
            }
            if (k+1 == kset_c[kc+1])
            {
                kc++;
            }
        }
        fprintf(fp, "\n");
    }
    ffclose(fp);
    do_view(oenv, fn_tc, "-nxy");

//.........这里部分代码省略.........

/* Allocate memory for the inputinfo struct: */
static void create_info(t_inputinfo *info)
{
    snew(info->fac, info->n_entries);
    snew(info->rcoulomb, info->n_entries);
    snew(info->rvdw, info->n_entries);
    snew(info->nkx, info->n_entries);
    snew(info->nky, info->n_entries);
    snew(info->nkz, info->n_entries);
    snew(info->fourier_sp, info->n_entries);
    snew(info->ewald_rtol, info->n_entries);
    snew(info->ewald_beta, info->n_entries);
    snew(info->pme_order, info->n_entries);
    snew(info->fn_out, info->n_entries);
    snew(info->e_dir, info->n_entries);
    snew(info->e_rec, info->n_entries);
}

int gmx_potential(int argc, char *argv[])
{
    const char        *desc[] = {
        "[THISMODULE] computes the electrostatical potential across the box. The potential is",
        "calculated by first summing the charges per slice and then integrating",
        "twice of this charge distribution. Periodic boundaries are not taken",
        "into account. Reference of potential is taken to be the left side of",
        "the box. It is also possible to calculate the potential in spherical",
        "coordinates as function of r by calculating a charge distribution in",
        "spherical slices and twice integrating them. epsilon_r is taken as 1,",
        "but 2 is more appropriate in many cases."
    };
    output_env_t       oenv;
    static int         axis       = 2;       /* normal to memb. default z  */
    static const char *axtitle    = "Z";
    static int         nslices    = 10;      /* nr of slices defined       */
    static int         ngrps      = 1;
    static gmx_bool    bSpherical = FALSE;   /* default is bilayer types   */
    static real        fudge_z    = 0;       /* translate coordinates      */
    static gmx_bool    bCorrect   = 0;
    t_pargs            pa []      = {
        {   "-d",   FALSE, etSTR, {&axtitle},
            "Take the normal on the membrane in direction X, Y or Z."
        },
        {   "-sl",  FALSE, etINT, {&nslices},
            "Calculate potential as function of boxlength, dividing the box"
            " in this number of slices."
        },
        {   "-cb",  FALSE, etINT, {&cb},
            "Discard this number of  first slices of box for integration"
        },
        {   "-ce",  FALSE, etINT, {&ce},
            "Discard this number of last slices of box for integration"
        },
        {   "-tz",  FALSE, etREAL, {&fudge_z},
            "Translate all coordinates by this distance in the direction of the box"
        },
        {   "-spherical", FALSE, etBOOL, {&bSpherical},
            "Calculate spherical thingie"
        },
        {   "-ng",       FALSE, etINT, {&ngrps},
            "Number of groups to consider"
        },
        {   "-correct",  FALSE, etBOOL, {&bCorrect},
            "Assume net zero charge of groups to improve accuracy"
        }
    };
    const char        *bugs[] = {
        "Discarding slices for integration should not be necessary."
    };

    double           **potential,              /* potential per slice        */
                     **charge,                                  /* total charge per slice     */
                     **field,                                   /* field per slice            */
                     slWidth;                /* width of one slice         */
    char      **grpname;                       /* groupnames                 */
    int        *ngx;                           /* sizes of groups            */
    t_topology *top;                           /* topology        */
    int         ePBC;
    atom_id   **index;                         /* indices for all groups     */
    t_filenm    fnm[] = {                      /* files for g_order       */
        { efTRX, "-f", NULL,  ffREAD },        /* trajectory file             */
        { efNDX, NULL, NULL,  ffREAD },        /* index file          */
        { efTPX, NULL, NULL,  ffREAD },        /* topology file               */
        { efXVG, "-o", "potential", ffWRITE }, /* xvgr output file    */
        { efXVG, "-oc", "charge", ffWRITE },   /* xvgr output file    */
        { efXVG, "-of", "field", ffWRITE },    /* xvgr output file    */
    };

#define NFILE asize(fnm)

    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
                           &oenv))
    {
        return 0;
    }

    /* Calculate axis */
    axis = toupper(axtitle[0]) - 'X';

    top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */

    snew(grpname, ngrps);
    snew(index, ngrps);
    snew(ngx, ngrps);

    rd_index(ftp2fn(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);


    calc_potential(ftp2fn(efTRX, NFILE, fnm), index, ngx,
                   &potential, &charge, &field,
                   &nslices, top, ePBC, axis, ngrps, &slWidth, fudge_z,
                   bSpherical, bCorrect, oenv);

    plot_potential(potential, charge, field, opt2fn("-o", NFILE, fnm),
                   opt2fn("-oc", NFILE, fnm), opt2fn("-of", NFILE, fnm),
                   nslices, ngrps, (const char**)grpname, slWidth, oenv);

    do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);  /* view xvgr file */
//.........这里部分代码省略.........

static void check_solvent(FILE *fp,t_topology *top,t_forcerec *fr,
			  t_mdatoms *md,t_nsborder *nsb)
{
  /* This routine finds out whether a charge group can be used as
   * solvent in the innerloops. The routine should be called once
   * at the beginning of the MD program.
   */
  t_block *cgs,*excl,*mols;
  atom_id *cgid;
  int     i,j,m,j0,j1,nj,k,aj,ak,tjA,tjB,nl_m,nl_n,nl_o;
  int     warncount;
  bool    bOneCG;
  bool    *bAllExcl,bAE,bOrder;
  bool    *bHaveLJ,*bHaveCoul;
  
  cgs  = &(top->blocks[ebCGS]);
  excl = &(top->atoms.excl);
  mols = &(top->blocks[ebMOLS]);

  if (fp)
    fprintf(fp,"Going to determine what solvent types we have.\n");
  snew(fr->solvent_type,cgs->nr+1);
  snew(fr->mno_index,(cgs->nr+1)*3);
  
  /* Generate charge group number for all atoms */
  cgid = make_invblock(cgs,cgs->nra);
  
  warncount=0;

  /* Loop over molecules */
  if (fp)
    fprintf(fp,"There are %d molecules, %d charge groups and %d atoms\n",
	    mols->nr,cgs->nr,cgs->nra);
  for(i=0; (i<mols->nr); i++) {
    /* Set boolean that determines whether the molecules consists of one CG */
    bOneCG = TRUE;
    /* Set some counters */
    j0     = mols->index[i];
    j1     = mols->index[i+1];
    nj     = j1-j0;
    for(j=j0+1; (j<j1); j++) {
      bOneCG = bOneCG && (cgid[mols->a[j]] == cgid[mols->a[j-1]]);
    }
    if (fr->bSolvOpt && bOneCG && nj>1) {
      /* Check whether everything is excluded */
      snew(bAllExcl,nj);
      bAE = TRUE;
      /* Loop over all atoms in molecule */
      for(j=j0; (j<j1) && bAE; j++) {
	/* Set a flag for each atom in the molecule that determines whether
	 * it is excluded or not 
	 */
	for(k=0; (k<nj); k++)
	  bAllExcl[k] = FALSE;
	/* Now check all the exclusions of this atom */
	for(k=excl->index[j]; (k<excl->index[j+1]); k++) {
	  ak = excl->a[k];
	  /* Consistency and range check */
	  if ((ak < j0) || (ak >= j1)) 
	    fatal_error(0,"Exclusion outside molecule? ak = %d, j0 = %d, j1 = 5d, mol is %d",ak,j0,j1,i);
	  bAllExcl[ak-j0] = TRUE;
	}
	/* Now sum up the booleans */
	for(k=0; (k<nj); k++)
	  bAE = bAE && bAllExcl[k];
      }
      if (bAE) {
	snew(bHaveCoul,nj);
	snew(bHaveLJ,nj);
	for(j=j0; (j<j1); j++) {
	  /* Check for coulomb */
	  aj = mols->a[j];
	  bHaveCoul[j-j0] = ( (fabs(top->atoms.atom[aj].q ) > GMX_REAL_MIN) ||
			      (fabs(top->atoms.atom[aj].qB) > GMX_REAL_MIN));
	  /* Check for LJ. */
	  tjA = top->atoms.atom[aj].type;
	  tjB = top->atoms.atom[aj].typeB;
	  bHaveLJ[j-j0] = FALSE;
	  for(k=0; (k<fr->ntype); k++) {
	    if (fr->bBHAM) 
	      bHaveLJ[j-j0] = (bHaveLJ[j-j0] || 
			       (fabs(BHAMA(fr->nbfp,fr->ntype,tjA,k)) > GMX_REAL_MIN) ||
			       (fabs(BHAMB(fr->nbfp,fr->ntype,tjA,k)) > GMX_REAL_MIN) ||
			       (fabs(BHAMC(fr->nbfp,fr->ntype,tjA,k)) > GMX_REAL_MIN) ||
			       (fabs(BHAMA(fr->nbfp,fr->ntype,tjB,k)) > GMX_REAL_MIN) ||
			       (fabs(BHAMB(fr->nbfp,fr->ntype,tjB,k)) > GMX_REAL_MIN) ||
			       (fabs(BHAMC(fr->nbfp,fr->ntype,tjB,k)) > GMX_REAL_MIN));
	    else
	      bHaveLJ[j-j0] = (bHaveLJ[j-j0] || 
			       (fabs(C6(fr->nbfp,fr->ntype,tjA,k))  > GMX_REAL_MIN) ||
			       (fabs(C12(fr->nbfp,fr->ntype,tjA,k)) > GMX_REAL_MIN) ||
			       (fabs(C6(fr->nbfp,fr->ntype,tjB,k))  > GMX_REAL_MIN) ||
			       (fabs(C12(fr->nbfp,fr->ntype,tjB,k)) > GMX_REAL_MIN));
	  }
	}
	/* Now we have determined what particles have which interactions 
	 * In the case of water-like molecules we only check for the number
	 * of particles and the LJ, not for the Coulomb. Let's just assume
	 * that the water loops are faster than the MNO loops anyway. DvdS
	 */
//.........这里部分代码省略.........

int gmx_pme_error(int argc, char *argv[])
{
    const char     *desc[] = {
        "[THISMODULE] estimates the error of the electrostatic forces",
        "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
        "the splitting parameter such that the error is equally",
        "distributed over the real and reciprocal space part.",
        "The part of the error that stems from self interaction of the particles "
        "is computationally demanding. However, a good a approximation is to",
        "just use a fraction of the particles for this term which can be",
        "indicated by the flag [TT]-self[tt].[PAR]",
    };

    real            fs        = 0.0; /* 0 indicates: not set by the user */
    real            user_beta = -1.0;
    real            fracself  = 1.0;
    t_inputinfo     info;
    t_state         state;     /* The state from the tpr input file */
    gmx_mtop_t      mtop;      /* The topology from the tpr input file */
    t_inputrec     *ir = NULL; /* The inputrec from the tpr file */
    FILE           *fp = NULL;
    t_commrec      *cr;
    unsigned long   PCA_Flags;
    gmx_bool        bTUNE    = FALSE;
    gmx_bool        bVerbose = FALSE;
    int             seed     = 0;


    static t_filenm fnm[] = {
        { efTPR, "-s",     NULL,    ffREAD },
        { efOUT, "-o",    "error",  ffWRITE },
        { efTPR, "-so",   "tuned",  ffOPTWR }
    };

    output_env_t    oenv = NULL;

    t_pargs         pa[] = {
        { "-beta",     FALSE, etREAL, {&user_beta},
          "If positive, overwrite ewald_beta from [TT].tpr[tt] file with this value" },
        { "-tune",     FALSE, etBOOL, {&bTUNE},
          "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
        { "-self",     FALSE, etREAL, {&fracself},
          "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
        { "-seed",     FALSE, etINT,  {&seed},
          "Random number seed used for Monte Carlo algorithm when [TT]-self[tt] is set to a value between 0.0 and 1.0" },
        { "-v",        FALSE, etBOOL, {&bVerbose},
          "Be loud and noisy" }
    };


#define NFILE asize(fnm)

    cr = init_commrec();

    PCA_Flags  = PCA_NOEXIT_ON_ARGS;

    if (!parse_common_args(&argc, argv, PCA_Flags,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc,
                           0, NULL, &oenv))
    {
        return 0;
    }

    if (!bTUNE)
    {
        bTUNE = opt2bSet("-so", NFILE, fnm);
    }

    info.n_entries = 1;

    /* Allocate memory for the inputinfo struct: */
    create_info(&info);
    info.fourier_sp[0] = fs;

    /* Read in the tpr file and open logfile for reading */
    if (MASTER(cr))
    {
        snew(ir, 1);
        read_tpr_file(opt2fn("-s", NFILE, fnm), &info, &state, &mtop, ir, user_beta, fracself);

        fp = fopen(opt2fn("-o", NFILE, fnm), "w");
    }

    /* Check consistency if the user provided fourierspacing */
    if (fs > 0 && MASTER(cr))
    {
        /* Recalculate the grid dimensions using fourierspacing from user input */
        info.nkx[0] = 0;
        info.nky[0] = 0;
        info.nkz[0] = 0;
        calc_grid(stdout, state.box, info.fourier_sp[0], &(info.nkx[0]), &(info.nky[0]), &(info.nkz[0]));
        if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
        {
            gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
                      fs, ir->nkx, ir->nky, ir->nkz, info.nkx[0], info.nky[0], info.nkz[0]);
        }
    }

    /* Estimate (S)PME force error */

//.........这里部分代码省略.........

void set_histp(t_atoms *pdba, rvec *x, real angle, real dist)
{
    static const char *prot_acc[] = {
        "O", "OD1", "OD2", "OE1", "OE2", "OG", "OG1", "OH", "OW"
    };
#define NPA asize(prot_acc)
    static const char *prot_don[] = {
        "N", "NH1", "NH2", "NE", "ND1", "ND2", "NE2", "NZ", "OG", "OG1", "OH", "NE1", "OW"
    };
#define NPD asize(prot_don)

    gmx_bool *donor, *acceptor;
    gmx_bool *hbond, bHaveH = FALSE;
    gmx_bool  bHDd, bHEd;
    rvec      xh1, xh2;
    int       natom;
    int       i, j, nd, na, aj, hisind, his0, type = -1;
    int       nd1, ne2, cg, cd2, ce1;
    t_blocka *hb;
    real      d;
    char     *atomnm;

    natom = pdba->nr;

    i = 0;
    while (i < natom &&
           gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name, "HIS") != 0)
    {
        i++;
    }
    if (natom == i)
    {
        return;
    }

    /* A histidine residue exists that requires automated assignment, so
     * doing the analysis of donors and acceptors is worthwhile. */
    fprintf(stderr,
            "Analysing hydrogen-bonding network for automated assignment of histidine\n"
            " protonation.");

    snew(donor, natom);
    snew(acceptor, natom);
    snew(hbond, natom);
    snew(hb, 1);

    nd = na = 0;
    for (j = 0; (j < natom); j++)
    {
        if (in_strings(*pdba->atomname[j], NPA, prot_acc) != -1)
        {
            acceptor[j] = TRUE;
            na++;
        }
        if (in_strings(*pdba->atomname[j], NPD, prot_don) != -1)
        {
            donor[j] = TRUE;
            nd++;
        }
    }
    fprintf(stderr, " %d donors and %d acceptors were found.\n", nd, na);
    chk_allhb(pdba, x, hb, donor, acceptor, dist);
    if (debug)
    {
        pr_hbonds(debug, hb, pdba);
    }
    fprintf(stderr, "There are %d hydrogen bonds\n", hb->nra);

    /* Now do the HIS stuff */
    hisind = -1;
    while (i < natom)
    {
        if (gmx_strcasecmp(*pdba->resinfo[pdba->atom[i].resind].name, "HIS") != 0)
        {
            i++;
        }
        else
        {
            if (pdba->atom[i].resind != hisind)
            {
                hisind = pdba->atom[i].resind;

                /* Find the  atoms in the ring */
                nd1 = ne2 = cg = cd2 = ce1 = -1;
                while (i < natom && pdba->atom[i].resind == hisind)
                {
                    atomnm = *pdba->atomname[i];
                    if (strcmp(atomnm, "CD2") == 0)
                    {
                        cd2 = i;
                    }
                    else if (strcmp(atomnm, "CG") == 0)
                    {
                        cg  = i;
                    }
                    else if (strcmp(atomnm, "CE1") == 0)
                    {
                        ce1 = i;
                    }
                    else if (strcmp(atomnm, "ND1") == 0)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    { efNDX, NULL,  NULL,   ffREAD  },
    { efTRX, "-f",  NULL,   ffREAD  },
    { efG87, "-to", NULL,   ffOPTWR },
    { efSTO, "-cz", "zconf",ffWRITE },
    { efSTO, "-co", "waver",ffWRITE }
  };
#define NFILE asize(fnm)

  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
  
  bRange=(opt2parg_bSet("-ahxstart",asize(pa),pa) &&
	  opt2parg_bSet("-ahxend",asize(pa),pa));
		        
  top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
  
  natoms=read_first_x(oenv,&status,opt2fn("-f",NFILE,fnm),&t,&x,box);

  if (opt2bSet("-to",NFILE,fnm)) {
    otrj=opt2FILE("-to",NFILE,fnm,"w");
    strcpy(prop,ppp[0]);
    fprintf(otrj,"%s Weighted Trajectory: %d atoms, NO box\n",prop,natoms);
  }
  else
    otrj=NULL;
    
  if (natoms != top->atoms.nr)
    gmx_fatal(FARGS,"Sorry can only run when the number of atoms in the run input file (%d) is equal to the number in the trajectory (%d)",
	    top->atoms.nr,natoms);
	    
  bb=mkbbind(ftp2fn(efNDX,NFILE,fnm),&nres,&nbb,r0,&nall,&allindex,
	     top->atoms.atomname,top->atoms.atom,top->atoms.resinfo);
  snew(bbindex,natoms);
  snew(caindex,nres);
  
  fprintf(stderr,"nall=%d\n",nall);
    
  /* Open output files, default x-axis is time */
  for(i=0; (i<efhNR); i++) {
    sprintf(buf,"%s.xvg",xf[i].filenm);
    remove(buf);
    xf[i].fp=xvgropen(buf,xf[i].title,
                      xf[i].xaxis ? xf[i].xaxis : "Time (ps)",
		      xf[i].yaxis,oenv);
    if (xf[i].bfp2) {
      sprintf(buf,"%s.out",xf[i].filenm);
      remove(buf);
      xf[i].fp2=ffopen(buf,"w");
    }
  }

  /* Read reference frame from tpx file to compute helix length */
  snew(xref,top->atoms.nr);
  read_tpx(ftp2fn(efTPX,NFILE,fnm),
	   NULL,NULL,&natoms,xref,NULL,NULL,NULL);
  calc_hxprops(nres,bb,xref,box);
  do_start_end(nres,bb,xref,&nbb,bbindex,&nca,caindex,bRange,rStart,rEnd);
  sfree(xref);
  if (bDBG) {
    fprintf(stderr,"nca=%d, nbb=%d\n",nca,nbb);
    pr_bb(stdout,nres,bb);
  }
  
  gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);

static int scan_ene_files(char **fnms, int nfiles,
                          real *readtime, real *timestep, int *nremax)
{
    /* Check number of energy terms and start time of all files */
    int          f, nre, nremin = 0, nresav = 0;
    ener_file_t  in;
    real         t1, t2;
    char         inputstring[STRLEN];
    gmx_enxnm_t *enm;
    t_enxframe  *fr;

    snew(fr, 1);

    for (f = 0; f < nfiles; f++)
    {
        in  = open_enx(fnms[f], "r");
        enm = NULL;
        do_enxnms(in, &nre, &enm);

        if (f == 0)
        {
            nresav  = nre;
            nremin  = nre;
            *nremax = nre;
            do_enx(in, fr);
            t1 = fr->t;
            do_enx(in, fr);
            t2          = fr->t;
            *timestep   = t2-t1;
            readtime[f] = t1;
            close_enx(in);
        }
        else
        {
            nremin  = std::min(nremin, fr->nre);
            *nremax = std::max(*nremax, fr->nre);
            if (nre != nresav)
            {
                fprintf(stderr,
                        "Energy files don't match, different number of energies:\n"
                        " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
                fprintf(stderr,
                        "\nContinue conversion using only the first %d terms (n/y)?\n"
                        "(you should be sure that the energy terms match)\n", nremin);
                if (NULL == fgets(inputstring, STRLEN-1, stdin))
                {
                    gmx_fatal(FARGS, "Error reading user input");
                }
                if (inputstring[0] != 'y' && inputstring[0] != 'Y')
                {
                    fprintf(stderr, "Will not convert\n");
                    exit(0);
                }
                nresav = fr->nre;
            }
            do_enx(in, fr);
            readtime[f] = fr->t;
            close_enx(in);
        }
        fprintf(stderr, "\n");
        free_enxnms(nre, enm);
    }

    free_enxframe(fr);
    sfree(fr);

    return nremin;
}

int gmx_eneconv(int argc, char *argv[])
{
    const char       *desc[] = {
        "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[PAR]",
        "Concatenates several energy files in sorted order.",
        "In the case of double time frames, the one",
        "in the later file is used. By specifying [TT]-settime[tt] you will be",
        "asked for the start time of each file. The input files are taken",
        "from the command line,",
        "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
        "the trick. [PAR]",
        "With [IT]one file[it] specified for [TT]-f[tt]:[PAR]",
        "Reads one energy file and writes another, applying the [TT]-dt[tt],",
        "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
        "converting to a different format if necessary (indicated by file",
        "extentions).[PAR]",
        "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
        "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
    };
    const char       *bugs[] = {
        "When combining trajectories the sigma and E^2 (necessary for statistics) are not updated correctly. Only the actual energy is correct. One thus has to compute statistics in another way."
    };
    ener_file_t       in  = NULL, out = NULL;
    gmx_enxnm_t      *enm = NULL;
#if 0
    ener_file_t       in, out = NULL;
    gmx_enxnm_t      *enm = NULL;
#endif
    t_enxframe       *fr, *fro;
    gmx_int64_t       ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
    t_energy         *ee_sum;
    gmx_int64_t       lastfilestep, laststep, startstep_file = 0;
    int               noutfr;
    int               nre, nremax, this_nre, nfile, f, i, kkk, nset, *set = NULL;
    double            last_t;
    char            **fnms;
    real             *readtime, *settime, timestep, tadjust;
    char              buf[22], buf2[22];
    int              *cont_type;
    gmx_bool          bNewFile, bFirst, bNewOutput;
    gmx_output_env_t *oenv;
    gmx_bool          warned_about_dh = FALSE;
    t_enxblock       *blocks          = NULL;
    int               nblocks         = 0;
    int               nblocks_alloc   = 0;

    t_filenm          fnm[] = {
        { efEDR, "-f", NULL,    ffRDMULT },
        { efEDR, "-o", "fixed", ffWRITE  },
    };

#define NFILE asize(fnm)
    gmx_bool         bWrite;
    static real      delta_t   = 0.0, toffset = 0, scalefac = 1;
    static gmx_bool  bSetTime  = FALSE;
    static gmx_bool  bSort     = TRUE, bError = TRUE;
    static real      begin     = -1;
    static real      end       = -1;
    gmx_bool         remove_dh = FALSE;

    t_pargs          pa[] = {
        { "-b",        FALSE, etREAL, {&begin},
          "First time to use"},
        { "-e",        FALSE, etREAL, {&end},
          "Last time to use"},
        { "-dt",       FALSE, etREAL, {&delta_t},
          "Only write out frame when t MOD dt = offset" },
        { "-offset",   FALSE, etREAL, {&toffset},
          "Time offset for [TT]-dt[tt] option" },
        { "-settime",  FALSE, etBOOL, {&bSetTime},
          "Change starting time interactively" },
        { "-sort",     FALSE, etBOOL, {&bSort},
          "Sort energy files (not frames)"},
        { "-rmdh",     FALSE, etBOOL, {&remove_dh},
          "Remove free energy block data" },
        { "-scalefac", FALSE, etREAL, {&scalefac},
          "Multiply energy component by this factor" },
        { "-error",    FALSE, etBOOL, {&bError},
          "Stop on errors in the file" }
    };

    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa),
                           pa, asize(desc), desc, asize(bugs), bugs, &oenv))
    {
        return 0;
    }
    tadjust  = 0;
    nremax   = 0;
    nset     = 0;
    timestep = 0.0;
    snew(fnms, argc);
    lastfilestep = 0;
    laststep     = 0;

    nfile = opt2fns(&fnms, "-f", NFILE, fnm);

    if (!nfile)
    {
        gmx_fatal(FARGS, "No input files!");
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
          "Change the number of steps" },
        { "-time",          FALSE, etREAL, {&start_t},
          "Continue from frame at this time (ps) instead of the last frame" },
        { "-zeroq",         FALSE, etBOOL, {&bZeroQ},
          "Set the charges of a group (from the index) to zero" },
        { "-vel",           FALSE, etBOOL, {&bVel},
          "Require velocities from trajectory" },
        { "-cont",          FALSE, etBOOL, {&bContinuation},
          "For exact continuation, the constraints should not be applied before the first step" },
        { "-init_fep_state", FALSE, etINT, {&init_fep_state},
          "fep state to initialize from" },
    };
    int             nerror = 0;

    /* Parse the command line */
    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
                           asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }

    /* Convert int to gmx_int64_t */
    nsteps_req = nsteps_req_int;
    bNsteps    = opt2parg_bSet("-nsteps", asize(pa), pa);
    bExtend    = opt2parg_bSet("-extend", asize(pa), pa);
    bUntil     = opt2parg_bSet("-until", asize(pa), pa);
    bFepState  = opt2parg_bSet("-init_fep_state", asize(pa), pa);
    bTime      = opt2parg_bSet("-time", asize(pa), pa);
    bTraj      = (opt2bSet("-f", NFILE, fnm) || bTime);

    top_fn = ftp2fn(efTPX, NFILE, fnm);
    fprintf(stderr, "Reading toplogy and stuff from %s\n", top_fn);

    snew(ir, 1);
    read_tpx_state(top_fn, ir, &state, NULL, &mtop);
    run_step = ir->init_step;
    run_t    = ir->init_step*ir->delta_t + ir->init_t;

    if (!EI_STATE_VELOCITY(ir->eI))
    {
        bVel = FALSE;
    }

    if (bTraj)
    {
        fprintf(stderr, "\n"
                "NOTE: Reading the state from trajectory is an obsolete feature of gmx convert-tpr.\n"
                "      Continuation should be done by loading a checkpoint file with mdrun -cpi\n"
                "      This guarantees that all state variables are transferred.\n"
                "      gmx convert-tpr is now only useful for increasing nsteps,\n"
                "      but even that can often be avoided by using mdrun -maxh\n"
                "\n");

        if (ir->bContinuation != bContinuation)
        {
            fprintf(stderr, "Modifying ir->bContinuation to %s\n",
                    bool_names[bContinuation]);
        }
        ir->bContinuation = bContinuation;


        bNeedEner = (ir->epc == epcPARRINELLORAHMAN || ir->etc == etcNOSEHOOVER);
        bReadEner = (bNeedEner && ftp2bSet(efEDR, NFILE, fnm));
        bScanEner = (bReadEner && !bTime);

        if (ir->epc != epcNO || EI_SD(ir->eI) || ir->eI == eiBD)

static game_params *dup_params(game_params *params)
{
    game_params *ret = snew(game_params);
    *ret = *params; /* struct copy */
    return ret;
}

static int strip_dssp(char *dsspfile, int nres,
                      gmx_bool bPhobres[], real t,
                      real *acc, FILE *fTArea,
                      t_matrix *mat, int average_area[],
                      const output_env_t oenv)
{
    static gmx_bool bFirst = TRUE;
    static char    *ssbuf;
    FILE           *tapeout;
    static int      xsize, frame;
    char            buf[STRLEN+1];
    char            SSTP;
    int             i, nr, iacc, nresidues;
    int             naccf, naccb; /* Count hydrophobic and hydrophilic residues */
    real            iaccf, iaccb;
    t_xpmelmt       c;

    tapeout = gmx_ffopen(dsspfile, "r");

    /* Skip header */
    do
    {
        fgets2(buf, STRLEN, tapeout);
    }
    while (strstr(buf, "KAPPA") == NULL);
    if (bFirst)
    {
        /* Since we also have empty lines in the dssp output (temp) file,
         * and some line content is saved to the ssbuf variable,
         * we need more memory than just nres elements. To be shure,
         * we allocate 2*nres-1, since for each chain there is a
         * separating line in the temp file. (At most each residue
         * could have been defined as a separate chain.) */
        snew(ssbuf, 2*nres-1);
    }

    iaccb     = iaccf = 0;
    nresidues = 0;
    naccf     = 0;
    naccb     = 0;
    for (nr = 0; (fgets2(buf, STRLEN, tapeout) != NULL); nr++)
    {
        if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
        {
            SSTP = '=';     /* Chain separator sign '=' */
        }
        else
        {
            SSTP = buf[16] == ' ' ? '~' : buf[16];
        }
        ssbuf[nr] = SSTP;

        buf[39] = '\0';

        /* Only calculate solvent accessible area if needed */
        if ((NULL != acc) && (buf[13] != '!'))
        {
            sscanf(&(buf[34]), "%d", &iacc);
            acc[nr] = iacc;
            /* average_area and bPhobres are counted from 0...nres-1 */
            average_area[nresidues] += iacc;
            if (bPhobres[nresidues])
            {
                naccb++;
                iaccb += iacc;
            }
            else
            {
                naccf++;
                iaccf += iacc;
            }
            /* Keep track of the residue number (do not count chain separator lines '!') */
            nresidues++;
        }
    }
    ssbuf[nr] = '\0';

    if (bFirst)
    {
        if (0 != acc)
        {
            fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
        }

        sprintf(mat->title, "Secondary structure");
        mat->legend[0] = 0;
        sprintf(mat->label_x, "%s", output_env_get_time_label(oenv));
        sprintf(mat->label_y, "Residue");
        mat->bDiscrete = TRUE;
        mat->ny        = nr;
        snew(mat->axis_y, nr);
        for (i = 0; i < nr; i++)
        {
            mat->axis_y[i] = i+1;
        }
        mat->axis_x = NULL;
        mat->matrix = NULL;
        xsize       = 0;
        frame       = 0;
        bFirst      = FALSE;
//.........这里部分代码省略.........

static void do_sham(const char *fn, const char *ndx,
                    const char *xpmP, const char *xpm, const char *xpm2,
                    const char *xpm3, const char *xpm4, const char *pdb,
                    const char *logf,
                    int n, int neig, real **eig,
                    gmx_bool bGE, int nenerT, real **enerT,
                    int nmap, real *mapindex, real **map,
                    real Tref,
                    real pmax, real gmax,
                    real *emin, real *emax, int nlevels, real pmin,
                    const char *mname, gmx_bool bSham, int *idim, int *ibox,
                    gmx_bool bXmin, real *xmin, gmx_bool bXmax, real *xmax)
{
    FILE        *fp;
    real        *min_eig, *max_eig;
    real        *axis_x, *axis_y, *axis_z, *axis = NULL;
    double      *P;
    real       **PP, *W, *E, **WW, **EE, *S, **SS, *M, **MM, *bE;
    rvec         xxx;
    char        *buf;
    double      *bfac, efac, bref, Pmax, Wmin, Wmax, Winf, Emin, Emax, Einf, Smin, Smax, Sinf, Mmin, Mmax, Minf;
    real        *delta;
    int          i, j, k, imin, len, index, d, *nbin, *bindex, bi;
    int         *nxyz, maxbox;
    t_blocka    *b;
    gmx_bool     bOutside;
    unsigned int flags;
    t_rgb        rlo  = { 0, 0, 0 };
    t_rgb        rhi  = { 1, 1, 1 };

    /* Determine extremes for the eigenvectors */
    snew(min_eig, neig);
    snew(max_eig, neig);
    snew(nxyz, neig);
    snew(bfac, neig);
    snew(delta, neig);

    for (i = 0; (i < neig); i++)
    {
        /* Check for input constraints */
        min_eig[i] = max_eig[i] = eig[i][0];
        for (j = 0; (j < n); j++)
        {
            min_eig[i] = min(min_eig[i], eig[i][j]);
            max_eig[i] = max(max_eig[i], eig[i][j]);
            delta[i]   = (max_eig[i]-min_eig[i])/(2.0*ibox[i]);
        }
        /* Add some extra space, half a bin on each side, unless the
         * user has set the limits.
         */
        if (bXmax)
        {
            if (max_eig[i] > xmax[i])
            {
                gmx_warning("Your xmax[%d] value %f is smaller than the largest data point %f", i, xmax[i], max_eig[i]);
            }
            max_eig[i] = xmax[i];
        }
        else
        {
            max_eig[i] += delta[i];
        }

        if (bXmin)
        {
            if (min_eig[i] < xmin[i])
            {
                gmx_warning("Your xmin[%d] value %f is larger than the smallest data point %f", i, xmin[i], min_eig[i]);
            }
            min_eig[i] = xmin[i];
        }
        else
        {
            min_eig[i] -= delta[i];
        }
        bfac[i]     = ibox[i]/(max_eig[i]-min_eig[i]);
    }
    /* Do the binning */
    bref = 1/(BOLTZ*Tref);
    snew(bE, n);
    if (bGE || nenerT == 2)
    {
        Emin = 1e8;
        for (j = 0; (j < n); j++)
        {
            if (bGE)
            {
                bE[j] = bref*enerT[0][j];
            }
            else
            {
                bE[j] = (bref - 1/(BOLTZ*enerT[1][j]))*enerT[0][j];
            }
            Emin  = min(Emin, bE[j]);
        }
    }
    else
    {
        Emin = 0;
    }
//.........这里部分代码省略.........

int main(int argc,char *argv[])
{
  static char *desc[] = {
    "The [TT]pmetest[tt] program tests the scaling of the PME code. When only given",
    "a [TT].tpr[tt] file it will compute PME for one frame. When given a trajectory",
    "it will do so for all the frames in the trajectory. Before the PME",
    "routine is called the coordinates are sorted along the X-axis.[PAR]",
    "As an extra service to the public the program can also compute",
    "long-range Coulomb energies for components of the system. When the",
    "[TT]-groups[tt] flag is given to the program the energy groups",
    "from the [TT].tpr[tt] file will be read, and half an energy matrix computed."
  };
  t_commrec    *cr,*mcr;
  static t_filenm fnm[] = {
    { efTPX, NULL,      NULL,       ffREAD  },
    { efTRN, "-o",      NULL,       ffWRITE },
    { efLOG, "-g",      "pme",      ffWRITE },
    { efTRX, "-f",      NULL,       ffOPTRD },
    { efXVG, "-x",      "ener-pme", ffWRITE }
  };
#define NFILE asize(fnm)

  /* Command line options ! */
  static gmx_bool bVerbose=FALSE;
  static gmx_bool bOptFFT=FALSE;
  static gmx_bool bSort=FALSE;
  static int  ewald_geometry=eewg3D;
  static int  nnodes=1;
  static int  pme_order=0;
  static rvec grid = { -1, -1, -1 };
  static real rc   = 0.0;
  static real dtol = 0.0;
  static gmx_bool bGroups = FALSE;
  static t_pargs pa[] = {
    { "-np",      FALSE, etINT, {&nnodes},
      "Number of nodes, must be the same as used for [TT]grompp[tt]" },
    { "-v",       FALSE, etBOOL,{&bVerbose},  
      "Be loud and noisy" },
    { "-sort",    FALSE, etBOOL,{&bSort},  
      "Sort coordinates. Crucial for domain decomposition." },
    { "-grid",    FALSE, etRVEC,{&grid},
      "Number of grid cells in X, Y, Z dimension (if -1 use from [TT].tpr[tt])" },
    { "-order",   FALSE, etINT, {&pme_order},
      "Order of the PME spreading algorithm" },
    { "-groups",  FALSE, etBOOL, {&bGroups},
      "Compute half an energy matrix based on the energy groups in your [TT].tpr[tt] file" },
    { "-rc",      FALSE, etREAL, {&rc},
      "Rcoulomb for Ewald summation" },
    { "-tol",     FALSE, etREAL, {&dtol},
      "Toleranc