int gmx_genpr(int argc,char *argv[])
{
  const char *desc[] = {
    "[TT]genrestr[tt] produces an include file for a topology containing",
    "a list of atom numbers and three force constants for the",
    "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-direction. A single isotropic force constant may",
    "be given on the command line instead of three components.[PAR]",
    "WARNING: position restraints only work for the one molecule at a time.",
    "Position restraints are interactions within molecules, therefore",
    "they should be included within the correct [TT][ moleculetype ][tt]",
    "block in the topology. Since the atom numbers in every moleculetype",
    "in the topology start at 1 and the numbers in the input file for",
    "[TT]genrestr[tt] number consecutively from 1, [TT]genrestr[tt] will only",
    "produce a useful file for the first molecule.[PAR]",
    "The [TT]-of[tt] option produces an index file that can be used for",
    "freezing atoms. In this case, the input file must be a [TT].pdb[tt] file.[PAR]",
    "With the [TT]-disre[tt] option, half a matrix of distance restraints",
    "is generated instead of position restraints. With this matrix, that",
    "one typically would apply to C[GRK]alpha[grk] atoms in a protein, one can",
    "maintain the overall conformation of a protein without tieing it to",
    "a specific position (as with position restraints)."
  };
  static rvec    fc={1000.0,1000.0,1000.0};
  static real    freeze_level = 0.0;
  static real    disre_dist = 0.1;
  static real    disre_frac = 0.0;
  static real    disre_up2  = 1.0;
  static gmx_bool    bDisre=FALSE;
  static gmx_bool    bConstr=FALSE;
  static real    cutoff = -1.0;
	
  t_pargs pa[] = {
    { "-fc", FALSE, etRVEC, {fc}, 
      "force constants (kJ/mol nm^2)" },
    { "-freeze", FALSE, etREAL, {&freeze_level},
      "if the [TT]-of[tt] option or this one is given an index file will be written containing atom numbers of all atoms that have a B-factor less than the level given here" },
    { "-disre", FALSE, etBOOL, {&bDisre},
      "Generate a distance restraint matrix for all the atoms in index" },
    { "-disre_dist", FALSE, etREAL, {&disre_dist},
      "Distance range around the actual distance for generating distance restraints" },
    { "-disre_frac", FALSE, etREAL, {&disre_frac},
      "Fraction of distance to be used as interval rather than a fixed distance. If the fraction of the distance that you specify here is less than the distance given in the previous option, that one is used instead." },
    { "-disre_up2", FALSE, etREAL, {&disre_up2},
      "Distance between upper bound for distance restraints, and the distance at which the force becomes constant (see manual)" },
    { "-cutoff", FALSE, etREAL, {&cutoff},
      "Only generate distance restraints for atoms pairs within cutoff (nm)" },
    { "-constr", FALSE, etBOOL, {&bConstr},
      "Generate a constraint matrix rather than distance restraints. Constraints of type 2 will be generated that do generate exclusions." }
  };
#define npargs asize(pa)

  output_env_t oenv;
  t_atoms      *atoms=NULL;
  int          i,j,k;
  FILE         *out;
  int          igrp;
  real         d,dd,lo,hi;
  atom_id      *ind_grp;
  const char   *xfn,*nfn;
  char         *gn_grp;
  char         title[STRLEN];
  matrix       box;
  gmx_bool         bFreeze;
  rvec         dx,*x=NULL,*v=NULL;
  
  t_filenm fnm[] = {
    { efSTX, "-f",  NULL,    ffREAD },
    { efNDX, "-n",  NULL,    ffOPTRD },
    { efITP, "-o",  "posre", ffWRITE },
    { efNDX, "-of", "freeze",    ffOPTWR }
  };
#define NFILE asize(fnm)
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,0,NFILE,fnm,npargs,pa,
		    asize(desc),desc,0,NULL,&oenv);
  
  bFreeze = opt2bSet("-of",NFILE,fnm) || opt2parg_bSet("-freeze",asize(pa),pa);
  bDisre  = bDisre || opt2parg_bSet("-disre_dist",npargs,pa);
  xfn     = opt2fn_null("-f",NFILE,fnm);
  nfn     = opt2fn_null("-n",NFILE,fnm);
  
  if (( nfn == NULL ) && ( xfn == NULL))
    gmx_fatal(FARGS,"no index file and no structure file suplied");
      
  if ((disre_frac < 0) || (disre_frac >= 1))
    gmx_fatal(FARGS,"disre_frac should be between 0 and 1");
  if (disre_dist < 0)
    gmx_fatal(FARGS,"disre_dist should be >= 0");
    
  if (xfn != NULL) {
    snew(atoms,1);
    get_stx_coordnum(xfn,&(atoms->nr));
    init_t_atoms(atoms,atoms->nr,TRUE);
    snew(x,atoms->nr);
    snew(v,atoms->nr);
    fprintf(stderr,"\nReading structure file\n");
    read_stx_conf(xfn,title,atoms,x,v,NULL,box);
  }
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  real        **val=NULL,*t=NULL,dt_remd;
  int         n,nset;
  t_filenm fnm[] = {
      { efTRX, "-f",     NULL,      ffRDMULT },
      { efTRO, "-o",     NULL,      ffWRMULT },
      { efNDX, "-n",     "index",   ffOPTRD  },
      { efXVG, "-demux", "remd",    ffOPTRD  }
  };
  
#define NFILE asize(fnm)
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_BE_NICE|PCA_TIME_UNIT,
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,
		    0,NULL);

  bIndex = ftp2bSet(efNDX,NFILE,fnm);
  bDeMux = ftp2bSet(efXVG,NFILE,fnm);
  bSort  = bSort && !bDeMux;
  
  imax=NO_ATID;
  if (bIndex) {
    printf("Select group for output\n");
    rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
    /* scan index */
    imax=index[0];
    for(i=1; i<isize; i++)
      imax = max(imax, index[i]);
  }
  if (bDeMux) {
    nset    = 0;
    dt_remd = 0;
    val=read_xvg_time(opt2fn("-demux",NFILE,fnm),TRUE,
		      opt2parg_bSet("-b",npargs,pa),begin,
		      opt2parg_bSet("-e",npargs,pa),end,
		      1,&nset,&n,&dt_remd,&t);
    printf("Read %d sets of %d points, dt = %g\n\n",nset,n,dt_remd);
    if (debug) {
      fprintf(debug,"Dump of replica_index.xvg\n");
      for(i=0; (i<n); i++) {
	fprintf(debug,"%10g",t[i]);
	for(j=0; (j<nset); j++) {
	  fprintf(debug,"  %3d",gmx_nint(val[j][i]));
	}
	fprintf(debug,"\n");
      }
    }
  }
  /* prec is in nr of decimal places, xtcprec is a multiplication factor: */
  xtcpr=1;
  for (i=0; i<prec; i++)
    xtcpr*=10;
  
  nfile_in = opt2fns(&fnms,"-f",NFILE,fnm);
  if (!nfile_in)
    gmx_fatal(FARGS,"No input files!");
    
  if (bDeMux && (nfile_in != nset)) 
    gmx_fatal(FARGS,"You have specified %d files and %d entries in the demux table",nfile_in,nset);
    
  nfile_out = opt2fns(&fnms_out,"-o",NFILE,fnm);
  if (!nfile_out)
    gmx_fatal(FARGS,"No output files!");
  if ((nfile_out > 1) && !bDeMux) 
    gmx_fatal(FARGS,"Don't know what to do with more than 1 output file if  not demultiplexing");
  else if (bDeMux && (nfile_out != nset) && (nfile_out != 1))

//.........这里部分代码省略.........
  
  /* Parse command line on all processors, arguments are passed on in 
   * init_par (see above)
   */
  parse_common_args(&argc,argv,
		    PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_BE_NICE |
		    PCA_CAN_SET_DEFFNM | (MASTER(cr) ? 0 : PCA_QUIET),
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
  
#ifndef GMX_MPI
  if (nnodes > 1) 
    gmx_fatal(FARGS,"GROMACS compiled without MPI support - can't do parallel runs");
#endif

  /* Open log files on all processors */
  open_log(ftp2fn(efLOG,NFILE,fnm),cr);
  snew(ir,1);
  
  if (MASTER(cr)) {
    /* Read tpr file etc. */
    read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&tpx,FALSE,NULL,NULL);
    snew(x,tpx.natoms);
    read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,ir,
	     box,&natoms,x,NULL,NULL,&top);
    /* Charges */
    qtot = 0;
    snew(charge,natoms);
    for(i=0; (i<natoms); i++) {
      charge[i] = top.atoms.atom[i].q;
      qtot += charge[i];
    }
  
    /* Grid stuff */
    if (opt2parg_bSet("-grid",asize(pa),pa)) {
      ir->nkx = grid[XX];
      ir->nky = grid[YY];
      ir->nkz = grid[ZZ];
    }
    /* Check command line parameters for consistency */
    if ((ir->nkx <= 0) || (ir->nky <= 0) || (ir->nkz <= 0))
      gmx_fatal(FARGS,"PME grid = %d %d %d",ir->nkx,ir->nky,ir->nkz);
    if (opt2parg_bSet("-rc",asize(pa),pa)) 
      ir->rcoulomb = rc;
    if (ir->rcoulomb <= 0)
      gmx_fatal(FARGS,"rcoulomb should be > 0 (not %f)",ir->rcoulomb);
    if (opt2parg_bSet("-order",asize(pa),pa)) 
      ir->pme_order = pme_order;
    if (ir->pme_order <= 0)
      gmx_fatal(FARGS,"pme_order should be > 0 (not %d)",ir->pme_order);
    if (opt2parg_bSet("-tol",asize(pa),pa))
      ir->ewald_rtol = dtol;
    if (ir->ewald_rtol <= 0)
      gmx_fatal(FARGS,"ewald_tol should be > 0 (not %f)",ir->ewald_rtol);
  }
  else {
    init_top(&top);
  }

  /* Add parallellization code here */
  snew(nsb,1);
  if (MASTER(cr)) {
    ncg = top.blocks[ebCGS].multinr[0];
    for(i=0; (i<cr->nnodes-1); i++)
      top.blocks[ebCGS].multinr[i] = min(ncg,(ncg*(i+1))/cr->nnodes);
    for( ; (i<MAXNODES); i++)
      top.blocks[ebCGS].multinr[i] = ncg;

int cmain (int argc, char *argv[])
{
    const char       *desc[] = {
        "tpbconv can edit run input files in four ways.[PAR]",
        "[BB]1.[bb] by modifying the number of steps in a run input file",
        "with options [TT]-extend[tt], [TT]-until[tt] or [TT]-nsteps[tt]",
        "(nsteps=-1 means unlimited number of steps)[PAR]",
        "[BB]2.[bb] (OBSOLETE) by creating a run input file",
        "for a continuation run when your simulation has crashed due to e.g.",
        "a full disk, or by making a continuation run input file.",
        "This option is obsolete, since mdrun now writes and reads",
        "checkpoint files.",
        "[BB]Note[bb] that a frame with coordinates and velocities is needed.",
        "When pressure and/or Nose-Hoover temperature coupling is used",
        "an energy file can be supplied to get an exact continuation",
        "of the original run.[PAR]",
        "[BB]3.[bb] by creating a [TT].tpx[tt] file for a subset of your original",
        "tpx file, which is useful when you want to remove the solvent from",
        "your [TT].tpx[tt] file, or when you want to make e.g. a pure C[GRK]alpha[grk] [TT].tpx[tt] file.",
        "Note that you may need to use [TT]-nsteps -1[tt] (or similar) to get",
        "this to work.",
        "[BB]WARNING: this [TT].tpx[tt] file is not fully functional[bb].[PAR]",
        "[BB]4.[bb] by setting the charges of a specified group",
        "to zero. This is useful when doing free energy estimates",
        "using the LIE (Linear Interaction Energy) method."
    };

    const char       *top_fn, *frame_fn;
    t_fileio         *fp;
    ener_file_t       fp_ener = NULL;
    t_trnheader       head;
    int               i;
    gmx_large_int_t   nsteps_req, run_step, frame;
    double            run_t, state_t;
    gmx_bool          bOK, bNsteps, bExtend, bUntil, bTime, bTraj;
    gmx_bool          bFrame, bUse, bSel, bNeedEner, bReadEner, bScanEner, bFepState;
    gmx_mtop_t        mtop;
    t_atoms           atoms;
    t_inputrec       *ir, *irnew = NULL;
    t_gromppopts     *gopts;
    t_state           state;
    rvec             *newx = NULL, *newv = NULL, *tmpx, *tmpv;
    matrix            newbox;
    int               gnx;
    char             *grpname;
    atom_id          *index = NULL;
    int               nre;
    gmx_enxnm_t      *enm     = NULL;
    t_enxframe       *fr_ener = NULL;
    char              buf[200], buf2[200];
    output_env_t      oenv;
    t_filenm          fnm[] = {
        { efTPX, NULL,  NULL,    ffREAD  },
        { efTRN, "-f",  NULL,    ffOPTRD },
        { efEDR, "-e",  NULL,    ffOPTRD },
        { efNDX, NULL,  NULL,    ffOPTRD },
        { efTPX, "-o",  "tpxout", ffWRITE }
    };
#define NFILE asize(fnm)

    /* Command line options */
    static int      nsteps_req_int = 0;
    static real     start_t        = -1.0, extend_t = 0.0, until_t = 0.0;
    static int      init_fep_state = 0;
    static gmx_bool bContinuation  = TRUE, bZeroQ = FALSE, bVel = TRUE;
    static t_pargs  pa[]           = {
        { "-extend",        FALSE, etREAL, {&extend_t},
          "Extend runtime by this amount (ps)" },
        { "-until",         FALSE, etREAL, {&until_t},
          "Extend runtime until this ending time (ps)" },
        { "-nsteps",        FALSE, etINT,  {&nsteps_req_int},
          "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;

    CopyRight(stderr, argv[0]);

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

    /* Convert int to gmx_large_int_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);

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

//.........这里部分代码省略.........
    { "-gradient",FALSE, etRVEC, {grad},
      "Re-scale colormap to a smooth gradient from white {1,1,1} to {r,g,b}" },
    { "-skip",    FALSE, etINT,  {&skip},
      "only write out every nr-th row and column" },
    { "-zeroline",FALSE, etBOOL, {&bZeroLine},
      "insert line in xpm matrix where axis label is zero"},
    { "-legoffset", FALSE, etINT, {&mapoffset},
      "Skip first N colors from xpm file for the legend" },
    { "-combine", FALSE, etENUM, {combine}, "Combine two matrices" },
    { "-cmin",    FALSE, etREAL, {&cmin}, "Minimum for combination output" },
    { "-cmax",    FALSE, etREAL, {&cmax}, "Maximum for combination output" }
  };
#define NPA asize(pa)
  t_filenm  fnm[] = {
    { efXPM, "-f",  NULL,      ffREAD },
    { efXPM, "-f2", "root2",   ffOPTRD },
    { efM2P, "-di", NULL,      ffLIBOPTRD },
    { efM2P, "-do", "out",     ffOPTWR },
    { efEPS, "-o",  NULL,      ffOPTWR },
    { efXPM, "-xpm",NULL,      ffOPTWR }
  };
#define NFILE asize(fnm)
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW,
		    NFILE,fnm,NPA,pa,
		    asize(desc),desc,0,NULL);

  etitle   = nenum(title);
  elegend  = nenum(legend);
  ediag    = nenum(diag);
  erainbow = nenum(rainbow);
  ecombine = nenum(combine);
  bGrad    = opt2parg_bSet("-gradient",NPA,pa);
  for(i=0; i<DIM; i++)
    if (grad[i] < 0 || grad[i] > 1)
      gmx_fatal(FARGS, "RGB value %g out of range (0.0-1.0)", grad[i]);
  if (!bFrame) {
    etitle = etNone;
    elegend = elNone;
  }

  epsfile=ftp2fn_null(efEPS,NFILE,fnm);
  xpmfile=opt2fn_null("-xpm",NFILE,fnm);
  if ( epsfile==NULL && xpmfile==NULL ) {
    if (ecombine!=ecHalves)
      xpmfile=opt2fn("-xpm",NFILE,fnm);
    else
      epsfile=ftp2fn(efEPS,NFILE,fnm);
  }
  if (ecombine!=ecHalves && epsfile) {
    fprintf(stderr,
	    "WARNING: can only write result of arithmetic combination "
	    "of two matrices to .xpm file\n"
	    "         file %s will not be written\n", epsfile);
    epsfile = NULL;
  }
  
  bDiag      = ediag!=edNone;
  bFirstDiag = ediag!=edSecond;
  
  fn=opt2fn("-f",NFILE,fnm);
  nmat=read_xpm_matrix(fn,&mat);
  fprintf(stderr,"There are %d matrices in %s\n",nmat,fn);
  fn=opt2fn_null("-f2",NFILE,fnm);
  if (fn) {

//.........这里部分代码省略.........
    { efSTO, "-3d",   "3dproj.pdb",  ffOPTWR },
    { efTRX, "-filt", "filtered",    ffOPTWR },
    { efTRX, "-extr", "extreme.pdb", ffOPTWR },
    { efXVG, "-over", "overlap",     ffOPTWR },
    { efXPM, "-inpr", "inprod",      ffOPTWR }
  }; 
#define NFILE asize(fnm) 

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

  indexfile=ftp2fn_null(efNDX,NFILE,fnm);

  VecFile         = opt2fn("-v",NFILE,fnm);
  Vec2File        = opt2fn_null("-v2",NFILE,fnm);
  topfile         = ftp2fn(efTPS,NFILE,fnm); 
  EigFile         = opt2fn_null("-eig",NFILE,fnm);
  Eig2File        = opt2fn_null("-eig2",NFILE,fnm);
  CompFile        = opt2fn_null("-comp",NFILE,fnm);
  RmsfFile        = opt2fn_null("-rmsf",NFILE,fnm);
  ProjOnVecFile   = opt2fn_null("-proj",NFILE,fnm);
  TwoDPlotFile    = opt2fn_null("-2d",NFILE,fnm);
  ThreeDPlotFile  = opt2fn_null("-3d",NFILE,fnm);
  FilterFile      = opt2fn_null("-filt",NFILE,fnm);
  ExtremeFile     = opt2fn_null("-extr",NFILE,fnm);
  OverlapFile     = opt2fn_null("-over",NFILE,fnm);
  InpMatFile      = ftp2fn_null(efXPM,NFILE,fnm);
  
  bTop   = fn2bTPX(topfile);
  bProj  = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile 
    || FilterFile || ExtremeFile;
  bFirstLastSet  = 
    opt2parg_bSet("-first",NPA,pa) && opt2parg_bSet("-last",NPA,pa);
  bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
    OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
  bVec2  = Vec2File || OverlapFile || InpMatFile;
  bM     = RmsfFile || bProj;
  bTraj  = ProjOnVecFile || FilterFile || (ExtremeFile && (max==0))
    || TwoDPlotFile || ThreeDPlotFile;
  bIndex = bM || bProj;
  bTPS   = ftp2bSet(efTPS,NFILE,fnm) || bM || bTraj ||
    FilterFile  || (bIndex && indexfile);
  bCompare = Vec2File || Eig2File;
  bPDB3D = fn2ftp(ThreeDPlotFile)==efPDB;
  
  read_eigenvectors(VecFile,&natoms,&bFit1,
		    &xref1,&bDMR1,&xav1,&bDMA1,
		    &nvec1,&eignr1,&eigvec1,&eigval1);
  neig1 = DIM*natoms;
  
  /* Overwrite eigenvalues from separate files if the user provides them */
  if (EigFile != NULL) {
    int neig_tmp = read_xvg(EigFile,&xvgdata,&i);
    if (neig_tmp != neig1)
      fprintf(stderr,"Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",neig1,natoms);
    neig1 = neig_tmp;
    srenew(eigval1,neig1);
    for(j=0;j<neig1;j++) {
      real tmp = eigval1[j];
      eigval1[j]=xvgdata[1][j];
      if (debug && (eigval1[j] != tmp))
	fprintf(debug,"Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
		j,tmp,eigval1[j]);
    }
    for(j=0;j<i;j++)

//.........这里部分代码省略.........
        { efXVG, "-ekt", "ektrans",   ffOPTWR },
        { efXVG, "-ekr", "ekrot",     ffOPTWR },
        { efXVG, "-vd",  "veldist",   ffOPTWR },
        { efPDB, "-cv",  "veloc",     ffOPTWR },
        { efPDB, "-cf",  "force",     ffOPTWR },
        { efXVG, "-av",  "all_veloc", ffOPTWR },
        { efXVG, "-af",  "all_force", ffOPTWR }
    };
#define NFILE asize(fnm)

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

    if (bMol)
    {
        fprintf(stderr, "Interpreting indexfile entries as molecules.\n"
                "Using center of mass.\n");
    }

    bOX  = opt2bSet("-ox", NFILE, fnm);
    bOXT = opt2bSet("-oxt", NFILE, fnm);
    bOV  = opt2bSet("-ov", NFILE, fnm);
    bOF  = opt2bSet("-of", NFILE, fnm);
    bOB  = opt2bSet("-ob", NFILE, fnm);
    bOT  = opt2bSet("-ot", NFILE, fnm);
    bEKT = opt2bSet("-ekt", NFILE, fnm);
    bEKR = opt2bSet("-ekr", NFILE, fnm);
    bCV  = opt2bSet("-cv", NFILE, fnm) || opt2bSet("-av", NFILE, fnm);
    bCF  = opt2bSet("-cf", NFILE, fnm) || opt2bSet("-af", NFILE, fnm);
    bVD  = opt2bSet("-vd", NFILE, fnm) || opt2parg_bSet("-bin", asize(pa), pa);
    if (bMol || bOT || bEKT || bEKR)
    {
        bCom = TRUE;
    }

    bDim[XX]  = bX;
    bDim[YY]  = bY;
    bDim[ZZ]  = bZ;
    bDim[DIM] = bNorm;

    if (bFP)
    {
        sprintf(sffmt, "\t%s", gmx_real_fullprecision_pfmt);
    }
    else
    {
        sprintf(sffmt, "\t%%g");
    }
    sprintf(sffmt6, "%s%s%s%s%s%s", sffmt, sffmt, sffmt, sffmt, sffmt, sffmt);

    bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC,
                         &xtop, NULL, topbox,
                         bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
    sfree(xtop);
    if ((bMol || bCV || bCF) && !bTop)
    {
        gmx_fatal(FARGS, "Need a run input file for option -mol, -cv or -cf");
    }

    if (bMol)
    {
        indexfn = ftp2fn(efNDX, NFILE, fnm);

//.........这里部分代码省略.........
        { efDAT, "-groups", "groups", ffREAD },
        { efDAT, "-eref",   "eref",   ffOPTRD },
        { efXPM, "-emat",   "emat",   ffWRITE },
        { efXVG, "-etot",   "energy", ffWRITE }
    };
#define NFILE asize(fnm)

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

    for (i = 0; (i < egNR+egSP); i++)
    {
        egrp_use[i] = FALSE;
    }
    egrp_use[egCOULSR] = bCoulSR;
    egrp_use[egLJSR]   = bLJSR;
    egrp_use[egBHAMSR] = bBhamSR;
    egrp_use[egCOUL14] = bCoul14;
    egrp_use[egLJ14]   = bLJ14;
    egrp_use[egTotal]  = TRUE;

    bRef = opt2bSet("-eref", NFILE, fnm);
    in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
    do_enxnms(in, &nre, &enm);

    if (nre == 0)
    {
        gmx_fatal(FARGS, "No energies!\n");
    }

    bCutmax = opt2parg_bSet("-max", asize(pa), pa);
    bCutmin = opt2parg_bSet("-min", asize(pa), pa);

    nenergy = 0;

    /* Read groupnames from input file and construct selection of
       energy groups from it*/

    fprintf(stderr, "Will read groupnames from inputfile\n");
    ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
    fprintf(stderr, "Read %d groups\n", ngroups);
    snew(set, static_cast<size_t>(gmx::square(ngroups)*egNR/2));
    n     = 0;
    prevk = 0;
    for (i = 0; (i < ngroups); i++)
    {
        fprintf(stderr, "\rgroup %d", i);
        for (j = i; (j < ngroups); j++)
        {
            for (m = 0; (m < egNR); m++)
            {
                if (egrp_use[m])
                {
                    sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
#ifdef DEBUG
                    fprintf(stderr, "\r%-15s %5d", groupname, n);
#endif
                    for (k = prevk; (k < prevk+nre); k++)
                    {
                        if (std::strcmp(enm[k%nre].name, groupname) == 0)
                        {
                            set[n++] = k;
                            break;

//.........这里部分代码省略.........
    /*to read topology file*/
    t_topology  top;
    int         ePBC;
    char        title[STRLEN];
    matrix      topbox;
    rvec       *xtop;
    gmx_bool    bTop, bFit1;

    t_filenm    fnm[] = {
        { efTRN, "-f",    "eigenvec",    ffREAD  },
        { efXVG, "-eig",  "eigenval",    ffOPTRD },
        { efTPS, NULL,    NULL,          ffREAD },
        { efNDX, NULL,    NULL,  ffOPTRD },
        { efSTX, "-tar", "target", ffOPTRD},
        { efSTX, "-ori", "origin", ffOPTRD},
        { efEDI, "-o", "sam", ffWRITE }
    };
#define NFILE asize(fnm)
    edi_params.outfrq = 100; edi_params.slope = 0.0; edi_params.maxedsteps = 0;
    if (!parse_common_args(&argc, argv, 0,
                           NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }

    indexfile       = ftp2fn_null(efNDX, NFILE, fnm);
    EdiFile         = ftp2fn(efEDI, NFILE, fnm);
    TargetFile      = opt2fn_null("-tar", NFILE, fnm);
    OriginFile      = opt2fn_null("-ori", NFILE, fnm);


    for (ev_class = 0; ev_class < evNr; ++ev_class)
    {
        if (opt2parg_bSet(evOptions[ev_class], NPA, pa))
        {
            /*get list of eigenvectors*/
            nvecs = sscan_list(&(listen[ev_class]), opt2parg_str(evOptions[ev_class], NPA, pa), evOptions[ev_class]);
            if (ev_class < evStepNr-2)
            {
                /*if apropriate get list of stepsizes for these eigenvectors*/
                if (opt2parg_bSet(evStepOptions[ev_class], NPA, pa))
                {
                    evStepList[ev_class] =
                        scan_vecparams(opt2parg_str(evStepOptions[ev_class], NPA, pa), evStepOptions[ev_class], nvecs);
                }
                else   /*if list is not given fill with zeros */
                {
                    snew(evStepList[ev_class], nvecs);
                    for (i = 0; i < nvecs; i++)
                    {
                        evStepList[ev_class][i] = 0.0;
                    }
                }
            }
            else if (ev_class == evRADFIX)
            {
                snew(evStepList[ev_class], nvecs);
                for (i = 0; i < nvecs; i++)
                {
                    evStepList[ev_class][i] = radstep;
                }
            }
            else if (ev_class == evFLOOD)
            {
                snew(evStepList[ev_class], nvecs);

int gmx_make_ndx(int argc, char *argv[])
{
    const char     *desc[] = {
        "Index groups are necessary for almost every GROMACS program.",
        "All these programs can generate default index groups. You ONLY",
        "have to use [THISMODULE] when you need SPECIAL index groups.",
        "There is a default index group for the whole system, 9 default",
        "index groups for proteins, and a default index group",
        "is generated for every other residue name.[PAR]",
        "When no index file is supplied, also [THISMODULE] will generate the",
        "default groups.",
        "With the index editor you can select on atom, residue and chain names",
        "and numbers.",
        "When a run input file is supplied you can also select on atom type.",
        "You can use NOT, AND and OR, you can split groups",
        "into chains, residues or atoms. You can delete and rename groups.[PAR]",
        "The atom numbering in the editor and the index file starts at 1.[PAR]",
        "The [TT]-twin[tt] switch duplicates all index groups with an offset of",
        "[TT]-natoms[tt], which is useful for Computational Electrophysiology",
        "double-layer membrane setups."
    };

    static int      natoms     = 0;
    static gmx_bool bVerbose   = FALSE;
    static gmx_bool bDuplicate = FALSE;
    t_pargs         pa[]       = {
        { "-natoms",  FALSE, etINT, {&natoms},
          "set number of atoms (default: read from coordinate or index file)" },
        { "-twin",     FALSE, etBOOL, {&bDuplicate},
          "Duplicate all index groups with an offset of -natoms" },
        { "-verbose", FALSE, etBOOL, {&bVerbose},
          "HIDDENVerbose output" }
    };
#define NPA asize(pa)

    output_env_t oenv;
    char         title[STRLEN];
    int          nndxin;
    const char  *stxfile;
    char       **ndxinfiles;
    const char  *ndxoutfile;
    gmx_bool     bNatoms;
    int          i, j;
    t_atoms     *atoms;
    rvec        *x, *v;
    int          ePBC;
    matrix       box;
    t_blocka    *block, *block2;
    char       **gnames, **gnames2;
    t_filenm     fnm[] = {
        { efSTX, "-f", NULL,     ffOPTRD  },
        { efNDX, "-n", NULL,     ffOPTRDMULT },
        { efNDX, "-o", NULL,     ffWRITE }
    };
#define NFILE asize(fnm)

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

    stxfile = ftp2fn_null(efSTX, NFILE, fnm);
    if (opt2bSet("-n", NFILE, fnm))
    {
        nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);
    }
    else
    {
        nndxin = 0;
    }
    ndxoutfile = opt2fn("-o", NFILE, fnm);
    bNatoms    = opt2parg_bSet("-natoms", NPA, pa);

    if (!stxfile && !nndxin)
    {
        gmx_fatal(FARGS, "No input files (structure or index)");
    }

    if (stxfile)
    {
        snew(atoms, 1);
        get_stx_coordnum(stxfile, &(atoms->nr));
        init_t_atoms(atoms, atoms->nr, TRUE);
        snew(x, atoms->nr);
        snew(v, atoms->nr);
        fprintf(stderr, "\nReading structure file\n");
        read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);
        natoms  = atoms->nr;
        bNatoms = TRUE;
    }
    else
    {
        atoms = NULL;
        x     = NULL;
    }

    /* read input file(s) */
    block  = new_blocka();
    gnames = NULL;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        { efNDX, "-n", "index", ffOPTRD },
        { efXVG, "-demux", "remd", ffOPTRD }
    };

#define NFILE asize(fnm)

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

    bIndex = ftp2bSet(efNDX, NFILE, fnm);
    bDeMux = ftp2bSet(efXVG, NFILE, fnm);
    bSort  = bSort && !bDeMux;

    imax = -1;
    if (bIndex)
    {
        printf("Select group for output\n");
        rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
        /* scan index */
        imax = index[0];
        for (i = 1; i < isize; i++)
        {
            imax = std::max(imax, index[i]);
        }
    }
    if (bDeMux)
    {
        nset    = 0;
        dt_remd = 0;
        val     = read_xvg_time(opt2fn("-demux", NFILE, fnm), TRUE,
                                opt2parg_bSet("-b", npargs, pa), begin,
                                opt2parg_bSet("-e", npargs, pa), end, 1, &nset, &n,
                                &dt_remd, &t);
        printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt_remd);
        if (debug)
        {
            fprintf(debug, "Dump of replica_index.xvg\n");
            for (i = 0; (i < n); i++)
            {
                fprintf(debug, "%10g", t[i]);
                for (j = 0; (j < nset); j++)
                {
                    fprintf(debug, "  %3d", static_cast<int>(std::round(val[j][i])));
                }
                fprintf(debug, "\n");
            }
        }
    }

    nfile_in = opt2fns(&fnms, "-f", NFILE, fnm);
    if (!nfile_in)
    {
        gmx_fatal(FARGS, "No input files!" );
    }

    if (bDeMux && (nfile_in != nset))
    {
        gmx_fatal(FARGS, "You have specified %d files and %d entries in the demux table", nfile_in, nset);
    }

    ftpin = fn2ftp(fnms[0]);

    for (i = 1; i < nfile_in; i++)

//.........这里部分代码省略.........
    real          **val, **et_val, *t, *e_t, e_dt, d_dt, dt, tot, error;
    real           *rmin, *rmax;
    double         *av, *sig, cum1, cum2, cum3, cum4, db;
    const char     *fn_ge, *fn_ene;
    output_env_t    oenv;
    gmx_int64_t     num_grid_points;

    t_filenm        fnm[] = {
        { efXVG, "-f",    "graph",    ffREAD   },
        { efXVG, "-ge",   "gibbs",    ffOPTRD  },
        { efXVG, "-ene",  "esham",    ffOPTRD  },
        { efXVG, "-dist", "ener",     ffOPTWR  },
        { efXVG, "-histo", "edist",    ffOPTWR  },
        { efNDX, "-bin",  "bindex",   ffOPTWR  },
        { efXPM, "-lp",   "prob",     ffOPTWR  },
        { efXPM, "-ls",   "gibbs",    ffOPTWR  },
        { efXPM, "-lsh",  "enthalpy", ffOPTWR  },
        { efXPM, "-lss",  "entropy",  ffOPTWR  },
        { efPDB, "-ls3",  "gibbs3",   ffOPTWR  },
        { efLOG, "-g",    "shamlog",  ffOPTWR  }
    };
#define NFILE asize(fnm)

    int     npargs;

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

    val = read_xvg_time(opt2fn("-f", NFILE, fnm), bHaveT,
                        opt2parg_bSet("-b", npargs, pa), tb-ttol,
                        opt2parg_bSet("-e", npargs, pa), te+ttol,
                        nsets_in, &nset, &n, &dt, &t);
    printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt);

    fn_ge  = opt2fn_null("-ge", NFILE, fnm);
    fn_ene = opt2fn_null("-ene", NFILE, fnm);

    if (fn_ge && fn_ene)
    {
        gmx_fatal(FARGS, "Can not do free energy and energy corrections at the same time");
    }

    if (fn_ge || fn_ene)
    {
        et_val = read_xvg_time(fn_ge ? fn_ge : fn_ene, bHaveT,
                               opt2parg_bSet("-b", npargs, pa), tb-ttol,
                               opt2parg_bSet("-e", npargs, pa), te+ttol,
                               1, &e_nset, &e_n, &e_dt, &e_t);
        if (fn_ge)
        {
            if (e_nset != 1)
            {
                gmx_fatal(FARGS, "Can only handle one free energy component in %s",
                          fn_ge);
            }
        }
        else
        {
            if (e_nset != 1 && e_nset != 2)
            {
                gmx_fatal(FARGS, "Can only handle one energy component or one energy and one T in %s",
                          fn_ene);

//.........这里部分代码省略.........
        { NULL, NULL, FALSE,  "JCaHa",   "J-Coupling Values",        "Residue", "Hz", 0.0 },
        { NULL, NULL, FALSE,  "helicity", "Helicity per Residue",     "Residue", "% of time", 0.0 }
    };

    output_env_t   oenv;
    char           buf[54];
    t_trxstatus   *status;
    int            natoms, nre, nres;
    t_bb          *bb;
    int            i, j, m, nall, nbb, nca, teller, nSel = 0;
    atom_id       *bbindex, *caindex, *allindex;
    t_topology    *top;
    int            ePBC;
    rvec          *x, *xref;
    real           t;
    real           rms;
    matrix         box;
    gmx_rmpbc_t    gpbc = NULL;
    gmx_bool       bRange;
    t_filenm       fnm[] = {
        { efTPX, NULL,  NULL,   ffREAD  },
        { efNDX, NULL,  NULL,   ffREAD  },
        { efTRX, "-f",  NULL,   ffREAD  },
        { efSTO, "-cz", "zconf", ffWRITE },
    };
#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, 0, NULL, &oenv))
    {
        return 0;
    }

    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 (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 = gmx_ffopen(buf, "w");

int gmx_dist(int argc,char *argv[])
{
  const char *desc[] = {
    "[TT]g_dist[tt] can calculate the distance between the centers of mass of two",
    "groups of atoms as a function of time. The total distance and its",
    "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-components are plotted.[PAR]",
    "Or when [TT]-dist[tt] is set, print all the atoms in group 2 that are",
    "closer than a certain distance to the center of mass of group 1.[PAR]",
    "With options [TT]-lt[tt] and [TT]-dist[tt] the number of contacts",
    "of all atoms in group 2 that are closer than a certain distance",
    "to the center of mass of group 1 are plotted as a function of the time",
    "that the contact was continuously present.[PAR]",
    "Other programs that calculate distances are [TT]g_mindist[tt]",
    "and [TT]g_bond[tt]."
  };
  
  t_topology *top=NULL;
  int  ePBC;
  real t,t0,cut2,dist2;
  rvec *x=NULL,*v=NULL,dx;
  matrix box;
  t_trxstatus *status;
  int natoms;

  int g,d,i,j,res,teller=0;
  atom_id aid;

  int     ngrps;     /* the number of index groups */
  atom_id **index,max;   /* the index for the atom numbers */
  int     *isize;    /* the size of each group */
  char    **grpname; /* the name of each group */
  rvec    *com;
  real    *mass;
  FILE    *fp=NULL,*fplt=NULL;
  gmx_bool    bCutoff,bPrintDist,bLifeTime;
  t_pbc   *pbc;
  int     *contact_time=NULL,*ccount=NULL,ccount_nalloc=0,sum;
  char    buf[STRLEN];
  output_env_t oenv;
  gmx_rmpbc_t  gpbc=NULL;
  
  const char *leg[4] = { "|d|","d\\sx\\N","d\\sy\\N","d\\sz\\N" };

  static real cut=0;
  
  static t_pargs pa[] = {
    { "-dist",      FALSE, etREAL, {&cut},
      "Print all atoms in group 2 closer than dist to the center of mass of group 1" }
  };
#define NPA asize(pa)

  t_filenm fnm[] = {
    { efTRX, "-f", NULL, ffREAD },
    { efTPX, NULL, NULL, ffREAD },
    { efNDX, NULL, NULL, ffOPTRD },
    { efXVG, NULL, "dist", ffOPTWR },
    { efXVG, "-lt", "lifetime", ffOPTWR },
  };
#define NFILE asize(fnm)


  CopyRight(stderr,argv[0]);

  parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
		    NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
  
  bCutoff = opt2parg_bSet("-dist",NPA,pa);
  cut2 = cut*cut;
  bLifeTime = opt2bSet("-lt",NFILE,fnm);
  bPrintDist = (bCutoff && !bLifeTime);
  
  top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
  
  /* read index files */
  ngrps = 2;
  snew(com,ngrps);
  snew(grpname,ngrps);
  snew(index,ngrps);
  snew(isize,ngrps);
  get_index(&top->atoms,ftp2fn(efNDX,NFILE,fnm),ngrps,isize,index,grpname);
  
  /* calculate mass */
  max=0;
  snew(mass,ngrps);
  for(g=0;(g<ngrps);g++) {
    mass[g]=0;
    for(i=0;(i<isize[g]);i++) {
      if (index[g][i]>max)
	max=index[g][i];
      if (index[g][i] >= top->atoms.nr)
	gmx_fatal(FARGS,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);
      mass[g]+=top->atoms.atom[index[g][i]].m;
    }
  }

  natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
  t0 = t;

  if (max>=natoms)
    gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);
//.........这里部分代码省略.........

int gmx_vanhove(int argc,char *argv[])
{
  const char *desc[] = {
    "g_vanhove computes the Van Hove correlation function.",
    "The Van Hove G(r,t) is the probability that a par