bool update_forcefield(FILE *fplog,
                       int nfile,t_filenm fnm[],t_forcerec *fr,
                       int natoms,rvec x[],matrix box)
{
    static int ntry,ntried;
    int    i,j;
    bool   bDone;

    /* First time around we have to read the parameters */
    if (nparm == 0) {
        range = read_range(ftp2fn(efDAT,nfile,fnm),&nparm);
        if (nparm == 0)
            gmx_fatal(FARGS,"No correct parameter info in %s",ftp2fn(efDAT,nfile,fnm));
        snew(param_val,nparm);

        if (opt2bSet("-ga",nfile,fnm)) {
            /* Genetic algorithm time */
            ga = init_ga(fplog,opt2fn("-ga",nfile,fnm),nparm,range);
        }
        else {
            /* Determine the grid size */
            ntry = 1;
            for(i=0; (i<nparm); i++)
                ntry *= range[i].np;
            ntried = 0;

            fprintf(fplog,"Going to try %d different combinations of %d parameters\n",
                    ntry,nparm);
        }
    }
    if (ga) {
        update_ga(fplog,range,ga);
    }
    else {
        /* Increment the counter
         * Non-trivial, since this is nparm nested loops in principle
         */
        for(i=0; (i<nparm); i++) {
            if (param_val[i] < (range[i].np-1)) {
                param_val[i]++;
                for(j=0; (j<i); j++)
                    param_val[j] = 0;
                ntried++;
                break;
            }
        }
        if (i == nparm) {
            fprintf(fplog,"Finished with %d out of %d iterations\n",ntried+1,ntry);
            return TRUE;
        }
    }

    /* Now do the real updating */
    update_ff(fr,nparm,range,param_val);

    /* Update box and coordinates if necessary */
    scale_box(natoms,x,box);

    return FALSE;
}

int cmain(int argc,char *argv[])
{
  int        step,natoms;
  real       t,lambda;
  t_inputrec ir;
  t_topology top;
  matrix     box;
  rvec       *x,*v,*f;
  t_filenm fnm[] = {
    { efTPX, NULL, NULL, ffREAD  },
    { efXML, "-r", NULL, ffREAD  },
    { efXML, "-o", NULL, ffWRITE }
  };  
#define NFILE asize(fnm)

  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,0,NFILE,fnm,0,NULL,0,NULL,0,NULL);
  
  init_top(&top);
  if (opt2bSet("-r",NFILE,fnm))
    read_xml(opt2fn("-r",NFILE,fnm),&step,&t,&lambda,&ir,
	     box,&natoms,&x,&v,&f,&top);
  else {
    t_tpxheader tpx;
    
    read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&tpx,FALSE);
    snew(x,tpx.natoms);
    snew(v,tpx.natoms);
    f = NULL;
    read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,&ir,
	     box,&natoms,x,v,f,&top);
  }
  /*write_xml(opt2fn("-o",NFILE,fnm),step,t,lambda,&ir,box,natoms,x,v,f,&top);*/
  
  return 0;
}

//.........这里部分代码省略.........
        xvgr_legend(fp, 2, leg, oenv);
        for (i = 0; (i < isize); i++)
        {
            if (!bRes || i+1 == isize ||
                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)
            {
                resind    = top.atoms.atom[index[i]].resind;
                pdb_bfac  = find_pdb_bfac(pdbatoms, &top.atoms.resinfo[resind],
                                          *(top.atoms.atomname[index[i]]));

                fprintf(fp, "%5d  %10.5f  %10.5f\n",
                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, rmsf[i]*bfac,
                        pdb_bfac);
            }
        }
        xvgrclose(fp);
    }
    else
    {
        fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS fluctuation", label, "(nm)", oenv);
        for (i = 0; i < isize; i++)
        {
            if (!bRes || i+1 == isize ||
                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)
            {
                fprintf(fp, "%5d %8.4f\n",
                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, std::sqrt(rmsf[i]));
            }
        }
        xvgrclose(fp);
    }

    for (i = 0; i < isize; i++)
    {
        pdbatoms->pdbinfo[index[i]].bfac = 800*M_PI*M_PI/3.0*rmsf[i];
    }

    if (devfn)
    {
        for (i = 0; i < isize; i++)
        {
            rmsf[i] = (rmsd_x[i][XX]+rmsd_x[i][YY]+rmsd_x[i][ZZ])/count;
        }
        if (bRes)
        {
            average_residues(rmsf, nullptr, 0, isize, index, w_rls, &top.atoms);
        }
        /* Write RMSD output */
        fp = xvgropen(devfn, "RMS Deviation", 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);
    }

    if (opt2bSet("-oq", NFILE, fnm))
    {
        /* Write a .pdb file with B-factors and optionally anisou records */
        for (i = 0; i < isize; i++)
        {
            rvec_inc(pdbx[index[i]], xcm);
        }
        write_sto_conf_indexed(opt2fn("-oq", NFILE, fnm), title, pdbatoms, pdbx,
                               nullptr, ePBC, pdbbox, isize, index);
    }
    if (opt2bSet("-ox", NFILE, fnm))
    {
        rvec *bFactorX;
        snew(bFactorX, top.atoms.nr);
        for (i = 0; i < isize; i++)
        {
            for (d = 0; d < DIM; d++)
            {
                bFactorX[index[i]][d] = xcm[d] + xav[i*DIM + d];
            }
        }
        /* Write a .pdb file with B-factors and optionally anisou records */
        write_sto_conf_indexed(opt2fn("-ox", NFILE, fnm), title, pdbatoms, bFactorX, nullptr,
                               ePBC, pdbbox, isize, index);
        sfree(bFactorX);
    }
    if (bAniso)
    {
        correlate_aniso(opt2fn("-oc", NFILE, fnm), refatoms, pdbatoms, oenv);
        do_view(oenv, opt2fn("-oc", NFILE, fnm), "-nxy");
    }
    do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
    if (devfn)
    {
        do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
    }

    return 0;
}

//.........这里部分代码省略.........
        while (qdelta != 0)
        {
            while (qdelta < 0)
            {
                p_num++;
                qdelta += p_q;
            }
            while (qdelta > 0)
            {
                n_num++;
                qdelta += n_q;
            }
        }
    }

    if ((p_num == 0) && (n_num == 0))
    {
        fprintf(stderr, "No ions to add, will just copy input configuration.\n");
    }
    else
    {
        printf("Will try to add %d %s ions and %d %s ions.\n",
               p_num, p_name, n_num, n_name);
        printf("Select a continuous group of solvent molecules\n");
        get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
        for (i = 1; i < nwa; i++)
        {
            if (index[i] != index[i-1]+1)
            {
                gmx_fatal(FARGS, "The solvent group %s is not continuous: "
                          "index[%d]=%d, index[%d]=%d",
                          grpname, i, index[i-1]+1, i+1, index[i]+1);
            }
        }
        nsa = 1;
        while ((nsa < nwa) &&
               (atoms.atom[index[nsa]].resind ==
                atoms.atom[index[nsa-1]].resind))
        {
            nsa++;
        }
        if (nwa % nsa)
        {
            gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
                      nwa, nsa);
        }
        nw = nwa/nsa;
        fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
        if (p_num+n_num > nw)
        {
            gmx_fatal(FARGS, "Not enough solvent for adding ions");
        }

        if (opt2bSet("-p", NFILE, fnm))
        {
            update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
        }

        snew(bSet, nw);
        snew(repl, nw);

        snew(v, atoms.nr);
        snew(atoms.pdbinfo, atoms.nr);

        set_pbc(&pbc, ePBC, box);

        if (seed == 0)
        {
            rng = gmx_rng_init(gmx_rng_make_seed());
        }
        else
        {
            rng = gmx_rng_init(seed);
        }
        /* Now loop over the ions that have to be placed */
        while (p_num-- > 0)
        {
            insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
                       1, p_q, p_name, &atoms, rmin, rng);
        }
        while (n_num-- > 0)
        {
            insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
                       -1, n_q, n_name, &atoms, rmin, rng);
        }
        gmx_rng_destroy(rng);
        fprintf(stderr, "\n");

        if (nw)
        {
            sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
        }

        sfree(atoms.pdbinfo);
        atoms.pdbinfo = NULL;
    }
    write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);

    return 0;
}

int gmx_pme_error(int argc,char *argv[])
{
    const char *desc[] = {
            "[TT]g_pme_error[tt] 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[] = {
      { efTPX, "-s",     NULL,    ffREAD },
      { efOUT, "-o",    "error",  ffWRITE },
      { efTPX, "-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_par(&argc,&argv);
    
#ifdef GMX_LIB_MPI
    MPI_Barrier(MPI_COMM_WORLD);
#endif

    if (MASTER(cr))
      CopyRight(stderr,argv[0]);
    
    PCA_Flags = PCA_NOEXIT_ON_ARGS;
    PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
    
    parse_common_args(&argc,argv,PCA_Flags,
                      NFILE,fnm,asize(pa),pa,asize(desc),desc,
                      0,NULL,&oenv);        

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

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

    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;
    }

int gmx_mdmat(int argc, char *argv[])
{
    const char     *desc[] = {
        "[TT]g_mdmat[tt] makes distance matrices consisting of the smallest distance",
        "between residue pairs. With [TT]-frames[tt], these distance matrices can be",
        "stored in order to see differences in tertiary structure as a",
        "function of time. If you choose your options unwisely, this may generate",
        "a large output file. By default, only an averaged matrix over the whole",
        "trajectory is output.",
        "Also a count of the number of different atomic contacts between",
        "residues over the whole trajectory can be made.",
        "The output can be processed with [TT]xpm2ps[tt] to make a PostScript (tm) plot."
    };
    static real     truncate = 1.5;
    static gmx_bool bAtom    = FALSE;
    static int      nlevels  = 40;
    t_pargs         pa[]     = {
        { "-t",   FALSE, etREAL, {&truncate},
          "trunc distance" },
        { "-nlevels",   FALSE, etINT,  {&nlevels},
          "Discretize distance in this number of levels" }
    };
    t_filenm        fnm[] = {
        { efTRX, "-f",  NULL, ffREAD },
        { efTPS, NULL,  NULL, ffREAD },
        { efNDX, NULL,  NULL, ffOPTRD },
        { efXPM, "-mean", "dm", ffWRITE },
        { efXPM, "-frames", "dmf", ffOPTWR },
        { efXVG, "-no", "num", ffOPTWR },
    };
#define NFILE asize(fnm)

    FILE          *out = NULL, *fp;
    t_topology     top;
    int            ePBC;
    t_atoms        useatoms;
    int            isize;
    atom_id       *index;
    char          *grpname;
    int           *rndx, *natm, prevres, newres;

    int            i, j, nres, natoms, nframes, it, trxnat;
    t_trxstatus   *status;
    int            nr0;
    gmx_bool       bCalcN, bFrames;
    real           t, ratio;
    char           title[256], label[234];
    t_rgb          rlo, rhi;
    rvec          *x;
    real         **mdmat, *resnr, **totmdmat;
    int          **nmat, **totnmat;
    real          *mean_n;
    int           *tot_n;
    matrix         box;
    output_env_t   oenv;
    gmx_rmpbc_t    gpbc = NULL;

    CopyRight(stderr, argv[0]);

    parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE, NFILE, fnm,
                      asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);

    fprintf(stderr, "Will truncate at %f nm\n", truncate);
    bCalcN  = opt2bSet("-no", NFILE, fnm);
    bFrames = opt2bSet("-frames", NFILE, fnm);
    if (bCalcN)
    {
        fprintf(stderr, "Will calculate number of different contacts\n");
    }

    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box, FALSE);

    fprintf(stderr, "Select group for analysis\n");
    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);

    natoms = isize;
    snew(useatoms.atom, natoms);
    snew(useatoms.atomname, natoms);

    useatoms.nres = 0;
    snew(useatoms.resinfo, natoms);

    prevres = top.atoms.atom[index[0]].resind;
    newres  = 0;
    for (i = 0; (i < isize); i++)
    {
        int ii = index[i];
        useatoms.atomname[i] = top.atoms.atomname[ii];
        if (top.atoms.atom[ii].resind != prevres)
        {
            prevres = top.atoms.atom[ii].resind;
            newres++;
            useatoms.resinfo[i] = top.atoms.resinfo[prevres];
            if (debug)
            {
                fprintf(debug, "New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
                        *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
                        *(top.atoms.atomname[ii]),
                        ii, i, newres);
            }
//.........这里部分代码省略.........

int main (int argc,char *argv[])
{
    const char *desc[] = {
        "[TT]g_protonate[tt] reads (a) conformation(s) and adds all missing",
        "hydrogens as defined in [TT]gmx2.ff/aminoacids.hdb[tt]. If only [TT]-s[tt] is",
        "specified, this conformation will be protonated, if also [TT]-f[tt]",
        "is specified, the conformation(s) will be read from this file, ",
        "which can be either a single conformation or a trajectory.",
        "[PAR]",
        "If a [TT].pdb[tt] file is supplied, residue names might not correspond to",
        "to the GROMACS naming conventions, in which case these residues will",
        "probably not be properly protonated.",
        "[PAR]",
        "If an index file is specified, please note that the atom numbers",
        "should correspond to the [BB]protonated[bb] state."
    };

    char        title[STRLEN+1];
    const char  *infile;
    char        *grpnm;
    t_topology  top;
    int         ePBC;
    t_atoms     *atoms,*iatoms;
    t_protonate protdata;
    atom_id     *index;
    t_trxstatus *status;
    t_trxstatus *out;
    t_trxframe  fr,frout;
    rvec        *x,*ix;
    int         nidx,natoms,natoms_out;
    matrix      box;
    int         i,frame,resind;
    gmx_bool        bReadMultiple;
    output_env_t oenv;

    const char *bugs[] = {
        "For the moment, only .pdb files are accepted to the -s flag"
    };

    t_filenm fnm[] = {
        { efTPS, NULL, NULL,         ffREAD  },
        { efTRX, "-f", NULL,         ffOPTRD },
        { efNDX, NULL, NULL,         ffOPTRD },
        { efTRO, "-o", "protonated", ffWRITE }
    };
#define NFILE asize(fnm)

    CopyRight(stderr,argv[0]);
    parse_common_args(&argc,argv,PCA_CAN_TIME,
                      NFILE,fnm,0,NULL,asize(desc),desc,asize(bugs),bugs,&oenv);

    infile=opt2fn("-s",NFILE,fnm);
    read_tps_conf(infile,title,&top,&ePBC,&x,NULL,box,FALSE);
    atoms=&(top.atoms);
    printf("Select group to process:\n");
    get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&nidx,&index,&grpnm);
    bReadMultiple = opt2bSet("-f",NFILE,fnm);
    if (bReadMultiple) {
        infile = opt2fn("-f",NFILE,fnm);
        if ( !read_first_frame(oenv,&status, infile, &fr, TRX_NEED_X ) ) {
            gmx_fatal(FARGS,"cannot read coordinate file %s",infile);
        }
        natoms = fr.natoms;
    } else {
        clear_trxframe(&fr,TRUE);
        fr.natoms = atoms->nr;
        fr.bTitle = TRUE;
        fr.title  = title;
        fr.bX     = TRUE;
        fr.x      = x;
        fr.bBox   = TRUE;
        copy_mat(box, fr.box);
        natoms = fr.natoms;
    }

    /* check input */
    if ( natoms == 0 ) {
        gmx_fatal(FARGS,"no atoms in coordinate file %s",infile);
    }

    if ( natoms > atoms->nr ) {
        gmx_fatal(FARGS,"topology with %d atoms does not match "
                  "coordinates with %d atoms",atoms->nr,natoms);
    }

    for(i=0; i<nidx; i++) {
        if (index[i] > natoms) {
            gmx_fatal(FARGS,"An atom number in group %s is larger than the number of "
                      "atoms (%d) in the coordinate file %s",grpnm,natoms,infile);
        }
    }

    /* get indexed copy of atoms */
    snew(iatoms,1);
    init_t_atoms(iatoms,nidx,FALSE);
    snew(iatoms->atom, iatoms->nr);
    resind = 0;
    for(i=0; i<nidx; i++) {
        iatoms->atom[i] = atoms->atom[index[i]];
        iatoms->atomname[i] = atoms->atomname[index[i]];
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    for (i = 0; (i < 2); i++)
    {
        for (j = 0; (j < isize[i]); j++)
        {
            if (index[i][j] > natoms)
            {
                gmx_fatal(FARGS, "An atom number in group %s is larger than the number of atoms in the trajectory");
            }
        }
    }

    if ((isize_sol % na) != 0)
    {
        gmx_fatal(FARGS, "Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
                  isize[1], na);
    }

    nwat = isize_sol/na;
    if (ref_a > na)
    {
        gmx_fatal(FARGS, "The reference atom can not be larger than the number of atoms in a molecule");
    }
    ref_a--;
    snew(xsol, nwat);
    snew(order, nwat);
    snew(swi, natoms);
    for (i = 0; (i < natoms); i++)
    {
        swi[i] = i;
    }

    out     = NULL;
    fp      = NULL;
    bNShell = ((opt2bSet("-nshell", NFILE, fnm)) ||
               (opt2parg_bSet("-r", asize(pa), pa)));
    bPDBout = FALSE;
    if (bNShell)
    {
        rcut2   = rcut*rcut;
        fp      = xvgropen(opt2fn("-nshell", NFILE, fnm), "Number of molecules",
                           "Time (ps)", "N", oenv);
        printf("Will compute the number of molecules within a radius of %g\n",
               rcut);
    }
    if (!bNShell || opt2bSet("-o", NFILE, fnm))
    {
        bPDBout = (fn2ftp(opt2fn("-o", NFILE, fnm)) == efPDB);
        if (bPDBout && !top.atoms.pdbinfo)
        {
            fprintf(stderr, "Creating pdbfino records\n");
            snew(top.atoms.pdbinfo, top.atoms.nr);
        }
        out = open_trx(opt2fn("-o", NFILE, fnm), "w");
    }
    gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
    do
    {
        gmx_rmpbc(gpbc, natoms, box, x);
        set_pbc(&pbc, ePBC, box);

        if (ref_a == -1)
        {
            /* Calculate the COM of all solvent molecules */
            for (i = 0; i < nwat; i++)
            {
                totmass = 0;

//.........这里部分代码省略.........
        { efMTX, "-f", "hessian",    ffREAD  },
        { efTPX, NULL, NULL,         ffREAD  },
        { efXVG, "-of", "eigenfreq", ffWRITE },
        { efXVG, "-ol", "eigenval",  ffWRITE },
        { efXVG, "-os", "spectrum",  ffOPTWR },
        { efXVG, "-qc", "quant_corr",  ffOPTWR },
        { efTRN, "-v", "eigenvec",  ffWRITE }
    };
#define NFILE asize(fnm)

    parse_common_args(&argc, argv, PCA_BE_NICE,
                      NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);

    /* Read tpr file for volume and number of harmonic terms */
    read_tpxheader(ftp2fn(efTPX, NFILE, fnm), &tpx, TRUE, &version, &generation);
    snew(top_x, tpx.natoms);

    read_tpx(ftp2fn(efTPX, NFILE, fnm), NULL, box, &natoms,
             top_x, NULL, NULL, &mtop);
    if (bCons)
    {
        nharm = get_nharm(&mtop, &nvsite);
    }
    else
    {
        nharm  = 0;
        nvsite = 0;
    }
    top = gmx_mtop_t_to_t_topology(&mtop);

    bM   = TRUE;
    ndim = DIM*natoms;

    if (opt2bSet("-qc", NFILE, fnm))
    {
        begin = 7+DIM*nvsite;
        end   = DIM*natoms;
    }
    if (begin < 1)
    {
        begin = 1;
    }
    if (end > ndim)
    {
        end = ndim;
    }
    printf("Using begin = %d and end = %d\n", begin, end);

    /*open Hessian matrix */
    gmx_mtxio_read(ftp2fn(efMTX, NFILE, fnm), &nrow, &ncol, &full_hessian, &sparse_hessian);

    /* Memory for eigenvalues and eigenvectors (begin..end) */
    snew(eigenvalues, nrow);
    snew(eigenvectors, nrow*(end-begin+1));

    /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
     * and they must start at the first one. If this is not valid we convert to full matrix
     * storage, but warn the user that we might run out of memory...
     */
    if ((sparse_hessian != NULL) && (begin != 1 || end == ndim))
    {
        if (begin != 1)
        {
            fprintf(stderr, "Cannot use sparse Hessian with first eigenvector != 1.\n");
        }
        else if (end == ndim)

int gmx_gyrate(int argc, char *argv[])
{
    const char       *desc[] = {
        "[THISMODULE] computes the radius of gyration of a molecule",
        "and the radii of gyration about the [IT]x[it]-, [IT]y[it]- and [IT]z[it]-axes,",
        "as a function of time. The atoms are explicitly mass weighted.[PAR]",
        "The axis components corresponds to the mass-weighted root-mean-square",
        "of the radii components orthogonal to each axis, for example:[PAR]",
        "Rg(x) = sqrt((sum_i m_i (R_i(y)^2 + R_i(z)^2))/(sum_i m_i)).[PAR]",
        "With the [TT]-nmol[tt] option the radius of gyration will be calculated",
        "for multiple molecules by splitting the analysis group in equally",
        "sized parts.[PAR]",
        "With the option [TT]-nz[tt] 2D radii of gyration in the [IT]x-y[it] plane",
        "of slices along the [IT]z[it]-axis are calculated."
    };
    static int        nmol = 1, nz = 0;
    static gmx_bool   bQ   = FALSE, bRot = FALSE, bMOI = FALSE;
    t_pargs           pa[] = {
        { "-nmol", FALSE, etINT, {&nmol},
          "The number of molecules to analyze" },
        { "-q", FALSE, etBOOL, {&bQ},
          "Use absolute value of the charge of an atom as weighting factor instead of mass" },
        { "-p", FALSE, etBOOL, {&bRot},
          "Calculate the radii of gyration about the principal axes." },
        { "-moi", FALSE, etBOOL, {&bMOI},
          "Calculate the moments of inertia (defined by the principal axes)." },
        { "-nz", FALSE, etINT, {&nz},
          "Calculate the 2D radii of gyration of this number of slices along the z-axis" },
    };
    FILE             *out;
    t_trxstatus      *status;
    t_topology        top;
    int               ePBC;
    rvec             *x, *x_s;
    rvec              xcm, gvec, gvec1;
    matrix            box, trans;
    gmx_bool          bACF;
    real            **moi_trans = NULL;
    int               max_moi   = 0, delta_moi = 100;
    rvec              d, d1; /* eigenvalues of inertia tensor */
    real              t, t0, tm, gyro;
    int               natoms;
    char             *grpname;
    int               j, m, gnx, nam, mol;
    atom_id          *index;
    gmx_output_env_t *oenv;
    gmx_rmpbc_t       gpbc   = NULL;
    const char       *leg[]  = { "Rg", "Rg\\sX\\N", "Rg\\sY\\N", "Rg\\sZ\\N" };
    const char       *legI[] = { "Itot", "I1", "I2", "I3" };
#define NLEG asize(leg)
    t_filenm          fnm[] = {
        { efTRX, "-f",   NULL,       ffREAD },
        { efTPS, NULL,   NULL,       ffREAD },
        { efNDX, NULL,   NULL,       ffOPTRD },
        { efXVG, NULL,   "gyrate",   ffWRITE },
        { efXVG, "-acf", "moi-acf",  ffOPTWR },
    };
#define NFILE asize(fnm)
    int               npargs;
    t_pargs          *ppa;

    npargs = asize(pa);
    ppa    = add_acf_pargs(&npargs, pa);

    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
                           NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }
    bACF = opt2bSet("-acf", NFILE, fnm);
    if (bACF && nmol != 1)
    {
        gmx_fatal(FARGS, "Can only do acf with nmol=1");
    }
    bRot = bRot || bMOI || bACF;
    /*
       if (nz > 0)
       bMOI = TRUE;
     */
    if (bRot)
    {
        printf("Will rotate system along principal axes\n");
        snew(moi_trans, DIM);
    }
    if (bMOI)
    {
        printf("Will print moments of inertia\n");
        bQ = FALSE;
    }
    if (bQ)
    {
        printf("Will print radius normalised by charge\n");
    }

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

int gmx_disre(int argc,char *argv[])
{
  const char *desc[] = {
    "g_disre computes violations of distance restraints.",
    "If necessary all protons can be added to a protein molecule ",
    "using the protonate program.[PAR]",
    "The program always",
    "computes the instantaneous violations rather than time-averaged,",
    "because this analysis is done from a trajectory file afterwards",
    "it does not make sense to use time averaging. However,",
    "the time averaged values per restraint are given in the log file.[PAR]",
    "An index file may be used to select specific restraints for",
    "printing.[PAR]",
    "When the optional[TT]-q[tt] flag is given a pdb file coloured by the",
    "amount of average violations.[PAR]",
    "When the [TT]-c[tt] option is given, an index file will be read",
    "containing the frames in your trajectory corresponding to the clusters",
    "(defined in another manner) that you want to analyze. For these clusters",
    "the program will compute average violations using the third power",
    "averaging algorithm and print them in the log file."
  };
  static int  ntop      = 0;
  static int  nlevels   = 20;
  static real max_dr    = 0;
  static gmx_bool bThird    = TRUE;
  t_pargs pa[] = {
    { "-ntop", FALSE, etINT,  {&ntop},
      "Number of large violations that are stored in the log file every step" },
    { "-maxdr", FALSE, etREAL, {&max_dr},
      "Maximum distance violation in matrix output. If less than or equal to 0 the maximum will be determined by the data." },
    { "-nlevels", FALSE, etINT, {&nlevels},
      "Number of levels in the matrix output" },
    { "-third", FALSE, etBOOL, {&bThird},
      "Use inverse third power averaging or linear for matrix output" }
  };
  
  FILE        *out=NULL,*aver=NULL,*numv=NULL,*maxxv=NULL,*xvg=NULL;
  t_tpxheader header;
  t_inputrec  ir;
  gmx_mtop_t  mtop;
  rvec        *xtop;
  gmx_localtop_t *top;
  t_atoms     *atoms=NULL;
  t_forcerec  *fr;
  t_fcdata    fcd;
  t_nrnb      nrnb;
  t_commrec   *cr;
  t_graph     *g;
  int         ntopatoms,natoms,i,j,kkk;
  t_trxstatus *status;
  real        t;
  rvec        *x,*f,*xav=NULL;
  matrix      box;
  gmx_bool        bPDB;
  int         isize;
  atom_id     *index=NULL,*ind_fit=NULL;
  char        *grpname;
  t_cluster_ndx *clust=NULL;
  t_dr_result dr,*dr_clust=NULL;
  char        **leg;
  real        *vvindex=NULL,*w_rls=NULL;
  t_mdatoms   *mdatoms;
  t_pbc       pbc,*pbc_null;
  int         my_clust;
  FILE        *fplog;
  output_env_t oenv;
  gmx_rmpbc_t  gpbc=NULL;
  
  t_filenm fnm[] = {
    { efTPX, NULL, NULL, ffREAD },
    { efTRX, "-f", NULL, ffREAD },
    { efXVG, "-ds", "drsum",  ffWRITE },
    { efXVG, "-da", "draver", ffWRITE },
    { efXVG, "-dn", "drnum",  ffWRITE },
    { efXVG, "-dm", "drmax",  ffWRITE },
    { efXVG, "-dr", "restr",  ffWRITE },
    { efLOG, "-l",  "disres", ffWRITE },
    { efNDX, NULL,  "viol",   ffOPTRD },
    { efPDB, "-q",  "viol",   ffOPTWR },
    { efNDX, "-c",  "clust",  ffOPTRD },
    { efXPM, "-x",  "matrix", ffOPTWR }
  };
#define NFILE asize(fnm)

  cr  = init_par(&argc,&argv);
  CopyRight(stderr,argv[0]);
  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);

  gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,FALSE,0,&fplog);
  
  if (ntop)
    init5(ntop);
  
  read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header,FALSE,NULL,NULL);
  snew(xtop,header.natoms);
  read_tpx(ftp2fn(efTPX,NFILE,fnm),&ir,box,&ntopatoms,xtop,NULL,NULL,&mtop);
  bPDB = opt2bSet("-q",NFILE,fnm);
  if (bPDB) {
    snew(xav,ntopatoms);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        skfnm = opt2fn("-Sk", NFILE, fnm);
        calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
                              opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm),
                              oenv);
        /* view xvgr files */
        do_view(oenv, opt2fn("-Sg", NFILE, fnm), NULL);
        do_view(oenv, opt2fn("-Sk", NFILE, fnm), NULL);
        if (nslices > 1)
        {
            do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), NULL);
            do_view(oenv, opt2fn("-Sksl", NFILE, fnm), NULL);
        }
    }
    else
    {
        /* tail order parameter */

        if (nslices > 1)
        {
            bSliced = TRUE;
            fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
        }

        if (bSzonly)
        {
            fprintf(stderr, "Only calculating Sz\n");
        }
        if (bUnsat)
        {
            fprintf(stderr, "Taking carbons as unsaturated!\n");
        }

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

        block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
        index = block->index;                   /* get indices from t_block block */
        a     = block->a;                       /* see block.h                    */
        ngrps = block->nr;

        if (permolecule)
        {
            nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
            fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
        }

        if (radial)
        {
            fprintf(stderr, "Calculating radial distances\n");
            if (!permolecule)
            {
                gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
            }
        }

        /* show atomtypes, to check if index file is correct */
        print_types(index, a, ngrps, grpname, top);

        calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order,
                   &slOrder, &slWidth, nslices, bSliced, bUnsat,
                   top, ePBC, ngrps, axis, permolecule, radial, distcalc, opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);

        if (radial)
        {
            ngrps--; /*don't print the last group--was used for
                               center-of-mass determination*/

        }
        order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
                   opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule, distvals, oenv);

        if (opt2bSet("-ob", NFILE, fnm))
        {
            if (!permolecule)
            {
                fprintf(stderr,
                        "Won't write B-factors with averaged order parameters; use -permolecule\n");
            }
            else
            {
                write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
            }
        }


        do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);  /* view xvgr file */
        do_view(oenv, opt2fn("-os", NFILE, fnm), NULL); /* view xvgr file */
        do_view(oenv, opt2fn("-od", NFILE, fnm), NULL); /* view xvgr file */
    }

    if (distvals != NULL)
    {
        for (i = 0; i < nslices; ++i)
        {
            sfree(distvals[i]);
        }
        sfree(distvals);
    }

    return 0;
}

//.........这里部分代码省略.........
    // further (or the working directory changed), so we can't check for their
    // existence during parsing.  It isn't useful to do any completion based on
    // file system contents, either.
    if (is_multisim_option_set(argc, argv))
    {
        PCA_Flags |= PCA_DISABLE_INPUT_FILE_CHECKING;
    }

    /* Comment this in to do fexist calls only on master
     * works not with rerun or tables at the moment
     * also comment out the version of init_forcerec in md.c
     * with NULL instead of opt2fn
     */
    /*
       if (!MASTER(cr))
       {
       PCA_Flags |= PCA_NOT_READ_NODE;
       }
     */

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


    dd_rank_order = nenum(ddrank_opt);

    hw_opt.thread_affinity = nenum(thread_aff_opt);

    /* now check the -multi and -multidir option */
    if (opt2bSet("-multidir", NFILE, fnm))
    {
        if (nmultisim > 0)
        {
            gmx_fatal(FARGS, "mdrun -multi and -multidir options are mutually exclusive.");
        }
        nmultisim = opt2fns(&multidir, "-multidir", NFILE, fnm);
    }


    if (repl_ex_nst != 0 && nmultisim < 2)
    {
        gmx_fatal(FARGS, "Need at least two replicas for replica exchange (option -multi)");
    }

    if (repl_ex_nex < 0)
    {
        gmx_fatal(FARGS, "Replica exchange number of exchanges needs to be positive");
    }

    if (nmultisim >= 1)
    {
#if !GMX_THREAD_MPI
        init_multisystem(cr, nmultisim, multidir, NFILE, fnm);
#else
        gmx_fatal(FARGS, "mdrun -multi or -multidir are not supported with the thread-MPI library. "
                  "Please compile GROMACS with a proper external MPI library.");
#endif
    }

    if (!opt2bSet("-cpi", NFILE, fnm))
    {
        // If we are not starting from a checkpoint we never allow files to be appended

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

        if (0 != system(dssp))
        {
            gmx_fatal(FARGS, "Failed to execute command: %s\n",
                      "Try specifying your dssp version with the -ver option.", dssp);
        }

        /* strip_dssp returns the number of lines found in the dssp file, i.e.
         * the number of residues plus the separator lines */

        if (bDoAccSurf)
        {
            accr_ptr = accr[nframe];
        }

        nres_plus_separators = strip_dssp(tmpfile, nres, b