void comp_trx(const output_env_t oenv,const char *fn1, const char *fn2, 
              gmx_bool bRMSD,real ftol,real abstol)
{
  int i;
  const char *fn[2];
  t_trxframe fr[2];
  t_trxstatus *status[2];
  gmx_bool b[2];
  
  fn[0]=fn1;
  fn[1]=fn2;
  fprintf(stderr,"Comparing trajectory files %s and %s\n",fn1,fn2);
  for (i=0; i<2; i++)
    b[i] = read_first_frame(oenv,&status[i],fn[i],&fr[i],TRX_READ_X|TRX_READ_V|TRX_READ_F);
  
  if (b[0] && b[1]) { 
    do {
      comp_frame(stdout, &(fr[0]), &(fr[1]), bRMSD, ftol, abstol);
      
      for (i=0; i<2; i++)
	b[i] = read_next_frame(oenv,status[i],&fr[i]);
    } while (b[0] && b[1]);
    
    for (i=0; i<2; i++) {
      if (b[i] && !b[1-i])
	fprintf(stdout,"\nEnd of file on %s but not on %s\n",fn[i],fn[1-i]);
      close_trj(status[i]);
    }
  }
  if (!b[0] && !b[1])
    fprintf(stdout,"\nBoth files read correctly\n");
}

void write_bfactors(t_filenm  *fnm, int nfile, atom_id *index, atom_id *a, int nslices, int ngrps, real **order, t_topology *top, real **distvals, output_env_t oenv)
{
    /*function to write order parameters as B factors in PDB file using
          first frame of trajectory*/
    t_trxstatus *status;
    int          natoms;
    t_trxframe   fr, frout;
    t_atoms      useatoms;
    int          i, j, ctr, nout;

    ngrps -= 2;  /*we don't have an order parameter for the first or
                       last atom in each chain*/
    nout   = nslices*ngrps;
    natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr,
                              TRX_NEED_X);
    close_trj(status);
    frout        = fr;
    frout.natoms = nout;
    frout.bF     = FALSE;
    frout.bV     = FALSE;
    frout.x      = 0;
    snew(frout.x, nout);

    init_t_atoms(&useatoms, nout, TRUE);
    useatoms.nr = nout;

    /*initialize PDBinfo*/
    for (i = 0; i < useatoms.nr; ++i)
    {
        useatoms.pdbinfo[i].type         = 0;
        useatoms.pdbinfo[i].occup        = 0.0;
        useatoms.pdbinfo[i].bfac         = 0.0;
        useatoms.pdbinfo[i].bAnisotropic = FALSE;
    }

    for (j = 0, ctr = 0; j < nslices; j++)
    {
        for (i = 0; i < ngrps; i++, ctr++)
        {
            /*iterate along each chain*/
            useatoms.pdbinfo[ctr].bfac = order[j][i+1];
            if (distvals)
            {
                useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
            }
            copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
            useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
            useatoms.atom[ctr]     = top->atoms.atom[a[index[i+1]+j]];
            useatoms.nres          = max(useatoms.nres, useatoms.atom[ctr].resind+1);
            useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
        }
    }

    write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL, frout.ePBC, frout.box);

    sfree(frout.x);
    free_t_atoms(&useatoms, FALSE);
}

int read_first_v(const output_env_t oenv, t_trxstatus **status,const char *fn,
                 real *t, rvec **v,matrix box)
{
  t_trxframe fr;

  read_first_frame(oenv,status,fn,&fr,TRX_NEED_V);
  *t = fr.time;
  clear_v(&fr);
  *v = fr.v;
  copy_mat(fr.box,box);
  
  return fr.natoms;
}

int read_first_x(const gmx_output_env_t *oenv, t_trxstatus **status, const char *fn,
                 real *t, rvec **x, matrix box)
{
    t_trxframe fr;

    read_first_frame(oenv, status, fn, &fr, TRX_NEED_X);

    snew((*status)->xframe, 1);
    (*(*status)->xframe) = fr;
    *t                   = (*status)->xframe->time;
    *x                   = (*status)->xframe->x;
    copy_mat((*status)->xframe->box, box);

    return (*status)->xframe->natoms;
}

static void scan_trj_files(char **fnms,int nfiles,
			   real *readtime,real *timestep,atom_id imax)
{
  /* Check start time of all files */
  int i,status,natoms=0;
  real t;
  t_trxframe fr;
  bool ok;
  
  for(i=0;i<nfiles;i++) {
    ok=read_first_frame(&status,fnms[i],&fr,FLAGS);
    
    if(!ok) 
      gmx_fatal(FARGS,"\nCouldn't read frame from file.");
    if(fr.bTime)
      readtime[i]=fr.time;
    else {
      readtime[i]=0;
      fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
    }
    
    if(i==0) {
      natoms=fr.natoms;
    }
    else {
      if (imax==NO_ATID) {
	if(natoms!=fr.natoms) 
	  gmx_fatal(FARGS,"\nDifferent numbers of atoms (%d/%d) in files",
		      natoms,fr.natoms);
      } else {
	if(fr.natoms <= imax)
	  gmx_fatal(FARGS,"\nNot enough atoms (%d) for index group (%d)",
		      fr.natoms,imax);
      }
    }
    ok=read_next_frame(status,&fr);
    if(ok && fr.bTime) {
      timestep[i] = fr.time - readtime[i];
    } else {
      timestep[i] = 0;
    }
    
    close_trj(status);
  }
  fprintf(stderr,"\n");
  
  sfree(fr.x);
}

bool
TrajectoryFrameReader::readNextFrame()
{
    if (haveProbedForNextFrame_)
    {
        if (nextFrameExists_)
        {
            GMX_THROW(APIError("This frame has already been probed for, it should be used before probing again."));
        }
        else
        {
            GMX_THROW(APIError("This frame has already been probed for, it doesn't exist, so there should not be subsequent attempts to probe for it."));
        }
    }
    haveProbedForNextFrame_ = true;
    // If there's a next frame, read it into trxframe_, and report the result.
    if (!haveReadFirstFrame_)
    {
        t_trxstatus *trajectoryFile;
        int          flags = TRX_READ_X | TRX_READ_V | TRX_READ_F;
        nextFrameExists_ = read_first_frame(oenvGuard_.get(),
                                            &trajectoryFile,
                                            filename_.c_str(),
                                            trxframeGuard_.get(),
                                            flags);
        if (!trajectoryFile)
        {
            GMX_THROW(FileIOError("Could not open trajectory file " + filename_ + " for reading"));
        }
        trajectoryFileGuard_.reset(trajectoryFile);
        haveReadFirstFrame_ = true;
    }
    else
    {
        nextFrameExists_ = read_next_frame(oenvGuard_.get(),
                                           trajectoryFileGuard_.get(),
                                           trxframeGuard_.get());
    }
    return nextFrameExists_;
}

//.........这里部分代码省略.........
    else
    {
        nkc = NKC0;
    }
    nk  = kset_c[nkc];
    ntc = nk*NPK;

    sprintf(title, "Velocity Autocorrelation Function for %s", grpname);

    sysmass = 0;
    for (i = 0; i < nk; i++)
    {
        if (iprod(v0[i], v1[i]) != 0)
        {
            gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
        }
        if (iprod(v0[i], v2[i]) != 0)
        {
            gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
        }
        if (iprod(v1[i], v2[i]) != 0)
        {
            gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
        }
        unitv(v1[i], v1[i]);
        unitv(v2[i], v2[i]);
    }
    snew(tc, ntc);
    for (i = 0; i < top.atoms.nr; i++)
    {
        sysmass += top.atoms.atom[i].m;
    }

    read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr,
                     TRX_NEED_X | TRX_NEED_V);
    t0 = fr.time;

    n_alloc = 0;
    nframes = 0;
    rho     = 0;

    do
    {

        if (nframes >= n_alloc)
        {
            n_alloc += 100;
            for (i = 0; i < ntc; i++)
            {
                srenew(tc[i], n_alloc);
            }
        }

        rho += 1/det(fr.box);
        for (k = 0; k < nk; k++)
        {
            for (d = 0; d < DIM; d++)
            {
                kfac[k][d] = 2*M_PI*v0[k][d]/fr.box[d][d];
            }
        }
        for (i = 0; i < ntc; i++)
        {
            tc[i][nframes] = 0;
        }

//.........这里部分代码省略.........
    if (bDump)
    {
        printf("Dumping reference figures. Thanks for your patience.\n");
        dump_fy(oenv, toler);
        dump_w(oenv, beta);
        exit(0);
    }

    fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
    fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
    please_cite(fplog, "Pascal2011a");
    please_cite(fplog, "Caleman2011b");

    read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
                  TRUE);
    V     = det(box);
    tmass = 0;
    for (i = 0; (i < top.atoms.nr); i++)
    {
        tmass += top.atoms.atom[i].m;
    }

    Natom = top.atoms.nr;
    Nmol  = top.mols.nr;
    gnx   = Natom*DIM;

    /* 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;
    nframes = 0;
    Vsum    = V2sum = 0;
    nV      = 0;
    do
    {
        if (fr.bBox)
        {
            V      = det(fr.box);
            V2sum += V*V;
            Vsum  += V;
            nV++;
        }
        if (nframes >= n_alloc)
        {
            n_alloc += 100;
            for (i = 0; i < gnx; i++)
            {
                srenew(c1[i], n_alloc);
            }
        }
        for (i = 0; i < gnx; i += DIM)
        {
            c1[i+XX][nframes] = fr.v[i/DIM][XX];
            c1[i+YY][nframes] = fr.v[i/DIM][YY];
            c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
        }

        t1 = fr.time;

//.........这里部分代码省略.........
        gmx_fatal(FARGS, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
    }

    printf("Select group with donor atom pairs defining the transition moment\n");
    get_index(NULL, ftp2fn_null(efNDX, NFILE, fnm), 1, &ndon, &donindex, &grpnm);

    printf("Select group with acceptor atom pairs defining the transition moment\n");
    get_index(NULL, ftp2fn_null(efNDX, NFILE, fnm), 1, &nacc, &accindex, &grpnm);

    /*check if groups are identical*/
    grident = TRUE;

    if (ndon == nacc)
    {
        for (i = 0; i < nacc; i++)
        {
            if (accindex[i] != donindex[i])
            {
                grident = FALSE;
                break;
            }
        }
    }

    if (grident)
    {
        gmx_fatal(FARGS, "Donor and acceptor group are identical. This makes no sense.");
    }

    printf("Reading first frame\n");
    /* open trx file for reading */
    flags           = 0;
    flags           = flags | TRX_READ_X;
    bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);

    if (bHaveFirstFrame)
    {
        printf("First frame is OK\n");
        natoms = fr.natoms;
        if ((ndon % 2 != 0) || (nacc % 2 != 0))
        {
            indexOK = FALSE;
        }
        else
        {
            for (i = 0; i < ndon; i++)
            {
                if (donindex[i] >= natoms)
                {
                    indexOK = FALSE;
                }
            }
            for (i = 0; i < nacc; i++)
            {
                if (accindex[i] >= natoms)
                {
                    indexOK = FALSE;
                }
            }
        }

        if (indexOK)
        {

            if (bDatout)
            {

//.........这里部分代码省略.........
        outekt    = xvgropen(opt2fn("-ekt", NFILE, fnm), "Center of mass translation",
                             output_env_get_xvgr_tlabel(oenv), "Energy (kJ mol\\S-1\\N)", oenv);
        make_legend(outekt, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
    }
    if (bEKR)
    {
        bDum[XX]  = FALSE;
        bDum[YY]  = FALSE;
        bDum[ZZ]  = FALSE;
        bDum[DIM] = TRUE;
        flags     = flags | TRX_READ_X | TRX_READ_V;
        outekr    = xvgropen(opt2fn("-ekr", NFILE, fnm), "Center of mass rotation",
                             output_env_get_xvgr_tlabel(oenv), "Energy (kJ mol\\S-1\\N)", oenv);
        make_legend(outekr, ngroups, isize[0], index[0], grpname, bCom, bMol, bDum, oenv);
    }
    if (bVD)
    {
        flags = flags | TRX_READ_V;
    }
    if (bCV)
    {
        flags = flags | TRX_READ_X | TRX_READ_V;
    }
    if (bCF)
    {
        flags = flags | TRX_READ_X | TRX_READ_F;
    }
    if ((flags == 0) && !bOB)
    {
        fprintf(stderr, "Please select one or more output file options\n");
        exit(0);
    }

    read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);


    if ((bOV || bOF) && fn2ftp(ftp2fn(efTRX, NFILE, fnm)) == efXTC)
    {
        gmx_fatal(FARGS, "Cannot extract velocities or forces since your input XTC file does not contain them.");
    }

    if (bCV || bCF)
    {
        snew(sumx, fr.natoms);
    }
    if (bCV)
    {
        snew(sumv, fr.natoms);
    }
    if (bCF)
    {
        snew(sumf, fr.natoms);
    }
    nr_xfr = 0;
    nr_vfr = 0;
    nr_ffr = 0;

    if (bCom && bPBC)
    {
        gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
    }

    do
    {
        time = output_env_conv_time(oenv, fr.time);

int main(int argc,char *argv[])
{
  static char *desc[] = {
    "[TT]g_anavel[tt] computes temperature profiles in a sample. The sample",
    "can be analysed radial, i.e. the temperature as a function of",
    "distance from the center, cylindrical, i.e. as a function of distance",
    "from the vector (0,0,1) through the center of the box, or otherwise",
    "(will be specified later)"
  };
  t_filenm fnm[] = {
    { efTRN,  "-f",  NULL, ffREAD },
    { efTPX,  "-s",  NULL, ffREAD },
    { efXPM,  "-o", "xcm", ffWRITE }
  };
#define NFILE asize(fnm)

  static int  mode = 0,   nlevels = 10;
  static real tmax = 300, xmax    = -1;
  t_pargs pa[] = {
    { "-mode",    FALSE, etINT,  {&mode},    "mode" },
    { "-nlevels", FALSE, etINT,  {&nlevels}, "number of levels" },
    { "-tmax",    FALSE, etREAL, {&tmax},    "max temperature in output" },
    { "-xmax",    FALSE, etREAL, {&xmax},    "max distance from center" }
  };
  
  FILE       *fp;
  int        *npts,nmax;
  int        status;
  int        i,j,idum,step,nframe=0,index;
  real       temp,rdum,hboxx,hboxy,scale,xnorm=0;
  real       **profile=NULL;
  real       *t_x=NULL,*t_y,hi=0;
  t_topology *top;
  int        d,m,n;
  matrix     box;
  atom_id    *sysindex;
  gmx_bool       bHaveV,bReadV;
  t_rgb      rgblo = { 0, 0, 1 },rgbhi = { 1, 0, 0 };
  int        flags = TRX_READ_X | TRX_READ_V;
  t_trxframe fr;

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

  top    = read_top(ftp2fn(efTPX,NFILE,fnm));

  read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,flags);
	
  if (xmax > 0) {
    scale  = 5;
    nmax   = xmax*scale;
  }
  else {
    scale  = 5;
    nmax   = (0.5*sqrt(sqr(box[XX][XX])+sqr(box[YY][YY])))*scale; 
  }
  snew(npts,nmax+1);
  snew(t_y,nmax+1);
  for(i=0; (i<=nmax); i++) {
    npts[i] = 0;
    t_y[i]  = i/scale;
  }
  do {
    srenew(profile,++nframe);
    snew(profile[nframe-1],nmax+1);
    srenew(t_x,nframe);
    t_x[nframe-1] = fr.time*1000;
    hboxx = box[XX][XX]/2;
    hboxy = box[YY][YY]/2;
    for(i=0; (i<fr.natoms); i++) {
      /* determine position dependent on mode */
      switch (mode) {
      case 0:
	xnorm = sqrt(sqr(fr.x[i][XX]-hboxx) + sqr(fr.x[i][YY]-hboxy));
	break;
      default:
	gmx_fatal(FARGS,"Unknown mode %d",mode);
      }
      index = xnorm*scale;
      if (index <= nmax) {
	temp = top->atoms.atom[i].m*iprod(fr.v[i],fr.v[i])/(2*BOLTZ);
	if (temp > hi)
	  hi = temp;
	npts[index]++;
	profile[nframe-1][index] += temp;
      }
    }
    for(i=0; (i<=nmax); i++) {
      if (npts[i] != 0) 
	profile[nframe-1][i] /= npts[i];
      npts[i] = 0;
    }
  } while (read_next_frame(status,&fr));
  close_trx(status);

  fp = ftp2FILE(efXPM,NFILE,fnm,"w");
  write_xpm(fp,0,"Temp. profile","T (a.u.)",
	    "t (fs)","R (nm)",
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
        bTPS = ftp2bSet(efTPS, NFILE, fnm) || !ftp2bSet(efNDX, NFILE, fnm);
    }

    if (bTPS)
    {
        bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &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;
        if (bMol)
        {
            for (i = 0; i < gnx; i++)
            {
                clear_rvec(mv_mol);
                k = top.mols.index[index[i]];
                l = top.mols.index[index[i]+1];
                for (j = k; j < l; j++)
                {
                    if (bMass)
                    {
                        mass = top.atoms.atom[j].m;
                    }
                    else
                    {
                        mass = normm[j];
                    }

//.........这里部分代码省略.........
     */
    out_file = fnms_out[0];
    n_append = -1;
    for(i=0; ((i<nfile_in) && (n_append==-1)); i++) {
      if (strcmp(fnms[i],out_file) == 0) {
	n_append = i;
      }
    }
    if (n_append == 0)
      fprintf(stderr,"Will append to %s rather than creating a new file\n",
	      out_file);
    else if (n_append != -1)
      gmx_fatal(FARGS,"Can only append to the first file which is %s (not %s)",
		fnms[0],out_file);
    
    earliersteps=0;    
    
    /* Not checking input format, could be dangerous :-) */
    /* Not checking output format, equally dangerous :-) */
    
    frame=-1;
    frame_out=-1;
    /* the default is not to change the time at all,
     * but this is overridden by the edit_files routine
     */
    t_corr=0;
    
    if (n_append == -1) {
      trxout = open_trx(out_file,"w");
      memset(&frout,0,sizeof(frout));
    }
    else {
      /* Read file to find what is the last frame in it */
      if (!read_first_frame(&status,out_file,&fr,FLAGS))
	gmx_fatal(FARGS,"Reading first frame from %s",out_file);
      while (read_next_frame(status,&fr))
	;
      close_trj(status);
      lasttime = fr.time;
      bKeepLast = TRUE;
      trxout = open_trx(out_file,"a");
      frout = fr;
    }
    /* Lets stitch up some files */
    timestep = timest[0];
    for(i=n_append+1; (i<nfile_in); i++) {
      /* Open next file */
      
      /* set the next time from the last frame in previous file */
      if (i > 0) {
	if (frame_out >= 0) {
	  if(cont_type[i]==TIME_CONTINUE) {
	    begin =frout.time;
	    begin += 0.5*timestep;
	    settime[i]=frout.time;
	    cont_type[i]=TIME_EXPLICIT;	  
	  }
	  else if(cont_type[i]==TIME_LAST) {
	    begin=frout.time;
	    begin += 0.5*timestep;
	  }
	  /* Or, if the time in the next part should be changed by the
	   * same amount, start at half a timestep from the last time
	   * so we dont repeat frames.
	   */
	  /* I don't understand the comment above, but for all the cases

static void scan_trj_files(char **fnms, int nfiles, real *readtime,
                           real *timestep, int imax,
                           const gmx_output_env_t *oenv)
{
    /* Check start time of all files */
    int          i, natoms = 0;
    t_trxstatus *status;
    t_trxframe   fr;
    gmx_bool     ok;

    for (i = 0; i < nfiles; i++)
    {
        ok = read_first_frame(oenv, &status, fnms[i], &fr, FLAGS);

        if (!ok)
        {
            gmx_fatal(FARGS, "\nCouldn't read frame from file." );
        }
        if (fr.bTime)
        {
            readtime[i] = fr.time;
        }
        else
        {
            readtime[i] = 0;
            fprintf(stderr, "\nWARNING: Couldn't find a time in the frame.\n");
        }

        if (i == 0)
        {
            natoms = fr.natoms;
        }
        else
        {
            if (imax == -1)
            {
                if (natoms != fr.natoms)
                {
                    gmx_fatal(FARGS, "\nDifferent numbers of atoms (%d/%d) in files",
                              natoms, fr.natoms);
                }
            }
            else
            {
                if (fr.natoms <= imax)
                {
                    gmx_fatal(FARGS, "\nNot enough atoms (%d) for index group (%d)",
                              fr.natoms, imax);
                }
            }
        }
        ok = read_next_frame(oenv, status, &fr);
        if (ok && fr.bTime)
        {
            timestep[i] = fr.time - readtime[i];
        }
        else
        {
            timestep[i] = 0;
        }

        close_trj(status);
        if (fr.bX)
        {
            sfree(fr.x);
        }
        if (fr.bV)
        {
            sfree(fr.v);
        }
        if (fr.bF)
        {
            sfree(fr.f);
        }
    }
    fprintf(stderr, "\n");

}

void
TrajectoryAnalysisRunnerCommon::initFirstFrame()
{
    // Return if we have already initialized the trajectory.
    if (impl_->fr)
    {
        return;
    }
    time_unit_t time_unit
        = static_cast<time_unit_t>(impl_->settings_.timeUnit() + 1);
    output_env_init(&impl_->oenv_, 0, NULL, time_unit, FALSE, exvgNONE, 0, 0);

    int frflags = impl_->settings_.frflags();
    frflags |= TRX_NEED_X;

    snew(impl_->fr, 1);

    const TopologyInformation &top = impl_->topInfo_;
    if (hasTrajectory())
    {
        if (!read_first_frame(impl_->oenv_, &impl_->status_,
                              impl_->trjfile_.c_str(), impl_->fr, frflags))
        {
            GMX_THROW(FileIOError("Could not read coordinates from trajectory"));
        }
        impl_->bTrajOpen_ = true;

        if (top.hasTopology() && impl_->fr->natoms > top.topology()->atoms.nr)
        {
            GMX_THROW(InconsistentInputError(formatString(
                                                     "Trajectory (%d atoms) does not match topology (%d atoms)",
                                                     impl_->fr->natoms, top.topology()->atoms.nr)));
        }
        // Check index groups if they have been initialized based on the topology.
        /*
           if (top)
           {
            for (int i = 0; i < impl_->sel->nr(); ++i)
            {
                gmx_ana_index_check(impl_->sel->sel(i)->indexGroup(),
                                    impl_->fr->natoms);
            }
           }
         */
    }
    else
    {
        // Prepare a frame from topology information.
        // TODO: Initialize more of the fields.
        if (frflags & (TRX_NEED_V))
        {
            GMX_THROW(NotImplementedError("Velocity reading from a topology not implemented"));
        }
        if (frflags & (TRX_NEED_F))
        {
            GMX_THROW(InvalidInputError("Forces cannot be read from a topology"));
        }
        impl_->fr->flags  = frflags;
        impl_->fr->natoms = top.topology()->atoms.nr;
        impl_->fr->bX     = TRUE;
        snew(impl_->fr->x, impl_->fr->natoms);
        memcpy(impl_->fr->x, top.xtop_,
               sizeof(*impl_->fr->x) * impl_->fr->natoms);
        impl_->fr->bBox   = TRUE;
        copy_mat(const_cast<rvec *>(top.boxtop_), impl_->fr->box);
    }

    set_trxframe_ePBC(impl_->fr, top.ePBC());
    if (top.hasTopology() && impl_->settings_.hasRmPBC())
    {
        impl_->gpbc_ = gmx_rmpbc_init(&top.topology()->idef, top.ePBC(),
                                      impl_->fr->natoms, impl_->fr->box);
    }
}

int main(int argn, char* args[])
{
	static char *desc[] = {
		"This script will read the xtc trajectory from an MD simulation",
		"performed by GROMACS and output all data related to the analysis",
		"of preferential interactions.",
		"[PAR]",
		"Generated: Thu Aug 20 16:49:21 EDT 2015.[BR]",
		"Last updated: Thu Aug 27 17:26:01 EDT 2015.",
		"[PAR]",
		"Future Developments:[BR]",
		"[BB]1.[bb] Partition calculation into groups (constituent amino acids).[BR]",
		"[BB]2.[bb] Calculate residence times.[BR]",
		"[BB]3.[bb] Include a measure for geometric orientation.[BR]",
		"[BB]4.[bb] Calculate molecular distributions.[BR]",
		"[BB]5.[bb] Parallelize it.",
		"[PAR]",
		"Remember to check the dependencie of the results on the number of block",
		" used by using [TT]-block[tt].",
		"Also, check that your simulation has been properly equilibrated by using [TT]-f0[tt].",
		"In terms of reducing the computational cost of the analysis, use the option",
		"[TT]-skip[tt] to only use every nr[TT]-th[tt] frame from the trajectory."
	};
	
	unsigned int protein_sequence=20; 	/* parse tpr.itp for info 			*/
        unsigned int N_co=4; 			/* look at topology for number of co-solvents 	*/
        unsigned int num_sol=7906; 		/* look at topology as well 			*/
	unsigned int skip_nr=1;
	unsigned int granularity=2*100;
	unsigned int f0=0;
	bool bIndex=false;

	t_pargs pa[] = {
		{"-start", FALSE, etINT, {&f0}, "Start after n-th frame."},
		{"-skip", FALSE, etINT, {&skip_nr}, "Only write every nr-th frame."},
		{"-grid", FALSE, etINT, {&granularity}, "Number points to calculate."},
		{"-seq", FALSE, etINT,{&protein_sequence}, "Number of amino acids in protein sequence."},
		{"-Nco", FALSE, etINT,{&N_co}, "Number of cosolvent molecules."},
		{"-Sol", FALSE, etINT,{&num_sol}, "Number of water molecules."}
	};	

	char 		title[STRLEN];	
	t_topology 	top;
	int		ePBC;
	rvec		*xtop;
	matrix		box;
	
	t_filenm fnm[] = { 			/* See filenm.h */
		{ efTOP, NULL,  NULL, ffREAD  },
		{ efTPS, NULL, NULL, ffREAD },	/*  topology 	*/
		{ efTRX, "-f", NULL, ffREAD },	/*  trajcetory  */
		{ efNDX, "-n", NULL, ffOPTRD }  /*  index file  */
	};

#define NFILE asize(fnm)

/* Interface. Adds default options. */	
CopyRight(stderr,args[0]);
parse_common_args(&argn,args,PCA_CAN_TIME | PCA_CAN_VIEW, NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
read_tps_conf( ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box, TRUE);
sfree(xtop);

/*******************************************************************************************************************/
/*
int iter;
for (iter=0;iter<top.atoms.nr/100;++iter)
{
printf("Atom name: %s\tAtom charge: %f\n", *(top.atoms.atomname[iter]), top.atoms.atom[iter].q );
printf("Atom type: %d\tAtomic Residue NUmber: %d\n", top.atoms.atom[iter].type, top.atoms.atom[iter].resnr );
printf("Chain Identifier: %u\tNr of residues names: %d\n", top.atoms.atom[iter].ptype, top.atoms.nres );
printf("Residue name: %s\n\n", *(top.atoms.resname[iter]) );
}
*/
/*******************************************************************************************************************/

int status;
t_trxframe fr;
int flags = TRX_READ_X;

/* Count Number of Frames */
        unsigned int bframes=0, frames=0;	
	int bwrite;
	read_first_frame(&status, ftp2fn(efTRX,NFILE,fnm), &fr, flags);
        do {
		bwrite = bframes % skip_nr; 
		++bframes; 
		if ((bframes>f0) && (bwrite==0)){ ++frames; }
	} while ( read_next_frame(status,&fr) );

/* Preparing Block Analysis */
	bIndex = ftp2bSet(efNDX,NFILE,fnm);
	std::vector<int> blocks_list;
	if (bIndex)
	{
		FILE *fp = ffopen(opt2fn("-n",NFILE,fnm),"r");
		int numblock;
		while(fscanf(fp,"%d",&numblock)!=EOF)
		{
			blocks_list.push_back(numblock);
		}
//.........这里部分代码省略.........

extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX,
                                    const char* fnXVG, const char *fnTRX,
                                    const char* fnDAT,
                                    real start_q, real end_q,
                                    real energy, int ng, const output_env_t oenv)
{
    int                     i, *isize, flags = TRX_READ_X, **index_atp;
    t_trxstatus            *status;
    char                  **grpname, title[STRLEN];
    atom_id               **index;
    t_topology              top;
    int                     ePBC;
    t_trxframe              fr;
    reduced_atom_t        **red;
    structure_factor       *sf;
    rvec                   *xtop;
    real                  **sf_table;
    int                     nsftable;
    matrix                  box;
    double                  r_tmp;

    gmx_structurefactors_t *gmx_sf;
    real                   *a, *b, c;
    int                     success;

    snew(a, 4);
    snew(b, 4);


    gmx_sf = gmx_structurefactors_init(fnDAT);

    success = gmx_structurefactors_get_sf(gmx_sf, 0, a, b, &c);

    snew (sf, 1);
    sf->energy = energy;

    /* Read the topology informations */
    read_tps_conf (fnTPS, title, &top, &ePBC, &xtop, NULL, box, TRUE);
    sfree (xtop);

    /* groups stuff... */
    snew (isize, ng);
    snew (index, ng);
    snew (grpname, ng);

    fprintf (stderr, "\nSelect %d group%s\n", ng,
             ng == 1 ? "" : "s");
    if (fnTPS)
    {
        get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
    }
    else
    {
        rd_index (fnNDX, ng, isize, index, grpname);
    }

    /* The first time we read data is a little special */
    read_first_frame (oenv, &status, fnTRX, &fr, flags);

    sf->total_n_atoms = fr.natoms;

    snew (red, ng);
    snew (index_atp, ng);

    r_tmp = max (box[XX][XX], box[YY][YY]);
    r_tmp = (double) max (box[ZZ][ZZ], r_tmp);

    sf->ref_k = (2.0 * M_PI) / (r_tmp);
    /* ref_k will be the reference momentum unit */
    sf->n_angles = (int) (end_q / sf->ref_k + 0.5);

    snew (sf->F, ng);
    for (i = 0; i < ng; i++)
    {
        snew (sf->F[i], sf->n_angles);
    }
    for (i = 0; i < ng; i++)
    {
        snew (red[i], isize[i]);
        rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE, gmx_sf);
        index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
    }

    sf_table = compute_scattering_factor_table (gmx_sf, (structure_factor_t *)sf);


    /* This is the main loop over frames */

    do
    {
        sf->nSteps++;
        for (i = 0; i < ng; i++)
        {
            rearrange_atoms (red[i], &fr, index[i], isize[i], &top, FALSE, gmx_sf);

            compute_structure_factor ((structure_factor_t *)sf, box, red[i], isize[i],
                                      start_q, end_q, i, sf_table);
        }
    }

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

static void clust_size(const char *ndx, const char *trx, const char *xpm,
                       const char *xpmw, const char *ncl, const char *acl,
                       const char *mcl, const char *histo, const char *tempf,
                       const char *mcn, gmx_bool bMol, gmx_bool bPBC, const char *tpr,
                       real cut, int nskip, int nlevels,
                       t_rgb rmid, t_rgb rhi, int ndf,
                       const output_env_t oenv)
{
    FILE                 *fp, *gp, *hp, *tp;
    atom_id              *index = NULL;
    int                   nindex, natoms;
    t_trxstatus          *status;
    rvec                 *x = NULL, *v = NULL, dx;
    t_pbc                 pbc;
    char                 *gname;
    char                  timebuf[32];
    gmx_bool              bSame, bTPRwarn = TRUE;
    /* Topology stuff */
    t_trxframe            fr;
    t_tpxheader           tpxh;
    gmx_mtop_t           *mtop = NULL;
    int                   ePBC = -1;
    t_block              *mols = NULL;
    gmx_mtop_atomlookup_t alook;
    t_atom               *atom;
    int                   version, generation, ii, jj;
    real                  temp, tfac;
    /* Cluster size distribution (matrix) */
    real                **cs_dist = NULL;
    real                  tf, dx2, cut2, *t_x = NULL, *t_y, cmid, cmax, cav, ekin;
    int                   i, j, k, ai, aj, ci, cj, nframe, nclust, n_x, max_size = 0;
    int                  *clust_index, *clust_size, max_clust_size, max_clust_ind, nav, nhisto;
    t_rgb                 rlo = { 1.0, 1.0, 1.0 };

    clear_trxframe(&fr, TRUE);
    sprintf(timebuf, "Time (%s)", output_env_get_time_unit(oenv));
    tf     = output_env_get_time_factor(oenv);
    fp     = xvgropen(ncl, "Number of clusters", timebuf, "N", oenv);
    gp     = xvgropen(acl, "Average cluster size", timebuf, "#molecules", oenv);
    hp     = xvgropen(mcl, "Max cluster size", timebuf, "#molecules", oenv);
    tp     = xvgropen(tempf, "Temperature of largest cluster", timebuf, "T (K)",
                      oenv);

    if (!read_first_frame(oenv, &status, trx, &fr, TRX_NEED_X | TRX_READ_V))
    {
        gmx_file(trx);
    }

    natoms = fr.natoms;
    x      = fr.x;

    if (tpr)
    {
        snew(mtop, 1);
        read_tpxheader(tpr, &tpxh, TRUE, &version, &generation);
        if (tpxh.natoms != natoms)
        {
            gmx_fatal(FARGS, "tpr (%d atoms) and trajectory (%d atoms) do not match!",
                      tpxh.natoms, natoms);
        }
        ePBC = read_tpx(tpr, NULL, NULL, &natoms, NULL, NULL, mtop);
    }
    if (ndf <= -1)
    {
        tfac = 1;
    }
    else
    {
        tfac = ndf/(3.0*natoms);
    }

    if (bMol)
    {
        if (ndx)
        {
            printf("Using molecules rather than atoms. Not reading index file %s\n",
                   ndx);
        }
        GMX_RELEASE_ASSERT(mtop != NULL, "Trying to access mtop->mols from NULL mtop pointer");
        mols = &(mtop->mols);

        /* Make dummy index */
        nindex = mols->nr;
        snew(index, nindex);
        for (i = 0; (i < nindex); i++)
        {
            index[i] = i;
        }
        gname = gmx_strdup("mols");
    }
    else
    {
        rd_index(ndx, 1, &nindex, &index, &gname);
    }

    alook = gmx_mtop_atomlookup_init(mtop);

    snew(clust_index, nindex);
    snew(clust_size, nindex);
    cut2   = cut*cut;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                sprintf(str, "f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
                        *(groups->grpname[groups->grps[egcENER].nm_ind[i]]));
                leg[e++] = strdup(str);
            }
            if (bRFExcl)
            {
                sprintf(str, "f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
                leg[e++] = strdup(str);
            }
            if (EEL_FULL(fr->eeltype))
            {
                sprintf(str, "f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
                leg[e++] = strdup(str);
            }
        }
        xvgr_legend(fp_tpi, 4+nener, (const char**)leg, oenv);
        for (i = 0; i < 4+nener; i++)
        {
            sfree(leg[i]);
        }
        sfree(leg);
    }
    clear_rvec(x_init);
    V_all     = 0;
    VembU_all = 0;

    invbinw = 10;
    nbin    = 10;
    snew(bin, nbin);

    /* Avoid frame step numbers <= -1 */
    frame_step_prev = -1;

    bNotLastFrame = read_first_frame(oenv, &status, opt2fn("-rerun", nfile, fnm),
                                     &rerun_fr, TRX_NEED_X);
    frame = 0;

    if (rerun_fr.natoms - (bCavity ? nat_cavity : 0) !=
        mdatoms->nr - (a_tp1 - a_tp0))
    {
        gmx_fatal(FARGS, "Number of atoms in trajectory (%d)%s "
                  "is not equal the number in the run input file (%d) "
                  "minus the number of atoms to insert (%d)\n",
                  rerun_fr.natoms, bCavity ? " minus one" : "",
                  mdatoms->nr, a_tp1-a_tp0);
    }

    refvolshift = log(det(rerun_fr.box));

    switch (inputrec->eI)
    {