//.........这里部分代码省略.........
    "would eliminate all intramolecular contributions to the rdf.",
    "Note that all atoms in the selected groups are used, also the ones",
    "that don't have Lennard-Jones interactions.[PAR]",
    "Option [TT]-cn[tt] produces the cumulative number rdf,",
    "i.e. the average number of particles within a distance r.[PAR]",
    "To bridge the gap between theory and experiment structure factors can",
    "be computed (option [TT]-sq[tt]). The algorithm uses FFT, the grid"
    "spacing of which is determined by option [TT]-grid[tt]."
  };
  static bool bCM=FALSE,bXY=FALSE,bPBC=TRUE,bNormalize=TRUE;
  static real cutoff=0,binwidth=0.002,grid=0.05,fade=0.0,lambda=0.1,distance=10;
  static int  npixel=256,nlevel=20,ngroups=1;
  static real start_q=0.0, end_q=60.0, energy=12.0;

  static char *rdft[]={ NULL, "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL };

  t_pargs pa[] = {
    { "-bin",      FALSE, etREAL, {&binwidth},
      "Binwidth (nm)" },
    { "-com",      FALSE, etBOOL, {&bCM},
      "RDF with respect to the center of mass of first group" },
    { "-rdf",   FALSE, etENUM, {rdft}, 
      "RDF type" },
    { "-pbc",      FALSE, etBOOL, {&bPBC},
      "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the larges box edge." },
    { "-norm",     FALSE, etBOOL, {&bNormalize},
      "Normalize for volume and density" },
    { "-xy",       FALSE, etBOOL, {&bXY},
      "Use only the x and y components of the distance" },
    { "-cut",      FALSE, etREAL, {&cutoff},
      "Shortest distance (nm) to be considered"},
    { "-ng",       FALSE, etINT, {&ngroups},
      "Number of secondary groups to compute RDFs around a central group" },
    { "-fade",     FALSE, etREAL, {&fade},
      "From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." },
    { "-grid",     FALSE, etREAL, {&grid},
      "[HIDDEN]Grid spacing (in nm) for FFTs when computing structure factors" },
    { "-npixel",   FALSE, etINT,  {&npixel},
      "[HIDDEN]# pixels per edge of the square detector plate" },
    { "-nlevel",   FALSE, etINT,  {&nlevel},
      "Number of different colors in the diffraction image" },
    { "-distance", FALSE, etREAL, {&distance},
      "[HIDDEN]Distance (in cm) from the sample to the detector" },
    { "-wave",     FALSE, etREAL, {&lambda},
      "[HIDDEN]Wavelength for X-rays/Neutrons for scattering. 0.1 nm corresponds to roughly 12 keV" },
    
    {"-startq", FALSE, etREAL, {&start_q},
     "Starting q (1/nm) "},
    {"-endq", FALSE, etREAL, {&end_q},
     "Ending q (1/nm)"},
    {"-energy", FALSE, etREAL, {&energy},
     "Energy of the incoming X-ray (keV) "}
  };
#define NPA asize(pa)
  char       *fnTPS,*fnNDX;
  bool       bSQ,bRDF;
  
  t_filenm   fnm[] = {
    { efTRX, "-f",  NULL,     ffREAD },
    { efTPS, NULL,  NULL,     ffOPTRD },
    { efNDX, NULL,  NULL,     ffOPTRD },
    { efXVG, "-o",  "rdf",    ffOPTWR },
    { efXVG, "-sq", "sq",     ffOPTWR },
    { efXVG, "-cn", "rdf_cn", ffOPTWR },
    { efXVG, "-hq", "hq",     ffOPTWR },
/*    { efXPM, "-image", "sq",  ffOPTWR }*/
  };
#define NFILE asize(fnm)
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
		    NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);

  bSQ   = opt2bSet("-sq",NFILE,fnm);
  bRDF  = opt2bSet("-o",NFILE,fnm) || !bSQ;
  if (bSQ || bCM || rdft[0][0]=='m' || rdft[0][6]=='m') {
    fnTPS = ftp2fn(efTPS,NFILE,fnm);
  } else {
    fnTPS = ftp2fn_null(efTPS,NFILE,fnm);
  }
  fnNDX = ftp2fn_null(efNDX,NFILE,fnm);

  if (!bSQ && (!fnTPS && !fnNDX))
    gmx_fatal(FARGS,"Neither index file nor topology file specified\n"
	      "Nothing to do!");
 
  if  (bSQ) 
   do_scattering_intensity(fnTPS,fnNDX,opt2fn("-sq",NFILE,fnm),ftp2fn(efTRX,NFILE,fnm),
		           start_q, end_q, energy, ngroups  );

  if (bRDF) 
    do_rdf(fnNDX,fnTPS,ftp2fn(efTRX,NFILE,fnm),
	   opt2fn("-o",NFILE,fnm),opt2fn_null("-cn",NFILE,fnm),
	   opt2fn_null("-hq",NFILE,fnm),
	   bCM,rdft,bXY,bPBC,bNormalize,cutoff,binwidth,fade,ngroups);

  thanx(stderr);
  
  return 0;
}

int gmx_sgangle(int argc,char *argv[])
{
  static char *desc[] = {
    "Compute the angle and distance between two groups. ",
    "The groups are defined by a number of atoms given in an index file and",
    "may be two or three atoms in size.",
    "If -one is set, only one group should be specified in the index",
    "file and the angle between this group at time 0 and t will be computed.",
    "The angles calculated depend on the order in which the atoms are ",
    "given. Giving for instance 5 6 will rotate the vector 5-6 with ",
    "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
    "the normal on the plane spanned by those three atoms will be",
    "calculated, using the formula  P1P2 x P1P3.",
    "The cos of the angle is calculated, using the inproduct of the two",
    "normalized vectors.[PAR]",
    "Here is what some of the file options do:[BR]",
    "-oa: Angle between the two groups specified in the index file. If a group contains three atoms the normal to the plane defined by those three atoms will be used. If a group contains two atoms, the vector defined by those two atoms will be used.[BR]",
    "-od: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
    "-od1: If one plane and one vector is given, the distances for each of the atoms from the center of the plane is given seperately.[BR]",
    "-od2: For two planes this option has no meaning."
  };

  char	    *fna, *fnd, *fnd1, *fnd2;		
  char      *grpname[2];          		/* name of the two groups */
  int       gnx[2];               		/* size of the two groups */
  t_topology *top;                		/* topology 		*/ 
  int       ePBC;
  atom_id   *index[2];            		
  static bool bOne = FALSE, bZ=FALSE;
  t_pargs pa[] = {
    { "-one", FALSE, etBOOL, {&bOne},
      "Only one group compute angle between vector at time zero and time t"},
    { "-z", FALSE, etBOOL, {&bZ},
        "Use the Z-axis as reference" }
  };
#define NPA asize(pa)

  t_filenm  fnm[] = {             		/* files for g_sgangle 	*/
    { efTRX, "-f", NULL,  ffREAD },    		/* trajectory file 	*/
    { efNDX, NULL, NULL,  ffREAD },    		/* index file 		*/
    { efTPX, NULL, NULL,  ffREAD },    		/* topology file 	*/
    { efXVG,"-oa","sg_angle",ffWRITE },		/* xvgr output file 	*/
    { efXVG, "-od","sg_dist",ffOPTWR }, 	/* xvgr output file 	*/
    { efXVG, "-od1", "sg_dist1",ffOPTWR }, 	/* xvgr output file 	*/
    { efXVG, "-od2", "sg_dist2",ffOPTWR } 	/* xvgr output file 	*/
  };

#define NFILE asize(fnm)

  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
		    NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
  

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

  fna = opt2fn("-oa",NFILE,fnm);
  fnd = opt2fn_null("-od",NFILE,fnm);
  fnd1 = opt2fn_null("-od1",NFILE,fnm);
  fnd2 = opt2fn_null("-od2",NFILE,fnm);
  
  /* read index file. */
  if(bOne) {
    rd_index(ftp2fn(efNDX,NFILE,fnm),1,gnx,index,grpname); 
    print_types(index[0],gnx[0],grpname[0],    
 		index[0],gnx[0],grpname[0],top); 

    sgangle_plot_single(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
			index[0],gnx[0],grpname[0],
			index[0],gnx[0],grpname[0],
			top,ePBC);
  }  else {
    rd_index(ftp2fn(efNDX,NFILE,fnm),bZ ? 1 : 2,gnx,index,grpname);
    if (!bZ)
      print_types(index[0],gnx[0],grpname[0],
		index[1],gnx[1],grpname[1],top); 
    else {
      gnx[1] = 0;
      grpname[1] = "Z-axis";
    }  
    sgangle_plot(ftp2fn(efTRX,NFILE,fnm), fna, fnd, fnd1, fnd2,
		 index[0],gnx[0],grpname[0],
		 index[1],gnx[1],grpname[1],
		 top,ePBC);
  }

  do_view(fna,"-nxy");     /* view xvgr file */
  do_view(fnd,"-nxy");     /* view xvgr file */
  do_view(fnd1,"-nxy");     /* view xvgr file */
  do_view(fnd2,"-nxy");     /* view xvgr file */

  thanx(stderr);
  return 0;
}

//.........这里部分代码省略.........
               "See the GROMACS manual for a description of the requirements that\n"
               "must be satisfied by descriptions of simulation cells.\n",
               check_box(epbcXYZ, box));
    }

    if (bDist && btype[0][0] == 't')
    {
        if (TRICLINIC(box))
        {
            printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
                   "distance from the solute to a box surface along the corresponding normal\n"
                   "vector might be somewhat smaller than your specified value %f.\n"
                   "You can check the actual value with g_mindist -pi\n", dist);
        }
        else if (!opt2parg_bSet("-bt", NPA, pa))
        {
            printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
                   "If the molecule rotates the actual distance will be smaller. You might want\n"
                   "to use a cubic box instead, or why not try a dodecahedron today?\n");
        }
    }
    if (bCONECT && (outftp == efPDB) && (inftp == efTPR))
    {
        conect = gmx_conect_generate(top);
    }
    else
    {
        conect = NULL;
    }

    if (bIndex)
    {
        fprintf(stderr, "\nSelect a group for output:\n");
        get_index(&atoms, opt2fn_null("-n", NFILE, fnm),
                  1, &isize, &index, &grpname);

        if (resnr_start >= 0)
        {
            renum_resnr(&atoms, isize, index, resnr_start);
        }

        if (opt2parg_bSet("-label", NPA, pa))
        {
            for (i = 0; (i < atoms.nr); i++)
            {
                atoms.resinfo[atoms.atom[i].resind].chainid = label[0];
            }
        }

        if (opt2bSet("-bf", NFILE, fnm) || bLegend)
        {
            gmx_fatal(FARGS, "Sorry, cannot do bfactors with an index group.");
        }

        if (outftp == efPDB)
        {
            out = ffopen(outfile, "w");
            write_pdbfile_indexed(out, title, &atoms, x, ePBC, box, ' ', 1, isize, index, conect, TRUE);
            ffclose(out);
        }
        else
        {
            write_sto_conf_indexed(outfile, title, &atoms, x, bHaveV ? v : NULL, ePBC, box, isize, index);
        }
    }
    else

int gmx_order(int argc, char *argv[])
{
    const char        *desc[] = {
        "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
        "vector i-1, i+1 is used together with an axis. ",
        "The index file should contain only the groups to be used for calculations,",
        "with each group of equivalent carbons along the relevant acyl chain in its own",
        "group. There should not be any generic groups (like System, Protein) in the index",
        "file to avoid confusing the program (this is not relevant to tetrahedral order",
        "parameters however, which only work for water anyway).[PAR]",
        "[THISMODULE] can also give all",
        "diagonal elements of the order tensor and even calculate the deuterium",
        "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
        "order tensor component (specified by the [TT]-d[tt] option) is given and the",
        "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
        "selected, all diagonal elements and the deuterium order parameter is",
        "given.[PAR]"
        "The tetrahedrality order parameters can be determined",
        "around an atom. Both angle an distance order parameters are calculated. See",
        "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
        "for more details."
    };

    static int         nslices       = 1;     /* nr of slices defined       */
    static gmx_bool    bSzonly       = FALSE; /* True if only Sz is wanted  */
    static gmx_bool    bUnsat        = FALSE; /* True if carbons are unsat. */
    static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
    static gmx_bool    permolecule   = FALSE; /*compute on a per-molecule basis */
    static gmx_bool    radial        = FALSE; /*compute a radial membrane normal */
    static gmx_bool    distcalc      = FALSE; /*calculate distance from a reference group */
    t_pargs            pa[]          = {
        { "-d",      FALSE, etENUM, {normal_axis},
          "Direction of the normal on the membrane" },
        { "-sl",     FALSE, etINT, {&nslices},
          "Calculate order parameter as function of box length, dividing the box"
          " into this number of slices." },
        { "-szonly", FALSE, etBOOL, {&bSzonly},
          "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
        { "-unsat",  FALSE, etBOOL, {&bUnsat},
          "Calculate order parameters for unsaturated carbons. Note that this can"
          "not be mixed with normal order parameters." },
        { "-permolecule", FALSE, etBOOL, {&permolecule},
          "Compute per-molecule Scd order parameters" },
        { "-radial", FALSE, etBOOL, {&radial},
          "Compute a radial membrane normal" },
        { "-calcdist", FALSE, etBOOL, {&distcalc},
          "Compute distance from a reference" },
    };

    rvec              *order;                         /* order par. for each atom   */
    real             **slOrder;                       /* same, per slice            */
    real               slWidth = 0.0;                 /* width of a slice           */
    char             **grpname;                       /* groupnames                 */
    int                ngrps,                         /* nr. of groups              */
                       i,
                       axis = 0;                      /* normal axis                */
    t_topology       *top;                            /* topology         */
    int               ePBC;
    atom_id          *index,                          /* indices for a              */
    *a;                                               /* atom numbers in each group */
    t_blocka         *block;                          /* data from index file       */
    t_filenm          fnm[] = {                       /* files for g_order    */
        { efTRX, "-f", NULL,  ffREAD },               /* trajectory file              */
        { efNDX, "-n", NULL,  ffREAD },               /* index file           */
        { efNDX, "-nr", NULL,  ffREAD },              /* index for radial axis calculation	  */
        { efTPR, NULL, NULL,  ffREAD },               /* topology file                */
        { efXVG, "-o", "order", ffWRITE },            /* xvgr output file     */
        { efXVG, "-od", "deuter", ffWRITE },          /* xvgr output file           */
        { efPDB, "-ob", NULL, ffWRITE },              /* write Scd as B factors to PDB if permolecule           */
        { efXVG, "-os", "sliced", ffWRITE },          /* xvgr output file           */
        { efXVG, "-Sg", "sg-ang", ffOPTWR },          /* xvgr output file           */
        { efXVG, "-Sk", "sk-dist", ffOPTWR },         /* xvgr output file           */
        { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR },  /* xvgr output file           */
        { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file           */
    };
    gmx_bool          bSliced = FALSE;                /* True if box is sliced      */
#define NFILE asize(fnm)
    real            **distvals = NULL;
    const char       *sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
    gmx_output_env_t *oenv;

    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;
    }
    if (nslices < 1)
    {
        gmx_fatal(FARGS, "Can not have nslices < 1");
    }
    sgfnm  = opt2fn_null("-Sg", NFILE, fnm);
    skfnm  = opt2fn_null("-Sk", NFILE, fnm);
    ndxfnm = opt2fn_null("-n", NFILE, fnm);
    tpsfnm = ftp2fn(efTPR, NFILE, fnm);
    trxfnm = ftp2fn(efTRX, NFILE, fnm);

    /* Calculate axis */
    GMX_RELEASE_ASSERT(normal_axis[0] != NULL, "Options inconsistency; normal_axis[0] is NULL");
    if (std::strcmp(normal_axis[0], "x") == 0)
    {
//.........这里部分代码省略.........

int gmx_do_dssp(int argc, char *argv[])
{
    const char        *desc[] = {
        "[THISMODULE] ",
        "reads a trajectory file and computes the secondary structure for",
        "each time frame ",
        "calling the dssp program. If you do not have the dssp program,",
        "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
        "that the dssp executable is located in ",
        "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
        "set an environment variable [TT]DSSP[tt] pointing to the dssp",
        "executable, e.g.: [PAR]",
        "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
        "Since version 2.0.0, dssp is invoked with a syntax that differs",
        "from earlier versions. If you have an older version of dssp,",
        "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
        "By default, do_dssp uses the syntax introduced with version 2.0.0.",
        "Even newer versions (which at the time of writing are not yet released)",
        "are assumed to have the same syntax as 2.0.0.[PAR]",
        "The structure assignment for each residue and time is written to an",
        "[TT].xpm[tt] matrix file. This file can be visualized with for instance",
        "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
        "Individual chains are separated by light grey lines in the [TT].xpm[tt] and",
        "postscript files.",
        "The number of residues with each secondary structure type and the",
        "total secondary structure ([TT]-sss[tt]) count as a function of",
        "time are also written to file ([TT]-sc[tt]).[PAR]",
        "Solvent accessible surface (SAS) per residue can be calculated, both in",
        "absolute values (A^2) and in fractions of the maximal accessible",
        "surface of a residue. The maximal accessible surface is defined as",
        "the accessible surface of a residue in a chain of glycines.",
        "[BB]Note[bb] that the program [gmx-sas] can also compute SAS",
        "and that is more efficient.[PAR]",
        "Finally, this program can dump the secondary structure in a special file",
        "[TT]ssdump.dat[tt] for usage in the program [gmx-chi]. Together",
        "these two programs can be used to analyze dihedral properties as a",
        "function of secondary structure type."
    };
    static gmx_bool    bVerbose;
    static const char *ss_string   = "HEBT";
    static int         dsspVersion = 2;
    t_pargs            pa[]        = {
        { "-v",  FALSE, etBOOL, {&bVerbose},
          "HIDDENGenerate miles of useless information" },
        { "-sss", FALSE, etSTR, {&ss_string},
          "Secondary structures for structure count"},
        { "-ver", FALSE, etINT, {&dsspVersion},
          "DSSP major version. Syntax changed with version 2"}
    };

    t_trxstatus       *status;
    FILE              *tapein;
    FILE              *ss, *acc, *fTArea, *tmpf;
    const char        *fnSCount, *fnArea, *fnTArea, *fnAArea;
    const char        *leg[] = { "Phobic", "Phylic" };
    t_topology         top;
    int                ePBC;
    t_atoms           *atoms;
    t_matrix           mat;
    int                nres, nr0, naccr, nres_plus_separators;
    gmx_bool          *bPhbres, bDoAccSurf;
    real               t;
    int                i, j, natoms, nframe = 0;
    matrix             box = {{0}};
    int                gnx;
    char              *grpnm, *ss_str;
    atom_id           *index;
    rvec              *xp, *x;
    int               *average_area;
    real             **accr, *accr_ptr = NULL, *av_area, *norm_av_area;
    char               pdbfile[32], tmpfile[32], title[256];
    char               dssp[256];
    const char        *dptr;
    output_env_t       oenv;
    gmx_rmpbc_t        gpbc = NULL;

    t_filenm           fnm[] = {
        { efTRX, "-f",   NULL,      ffREAD },
        { efTPS, NULL,   NULL,      ffREAD },
        { efNDX, NULL,   NULL,      ffOPTRD },
        { efDAT, "-ssdump", "ssdump", ffOPTWR },
        { efMAP, "-map", "ss",      ffLIBRD },
        { efXPM, "-o",   "ss",      ffWRITE },
        { efXVG, "-sc",  "scount",  ffWRITE },
        { efXPM, "-a",   "area",    ffOPTWR },
        { efXVG, "-ta",  "totarea", ffOPTWR },
        { efXVG, "-aa",  "averarea", ffOPTWR }
    };
#define NFILE asize(fnm)

    if (!parse_common_args(&argc, argv,
                           PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }
    fnSCount   = opt2fn("-sc", NFILE, fnm);
    fnArea     = opt2fn_null("-a", NFILE, fnm);
    fnTArea    = opt2fn_null("-ta", NFILE, fnm);
    fnAArea    = opt2fn_null("-aa", NFILE, fnm);
//.........这里部分代码省略.........

int gmx_check(int argc, char *argv[])
{
    const char       *desc[] = {
        "[THISMODULE] reads a trajectory ([REF].tng[ref], [REF].trr[ref] or ",
        "[REF].xtc[ref]), an energy file ([REF].edr[ref])",
        "or an index file ([REF].ndx[ref])",
        "and prints out useful information about them.[PAR]",
        "Option [TT]-c[tt] checks for presence of coordinates,",
        "velocities and box in the file, for close contacts (smaller than",
        "[TT]-vdwfac[tt] and not bonded, i.e. not between [TT]-bonlo[tt]",
        "and [TT]-bonhi[tt], all relative to the sum of both Van der Waals",
        "radii) and atoms outside the box (these may occur often and are",
        "no problem). If velocities are present, an estimated temperature",
        "will be calculated from them.[PAR]",
        "If an index file, is given its contents will be summarized.[PAR]",
        "If both a trajectory and a [REF].tpr[ref] file are given (with [TT]-s1[tt])",
        "the program will check whether the bond lengths defined in the tpr",
        "file are indeed correct in the trajectory. If not you may have",
        "non-matching files due to e.g. deshuffling or due to problems with",
        "virtual sites. With these flags, [TT]gmx check[tt] provides a quick check for such problems.[PAR]",
        "The program can compare two run input ([REF].tpr[ref])",
        "files",
        "when both [TT]-s1[tt] and [TT]-s2[tt] are supplied.",
        "Similarly a pair of trajectory files can be compared (using the [TT]-f2[tt]",
        "option), or a pair of energy files (using the [TT]-e2[tt] option).[PAR]",
        "For free energy simulations the A and B state topology from one",
        "run input file can be compared with options [TT]-s1[tt] and [TT]-ab[tt].[PAR]",
        "In case the [TT]-m[tt] flag is given a LaTeX file will be written",
        "consisting of a rough outline for a methods section for a paper."
    };
    t_filenm          fnm[] = {
        { efTRX, "-f",  NULL, ffOPTRD },
        { efTRX, "-f2",  NULL, ffOPTRD },
        { efTPR, "-s1", "top1", ffOPTRD },
        { efTPR, "-s2", "top2", ffOPTRD },
        { efTPS, "-c",  NULL, ffOPTRD },
        { efEDR, "-e",  NULL, ffOPTRD },
        { efEDR, "-e2", "ener2", ffOPTRD },
        { efNDX, "-n",  NULL, ffOPTRD },
        { efTEX, "-m",  NULL, ffOPTWR }
    };
#define NFILE asize(fnm)
    const char       *fn1 = NULL, *fn2 = NULL, *tex = NULL;

    gmx_output_env_t *oenv;
    static real       vdw_fac  = 0.8;
    static real       bon_lo   = 0.4;
    static real       bon_hi   = 0.7;
    static gmx_bool   bRMSD    = FALSE;
    static real       ftol     = 0.001;
    static real       abstol   = 0.001;
    static gmx_bool   bCompAB  = FALSE;
    static char      *lastener = NULL;
    static t_pargs    pa[]     = {
        { "-vdwfac", FALSE, etREAL, {&vdw_fac},
          "Fraction of sum of VdW radii used as warning cutoff" },
        { "-bonlo",  FALSE, etREAL, {&bon_lo},
          "Min. fract. of sum of VdW radii for bonded atoms" },
        { "-bonhi",  FALSE, etREAL, {&bon_hi},
          "Max. fract. of sum of VdW radii for bonded atoms" },
        { "-rmsd",   FALSE, etBOOL, {&bRMSD},
          "Print RMSD for x, v and f" },
        { "-tol",    FALSE, etREAL, {&ftol},
          "Relative tolerance for comparing real values defined as [MATH]2*(a-b)/([MAG]a[mag]+[MAG]b[mag])[math]" },
        { "-abstol",    FALSE, etREAL, {&abstol},
          "Absolute tolerance, useful when sums are close to zero." },
        { "-ab",     FALSE, etBOOL, {&bCompAB},
          "Compare the A and B topology from one file" },
        { "-lastener", FALSE, etSTR,  {&lastener},
          "Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure." }
    };

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

    fn1 = opt2fn_null("-f", NFILE, fnm);
    fn2 = opt2fn_null("-f2", NFILE, fnm);
    tex = opt2fn_null("-m", NFILE, fnm);
    if (fn1 && fn2)
    {
        comp_trx(oenv, fn1, fn2, bRMSD, ftol, abstol);
    }
    else if (fn1)
    {
        chk_trj(oenv, fn1, opt2fn_null("-s1", NFILE, fnm), ftol);
    }
    else if (fn2)
    {
        fprintf(stderr, "Please give me TWO trajectory (.xtc/.trr/.tng) files!\n");
    }

    fn1 = opt2fn_null("-s1", NFILE, fnm);
    fn2 = opt2fn_null("-s2", NFILE, fnm);
    if ((fn1 && fn2) || bCompAB)
    {
        if (bCompAB)
        {
//.........这里部分代码省略.........

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

    t_filenm           fnm[] = {
        { efTRX, NULL, NULL,  ffREAD },
        { efTPS, NULL, NULL,  ffREAD },
        { efNDX, NULL, NULL,  ffOPTRD },
        { efXVG, NULL, "msd", ffWRITE },
        { efXVG, "-mol", "diff_mol", ffOPTWR },
        { efPDB, "-pdb", "diff_mol", ffOPTWR }
    };
#define NFILE asize(fnm)

    t_topology      top;
    int             ePBC;
    matrix          box;
    char            title[256];
    const char     *trx_file, *tps_file, *ndx_file, *msd_file, *mol_file, *pdb_file;
    rvec           *xdum;
    gmx_bool        bTop;
    int             axis, type;
    real            dim_factor;
    output_env_t    oenv;

    if (!parse_common_args(&argc, argv,
                           PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_TIME_UNIT,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }
    trx_file = ftp2fn_null(efTRX, NFILE, fnm);
    tps_file = ftp2fn_null(efTPS, NFILE, fnm);
    ndx_file = ftp2fn_null(efNDX, NFILE, fnm);
    msd_file = ftp2fn_null(efXVG, NFILE, fnm);
    pdb_file = opt2fn_null("-pdb", NFILE, fnm);
    if (pdb_file)
    {
        mol_file = opt2fn("-mol", NFILE, fnm);
    }
    else
    {
        mol_file = opt2fn_null("-mol", NFILE, fnm);
    }

    if (ngroup < 1)
    {
        gmx_fatal(FARGS, "Must have at least 1 group (now %d)", ngroup);
    }
    if (mol_file && ngroup > 1)
    {
        gmx_fatal(FARGS, "With molecular msd can only have 1 group (now %d)",
                  ngroup);
    }


    if (mol_file)
    {
        bMW  = TRUE;
        fprintf(stderr, "Calculating diffusion coefficients for molecules.\n");
    }

    GMX_RELEASE_ASSERT(normtype[0] != 0, "Options inconsistency; normtype[0] is NULL");
    GMX_RELEASE_ASSERT(axtitle[0] != 0, "Options inconsistency; axtitle[0] is NULL");

    if (normtype[0][0] != 'n')
    {
        type       = normtype[0][0] - 'x' + X; /* See defines above */

//.........这里部分代码省略.........
    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);
            }
        }
        if (e_n != n)
        {
            gmx_fatal(FARGS, "Number of energies (%d) does not match number of entries (%d) in %s", e_n, n, opt2fn("-f", NFILE, fnm));

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

//.........这里部分代码省略.........
        }
        /*  realy calc p(r) */
        prframecurrent = calc_radial_distribution_histogram(gsans, x, box, index, isize, binwidth, bMC, bNORM, mcover, seed);
        /* copy prframecurrent -> pr and summ up pr->gr[i] */
        /* allocate and/or resize memory for pr->gr[i] and pr->r[i] */
        if (pr->gr == NULL)
        {
            /* check if we use pr->gr first time */
            snew(pr->gr, prframecurrent->grn);
            snew(pr->r, prframecurrent->grn);
        }
        else
        {
            /* resize pr->gr and pr->r if needed to preven overruns */
            if (prframecurrent->grn > pr->grn)
            {
                srenew(pr->gr, prframecurrent->grn);
                srenew(pr->r, prframecurrent->grn);
            }
        }
        pr->grn      = prframecurrent->grn;
        pr->binwidth = prframecurrent->binwidth;
        /* summ up gr and fill r */
        for (i = 0; i < prframecurrent->grn; i++)
        {
            pr->gr[i] += prframecurrent->gr[i];
            pr->r[i]   = prframecurrent->r[i];
        }
        /* normalize histo */
        normalize_probability(prframecurrent->grn, prframecurrent->gr);
        /* convert p(r) to sq */
        sqframecurrent = convert_histogram_to_intensity_curve(prframecurrent, start_q, end_q, q_step);
        /* print frame data if needed */
        if (opt2fn_null("-prframe", NFILE, fnm))
        {
            snew(hdr, 25);
            snew(suffix, GMX_PATH_MAX);
            /* prepare header */
            sprintf(hdr, "g(r), t = %f", t);
            /* prepare output filename */
            fnmdup = dup_tfn(NFILE, fnm);
            sprintf(suffix, "-t%.2f", t);
            add_suffix_to_output_names(fnmdup, NFILE, suffix);
            fp = xvgropen(opt2fn_null("-prframe", NFILE, fnmdup), hdr, "Distance (nm)", "Probability", oenv);
            for (i = 0; i < prframecurrent->grn; i++)
            {
                fprintf(fp, "%10.6f%10.6f\n", prframecurrent->r[i], prframecurrent->gr[i]);
            }
            done_filenms(NFILE, fnmdup);
            xvgrclose(fp);
            sfree(hdr);
            sfree(suffix);
            sfree(fnmdup);
        }
        if (opt2fn_null("-sqframe", NFILE, fnm))
        {
            snew(hdr, 25);
            snew(suffix, GMX_PATH_MAX);
            /* prepare header */
            sprintf(hdr, "I(q), t = %f", t);
            /* prepare output filename */
            fnmdup = dup_tfn(NFILE, fnm);
            sprintf(suffix, "-t%.2f", t);
            add_suffix_to_output_names(fnmdup, NFILE, suffix);
            fp = xvgropen(opt2fn_null("-sqframe", NFILE, fnmdup), hdr, "q (nm^-1)", "s(q)/s(0)", oenv);
            for (i = 0; i < sqframecurrent->qn; i++)

int main(int argc,char *argv[])
{
    const char *desc[] = {
        "[TT]do_multiprot[tt] ",
        "reads a trajectory file and aligns it to a reference structure  ",
        "each time frame",
        "calling the multiprot program. This allows you to use a reference",
        "structure whose sequence is different than that of the protein in the ",
        "trajectory, since the alignment is based on the geometry, not sequence.",
        "The output of [TT]do_multiprot[tt] includes the rmsd and the number of residues",
        "on which it was calculated.",
        "[PAR]",
        "An aligned trajectory file is generated with the [TT]-ox[tt] option.[PAR]",
        "With the [TT]-cr[tt] option, the number of hits in the alignment is given",
        "per residue. This number can be between 0 and the number of frames, and",
        "indicates the structural conservation of this residue.[PAR]",
        "If you do not have the [TT]multiprot[tt] program, get it. [TT]do_multiprot[tt] assumes",
        "that the [TT]multiprot[tt] executable is [TT]/usr/local/bin/multiprot[tt]. If this is ",
        "not the case, then you should set an environment variable [BB]MULTIPROT[bb]",
        "pointing to the [TT]multiprot[tt] executable, e.g.: [PAR]",
        "[TT]setenv MULTIPROT /usr/MultiProtInstall/multiprot.Linux[tt][PAR]",
        "Note that at the current implementation only binary alignment (your",
        "molecule to a reference) is supported. In addition, note that the ",
        "by default [TT]multiprot[tt] aligns the two proteins on their C-alpha carbons.",
        "and that this depends on the [TT]multiprot[tt] parameters which are not dealt ",
        "with here. Thus, the C-alpha carbons is expected to give the same "
        "results as choosing the whole protein and will be slightly faster.[PAR]",
        "For information about [TT]multiprot[tt], see:",
        "http://bioinfo3d.cs.tau.ac.il/MultiProt/.[PAR]"
    };
    static bool bVerbose;
    t_pargs pa[] = {
        {   "-v",  FALSE, etBOOL, {&bVerbose},
            "HIDDENGenerate miles of useless information"
        }
    };

    const char *bugs[] = {
        "The program is very slow, since multiprot is run externally"
    };

    t_trxstatus *status;
    t_trxstatus *trxout=NULL;
    FILE        *tapein,*fo,*frc,*tmpf,*out=NULL,*fres=NULL;
    const char  *fnRef;
    const char  *fn="2_sol.res";
    t_topology  top;
    int         ePBC;
    t_atoms     *atoms,ratoms,useatoms;
    t_trxframe  fr;
    t_pdbinfo   p;
    int         nres,nres2,nr0;
    real        t;
    int         i,j,natoms,nratoms,nframe=0,model_nr=-1;
    int         cur_res,prev_res;
    int         nout;
    t_countres  *countres=NULL;
    matrix      box,rbox;
    int         gnx;
    char        *grpnm,*ss_str;
    atom_id     *index;
    rvec        *xp,*x,*xr;
    char        pdbfile[32],refpdb[256],title[256],rtitle[256],filemode[5];
    char        out_title[256];
    char        multiprot[256],*mptr;
    int         ftp;
    int         outftp=-1;
    real        rmsd;
    bool        bTrjout,bCountres;
    const char  *TrjoutFile=NULL;
    output_env_t oenv;
    static rvec translation= {0,0,0},rotangles= {0,0,0};
    gmx_rmpbc_t gpbc=NULL;

    t_filenm   fnm[] = {
        { efTRX, "-f",   NULL,      ffREAD },
        { efTPS, NULL,   NULL,      ffREAD },
        { efNDX, NULL,   NULL,      ffOPTRD },
        { efSTX, "-r",   NULL     , ffREAD },
        { efXVG, "-o",  "rmss",     ffWRITE },
        { efXVG, "-rc", "rescount", ffWRITE},
        { efXVG, "-cr", "countres", ffOPTWR},
        { efTRX, "-ox", "aligned",  ffOPTWR }
    };
#define NFILE asize(fnm)

    CopyRight(stderr,argv[0]);
    parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE ,
                      NFILE,fnm, asize(pa),pa, asize(desc),desc,
                      asize(bugs),bugs,&oenv
                     );
    fnRef=opt2fn("-r",NFILE,fnm);
    bTrjout = opt2bSet("-ox",NFILE,fnm);
    bCountres=  opt2bSet("-cr",NFILE,fnm);

    if (bTrjout) {
        TrjoutFile = opt2fn_null("-ox",NFILE,fnm);
    }

    read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xp,NULL,box,FALSE);
//.........这里部分代码省略.........

int	gmx_energy2bfac (int argc, char *argv[])		{
  const char *desc[] =
    { "It maps the binding energy contribution of each residue on the structure.",
      "The energy will be written in the B-factor field of the output PDB file/s",
      "This PDB file can be used with any molecular visualizer and ",
      "residues can be viewed in color a