bool new_data(t_xrama *xr)
{
  if (!read_next_x(xr->traj,&xr->t,xr->natoms,xr->x,xr->box))
    return FALSE;

  calc_dihs(xr);

  return TRUE;
}

gmx_bool new_data(t_xrama *xr)
{
    if (!read_next_x(xr->oenv, xr->traj, &xr->t, xr->x, xr->box))
    {
        return FALSE;
    }

    calc_dihs(xr);

    return TRUE;
}

static bool step_man(t_manager *man, int *nat)
{
    static int      ncount = 0;
    static bool     bWarn  = false;
    bool            bEof;
    const char     *warn;

    if (!man->natom)
    {
        fprintf(stderr, "Not initiated yet!");
        exit(1);
    }
    bEof = read_next_x(man->oenv, man->status, &man->time, man->x, man->box);
    *nat = man->natom;
    if (ncount == man->nSkip)
    {
        switch (man->molw->boxtype)
        {
            case esbTri:
                put_atoms_in_triclinic_unitcell(ecenterDEF, man->box, man->natom, man->x);
                break;
            case esbTrunc:
                warn = put_atoms_in_compact_unitcell(man->molw->ePBC, ecenterDEF, man->box,
                                                     man->natom, man->x);
                if (warn && !bWarn)
                {
                    fprintf(stderr, "\n%s\n", warn);
                    bWarn = true;
                }
                break;
            case esbRect:
            case esbNone:
            default:
                break;
        }
        if (man->bPbc)
        {
            gmx_rmpbc(man->gpbc, man->natom, man->box, man->x);
            reset_mols(&(man->top.mols), man->box, man->x);
        }
        ncount = 0;
    }
    else
    {
        if (man->nSkip > 0)
        {
            ncount++;
            return step_man(man, nat);
        }
    }

    return bEof;
}

static bool step_man(t_manager *man, int *nat)
{
    static int      ncount = 0;
    bool            bEof;

    if (!man->natom)
    {
        std::fprintf(stderr, "Not initiated yet!");
        std::exit(1);
    }
    bEof = read_next_x(man->oenv, man->status, &man->time, man->x, man->box);
    *nat = man->natom;
    if (ncount == man->nSkip)
    {
        auto atomsArrayRef = gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(man->x), man->natom);
        switch (man->molw->boxtype)
        {
            case esbTri:
                put_atoms_in_triclinic_unitcell(ecenterDEF, man->box, atomsArrayRef);
                break;
            case esbTrunc:
                put_atoms_in_compact_unitcell(man->molw->ePBC, ecenterDEF, man->box,
                                              atomsArrayRef);
                break;
            case esbRect:
            case esbNone:
            default:
                break;
        }
        if (man->bPbc)
        {
            gmx_rmpbc(man->gpbc, man->natom, man->box, man->x);
            reset_mols(&(man->top.mols), man->box, man->x);
        }
        ncount = 0;
    }
    else
    {
        if (man->nSkip > 0)
        {
            ncount++;
            return step_man(man, nat);
        }
    }

    return bEof;
}

//.........这里部分代码省略.........
    { efNDX, NULL, NULL,  ffREAD  },
    { efXVG, "-o", "rotacf",  ffWRITE }
  };
#define NFILE asize(fnm)
  int     npargs;
  t_pargs *ppa;
  
  CopyRight(stderr,argv[0]);
  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);
  
  rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
  
  if (bVec) 
    nvec = isize/2;
  else
    nvec = isize/3;
  
  if (((isize % 3) != 0) && !bVec)
    gmx_fatal(FARGS,"number of index elements not multiple of 3, "
		"these can not be atom triplets\n");
  if (((isize % 2) != 0) && bVec)
    gmx_fatal(FARGS,"number of index elements not multiple of 2, "
		"these can not be atom doublets\n");
  
  top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
  
  snew(c1,nvec);
  for (i=0; (i<nvec); i++)
    c1[i]=NULL;
  n_alloc=0;

  natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
  snew(x_s,natoms);
  
  /* Start the loop over frames */
  t1 = t0 = t;
  teller  = 0;
  do {
    if (teller >= n_alloc) {
      n_alloc+=100;
      for (i=0; (i<nvec); i++)
	srenew(c1[i],DIM*n_alloc);
    }
    t1 = t;
    
    /* Remove periodicity */
    rm_pbc(&(top->idef),ePBC,natoms,box,x,x_s);
    
    /* Compute crossproducts for all vectors, if triplets.
     * else, just get the vectors in case of doublets.
     */
    if (bVec == FALSE) {
      for (i=0; (i<nvec); i++) {
	ai=index[3*i];
	aj=index[3*i+1];
	ak=index[3*i+2];
	rvec_sub(x_s[ai],x_s[aj],xij);
	rvec_sub(x_s[aj],x_s[ak],xjk);
	cprod(xij,xjk,n);
	for(m=0; (m<DIM); m++)
	  c1[i][DIM*teller+m]=n[m];
      }
    }
    else {
      for (i=0; (i<nvec); i++) {
	ai=index[2*i];
	aj=index[2*i+1];
	rvec_sub(x_s[ai],x_s[aj],n);
	for(m=0; (m<DIM); m++)
	  c1[i][DIM*teller+m]=n[m];
      }
    }
    /* Increment loop counter */
    teller++;
  } while (read_next_x(status,&t,natoms,x,box));  
  close_trj(status); 
  fprintf(stderr,"\nDone with trajectory\n");
  
  /* Autocorrelation function */
  if (teller < 2)
    fprintf(stderr,"Not enough frames for correlation function\n");
  else {
    dt=(t1 - t0)/(teller-1);
    
    mode = eacVector;
    
    do_autocorr(ftp2fn(efXVG,NFILE,fnm),"Rotational Correlation Function",
		teller,nvec,c1,dt,mode,bAver);
  }

  do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
    
  thanx(stderr);
    
  return 0;
}

//.........这里部分代码省略.........
                {
                    real dd = angles[cur][i];
                    angles[cur][i] = atan2(sin(dd), cos(dd));
                }
            }
            else
            {
                if (teller > 1)
                {
                    for (i = 0; (i < nangles); i++)
                    {
                        while (angles[cur][i] <= angles[prev][i] - M_PI)
                        {
                            angles[cur][i] += 2*M_PI;
                        }
                        while (angles[cur][i] > angles[prev][i] + M_PI)
                        {
                            angles[cur][i] -= 2*M_PI;
                        }
                    }
                }
            }
        }

        /* Average angles */
        aa = 0;
        for (i = 0; (i < nangles); i++)
        {
            aa = aa+angles[cur][i];

            /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
               even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
               (angle) Basically: translate the x-axis by Pi. Translate it back by
               -Pi when plotting.
             */

            angle = angles[cur][i];
            if (!bAngles)
            {
                while (angle < -M_PI)
                {
                    angle += 2*M_PI;
                }
                while (angle >= M_PI)
                {
                    angle -= 2*M_PI;
                }

                angle += M_PI;
            }

            /* Update the distribution histogram */
            angind = (int) ((angle*maxangstat)/pifac + 0.5);
            if (angind == maxangstat)
            {
                angind = 0;
            }
            if ( (angind < 0) || (angind >= maxangstat) )
            {
                /* this will never happen */
                gmx_fatal(FARGS, "angle (%f) index out of range (0..%d) : %d\n",
                          angle, maxangstat, angind);
            }

            angstat[angind]++;
            if (angind == maxangstat)
            {
                fprintf(stderr, "angle %d fr %d = %g\n", i, cur, angle);
            }

            total++;
        }

        /* average over all angles */
        (*aver_angle)[teller] = (aa/nangles);

        /* this copies all current dih. angles to dih[i], teller is frame */
        if (bSaveAll)
        {
            for (i = 0; i < nangles; i++)
            {
                dih[i][teller] = angles[cur][i];
            }
        }

        /* Swap buffers */
        cur = prev;

        /* Increment loop counter */
        teller++;
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);

    sfree(x);
    sfree(angles[cur]);
    sfree(angles[prev]);

    *nframes = teller;
}

//.........这里部分代码省略.........
    trxnat = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);

    nframes = 0;

    rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;
    rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;

    gpbc = gmx_rmpbc_init(&top.idef, ePBC, trxnat, box);

    if (bFrames)
    {
        out = opt2FILE("-frames", NFILE, fnm, "w");
    }
    do
    {
        gmx_rmpbc(gpbc, trxnat, box, x);
        nframes++;
        calc_mat(nres, natoms, rndx, x, index, truncate, mdmat, nmat, ePBC, box);
        for (i = 0; (i < nres); i++)
        {
            for (j = 0; (j < natoms); j++)
            {
                if (nmat[i][j])
                {
                    totnmat[i][j]++;
                }
            }
        }
        for (i = 0; (i < nres); i++)
        {
            for (j = 0; (j < nres); j++)
            {
                totmdmat[i][j] += mdmat[i][j];
            }
        }
        if (bFrames)
        {
            sprintf(label, "t=%.0f ps", t);
            write_xpm(out, 0, label, "Distance (nm)", "Residue Index", "Residue Index",
                      nres, nres, resnr, resnr, mdmat, 0, truncate, rlo, rhi, &nlevels);
        }
    }
    while (read_next_x(oenv, status, &t, trxnat, x, box));
    fprintf(stderr, "\n");
    close_trj(status);
    gmx_rmpbc_done(gpbc);
    if (bFrames)
    {
        ffclose(out);
    }

    fprintf(stderr, "Processed %d frames\n", nframes);

    for (i = 0; (i < nres); i++)
    {
        for (j = 0; (j < nres); j++)
        {
            totmdmat[i][j] /= nframes;
        }
    }
    write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean smallest distance",
              "Distance (nm)", "Residue Index", "Residue Index",
              nres, nres, resnr, resnr, totmdmat, 0, truncate, rlo, rhi, &nlevels);

    if (bCalcN)
    {
        tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);
        fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
                      "Increase in number of contacts", "Residue", "Ratio", oenv);
        fprintf(fp, "@ legend on\n");
        fprintf(fp, "@ legend box on\n");
        fprintf(fp, "@ legend loctype view\n");
        fprintf(fp, "@ legend 0.75, 0.8\n");
        fprintf(fp, "@ legend string 0 \"Total/mean\"\n");
        fprintf(fp, "@ legend string 1 \"Total\"\n");
        fprintf(fp, "@ legend string 2 \"Mean\"\n");
        fprintf(fp, "@ legend string 3 \"# atoms\"\n");
        fprintf(fp, "@ legend string 4 \"Mean/# atoms\"\n");
        fprintf(fp, "#%3s %8s  %3s  %8s  %3s  %8s\n",
                "res", "ratio", "tot", "mean", "natm", "mean/atm");
        for (i = 0; (i < nres); i++)
        {
            if (mean_n[i] == 0)
            {
                ratio = 1;
            }
            else
            {
                ratio = tot_n[i]/mean_n[i];
            }
            fprintf(fp, "%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
                    i+1, ratio, tot_n[i], mean_n[i], natm[i], mean_n[i]/natm[i]);
        }
        ffclose(fp);
    }

    thanx(stderr);

    return 0;
}

/* scan trajectory and accumulate data */
static int dotrj(hist_t *hsrms, bb_t *bb, int rescnt, FILE *fplog,
    trj_t *trj, const char *fxtc, rvec *xrefbb,
    int nat, matrix box, double *wtot)
{
  int nfr, match, natoms;
  int fpxtc = 0; /* gromacs file handle for trj */
  rvec *x = NULL, *xbb;
  double t1, w, rmsd;
  real t = 0.f;

  xnew(xbb, rescnt * 3);
  /* loop over xtc frames */
  for (nfr = 0; ; nfr++) {
    if (nfr == 0) {
      /* first frame, get the file handle fpxtc, allocate x
       * and check the number of atoms */
      natoms = read_first_x(&fpxtc, fxtc, &t, &x, box);
      if (natoms != nat) {
        fprintf(stderr, "natoms mismatch! top %d, xtc %d\n", nat, natoms);
        return -1;
      }
    } else {
      if (!read_next_x(fpxtc, &t, natoms, x, box))
        break;
    }

    if (ctmd) { /* each frame has a weight */
      w = 1.;
    } else {
      /* for each xtc frame, we search a corresponding
       * frame in trj at the same time */
      for (match = 0; trj->pos < trj->n; trj->pos++) {
        t1 = trj->fr[trj->pos].t;
        if (fabs(t1 - t) < 0.5) {
          match = 1;
          break;
        } else if (t1 > t + 0.5) {
          fprintf(stderr, "time difference %g (trace), %g (xtc)\n", t1, t);
          return -1;
        }
      }
      if (!match) {
        fprintf(stderr, "cannot find a frame for t = %g\n", t);
        break;
      }
      if (!trj->fr[trj->pos].in) {
        w = 0.0;
      } else {
        w = exp(trj->fr[trj->pos].lnw);
      }
    }

    /* skip trivial frames */
    if (fplog == NULL && w < 1e-6) continue;

    /* extract backbone indices */
    xgetbb(xbb, x, bb, rescnt);

    /* fit x to xref */
    rmsd = rv3_rmsd(xbb, NULL, xrefbb, vmass, rescnt*3, 0, NULL, NULL);

    if (fplog) {
      fprintf(fplog, "%10.0f %8.4f %g %d\n", t, rmsd, w, trj->fr[trj->pos].dataid);
    }
    hs_add(hsrms, &rmsd, w, HIST_VERBOSE);
    *wtot += w;
  }
  free(xbb);
  return 0;
}

//.........这里部分代码省略.........
  
  /* Prepare reference frame */
  if (bPBC) {
    gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,box);
    gmx_rmpbc(gpbc,atoms->nr,box,xref);
  }
  if (bFit)
    reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);

  snew(x,natoms);
  snew(xav,natoms);
  ndim=natoms*DIM;
  if (sqrt(GMX_LARGE_INT_MAX)<ndim) {
    gmx_fatal(FARGS,"Number of degrees of freedoms to large for matrix.\n");
  }
  snew(mat,ndim*ndim);

  fprintf(stderr,"Calculating the average structure ...\n");
  nframes0 = 0;
  nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
  if (nat != atoms->nr)
    fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat);
  do {
    nframes0++;
    /* calculate x: a fitted struture of the selected atoms */
    if (bPBC)
      gmx_rmpbc(gpbc,nat,box,xread);
    if (bFit) {
      reset_x(nfit,ifit,nat,NULL,xread,w_rls);
      do_fit(nat,w_rls,xref,xread);
    }
    for (i=0; i<natoms; i++)
      rvec_inc(xav[i],xread[index[i]]);
  } while (read_next_x(oenv,status,&t,nat,xread,box));
  close_trj(status);
  
  inv_nframes = 1.0/nframes0;
  for(i=0; i<natoms; i++)
    for(d=0; d<DIM; d++) {
      xav[i][d] *= inv_nframes;
      xread[index[i]][d] = xav[i][d];
    }
  write_sto_conf_indexed(opt2fn("-av",NFILE,fnm),"Average structure",
			 atoms,xread,NULL,epbcNONE,zerobox,natoms,index);
  sfree(xread);

  fprintf(stderr,"Constructing covariance matrix (%dx%d) ...\n",(int)ndim,(int)ndim);
  nframes=0;
  nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
  tstart = t;
  do {
    nframes++;
    tend = t;
    /* calculate x: a (fitted) structure of the selected atoms */
    if (bPBC)
      gmx_rmpbc(gpbc,nat,box,xread);
    if (bFit) {
      reset_x(nfit,ifit,nat,NULL,xread,w_rls);
      do_fit(nat,w_rls,xref,xread);
    }
    if (bRef)
      for (i=0; i<natoms; i++)
	rvec_sub(xread[index[i]],xref[index[i]],x[i]);
    else
      for (i=0; i<natoms; i++)
	rvec_sub(xread[index[i]],xav[i],x[i]);

//.........这里部分代码省略.........
    pbc = NULL;
    
  gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
  do {
    /* initialisation for correct distance calculations */
    if (pbc) {
      set_pbc(pbc,ePBC,box);
      /* make molecules whole again */
      gmx_rmpbc(gpbc,natoms,box,x);
    }
    /* calculate center of masses */
    for(g=0;(g<ngrps);g++) {
      if (isize[g] == 1) {
	copy_rvec(x[index[g][0]],com[g]);
      } else {
	for(d=0;(d<DIM);d++) {
	  com[g][d]=0;
	  for(i=0;(i<isize[g]);i++) {
	    com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
	  }
	  com[g][d] /= mass[g];
	}
      }
    }
    
    if (!bCutoff) {
      /* write to output */
      fprintf(fp,"%12.7f ",t);
      for(g=0;(g<ngrps/2);g++) {
	if (pbc)
	  pbc_dx(pbc,com[2*g],com[2*g+1],dx);
	else
	  rvec_sub(com[2*g],com[2*g+1],dx);
	
	fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
		norm(dx),dx[XX],dx[YY],dx[ZZ]);
      }
      fprintf(fp,"\n");
    } else {
      for(i=0;(i<isize[1]);i++) { 
	j=index[1][i];
	if (pbc)
	  pbc_dx(pbc,x[j],com[0],dx);
	else
	  rvec_sub(x[j],com[0],dx);
	
	dist2 = norm2(dx);
	if (dist2<cut2) {
	  if (bPrintDist) {
	    res=top->atoms.atom[j].resind;
	    fprintf(stdout,"\rt: %g  %d %s %d %s  %g (nm)\n",
		    t,top->atoms.resinfo[res].nr,*top->atoms.resinfo[res].name,
		    j+1,*top->atoms.atomname[j],sqrt(dist2));
	  }
	  if (bLifeTime)
	    contact_time[i]++;
	} else {
	  if (bLifeTime) {
	    if (contact_time[i]) {
	      add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
	      contact_time[i] = 0;
	    }
	  }
	}
      }
    }
    
    teller++;
  } while (read_next_x(oenv,status,&t,natoms,x,box));
  gmx_rmpbc_done(gpbc);

  if (!bCutoff)
    ffclose(fp);

  close_trj(status);
  
  if (bCutoff && bLifeTime) {
    /* Add the contacts still present in the last frame */
    for(i=0; i<isize[1]; i++)
      if (contact_time[i])
	add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);

    sprintf(buf,"%s - %s within %g nm",
	    grpname[0],grpname[1],cut);
    fp = xvgropen(opt2fn("-lt",NFILE,fnm),
		  buf,"Time (ps)","Number of contacts",oenv);
    for(i=0; i<min(ccount_nalloc,teller-1); i++) {
      /* Account for all subintervals of longer intervals */
      sum = 0;
      for(j=i; j<ccount_nalloc; j++)
	sum += (j-i+1)*ccount[j];

      fprintf(fp,"%10.3f %10.3f\n",i*(t-t0)/(teller-1),sum/(double)(teller-i));
    }
    ffclose(fp);
  }
  
  thanx(stderr);
  return 0;
}

//.........这里部分代码省略.........
        xvgr_legend(out, NLEG, legI, oenv);
    }
    else
    {
        if (bRot)
        {
            if (output_env_get_print_xvgr_codes(oenv))
            {
                fprintf(out, "@ subtitle \"Axes are principal component axes\"\n");
            }
        }
        xvgr_legend(out, NLEG, leg, oenv);
    }
    if (nz == 0)
    {
        gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
    }
    do
    {
        if (nz == 0)
        {
            gmx_rmpbc_copy(gpbc, natoms, box, x, x_s);
        }
        gyro = 0;
        clear_rvec(gvec);
        clear_rvec(gvec1);
        clear_rvec(d);
        clear_rvec(d1);
        for (mol = 0; mol < nmol; mol++)
        {
            tm    = sub_xcm(nz == 0 ? x_s : x, nam, index+mol*nam, top.atoms.atom, xcm, bQ);
            if (nz == 0)
            {
                gyro += calc_gyro(x_s, nam, index+mol*nam, top.atoms.atom,
                                  tm, gvec1, d1, bQ, bRot, bMOI, trans);
            }
            else
            {
                calc_gyro_z(x, box, nam, index+mol*nam, top.atoms.atom, nz, t, out);
            }
            rvec_inc(gvec, gvec1);
            rvec_inc(d, d1);
        }
        if (nmol > 0)
        {
            gyro /= nmol;
            svmul(1.0/nmol, gvec, gvec);
            svmul(1.0/nmol, d, d);
        }

        if (nz == 0)
        {
            if (bRot)
            {
                if (j >= max_moi)
                {
                    max_moi += delta_moi;
                    for (m = 0; (m < DIM); m++)
                    {
                        srenew(moi_trans[m], max_moi*DIM);
                    }
                }
                for (m = 0; (m < DIM); m++)
                {
                    copy_rvec(trans[m], moi_trans[m]+DIM*j);
                }
                fprintf(out, "%10g  %10g  %10g  %10g  %10g\n",
                        t, gyro, d[XX], d[YY], d[ZZ]);
            }
            else
            {
                fprintf(out, "%10g  %10g  %10g  %10g  %10g\n",
                        t, gyro, gvec[XX], gvec[YY], gvec[ZZ]);
            }
        }
        j++;
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);
    if (nz == 0)
    {
        gmx_rmpbc_done(gpbc);
    }

    xvgrclose(out);

    if (bACF)
    {
        int mode = eacVector;

        do_autocorr(opt2fn("-acf", NFILE, fnm), oenv,
                    "Moment of inertia vector ACF",
                    j, 3, moi_trans, (t-t0)/j, mode, FALSE);
        do_view(oenv, opt2fn("-acf", NFILE, fnm), "-nxy");
    }

    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");

    return 0;
}

//.........这里部分代码省略.........
  mdatoms = init_mdatoms(fplog,&mtop,ir.efep!=efepNO);
  atoms2md(&mtop,&ir,0,NULL,0,mtop.natoms,mdatoms);
  update_mdatoms(mdatoms,ir.init_lambda);
  fr      = mk_forcerec();
  fprintf(fplog,"Made forcerec\n");
  init_forcerec(fplog,oenv,fr,NULL,&ir,&mtop,cr,box,FALSE,NULL,NULL,NULL,
                FALSE,-1);
  init_nrnb(&nrnb);
  if (ir.ePBC != epbcNONE)
    gpbc = gmx_rmpbc_init(&top->idef,ir.ePBC,natoms,box);
  
  j=0;
  do {
    if (ir.ePBC != epbcNONE) {
      if (ir.bPeriodicMols)
	set_pbc(&pbc,ir.ePBC,box);
      else
	gmx_rmpbc(gpbc,natoms,box,x);
    }
    
    if (clust) {
      if (j > clust->maxframe)
	gmx_fatal(FARGS,"There are more frames in the trajectory than in the cluster index file. t = %8f\n",t);
      my_clust = clust->inv_clust[j];
      range_check(my_clust,0,clust->clust->nr);
      check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
		 top->idef.iparams,top->idef.functype,
		 x,f,fr,pbc_null,g,dr_clust,my_clust,isize,index,vvindex,&fcd);
    }
    else
      check_viol(fplog,cr,&(top->idef.il[F_DISRES]),
		 top->idef.iparams,top->idef.functype,
		 x,f,fr,pbc_null,g,&dr,0,isize,index,vvindex,&fcd);
    if (bPDB) {
      reset_x(atoms->nr,ind_fit,atoms->nr,NULL,x,w_rls);
      do_fit(atoms->nr,w_rls,x,x);
      if (j == 0) {
	/* Store the first frame of the trajectory as 'characteristic'
	 * for colouring with violations.
	 */
	for(kkk=0; (kkk<atoms->nr); kkk++)
	  copy_rvec(x[kkk],xav[kkk]);
      }
    }
    if (!clust) {
      if (isize > 0) {
	fprintf(xvg,"%10g",t);
	for(i=0; (i<isize); i++)
	  fprintf(xvg,"  %10g",vvindex[i]);
	fprintf(xvg,"\n");
      }    
      fprintf(out,  "%10g  %10g\n",t,dr.sumv);
      fprintf(aver, "%10g  %10g\n",t,dr.averv);
      fprintf(maxxv,"%10g  %10g\n",t,dr.maxv);
      fprintf(numv, "%10g  %10d\n",t,dr.nv);
    }
    j++;
  } while (read_next_x(oenv,status,&t,natoms,x,box));
  close_trj(status);
  if (ir.ePBC != epbcNONE)
    gmx_rmpbc_done(gpbc);

  if (clust) {
    dump_clust_stats(fplog,fcd.disres.nres,&(top->idef.il[F_DISRES]),
		     top->idef.iparams,clust->clust,dr_clust,
		     clust->grpname,isize,index);
  }
  else {
    dump_stats(fplog,j,fcd.disres.nres,&(top->idef.il[F_DISRES]),
	       top->idef.iparams,&dr,isize,index,
	       bPDB ? atoms : NULL);
    if (bPDB) {
      write_sto_conf(opt2fn("-q",NFILE,fnm),
		     "Coloured by average violation in Angstrom",
		     atoms,xav,NULL,ir.ePBC,box);
    }
    dump_disre_matrix(opt2fn_null("-x",NFILE,fnm),&dr,fcd.disres.nres,
		      j,&top->idef,&mtop,max_dr,nlevels,bThird);
    ffclose(out);
    ffclose(aver);
    ffclose(numv);
    ffclose(maxxv);
    if (isize > 0) {
      ffclose(xvg);
      do_view(oenv,opt2fn("-dr",NFILE,fnm),"-nxy");
    }
    do_view(oenv,opt2fn("-dn",NFILE,fnm),"-nxy");
    do_view(oenv,opt2fn("-da",NFILE,fnm),"-nxy");
    do_view(oenv,opt2fn("-ds",NFILE,fnm),"-nxy");
    do_view(oenv,opt2fn("-dm",NFILE,fnm),"-nxy");
  }
  thanx(stderr);

  if (gmx_parallel_env_initialized())
    gmx_finalize();

  gmx_log_close(fplog);
  
  return 0;
}

//.........这里部分代码省略.........
                {
                    /*  store per-molecule order parameter
                     *  To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
                     *  following is for Scd order: */
                    (*slOrder)[j][i] += -1* (1.0/3.0 * (3 * cossum[XX] - 1) + 1.0/3.0 * 0.5 * (3.0 * cossum[YY] - 1));
                }
                if (distcalc)
                {
                    if (radial)
                    {
                        /* bin order parameter by arc distance from reference group*/
                        arcdist            = gmx_angle(dvec, direction);
                        (*distvals)[j][i] += arcdist;
                    }
                    else if (i == 1)
                    {
                        /* Want minimum lateral distance to first group calculated */
                        tmpdist = trace(box);  /* should be max value */
                        for (k = 0; k < distsize; k++)
                        {
                            pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
                            /* at the moment, just remove dvec[axis] */
                            dvec[axis] = 0;
                            tmpdist    = std::min(tmpdist, norm2(dvec));
                        }
                        //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
                        (*distvals)[j][i] += std::sqrt(tmpdist);
                    }
                }
            } /* end loop j, over all atoms in group */

            for (m = 0; m < DIM; m++)
            {
                (*order)[i][m] += (frameorder[m]/size);
            }

            if (!permolecule)
            {   /*Skip following if doing per-molecule*/
                for (k = 0; k < nslices; k++)
                {
                    if (slCount[k]) /* if no elements, nothing has to be added */
                    {
                        (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
                        slFrameorder[k]   = 0; slCount[k] = 0;
                    }
                }
            } /* end loop i, over all groups in indexfile */
        }
        nr_frames++;

    }
    while (read_next_x(oenv, status, &t, x0, box));
    /*********** done with status file **********/

    fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
    gmx_rmpbc_done(gpbc);

    /* average over frames */
    for (i = 1; i < ngrps - 1; i++)
    {
        svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
        fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX],
                (*order)[i][YY], (*order)[i][ZZ]);
        if (bSliced || permolecule)
        {
            for (k = 0; k < nslices; k++)
            {
                (*slOrder)[k][i] /= nr_frames;
            }
        }
        if (distcalc)
        {
            for (k = 0; k < nslices; k++)
            {
                (*distvals)[k][i] /= nr_frames;
            }
        }
    }

    if (bUnsat)
    {
        fprintf(stderr, "Average angle between double bond and normal: %f\n",
                180*sdbangle/(nr_frames * size*M_PI));
    }

    sfree(x0); /* free memory used by coordinate arrays */
    sfree(x1);
    if (comidx != NULL)
    {
        sfree(comidx);
    }
    if (distidx != NULL)
    {
        sfree(distidx);
    }
    if (grpname != NULL)
    {
        sfree(grpname);
    }
}

static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
                                  const char *fnTRX, const char *sgfn,
                                  const char *skfn,
                                  int nslice, int slice_dim,
                                  const char *sgslfn, const char *skslfn,
                                  const gmx_output_env_t *oenv)
{
    FILE        *fpsg = NULL, *fpsk = NULL;
    t_topology   top;
    int          ePBC;
    t_trxstatus *status;
    int          natoms;
    real         t;
    rvec        *xtop, *x;
    matrix       box;
    real         sg, sk;
    atom_id    **index;
    char       **grpname;
    int          i, *isize, ng, nframes;
    real        *sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
    gmx_rmpbc_t  gpbc = NULL;


    read_tps_conf(fnTPS, &top, &ePBC, &xtop, NULL, box, FALSE);

    snew(sg_slice, nslice);
    snew(sk_slice, nslice);
    snew(sg_slice_tot, nslice);
    snew(sk_slice_tot, nslice);
    ng = 1;
    /* get index groups */
    printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
    snew(grpname, ng);
    snew(index, ng);
    snew(isize, ng);
    get_index(&top.atoms, fnNDX, ng, isize, index, grpname);

    /* Analyze trajectory */
    natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
    if (natoms > top.atoms.nr)
    {
        gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
                  top.atoms.nr, natoms);
    }
    check_index(NULL, ng, index[0], NULL, natoms);

    fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
                    oenv);
    fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
                    oenv);

    /* loop over frames */
    gpbc    = gmx_rmpbc_init(&top.idef, ePBC, natoms);
    nframes = 0;
    do
    {
        find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0],
                                &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
        for (i = 0; (i < nslice); i++)
        {
            sg_slice_tot[i] += sg_slice[i];
            sk_slice_tot[i] += sk_slice[i];
        }
        fprintf(fpsg, "%f %f\n", t, sg);
        fprintf(fpsk, "%f %f\n", t, sk);
        nframes++;
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);
    gmx_rmpbc_done(gpbc);

    sfree(grpname);
    sfree(index);
    sfree(isize);

    xvgrclose(fpsg);
    xvgrclose(fpsk);

    fpsg = xvgropen(sgslfn,
                    "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
                    oenv);
    fpsk = xvgropen(skslfn,
                    "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
                    oenv);
    for (i = 0; (i < nslice); i++)
    {
        fprintf(fpsg, "%10g  %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
                sg_slice_tot[i]/nframes);
        fprintf(fpsk, "%10g  %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
                sk_slice_tot[i]/nframes);
    }
    xvgrclose(fpsg);
    xvgrclose(fpsk);
}

//.........这里部分代码省略.........
  }
  
  read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,boxtop,
		FALSE); 
  get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
  
  nalloc = 0;
  time = NULL;
  sbox = NULL;
  sx   = NULL;
  clear_mat(avbox);

  natom=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
  nfr = 0;
  do {
    if (nfr >= nalloc) {
      nalloc += 100;
      srenew(time,nalloc);
      srenew(sbox,nalloc);
      srenew(sx,nalloc);
    }
    
    time[nfr] = t;
    copy_mat(box,sbox[nfr]);
    /* This assumes that the off-diagonal box elements
     * are not affected by jumps across the periodic boundaries.
     */
    m_add(avbox,box,avbox);
    snew(sx[nfr],isize);
    for(i=0; i<isize; i++)
     copy_rvec(x[index[i]],sx[nfr][i]);
    
    nfr++;
  } while (read_next_x(oenv,status,&t,natom,x,box));

  /* clean up */
  sfree(x);
  close_trj(status);
  
  fprintf(stderr,"Read %d frames\n",nfr);

  dt = (time[nfr-1] - time[0])/(nfr - 1);
  /* Some ugly rounding to get nice nice times in the output */
  dt = (int)(10000.0*dt + 0.5)/10000.0;

  invbin = 1.0/rbin;

  if (matfile) {
    if (fmmax <= 0 || fmmax >= nfr)
      fmmax = nfr - 1;
    snew(mcount,fmmax);
    nbin = (int)(rmax*invbin + 0.5);
    if (sbin == 0) {
      mat_nx = fmmax + 1;
    } else {
      invsbin = 1.0/sbin;
      mat_nx = sqrt(fmmax*dt)*invsbin + 1;
    }
    snew(mat,mat_nx);
    for(f=0; f<mat_nx; f++)
      snew(mat[f],nbin);
    rmax2 = sqr(nbin*rbin);
    /* Initialize time zero */
    mat[0][0] = nfr*isize;
    mcount[0] += nfr;
  } else {

//.........这里部分代码省略.........
        }
        pr->grn      = prframecurrent->grn;
        pr->binwidth = prframecurrent->binwidth;
        /* summ up gr and fill r */
        for (i = 0; i < prframecurrent->grn; i++)
        {
            pr->gr[i] += prframecurrent->gr[i];
            pr->r[i]   = prframecurrent->r[i];
        }
        /* normalize histo */
        normalize_probability(prframecurrent->grn, prframecurrent->gr);
        /* convert p(r) to sq */
        sqframecurrent = convert_histogram_to_intensity_curve(prframecurrent, start_q, end_q, q_step);
        /* print frame data if needed */
        if (opt2fn_null("-prframe", NFILE, fnm))
        {
            snew(hdr, 25);
            snew(suffix, GMX_PATH_MAX);
            /* prepare header */
            sprintf(hdr, "g(r), t = %f", t);
            /* prepare output filename */
            fnmdup = dup_tfn(NFILE, fnm);
            sprintf(suffix, "-t%.2f", t);
            add_suffix_to_output_names(fnmdup, NFILE, suffix);
            fp = xvgropen(opt2fn_null("-prframe", NFILE, fnmdup), hdr, "Distance (nm)", "Probability", oenv);
            for (i = 0; i < prframecurrent->grn; i++)
            {
                fprintf(fp, "%10.6f%10.6f\n", prframecurrent->r[i], prframecurrent->gr[i]);
            }
            done_filenms(NFILE, fnmdup);
            fclose(fp);
            sfree(hdr);
            sfree(suffix);
            sfree(fnmdup);
        }
        if (opt2fn_null("-sqframe", NFILE, fnm))
        {
            snew(hdr, 25);
            snew(suffix, GMX_PATH_MAX);
            /* prepare header */
            sprintf(hdr, "I(q), t = %f", t);
            /* prepare output filename */
            fnmdup = dup_tfn(NFILE, fnm);
            sprintf(suffix, "-t%.2f", t);
            add_suffix_to_output_names(fnmdup, NFILE, suffix);
            fp = xvgropen(opt2fn_null("-sqframe", NFILE, fnmdup), hdr, "q (nm^-1)", "s(q)/s(0)", oenv);
            for (i = 0; i < sqframecurrent->qn; i++)
            {
                fprintf(fp, "%10.6f%10.6f\n", sqframecurrent->q[i], sqframecurrent->s[i]);
            }
            done_filenms(NFILE, fnmdup);
            fclose(fp);
            sfree(hdr);
            sfree(suffix);
            sfree(fnmdup);
        }
        /* free pr structure */
        sfree(prframecurrent->gr);
        sfree(prframecurrent->r);
        sfree(prframecurrent);
        /* free sq structure */
        sfree(sqframecurrent->q);
        sfree(sqframecurrent->s);
        sfree(sqframecurrent);
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);

    /* normalize histo */
    normalize_probability(pr->grn, pr->gr);
    sq = convert_histogram_to_intensity_curve(pr, start_q, end_q, q_step);
    /* prepare pr.xvg */
    fp = xvgropen(opt2fn_null("-pr", NFILE, fnm), "G(r)", "Distance (nm)", "Probability", oenv);
    for (i = 0; i < pr->grn; i++)
    {
        fprintf(fp, "%10.6f%10.6f\n", pr->r[i], pr->gr[i]);
    }
    xvgrclose(fp);

    /* prepare sq.xvg */
    fp = xvgropen(opt2fn_null("-sq", NFILE, fnm), "I(q)", "q (nm^-1)", "s(q)/s(0)", oenv);
    for (i = 0; i < sq->qn; i++)
    {
        fprintf(fp, "%10.6f%10.6f\n", sq->q[i], sq->s[i]);
    }
    xvgrclose(fp);
    /*
     * Clean up memory
     */
    sfree(pr->gr);
    sfree(pr->r);
    sfree(pr);
    sfree(sq->q);
    sfree(sq->s);
    sfree(sq);

    please_cite(stdout, "Garmay2012");

    return 0;
}

//.........这里部分代码省略.........
                        if (bRandom)
                        {
                            x[ndx+l][m] = xrot[l][m];
                            v[ndx+l][m] = vrot[l][m];
                        }
                        else
                        {
                            x[ndx+l][m] = xx[l][m];
                            v[ndx+l][m] = v[l][m];
                        }
                    }
                    if (ePBC == epbcSCREW && i % 2 == 1)
                    {
                        /* Rotate around x axis */
                        for (m = YY; m <= ZZ; m++)
                        {
                            x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
                            v[ndx+l][m] = -v[ndx+l][m];
                        }
                    }
                    for (m = 0; (m < DIM); m++)
                    {
                        x[ndx+l][m] += shift[m];
                    }
                    atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
                    atoms->atomname[ndx+l]    = atoms->atomname[l];
                }

                for (l = 0; (l < nres); l++)
                {
                    atoms->resinfo[nrdx+l] = atoms->resinfo[l];
                    if (bRenum)
                    {
                        atoms->resinfo[nrdx+l].nr += nrdx;
                    }
                }
                if (bTRX)
                {
                    if (!read_next_x(oenv, status, &t, natoms, xx, boxx) &&
                        ((i+1)*(j+1)*(k+1) < vol))
                    {
                        gmx_fatal(FARGS, "Not enough frames in trajectory");
                    }
                }
            }
        }
    }
    if (bTRX)
    {
        close_trj(status);
    }

    /* make box bigger */
    for (m = 0; (m < DIM); m++)
    {
        box[m][m] += dist[m];
    }
    svmul(nx, box[XX], box[XX]);
    svmul(ny, box[YY], box[YY]);
    svmul(nz, box[ZZ], box[ZZ]);
    if (ePBC == epbcSCREW && nx % 2 == 0)
    {
        /* With an even number of boxes in x we can forgot about the screw */
        ePBC = epbcXYZ;
    }

    /* move_x(natoms*vol,x,box); */          /* put atoms in box? */

    atoms->nr   *= vol;
    atoms->nres *= vol;

    /*depending on how you look at it, this is either a nasty hack or the way it should work*/
    if (bRenum)
    {
        for (i = 0; i < atoms->nres; i++)
        {
            atoms->resinfo[i].nr = i+1;
        }
    }


    if (bShuffle)
    {
        randwater(0, atoms->nr/nmolat, nmolat, x, v, &seed);
    }
    else if (bSort)
    {
        sortwater(0, atoms->nr/nmolat, nmolat, x, v);
    }
    else if (opt2parg_bSet("-block", asize(pa), pa))
    {
        mkcompact(0, atoms->nr/nmolat, nmolat, x, v, nblock, box);
    }

    write_sto_conf(opt2fn("-o", NFILE, fnm), title, atoms, x, v, ePBC, box);

    thanx(stderr);