C++ parse_common_args函数代码示例

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.",
        "",
        "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 boolean operations, you can split groups",
        "into chains, residues or atoms. You can delete and rename groups.",
        "Type 'h' in the editor for more details.",
        "",
        "The atom numbering in the editor and the index file starts at 1.",
        "",
        "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.",
        "",
        "See also [gmx-select] [TT]-on[tt], which provides an alternative way",
        "for constructing index groups.  It covers nearly all of [THISMODULE]",
        "functionality, and in many cases much more."
    };

    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)

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

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

    stxfile = ftp2fn_null(efSTX, NFILE, fnm);
    gmx::ArrayRef<const std::string> ndxInFiles = opt2fnsIfOptionSet("-n", NFILE, fnm);
    ndxoutfile = opt2fn("-o", NFILE, fnm);
    bNatoms    = opt2parg_bSet("-natoms", NPA, pa);

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

    if (stxfile)
    {
        t_topology *top;
        snew(top, 1);
        fprintf(stderr, "\nReading structure file\n");
        read_tps_conf(stxfile, top, &ePBC, &x, &v, box, FALSE);
        atoms = &top->atoms;
        if (atoms->pdbinfo == nullptr)
        {
            snew(atoms->pdbinfo, atoms->nr);
        }
        natoms  = atoms->nr;
        bNatoms = TRUE;
    }
    else
    {
        atoms = nullptr;
        x     = nullptr;
    }

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

//.........这里部分代码省略.........
    t_trxframe        fr;
    rvec             *xtop;
    matrix            box, box_pbc;
    t_trxstatus      *status;
    int               flags = TRX_READ_X;
    t_pbc             pbc;
    t_atoms          *atoms;
    int               natoms;
    char             *grpnm, *grpnmp;
    int              *index, *indexp;
    int               i, nidx, nidxp;
    int               v;
    int               j, k;
    int            ***bin = NULL;
    int               nbin[3];
    FILE             *flp;
    int               x, y, z, minx, miny, minz, maxx, maxy, maxz;
    int               numfr, numcu;
    int               tot, maxval, minval;
    double            norm;
    gmx_output_env_t *oenv;
    gmx_rmpbc_t       gpbc = NULL;

    t_filenm          fnm[] = {
        { efTPS,  NULL,  NULL, ffREAD }, /* this is for the topology */
        { efTRX, "-f", NULL, ffREAD },   /* and this for the trajectory */
        { efNDX, NULL, NULL, ffOPTRD }
    };

#define NFILE asize(fnm)

    /* This is the routine responsible for adding default options,
     * calling the X/motif interface, etc. */
    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
    {
        return 0;
    }

    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, NULL, box, TRUE);
    sfree(xtop);

    atoms = &(top.atoms);
    printf("Select group to generate SDF:\n");
    get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidx, &index, &grpnm);
    printf("Select group to output coords (e.g. solute):\n");
    get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidxp, &indexp, &grpnmp);

    /* The first time we read data is a little special */
    natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);

    /* Memory Allocation */
    MINBIN[XX] = MAXBIN[XX] = fr.x[0][XX];
    MINBIN[YY] = MAXBIN[YY] = fr.x[0][YY];
    MINBIN[ZZ] = MAXBIN[ZZ] = fr.x[0][ZZ];
    for (i = 1; i < top.atoms.nr; ++i)
    {
        if (fr.x[i][XX] < MINBIN[XX])
        {
            MINBIN[XX] = fr.x[i][XX];
        }
        if (fr.x[i][XX] > MAXBIN[XX])
        {
            MAXBIN[XX] = fr.x[i][XX];
        }
        if (fr.x[i][YY] < MINBIN[YY])

int gmx_pme_error(int argc,char *argv[])
{
    const char *desc[] = {
            "g_pme_error 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.",
            "As a part of the error stems from self interaction of the particles "
            "and is computationally very demanding a good a approximation is possible",
            "if just a fraction of the particles is used to calculate the average",
            "of this error by using the flag [TT]-self[tt].[PAR]",
    };

    int        repeats=2;
    real       fs=0.0;             /* 0 indicates: not set by the user */

    real        user_beta=-1.0;
    real        fracself=-1.0;
    
    
    t_perf      **perfdata;
    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;


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


    output_env_t oenv=NULL;

    t_pargs pa[] = {
      /***********************/
      /* g_tune_pme options: */
      /***********************/
        { "-beta",     FALSE, etREAL, {&user_beta},
            "If positive, overwrite ewald_beta from tpr file with this value" },
        { "-tune",     FALSE, etBOOL, {&bTUNE},
            "If flag is set the splitting parameter will be tuned to distribute the error equally in real and rec. space" },
        { "-self",     FALSE, etREAL, {&fracself},
            "If positive, determine selfinteraction error just over this fraction (default=1.0)" }
    };

    
#define NFILE asize(fnm)
    
    cr = init_par(&argc,&argv);
    
    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 */
//.........这里部分代码省略.........

int gmx_sorient(int argc, char *argv[])
{
    t_topology        top;
    int               ePBC = -1;
    t_trxstatus      *status;
    int               natoms;
    real              t;
    rvec             *xtop, *x;
    matrix            box;

    FILE             *fp;
    int               i, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin;
    real             *histi1, *histi2, invbw, invrbw;
    double            sum1, sum2;
    int              *isize, nrefgrp, nrefat;
    int             **index;
    char            **grpname;
    real              inp, outp, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r;
    real              c1, c2;
    char              str[STRLEN];
    gmx_bool          bTPS;
    rvec              xref, dx, dxh1, dxh2, outer;
    gmx_rmpbc_t       gpbc = NULL;
    t_pbc             pbc;
    const char       *legr[] = {
        "<cos(\\8q\\4\\s1\\N)>",
        "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
    };
    const char       *legc[] = {
        "cos(\\8q\\4\\s1\\N)",
        "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
    };

    const char       *desc[] = {
        "[THISMODULE] analyzes solvent orientation around solutes.",
        "It calculates two angles between the vector from one or more",
        "reference positions to the first atom of each solvent molecule:",
        "",
        " * [GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
        "   molecule to the midpoint between atoms 2 and 3.",
        " * [GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
        "   same three atoms, or, when the option [TT]-v23[tt] is set, ",
        "   the angle with the vector between atoms 2 and 3.",
        "",
        "The reference can be a set of atoms or",
        "the center of mass of a set of atoms. The group of solvent atoms should",
        "consist of 3 atoms per solvent molecule.",
        "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
        "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
        "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
        "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
        "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a function of the",
        "distance.[PAR]",
        "[TT]-co[tt]: the sum over all solvent molecules within distance r",
        "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and [MATH]3[COS]^2([GRK]theta[grk][SUB]2[sub])-1[cos][math] as a function of r.[PAR]",
        "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
    };

    gmx_output_env_t *oenv;
    static gmx_bool   bCom = FALSE, bVec23 = FALSE, bPBC = FALSE;
    static real       rmin = 0.0, rmax = 0.5, binwidth = 0.02, rbinw = 0.02;
    t_pargs           pa[] = {
        {   "-com",  FALSE, etBOOL,  {&bCom},
            "Use the center of mass as the reference postion"
        },
        {   "-v23",  FALSE, etBOOL,  {&bVec23},
            "Use the vector between atoms 2 and 3"
        },
        { "-rmin",  FALSE, etREAL, {&rmin}, "Minimum distance (nm)" },
        { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
        { "-cbin",  FALSE, etREAL, {&binwidth}, "Binwidth for the cosine" },
        { "-rbin",  FALSE, etREAL, {&rbinw}, "Binwidth for r (nm)" },
        { "-pbc",   FALSE, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
    };

    t_filenm          fnm[] = {
        { efTRX, NULL,  NULL,  ffREAD },
        { efTPS, NULL,  NULL,  ffREAD },
        { efNDX, NULL,  NULL,  ffOPTRD },
        { efXVG, NULL,  "sori",   ffWRITE },
        { efXVG, "-no", "snor",   ffWRITE },
        { efXVG, "-ro", "sord",   ffWRITE },
        { efXVG, "-co", "scum",   ffWRITE },
        { efXVG, "-rc", "scount", ffWRITE }
    };
#define NFILE asize(fnm)

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

    bTPS = (opt2bSet("-s", NFILE, fnm) || !opt2bSet("-n", NFILE, fnm) || bCom);
    if (bTPS)
    {
        read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, NULL, box,
                      bCom);
    }

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

int gmx_dielectric(int argc, char *argv[])
{
    const char  *desc[] = {
        "[THISMODULE] calculates frequency dependent dielectric constants",
        "from the autocorrelation function of the total dipole moment in",
        "your simulation. This ACF can be generated by [gmx-dipoles].",
        "The functional forms of the available functions are:[PAR]",
        "One parameter:    y = [EXP]-a[SUB]1[sub] x[exp],[BR]",
        "Two parameters:   y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp],[BR]",
        "Three parameters: y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp] + (1 - a[SUB]2[sub]) [EXP]-a[SUB]3[sub] x[exp].[BR]",
        "Start values for the fit procedure can be given on the command line.",
        "It is also possible to fix parameters at their start value, use [TT]-fix[tt]",
        "with the number of the parameter you want to fix.",
        "[PAR]",
        "Three output files are generated, the first contains the ACF,",
        "an exponential fit to it with 1, 2 or 3 parameters, and the",
        "numerical derivative of the combination data/fit.",
        "The second file contains the real and imaginary parts of the",
        "frequency-dependent dielectric constant, the last gives a plot",
        "known as the Cole-Cole plot, in which the imaginary",
        "component is plotted as a function of the real component.",
        "For a pure exponential relaxation (Debye relaxation) the latter",
        "plot should be one half of a circle."
    };
    t_filenm     fnm[] = {
        { efXVG, "-f", "dipcorr", ffREAD  },
        { efXVG, "-d", "deriv",  ffWRITE },
        { efXVG, "-o", "epsw",   ffWRITE },
        { efXVG, "-c", "cole",   ffWRITE }
    };
#define NFILE asize(fnm)
    output_env_t oenv;
    int          i, j, nx, ny, nxtail, eFitFn, nfitparm;
    real         dt, integral, fitintegral, fac, rffac;
    double      *fitparms;
    double     **yd;
    real       **y;
    const char  *legend[] = { "Correlation", "Std. Dev.", "Fit", "Combined", "Derivative" };
    static int   fix      = 0, bX = 1, nsmooth = 3;
    static real  tendInt  = 5.0, tbegin = 5.0, tend = 500.0;
    static real  A        = 0.5, tau1 = 10.0, tau2 = 1.0, eps0 = 80, epsRF = 78.5, tail = 500.0;
    real         lambda;
    t_pargs      pa[] = {
        { "-x1",  FALSE, etBOOL, {&bX},
          "use first column as [IT]x[it]-axis rather than first data set" },
        { "-eint", FALSE, etREAL, {&tendInt},
          "Time to end the integration of the data and start to use the fit"},
        { "-bfit", FALSE, etREAL, {&tbegin},
          "Begin time of fit" },
        { "-efit", FALSE, etREAL, {&tend},
          "End time of fit" },
        { "-tail", FALSE, etREAL, {&tail},
          "Length of function including data and tail from fit" },
        { "-A", FALSE, etREAL, {&A},
          "Start value for fit parameter A" },
        { "-tau1", FALSE, etREAL, {&tau1},
          "Start value for fit parameter [GRK]tau[grk]1" },
        { "-tau2", FALSE, etREAL, {&tau2},
          "Start value for fit parameter [GRK]tau[grk]2" },
        { "-eps0", FALSE, etREAL, {&eps0},
          "[GRK]epsilon[grk]0 of your liquid" },
        { "-epsRF", FALSE, etREAL, {&epsRF},
          "[GRK]epsilon[grk] of the reaction field used in your simulation. A value of 0 means infinity." },
        { "-fix", FALSE, etINT,  {&fix},
          "Fix parameters at their start values, A (2), tau1 (1), or tau2 (4)" },
        { "-ffn",    FALSE, etENUM, {s_ffn},
          "Fit function" },
        { "-nsmooth", FALSE, etINT, {&nsmooth},
          "Number of points for smoothing" }
    };

    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }
    please_cite(stdout, "Spoel98a");
    printf("WARNING: non-polarizable models can never yield an infinite\n"
           "dielectric constant that is different from 1. This is incorrect\n"
           "in the reference given above (Spoel98a).\n\n");


    nx     = read_xvg(opt2fn("-f", NFILE, fnm), &yd, &ny);
    dt     = yd[0][1] - yd[0][0];
    nxtail = min(tail/dt, nx);

    printf("Read data set containing %d colums and %d rows\n", ny, nx);
    printf("Assuming (from data) that timestep is %g, nxtail = %d\n",
           dt, nxtail);
    snew(y, 6);
    for (i = 0; (i < ny); i++)
    {
        snew(y[i], max(nx, nxtail));
    }
    for (i = 0; (i < nx); i++)
    {
        y[0][i] = yd[0][i];
        for (j = 1; (j < ny); j++)
        {
            y[j][i] = yd[j][i];
//.........这里部分代码省略.........

int gmx_mindist(int argc, char *argv[])
{
    const char     *desc[] = {
        "[THISMODULE] computes the distance between one group and a number of",
        "other groups. Both the minimum distance",
        "(between any pair of atoms from the respective groups)",
        "and the number of contacts within a given",
        "distance are written to two separate output files.",
        "With the [TT]-group[tt] option a contact of an atom in another group",
        "with multiple atoms in the first group is counted as one contact",
        "instead of as multiple contacts.",
        "With [TT]-or[tt], minimum distances to each residue in the first",
        "group are determined and plotted as a function of residue number.[PAR]",
        "With option [TT]-pi[tt] the minimum distance of a group to its",
        "periodic image is plotted. This is useful for checking if a protein",
        "has seen its periodic image during a simulation. Only one shift in",
        "each direction is considered, giving a total of 26 shifts.",
        "It also plots the maximum distance within the group and the lengths",
        "of the three box vectors.[PAR]",
        "Also [gmx-distance] calculates distances."
    };
    const char     *bugs[] = {
        "The [TT]-pi[tt] option is very slow."
    };

    static gmx_bool bMat             = FALSE, bPI = FALSE, bSplit = FALSE, bMax = FALSE, bPBC = TRUE;
    static gmx_bool bGroup           = FALSE;
    static real     rcutoff          = 0.6;
    static int      ng               = 1;
    static gmx_bool bEachResEachTime = FALSE, bPrintResName = FALSE;
    t_pargs         pa[]             = {
        { "-matrix", FALSE, etBOOL, {&bMat},
          "Calculate half a matrix of group-group distances" },
        { "-max",    FALSE, etBOOL, {&bMax},
          "Calculate *maximum* distance instead of minimum" },
        { "-d",      FALSE, etREAL, {&rcutoff},
          "Distance for contacts" },
        { "-group",      FALSE, etBOOL, {&bGroup},
          "Count contacts with multiple atoms in the first group as one" },
        { "-pi",     FALSE, etBOOL, {&bPI},
          "Calculate minimum distance with periodic images" },
        { "-split",  FALSE, etBOOL, {&bSplit},
          "Split graph where time is zero" },
        { "-ng",       FALSE, etINT, {&ng},
          "Number of secondary groups to compute distance to a central group" },
        { "-pbc",    FALSE, etBOOL, {&bPBC},
          "Take periodic boundary conditions into account" },
        { "-respertime",  FALSE, etBOOL, {&bEachResEachTime},
          "When writing per-residue distances, write distance for each time point" },
        { "-printresname",  FALSE, etBOOL, {&bPrintResName},
          "Write residue names" }
    };
    output_env_t    oenv;
    t_topology     *top  = NULL;
    int             ePBC = -1;
    char            title[256];
    real            t;
    rvec           *x;
    matrix          box;
    gmx_bool        bTop = FALSE;

    FILE           *atm;
    int             i, j, nres = 0;
    const char     *trxfnm, *tpsfnm, *ndxfnm, *distfnm, *numfnm, *atmfnm, *oxfnm, *resfnm;
    char          **grpname;
    int            *gnx;
    atom_id       **index, *residues = NULL;
    t_filenm        fnm[] = {
        { efTRX, "-f",  NULL,      ffREAD },
        { efTPS,  NULL, NULL,      ffOPTRD },
        { efNDX,  NULL, NULL,      ffOPTRD },
        { efXVG, "-od", "mindist",  ffWRITE },
        { efXVG, "-on", "numcont",  ffOPTWR },
        { efOUT, "-o", "atm-pair", ffOPTWR },
        { efTRO, "-ox", "mindist",  ffOPTWR },
        { efXVG, "-or", "mindistres", ffOPTWR }
    };
#define NFILE asize(fnm)

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

    trxfnm  = ftp2fn(efTRX, NFILE, fnm);
    ndxfnm  = ftp2fn_null(efNDX, NFILE, fnm);
    distfnm = opt2fn("-od", NFILE, fnm);
    numfnm  = opt2fn_null("-on", NFILE, fnm);
    atmfnm  = ftp2fn_null(efOUT, NFILE, fnm);
    oxfnm   = opt2fn_null("-ox", NFILE, fnm);
    resfnm  = opt2fn_null("-or", NFILE, fnm);
    if (bPI || resfnm != NULL)
    {
        /* We need a tps file */
        tpsfnm = ftp2fn(efTPS, NFILE, fnm);
    }
    else
    {
//.........这里部分代码省略.........

int gmx_spol(int argc, char *argv[])
{
    t_topology  *top;
    t_inputrec  *ir;
    t_atom      *atom;
    char         title[STRLEN];
    t_trxstatus *status;
    int          nrefat, natoms, nf, ntot;
    real         t;
    rvec        *xtop, *x, xref, trial, dx = {0}, dip, dir;
    matrix       box;

    FILE        *fp;
    int         *isize, nrefgrp;
    atom_id    **index, *molindex;
    char       **grpname;
    real         rmin2, rmax2, rcut, rcut2, rdx2 = 0, rtry2, qav, q, dip2, invbw;
    int          nbin, i, m, mol, a0, a1, a, d;
    double       sdip, sdip2, sinp, sdinp, nmol;
    int         *hist;
    t_pbc        pbc;
    gmx_rmpbc_t  gpbc = NULL;


    const char     *desc[] = {
        "[THISMODULE] analyzes dipoles around a solute; it is especially useful",
        "for polarizable water. A group of reference atoms, or a center",
        "of mass reference (option [TT]-com[tt]) and a group of solvent",
        "atoms is required. The program splits the group of solvent atoms",
        "into molecules. For each solvent molecule the distance to the",
        "closest atom in reference group or to the COM is determined.",
        "A cumulative distribution of these distances is plotted.",
        "For each distance between [TT]-rmin[tt] and [TT]-rmax[tt]",
        "the inner product of the distance vector",
        "and the dipole of the solvent molecule is determined.",
        "For solvent molecules with net charge (ions), the net charge of the ion",
        "is subtracted evenly from all atoms in the selection of each ion.",
        "The average of these dipole components is printed.",
        "The same is done for the polarization, where the average dipole is",
        "subtracted from the instantaneous dipole. The magnitude of the average",
        "dipole is set with the option [TT]-dip[tt], the direction is defined",
        "by the vector from the first atom in the selected solvent group",
        "to the midpoint between the second and the third atom."
    };

    output_env_t    oenv;
    static gmx_bool bCom   = FALSE, bPBC = FALSE;
    static int      srefat = 1;
    static real     rmin   = 0.0, rmax = 0.32, refdip = 0, bw = 0.01;
    t_pargs         pa[]   = {
        { "-com",  FALSE, etBOOL,  {&bCom},
          "Use the center of mass as the reference postion" },
        { "-refat",  FALSE, etINT, {&srefat},
          "The reference atom of the solvent molecule" },
        { "-rmin",  FALSE, etREAL, {&rmin}, "Maximum distance (nm)" },
        { "-rmax",  FALSE, etREAL, {&rmax}, "Maximum distance (nm)" },
        { "-dip",   FALSE, etREAL, {&refdip}, "The average dipole (D)" },
        { "-bw",    FALSE, etREAL, {&bw}, "The bin width" }
    };

    t_filenm        fnm[] = {
        { efTRX, NULL,  NULL,  ffREAD },
        { efTPR, NULL,  NULL,  ffREAD },
        { efNDX, NULL,  NULL,  ffOPTRD },
        { efXVG, NULL,  "scdist",  ffWRITE }
    };
#define NFILE asize(fnm)

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

    snew(top, 1);
    snew(ir, 1);
    read_tpx_top(ftp2fn(efTPR, NFILE, fnm),
                 ir, box, &natoms, NULL, NULL, NULL, top);

    /* get index groups */
    printf("Select a group of reference particles and a solvent group:\n");
    snew(grpname, 2);
    snew(index, 2);
    snew(isize, 2);
    get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), 2, isize, index, grpname);

    if (bCom)
    {
        nrefgrp = 1;
        nrefat  = isize[0];
    }
    else
    {
        nrefgrp = isize[0];
        nrefat  = 1;
    }

    spol_atom2molindex(&(isize[1]), index[1], &(top->mols));
    srefat--;

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

int gmx_trjorder(int argc, char *argv[])
{
    const char       *desc[] = {
        "[THISMODULE] orders molecules according to the smallest distance",
        "to atoms in a reference group",
        "or on z-coordinate (with option [TT]-z[tt]).",
        "With distance ordering, it will ask for a group of reference",
        "atoms and a group of molecules. For each frame of the trajectory",
        "the selected molecules will be reordered according to the shortest",
        "distance between atom number [TT]-da[tt] in the molecule and all the",
        "atoms in the reference group. The center of mass of the molecules can",
        "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
        "All atoms in the trajectory are written",
        "to the output trajectory.[PAR]",
        "[THISMODULE] can be useful for e.g. analyzing the n waters closest to a",
        "protein.",
        "In that case the reference group would be the protein and the group",
        "of molecules would consist of all the water atoms. When an index group",
        "of the first n waters is made, the ordered trajectory can be used",
        "with any GROMACS program to analyze the n closest waters.",
        "[PAR]",
        "If the output file is a [REF].pdb[ref] file, the distance to the reference target",
        "will be stored in the B-factor field in order to color with e.g. Rasmol.",
        "[PAR]",
        "With option [TT]-nshell[tt] the number of molecules within a shell",
        "of radius [TT]-r[tt] around the reference group are printed."
    };
    static int        na   = 3, ref_a = 1;
    static real       rcut = 0;
    static gmx_bool   bCOM = FALSE, bZ = FALSE;
    t_pargs           pa[] = {
        { "-na", FALSE, etINT,  {&na},
          "Number of atoms in a molecule" },
        { "-da", FALSE, etINT,  {&ref_a},
          "Atom used for the distance calculation, 0 is COM" },
        { "-com", FALSE, etBOOL, {&bCOM},
          "Use the distance to the center of mass of the reference group" },
        { "-r",  FALSE, etREAL, {&rcut},
          "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
        { "-z", FALSE, etBOOL, {&bZ},
          "Order molecules on z-coordinate" }
    };
    FILE             *fp;
    t_trxstatus      *out;
    t_trxstatus      *status;
    gmx_bool          bNShell, bPDBout;
    t_topology        top;
    int               ePBC;
    rvec             *x, *xsol, xcom, dx;
    matrix            box;
    t_pbc             pbc;
    gmx_rmpbc_t       gpbc;
    real              t, totmass, mass, rcut2 = 0, n2;
    int               natoms, nwat, ncut;
    char            **grpname;
    int               i, j, d, *isize, isize_ref = 0, isize_sol;
    int               sa, sr, *swi, **index, *ind_ref = NULL, *ind_sol;
    gmx_output_env_t *oenv;
    t_filenm          fnm[] = {
        { efTRX, "-f", NULL, ffREAD  },
        { efTPS, NULL, NULL, ffREAD  },
        { efNDX, NULL, NULL, ffOPTRD },
        { efTRO, "-o", "ordered", ffOPTWR },
        { efXVG, "-nshell", "nshell", ffOPTWR }
    };
#define NFILE asize(fnm)

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

    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, NULL, box, TRUE);
    sfree(x);

    /* get index groups */
    printf("Select %sa group of molecules to be ordered:\n",
           bZ ? "" : "a group of reference atoms and ");
    snew(grpname, 2);
    snew(index, 2);
    snew(isize, 2);
    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), bZ ? 1 : 2,
              isize, index, grpname);

    if (!bZ)
    {
        isize_ref = isize[0];
        isize_sol = isize[1];
        ind_ref   = index[0];
        ind_sol   = index[1];
    }
    else
    {
        isize_sol = isize[0];
        ind_sol   = index[0];
    }

    natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
    if (natoms > top.atoms.nr)
//.........这里部分代码省略.........

int gmx_velacc(int argc,char *argv[])
{
    const char *desc[] = {
        "[TT]g_velacc[tt] computes the velocity autocorrelation function.",
        "When the [TT]-m[tt] option is used, the momentum autocorrelation",
        "function is calculated.[PAR]",
        "With option [TT]-mol[tt] the velocity autocorrelation function of",
        "molecules is calculated. In this case the index group should consist",
        "of molecule numbers instead of atom numbers.[PAR]",
        "Be sure that your trajectory contains frames with velocity information",
        "(i.e. [TT]nstvout[tt] was set in your original [TT].mdp[tt] file),",
        "and that the time interval between data collection points is",
        "much shorter than the time scale of the autocorrelation."
    };
  
    static gmx_bool bMass=FALSE,bMol=FALSE,bRecip=TRUE;
    t_pargs pa[] = {
        { "-m", FALSE, etBOOL, {&bMass},
          "Calculate the momentum autocorrelation function" },
        { "-recip", FALSE, etBOOL, {&bRecip},
          "Use cm^-1 on X-axis instead of 1/ps for spectra." },
        { "-mol", FALSE, etBOOL, {&bMol},
          "Calculate the velocity acf of molecules" }
    };

    t_topology top;
    int        ePBC=-1;
    t_trxframe fr;
    matrix     box;
    gmx_bool       bTPS=FALSE,bTop=FALSE;
    int        gnx;
    atom_id    *index;
    char       *grpname;
    char       title[256];
    /* t0, t1 are the beginning and end time respectively. 
     * dt is the time step, mass is temp variable for atomic mass.
     */
    real       t0,t1,dt,mass;
    t_trxstatus *status;
    int        counter,n_alloc,i,j,counter_dim,k,l;
    rvec       mv_mol;
    /* Array for the correlation function */
    real       **c1;
    real             *normm=NULL;
    output_env_t oenv;
  
#define NHISTO 360
    
    t_filenm  fnm[] = {
        { efTRN, "-f",    NULL,   ffREAD  },
        { efTPS, NULL,    NULL,   ffOPTRD }, 
        { efNDX, NULL,    NULL,   ffOPTRD },
        { efXVG, "-o",    "vac",  ffWRITE },
        { efXVG, "-os",   "spectrum", ffOPTWR }
    };
#define NFILE asize(fnm)
    int     npargs;
    t_pargs *ppa;

    npargs = asize(pa);
    ppa    = add_acf_pargs(&npargs,pa);
    parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
                      NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);

    if (bMol || bMass) {
        bTPS = ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm);
    }

    if (bTPS) {
        bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
                           TRUE);
        get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
    } else
        rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);

    if (bMol) {
        if (!bTop)
            gmx_fatal(FARGS,"Need a topology to determine the molecules");
        snew(normm,top.atoms.nr);
        precalc(top,normm);
        index_atom2mol(&gnx,index,&top.mols);
    }
  
    /* Correlation stuff */
    snew(c1,gnx);
    for(i=0; (i<gnx); i++)
        c1[i]=NULL;
  
    read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,TRX_NEED_V);
    t0=fr.time;
      
    n_alloc=0;
    counter=0;
    do {
        if (counter >= n_alloc) {
            n_alloc+=100;
            for(i=0; i<gnx; i++)
                srenew(c1[i],DIM*n_alloc);
        }
        counter_dim=DIM*counter;
//.........这里部分代码省略.........

int gmx_dih(int argc,char *argv[])
{
  const char *desc[] = {
    "g_dih can do two things. The default is to analyze dihedral transitions",
    "by merely computing all the dihedral angles defined in your topology",
    "for the whole trajectory. When a dihedral flips over to another minimum",
    "an angle/time plot is made.[PAR]",
    "The opther option is to discretize the dihedral space into a number of",
    "bins, and group each conformation in dihedral space in the",
    "appropriate bin. The output is then given as a number of dihedral",
    "conformations sorted according to occupancy."
  };
  static int  mult = -1;
  static gmx_bool bSA  = FALSE;
  t_pargs pa[] = {
    { "-sa", FALSE, etBOOL, {&bSA},
      "Perform cluster analysis in dihedral space instead of analysing dihedral transitions." },
    { "-mult", FALSE, etINT, {&mult},
      "mulitiplicity for dihedral angles (by default read from topology)" }
  };
  FILE       *out;
  t_xrama    *xr;
  t_topology *top;
  real       **dih,*time;
  real       dd;
  int        i,nframes,maxframes=1000;
  output_env_t oenv;
  t_filenm   fnm[] = {
    { efTRX, "-f", NULL, ffREAD },
    { efTPX, NULL, NULL, ffREAD },
    { efOUT, NULL, NULL, ffWRITE }
  };
#define NFILE asize(fnm)

  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
  
  if (mult != -1)
    fprintf(stderr,"Using %d for dihedral multiplicity rather than topology values\n",mult);
    
  snew(xr,1);
  init_rama(oenv,ftp2fn(efTRX,NFILE,fnm),
	    ftp2fn(efTPX,NFILE,fnm),xr,3);
  top=read_top(ftp2fn(efTPX,NFILE,fnm),NULL);
	       
  /* Brute force malloc, may be too big... */
  snew(dih,xr->ndih);
  for(i=0; (i<xr->ndih); i++)
    snew(dih[i],maxframes);
  snew(time,maxframes);

  fprintf(stderr,"\n");
  nframes = 0;
  while (new_data(xr)) {
    for(i=0; (i<xr->ndih); i++) {
      dd=xr->dih[i].ang*RAD2DEG;
      while (dd < 0)
	dd+=360;
      while (dd > 360)
	dd-=360;
      dih[i][nframes]=dd;
    }
    time[nframes]=xr->t;
    nframes++;
    if (nframes > maxframes) {
      maxframes += 1000;
      for(i=0; (i<xr->ndih); i++)
	srenew(dih[i],maxframes);
      srenew(time,maxframes);
    }
  } 

  fprintf(stderr,"\nCalculated all dihedrals, now analysing...\n");

  out=ftp2FILE(efOUT,NFILE,fnm,"w");

  if (bSA) {
    /* Cluster and structure analysis */
    ana_cluster(out,xr,dih,time,top,nframes,mult);
  }
  else {
    /* Analyse transitions... */
    ana_trans(out,xr,dih,time,top,nframes,oenv);
  }
  ffclose(out);
    
  thanx(stderr);
    
  return 0;
}

//.........这里部分代码省略.........
        { "-yonce",   FALSE, etBOOL, {&bYonce}, "Show y-label only once" },
        { "-legend",  FALSE, etENUM, {legend},  "Show legend" },
        { "-diag",    FALSE, etENUM, {diag},    "Diagonal" },
        { "-size",    FALSE, etREAL, {&size},
          "Horizontal size of the matrix in ps units" },
        { "-bx",      FALSE, etREAL, {&boxx},
          "Element x-size, overrides [TT]-size[tt] (also y-size when [TT]-by[tt] is not set)" },
        { "-by",      FALSE, etREAL, {&boxy},   "Element y-size" },
        { "-rainbow", FALSE, etENUM, {rainbow},
          "Rainbow colors, convert white to" },
        { "-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 [TT].xpm[tt] matrix where axis label is zero"},
        { "-legoffset", FALSE, etINT, {&mapoffset},
          "Skip first N colors from [TT].xpm[tt] 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)

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

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

int main(int argc, char *argv[]) {
    const char *desc[] = {
        "g_ensemble_comp evaluates the difference between two conformational ensembles, R and R'.",
        "Quanitification is in terms of a true metric that satisfies the conditions set forth by the zeroth law of thermodynamics.",
        "The quantification metric eta=1-|Overlap|=|R'|-|Overlap|=DeltaR is normalized, that is, 0<=eta<1,",
        "and takes up a value closer to unity as the difference between the ensembles increases.",
        "For two Gaussian distributions with identical standard deviations of 0.5 A,",
        "eta=0.68 represents a geometric center deviation of 1 A.\n",
        "The two ensembles are provided as two trajectory files specified by the -f1 and -f2 options (supported formats=xtc,trr,pdb).",
        "We recommend that frames numbers in the trajectory files are in the range 2500-5000.",
        "While the speed of the algorithm decreases with increase in ensemble size,", 
        "the numerical accuracy of the calculation reduces with decrease in ensemble size,",
        "and a small number of frames may not provide a good representation of the ensemble.",
        "Note that if you are not interested in eta to reflect changes in whole molecule translation/rotation,",
        "then these degress of freedom need to be removed prior to ensemble comparison.",
        "For most cases, this can be accomplished by fitting all conformations",
        "in the two ensembles on to one single representative structure, such as the x-ray structure",
        "For this, we recommend the use of trjconv with the -fit rot+trans option.\n"
        "By default, differences (eta) are estimated for all atoms, but comparisons can be done for a smaller specific group of atoms,",
        "which can be selected from index files -n1 and -n2. ",
        "The average eta can also be calculated for each residue specified in a structure file given by -res (.pdb and .gro are supported).\n",
        "Overlaps are estimated by training a support vector machine in a pre-defined Hilbert space specified by",
        "the width of the RDF Kernel (gamma=0.4) and",
        "the maximum value that can be taken up by the Lagrange multiplier (C=100.0).", 
        "The values of C and gamma can be changed with -c and -g,",
        "but such changes may increase mean absolute error (MAE=3.26%) of the method.\n",
        "If g_ensemble_comp was built with OPENMP, you can set the number of threads to use with -nthreads X,\n",
        "where X is the number of threads to use. The default behavior is to use the maximum number of cores available.\n\n",
        "Methodoligical details and example applications can be found in\n",
        "Leighty and Varma, JCTC, 2013, 9: 868-875.\n",
        "Varma, Botlani and Leighty, Proteins, 2014, 82: 3241-3254.\n",
        "Dutta, Botlani and Varma, JPC B, 2014, 118: 14795-14807.\n"
    };
    const char *fnames[eNUMFILES];
    output_env_t oenv = NULL;
    real gamma = GAMMA, c = COST;
    int nthreads = -1;

    /* eta values */
    real *eta;
    int natoms;

    init_log("eta.log", argv[0]);

    t_filenm fnm[] = {
        {efTRX, "-f1", "traj1.xtc", ffREAD},
        {efTRX, "-f2", "traj2.xtc", ffREAD},
        {efNDX, "-n1", "index1.ndx", ffOPTRD},
        {efNDX, "-n2", "index2.ndx", ffOPTRD},
        {efSTX, "-res", "res.pdb", ffOPTRD},
        {efDAT, "-eta_atom", "eta_atom.dat", ffWRITE},
        {efDAT, "-eta_res", "eta_res.dat", ffWRITE}
    };

    t_pargs pa[] = {
        {"-g", FALSE, etREAL, {&gamma}, "RBD Kernel width (default=0.4)"},
        {"-c", FALSE, etREAL, {&c}, "Max value of Lagrange multiplier (default=100)"},
        {"-nthreads", FALSE, etINT, {&nthreads}, "set the number of parallel threads to use (default is max available)"}
    };

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

    fnames[eTRAJ1] = opt2fn("-f1", eNUMFILES, fnm);
    fnames[eTRAJ2] = opt2fn("-f2", eNUMFILES, fnm);
    fnames[eNDX1] = opt2fn_null("-n1", eNUMFILES, fnm);
    fnames[eNDX2] = opt2fn_null("-n2", eNUMFILES, fnm);
    fnames[eRES1] = opt2fn_null("-res", eNUMFILES, fnm);
    fnames[eETA_ATOM] = opt2fn("-eta_atom", eNUMFILES, fnm);
    fnames[eETA_RES] = opt2fn("-eta_res", eNUMFILES, fnm);

    ensemble_comp(fnames, gamma, c, &eta, &natoms, nthreads, &oenv);

    save_eta(eta, natoms, fnames[eETA_ATOM]);

    if(fnames[eRES1] != NULL) {
        eta_res_t eta_res;

        calc_eta_res(fnames[eRES1], eta, natoms, &eta_res);
        save_eta_res(&eta_res, fnames[eETA_RES]);

        free_eta_res(&eta_res);
    }

    sfree(eta);

    print_log("%s completed successfully.\n", argv[0]);
    close_log();

    return 0;
}

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 continously 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_morph(int argc,char *argv[])
{
  char *desc[] = {
    "g_morph does a linear interpolation of conformations in order to",
    "create intermediates. Of course these are completely unphysical, but",
    "that you may try to justify yourself. Output is in the form of a ",
    "generic trajectory. The number of intermediates can be controlled with",
    "the -ninterm flag. The first and last flag correspond to the way of",
    "interpolating: 0 corresponds to input structure 1 while",
    "1 corresponds to input strucutre 2.",
    "If you specify first < 0 or last > 1 extrapolation will be",
    "on the path from input structure x1 to x2. In general the coordinates",
    "of the intermediate x(i) out of N total intermidates correspond to:[PAR]",
    "x(i) = x1 + (first+(i/(N-1))*(last-first))*(x2-x1)[PAR]",
    "Finally the RMSD with respect to both input structures can be computed",
    "if explicitly selected (-or option). In that case an index file may be",
    "read to select what group RMS is computed from."
  };
  t_filenm fnm[] = {
    { efSTX, "-f1", "conf1",  ffREAD },
    { efSTX, "-f2", "conf2",  ffREAD },
    { efTRO, "-o",  "interm", ffWRITE },
    { efXVG, "-or", "rms-interm", ffOPTWR },
    { efNDX, "-n",  "index",  ffOPTRD }
  };
#define NFILE asize(fnm)
  static  int  ninterm = 11;
  static  real first   = 0.0;
  static  real last    = 1.0;
  static  bool bFit    = TRUE;
  t_pargs pa [] = {
    { "-ninterm", FALSE, etINT,  {&ninterm},
      "Number of intermediates" },
    { "-first",   FALSE, etREAL, {&first},
      "Corresponds to first generated structure (0 is input x0, see above)" },
    { "-last",    FALSE, etREAL, {&last},
      "Corresponds to last generated structure (1 is input x1, see above)" },
    { "-fit",     FALSE, etBOOL, {&bFit},
      "Do a least squares fit of the second to the first structure before interpolating" }
  };
  char     *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
  FILE     *fp=NULL;
  int      i,isize,is_lsq,status,nat1,nat2;
  atom_id  *index,*index_lsq,*index_all,*dummy;
  t_atoms  atoms;
  rvec     *x1,*x2,*xx,*v;
  matrix   box;
  real     rms1,rms2,fac,*mass;
  char     title[STRLEN],*grpname;
  bool     bRMS;
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW,
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,
		    0,NULL);
  get_stx_coordnum (opt2fn("-f1",NFILE,fnm),&nat1);
  get_stx_coordnum (opt2fn("-f2",NFILE,fnm),&nat2);
  if (nat1 != nat2)
    gmx_fatal(FARGS,"Number of atoms in first structure is %d, in second %d",
		nat1,nat2);
  
  init_t_atoms(&atoms,nat1,TRUE);
  snew(x1,nat1);
  snew(x2,nat1);
  snew(xx,nat1);
  snew(v,nat1);
  
  read_stx_conf(opt2fn("-f1",NFILE,fnm),title,&atoms,x1,v,NULL,box);
  read_stx_conf(opt2fn("-f2",NFILE,fnm),title,&atoms,x2,v,NULL,box);

  snew(mass,nat1);
  snew(index_all,nat1);
  for(i=0; (i<nat1); i++) {
    mass[i] = 1;
    index_all[i] = i;
  }
  if (bFit) {
    printf("Select group for LSQ superposition:\n");
    get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&is_lsq,&index_lsq,
	      &grpname);
    reset_x(is_lsq,index_lsq,nat1,index_all,x1,mass);
    reset_x(is_lsq,index_lsq,nat1,index_all,x2,mass);
    do_fit(nat1,mass,x1,x2);
  }
  
  bRMS = opt2bSet("-or",NFILE,fnm);
  if (bRMS) {
    fp = xvgropen(opt2fn("-or",NFILE,fnm),"RMSD","Conf","(nm)");
    xvgr_legend(fp,asize(leg),leg);
    printf("Select group for RMSD calculation:\n");
    get_index(&atoms,opt2fn_null("-n",NFILE,fnm),1,&isize,&index,&grpname);
    printf("You selected group %s, containing %d atoms\n",grpname,isize);
    rms1 = rmsdev_ind(isize,index,mass,x1,x2);  
    fprintf(stderr,"RMSD between input conformations is %g nm\n",rms1);
  }
  
  snew(dummy,nat1);
  for(i=0; (i<nat1); i++)
    dummy[i] = i;
  status = open_trx(ftp2fn(efTRX,NFILE,fnm),"w");
//.........这里部分代码省略.........

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.[BR]",
        ""
    };

    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	  */
        { efTPX, 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;
    output_env_t  oenv;

    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;
    }
    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(efTPX, NFILE, fnm);
    trxfnm = ftp2fn(efTRX, NFILE, fnm);

    /* Calculate axis */
    if (strcmp(normal_axis[0], "x") == 0)
    {
//.........这里部分代码省略.........