void
TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
{
    if (bTrajOpen_)
    {
        close_trx(status_);
        bTrajOpen_ = false;
    }
    if (gpbc_ != NULL)
    {
        gmx_rmpbc_done(gpbc_);
        gpbc_ = NULL;
    }
}

static void init_gmx(t_x11 *x11, char *program, int nfile, t_filenm fnm[],
                     const output_env_t oenv)
{
    Pixmap                pm;
    t_gmx                *gmx;
    XSizeHints            hints;
    int                   w0, h0;
    int                   natom, natom_trx;
    t_topology            top;
    int                   ePBC;
    matrix                box;
    t_trxframe            fr;
    t_trxstatus          *status;
    char                  quote[256];

    snew(gmx, 1);
    snew(gmx->wd, 1);

    ePBC = read_tpx_top(ftp2fn(efTPR, nfile, fnm),
                        NULL, box, &natom, NULL, NULL, &top);

    read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_DONT_SKIP);
    close_trx(status);
    natom_trx = fr.natoms;

    /* Creates a simple window */
    w0 = DisplayWidth(x11->disp, x11->screen)-132;
    h0 = DisplayHeight(x11->disp, x11->screen)-140;
    bromacs(quote, 255);
    InitWin(gmx->wd, 0, 0, w0, h0, 3, quote);
    gmx->wd->self = XCreateSimpleWindow(x11->disp, x11->root,
                                        gmx->wd->x, gmx->wd->y,
                                        gmx->wd->width, gmx->wd->height,
                                        gmx->wd->bwidth, WHITE, BLACK);
    pm = XCreatePixmapFromBitmapData(x11->disp, x11->root,
                                     (char *)gromacs_bits, gromacs_width,
                                     gromacs_height,
                                     WHITE, BLACK, 1);
    hints.flags      = PMinSize;
    hints.min_width  = 2*EWIDTH+40;
    hints.min_height = EHEIGHT+LDHEIGHT+LEGHEIGHT+40;
    XSetStandardProperties(x11->disp, gmx->wd->self, gmx->wd->text, program,
                           pm, NULL, 0, &hints);

    x11->RegisterCallback(x11, gmx->wd->self, x11->root, MainCallBack, gmx);
    x11->SetInputMask(x11, gmx->wd->self,
                      ButtonPressMask     | ButtonReleaseMask |
                      OwnerGrabButtonMask | ExposureMask      |
                      StructureNotifyMask);

    /* The order of creating windows is important here! */
    /* Manager */
    gmx->man  = init_man(x11, gmx->wd->self, 0, 0, 1, 1, WHITE, BLACK, ePBC, box, oenv);
    gmx->logo = init_logo(x11, gmx->wd->self, false);

    /* Now put all windows in the proper place */
    move_gmx(x11, gmx, w0, h0, false);

    XMapWindow(x11->disp, gmx->wd->self);
    map_man(x11, gmx->man);

    /* Pull Down menu */
    gmx->pd = init_pd(x11, gmx->wd->self, gmx->wd->width,
                      x11->fg, x11->bg,
                      MSIZE, gmx_pd_size, gmx_pd, MenuTitle);

    /* Dialogs & Filters */

    gmx->filter = init_filter(&(top.atoms), ftp2fn_null(efNDX, nfile, fnm),
                              natom_trx);

    init_dlgs(x11, gmx);

    /* Now do file operations */
    set_file(x11, gmx->man, ftp2fn(efTRX, nfile, fnm), ftp2fn(efTPR, nfile, fnm));

    ShowDlg(gmx->dlgs[edFilter]);
}

//.........这里部分代码省略.........
	   * I tried the code seems to work properly. B. Hess 2008-4-2.
	   */
	}
	/* Or, if time is set explicitly, we check for overlap/gap */
	if(cont_type[i]==TIME_EXPLICIT) 
	  if( ( i < nfile_in ) &&
	      ( frout.time < settime[i]-1.5*timestep ) ) 
	    fprintf(stderr, "WARNING: Frames around t=%f %s have a different "
		    "spacing than the rest,\n"
		    "might be a gap or overlap that couldn't be corrected "
		    "automatically.\n",convert_time(frout.time),time_unit());
      }
      
      /* if we don't have a timestep in the current file, use the old one */
      if ( timest[i] != 0 )
	timestep = timest[i];
      
      read_first_frame(&status,fnms[i],&fr,FLAGS);
      if(!fr.bTime) {
	fr.time=0;
	fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
      }
      
      if(cont_type[i]==TIME_EXPLICIT)
	t_corr=settime[i]-fr.time;
      /* t_corr is the amount we want to change the time.
       * If the user has chosen not to change the time for
       * this part of the trajectory t_corr remains at 
       * the value it had in the last part, changing this
       * by the same amount.
       * If no value was given for the first trajectory part
       * we let the time start at zero, see the edit_files routine.
       */
	
      bNewFile=TRUE;
      
      printf("\n");
      if (lasttime != NOTSET)
	printf("lasttime %g\n", lasttime);
      
      do {
	/* copy the input frame to the output frame */
	frout=fr;
	/* set the new time by adding the correct calculated above */
	frout.time += t_corr; 
	/* quit if we have reached the end of what should be written */
	if((end > 0) && (frout.time > end+GMX_REAL_EPS)) {
	  i=nfile_in;
	  break;
	}
	
	/* determine if we should write this frame (dt is handled elsewhere) */
	if (bCat) /* write all frames of all files */ 
	  bWrite = TRUE;
	else if ( bKeepLast ) /* write till last frame of this traj
				 and skip first frame(s) of next traj */
	  bWrite = ( frout.time > lasttime+0.5*timestep );
	else /* write till first frame of next traj */
	  bWrite = ( frout.time < settime[i+1]-0.5*timestep );
	
	if( bWrite && (frout.time >= begin) ) {
	  frame++;
	  if (frame_out == -1)
	    first_time = frout.time;
	  lasttime = frout.time;
	  if (dt==0 || bRmod(frout.time,first_time,dt)) {
	    frame_out++;
	    last_ok_t=frout.time;
	    if(bNewFile) {
	      fprintf(stderr,"\nContinue writing frames from %s t=%g %s, "
		      "frame=%d      \n",
		      fnms[i],convert_time(frout.time),time_unit(),frame);
	      bNewFile=FALSE;
	    }
	    
	    if (bIndex)
	      write_trxframe_indexed(trxout,&frout,isize,index);
	    else
	      write_trxframe(trxout,&frout);
	    if ( ((frame % 10) == 0) || (frame < 10) )
	      fprintf(stderr," ->  frame %6d time %8.3f %s     \r",
		      frame_out,convert_time(frout.time),time_unit());
	  }
	}
      } while( read_next_frame(status,&fr));
      
      close_trj(status);
      
      earliersteps+=step;	  
    }  
    if (trxout >= 0)
      close_trx(trxout);
     
    fprintf(stderr,"\nLast frame written was %d, time %f %s\n",
	    frame,convert_time(last_ok_t),time_unit()); 
  }
  thanx(stderr);
  
  return 0;
}

//.........这里部分代码省略.........
        sum_gyro_tot += sum_gyro;

        if (outp)
        {
            i = -1;
            for (j = 0; j < nat_min/2; j += 2)
            {
                sum_inp[j] /= ninp[j];
                if (i == -1 && sum_inp[j] <= std::exp(-1.0))
                {
                    i = j;
                }
            }
            if (i == -1)
            {
                pers = j;
            }
            else
            {
                /* Do linear interpolation on a log scale */
                pers = i - 2.0
                    + 2.0*(std::log(sum_inp[i-2]) + 1.0)/(std::log(sum_inp[i-2]) - std::log(sum_inp[i]));
            }
            fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
            sum_pers_tot += pers;
        }

        frame++;
    }
    while (read_next_x(oenv, status, &t, x, box));

    gmx_rmpbc_done(gpbc);

    close_trx(status);

    xvgrclose(out);
    if (outv)
    {
        xvgrclose(outv);
    }
    if (outp)
    {
        xvgrclose(outp);
    }

    sum_eed2_tot /= frame;
    sum_gyro_tot /= frame;
    sum_pers_tot /= frame;
    fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
            std::sqrt(sum_eed2_tot));
    fprintf(stdout, "\nAverage radius of gyration:  %.3f (nm)\n",
            std::sqrt(sum_gyro_tot));
    if (opt2bSet("-p", NFILE, fnm))
    {
        fprintf(stdout, "\nAverage persistence length:  %.2f bonds\n",
                sum_pers_tot);
    }

    /* Handle printing of internal distances. */
    if (outi)
    {
        if (output_env_get_print_xvgr_codes(oenv))
        {
            fprintf(outi, "@    xaxes scale Logarithmic\n");
        }
        ymax = -1;

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

//.........这里部分代码省略.........
        printf("Select eigenvectors for output, end your selection with 0\n");
        nout = -1;
        iout = NULL;
        do
        {
            nout++;
            srenew(iout, nout+1);
            if (1 != scanf("%d", &iout[nout]))
            {
                gmx_fatal(FARGS, "Error reading user input");
            }
            iout[nout]--;
        }
        while (iout[nout] >= 0);
        printf("\n");
    }

    /* make an index of the eigenvectors which are present */
    snew(outvec, nout);
    noutvec = 0;
    for (i = 0; i < nout; i++)
    {
        j = 0;
        while ((j < nvec) && (eignr[j] != iout[i]))
        {
            j++;
        }
        if ((j < nvec) && (eignr[j] == iout[i]))
        {
            outvec[noutvec] = j;
            iout[noutvec]   = iout[i];
            noutvec++;
        }
    }

    fprintf(stderr, "%d eigenvectors selected for output\n", noutvec);

    if (seed == -1)
    {
        seed = (int)gmx_rng_make_seed();
        rng  = gmx_rng_init(seed);
    }
    else
    {
        rng = gmx_rng_init(seed);
    }
    fprintf(stderr, "Using seed %d and a temperature of %g K\n", seed, temp);

    snew(xout1, natoms);
    snew(xout2, atoms->nr);
    out  = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
    jran = (unsigned long)((real)im*gmx_rng_uniform_real(rng));
    gmx_rng_destroy(rng);
    for (s = 0; s < nstruct; s++)
    {
        for (i = 0; i < natoms; i++)
        {
            copy_rvec(xav[i], xout1[i]);
        }
        for (j = 0; j < noutvec; j++)
        {
            v = outvec[j];
            /* (r-0.5) n times:  var_n = n * var_1 = n/12
               n=4:  var_n = 1/3, so multiply with 3 */

            rfac  = sqrt(3.0 * BOLTZ*temp/eigval[iout[j]]);
            rhalf = 2.0*rfac;
            rfac  = rfac/(real)im;

            jran = (jran*ia+ic) & im;
            jr   = (real)jran;
            jran = (jran*ia+ic) & im;
            jr  += (real)jran;
            jran = (jran*ia+ic) & im;
            jr  += (real)jran;
            jran = (jran*ia+ic) & im;
            jr  += (real)jran;
            disp = rfac * jr - rhalf;

            for (i = 0; i < natoms; i++)
            {
                for (d = 0; d < DIM; d++)
                {
                    xout1[i][d] += disp*eigvec[v][i][d]*invsqrtm[i];
                }
            }
        }
        for (i = 0; i < natoms; i++)
        {
            copy_rvec(xout1[i], xout2[index[i]]);
        }
        t = s+1;
        write_trx(out, natoms, index, atoms, 0, t, box, xout2, NULL, NULL);
        fprintf(stderr, "\rGenerated %d structures", s+1);
    }
    fprintf(stderr, "\n");
    close_trx(out);

    return 0;
}

//.........这里部分代码省略.........
            }
        }
        else if (bCOM)
        {
            totmass = 0;
            clear_rvec(xcom);
            for (i = 0; i < isize_ref; i++)
            {
                mass     = top.atoms.atom[ind_ref[i]].m;
                totmass += mass;
                for (j = 0; j < DIM; j++)
                {
                    xcom[j] += mass*x[ind_ref[i]][j];
                }
            }
            svmul(1/totmass, xcom, xcom);
            for (i = 0; (i < nwat); i++)
            {
                sa = ind_sol[na*i];
                pbc_dx(&pbc, xcom, xsol[i], dx);
                order[i].i   = sa;
                order[i].d2  = norm2(dx);
            }
        }
        else
        {
            /* Set distance to first atom */
            for (i = 0; (i < nwat); i++)
            {
                sa = ind_sol[na*i];
                pbc_dx(&pbc, x[ind_ref[0]], xsol[i], dx);
                order[i].i   = sa;
                order[i].d2  = norm2(dx);
            }
            for (j = 1; (j < isize_ref); j++)
            {
                sr = ind_ref[j];
                for (i = 0; (i < nwat); i++)
                {
                    pbc_dx(&pbc, x[sr], xsol[i], dx);
                    n2 = norm2(dx);
                    if (n2 < order[i].d2)
                    {
                        order[i].d2  = n2;
                    }
                }
            }
        }

        if (bNShell)
        {
            ncut = 0;
            for (i = 0; (i < nwat); i++)
            {
                if (order[i].d2 <= rcut2)
                {
                    ncut++;
                }
            }
            fprintf(fp, "%10.3f  %8d\n", t, ncut);
        }
        if (out)
        {
            qsort(order, nwat, sizeof(*order), ocomp);
            for (i = 0; (i < nwat); i++)
            {
                for (j = 0; (j < na); j++)
                {
                    swi[ind_sol[na*i]+j] = order[i].i+j;
                }
            }

            /* Store the distance as the B-factor */
            if (bPDBout)
            {
                for (i = 0; (i < nwat); i++)
                {
                    for (j = 0; (j < na); j++)
                    {
                        top.atoms.pdbinfo[order[i].i+j].bfac = std::sqrt(order[i].d2);
                    }
                }
            }
            write_trx(out, natoms, swi, &top.atoms, 0, t, box, x, NULL, NULL);
        }
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);
    if (out)
    {
        close_trx(out);
    }
    if (fp)
    {
        xvgrclose(fp);
    }
    gmx_rmpbc_done(gpbc);

    return 0;
}

//.........这里部分代码省略.........
                       the sphere
                       if (slice > (*nslices))
                       {
                       fprintf(stderr,"Warning: slice = %d\n",slice);
                       }
                     */
                    (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
                }
                else
                {
                    z = x0[index[n][i]][axis];
                    z = z + fudge_z;
                    if (z < 0)
                    {
                        z += box[axis][axis];
                    }
                    if (z > box[axis][axis])
                    {
                        z -= box[axis][axis];
                    }
                    /* determine which slice atom is in */
                    slice                  = static_cast<int>((z / (*slWidth)));
                    (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
                }
            }
        }
        nr_frames++;
    }
    while (read_next_x(oenv, status, &t, x0, box));

    gmx_rmpbc_done(gpbc);

    /*********** done with status file **********/
    close_trx(status);

    /* slCharge now contains the total charge per slice, summed over all
       frames. Now divide by nr_frames and integrate twice
     */


    if (bSpherical)
    {
        fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential"
                "in spherical coordinates\n", nr_frames);
    }
    else
    {
        fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential\n",
                nr_frames);
    }

    for (n = 0; n < nr_grps; n++)
    {
        for (i = 0; i < *nslices; i++)
        {
            if (bSpherical)
            {
                /* charge per volume is now the summed charge, divided by the nr
                   of frames and by the volume of the slice it's in, 4pi r^2 dr
                 */
                slVolume = 4*M_PI * gmx::square(i) * gmx::square(*slWidth) * *slWidth;
                if (slVolume == 0)
                {
                    (*slCharge)[n][i] = 0;
                }
                else

//.........这里部分代码省略.........
    }
    else
    {
        for (i = 0; (i < natoms); i++)
        {
            invsqrtm[i] = 1.0;
        }
    }

    snew(xout, natoms);
    snew(amplitude, nmodes);

    printf("mode phases: %g %g\n", phases[0], phases[1]);

    for (i = 0; i < nmodes; i++)
    {
        kmode       = out_eigidx[i];
        this_eigvec = eigvec[kmode];

        if ( (kmode >= 6) && (eigval[kmode] > 0))
        {
            /* Derive amplitude from temperature and eigenvalue if we can */

            /* Convert eigenvalue to angular frequency, in units s^(-1) */
            omega = std::sqrt(eigval[kmode]*1.0E21/(AVOGADRO*AMU));
            /* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
             * The velocity is thus:
             *
             * v = A*omega*cos(omega*t)*eigenvec.
             *
             * And the average kinetic energy the integral of mass*v*v/2 over a
             * period:
             *
             * (1/4)*mass*A*omega*eigenvec
             *
             * For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
             * The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
             * and the average over a period half of this.
             */

            Ekin = 0;
            for (k = 0; k < natoms; k++)
            {
                m = atoms->atom[k].m;
                for (d = 0; d < DIM; d++)
                {
                    vel   = omega*this_eigvec[k][d];
                    Ekin += 0.5*0.5*m*vel*vel;
                }
            }

            /* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
             * This will also be proportional to A^2
             */
            Ekin *= AMU*1E-18;

            /* Set the amplitude so the energy is kT/2 */
            amplitude[i] = std::sqrt(0.5*BOLTZMANN*temp/Ekin);
        }
        else
        {
            amplitude[i] = refamplitude;
        }
    }

    out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");

    /* Write a sine oscillation around the average structure,
     * modulated by the eigenvector with selected amplitude.
     */

    for (i = 0; i < nframes; i++)
    {
        fraction = static_cast<real>(i)/nframes;
        for (j = 0; j < natoms; j++)
        {
            copy_rvec(xav[j], xout[j]);
        }

        for (k = 0; k < nmodes; k++)
        {
            kmode       = out_eigidx[k];
            this_eigvec = eigvec[kmode];

            for (j = 0; j < natoms; j++)
            {
                for (d = 0; d < DIM; d++)
                {
                    xout[j][d] += amplitude[k]*std::sin(2*M_PI*(fraction+phases[k]/360.0))*this_eigvec[j][d];
                }
            }
        }
        write_trx(out, natoms, dummy, atoms, i, static_cast<real>(i)/nframes, box, xout, NULL, NULL);
    }

    fprintf(stderr, "\n");
    close_trx(out);

    return 0;
}

//.........这里部分代码省略.........
        out_gro = ffopen(out_file,"w");
    }
    else if (ftp == efXTC) {
        trxout = open_trx(out_file,"w");
    }
    
    /* read file and loop through frames */
    int frameN = 0;
    do {
        if (ftp == efGRO) {
            fprintf(out_gro,"Dummy atom inserted into %s, FRAME %i\n",traj_file,frameN);
            fprintf(out_gro,"%i\n",top.atoms.nr+1);
            float CD[3] = { fr.x[a1][0], fr.x[a1][1], fr.x[a1][2] };
            float NE[3] = { fr.x[a2][0], fr.x[a2][1], fr.x[a2][2] };
            float MP[3] ;
            for (int i=0;i<3;i++) {
                MP[i] = (CD[i]+NE[i])*0.5;
            }
            // GRO format:
            // RESID, RESNAME, ATOM, INDEX, X, Y, Z, vX, vY, vZ
            //"%5d%-5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f"
            int index = 1;
			if (insert_at <= 0) {
            	fprintf(out_gro,"%5d%-5s%5s%5d%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n",
                    0,"TCHG","TCHG",0,MP[0],MP[1],MP[2],0.0f,0.0f,0.0f);
					index++;
			}

            /* Loop over atoms */
			int i;
			int resid_offset = 0;
            for (i=0;i<top.atoms.nr;i++){
				if (insert_at == index) {
            		fprintf(out_gro,"%5d%-5s%5s%5d%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n",
                    	top.atoms.atom[i-1].resind+2,"TCHG","TCHG",index,MP[0],MP[1],MP[2],0.0f,0.0f,0.0f);
					index++;
					resid_offset++;
				}
                // Ignoring velocities since I'm using this with mdrun -rerun for
                // force calculations only, which don't care about velocities
                fprintf(out_gro,"%5d%-5s%5s%5d%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n",
                        top.atoms.atom[i].resind+1 + resid_offset,
                        *top.atoms.resinfo[top.atoms.atom[i].resind].name,
                        *top.atoms.atomname[i],
                        index,
                        fr.x[i][0], fr.x[i][1], fr.x[i][2],
                        0.0f,0.0f,0.0f);
                index++;
                if (index > 99999) {
                    index = 0;
                }
            }
            /* Get box information */
            write_hconf_box(out_gro,1,box);
        }
        else if (ftp == efXTC) {
            float CD[3] = { fr.x[a1][0], fr.x[a1][1], fr.x[a1][2] };
            float NE[3] = { fr.x[a2][0], fr.x[a2][1], fr.x[a2][2] };
            rvec MP ;
            for (int i=0;i<3;i++) {
                MP[i] = (CD[i]+NE[i])*0.5;
            }
            rvec * newX = new rvec [top.atoms.nr+1];
			int i = 0;
			int offset = 0;
			if (insert_at <= 0) {
            	for (int j=0;j<3;j++) {
                	newX[i][j] = MP[j];
            	}
				offset++;
			}
            for (i=0; i<top.atoms.nr; i++) {
				if (insert_at == i) {
            		for (int j=0;j<3;j++) {
                		newX[i][j] = MP[j];
					}
					offset++;
				}
                for (int j=0;j<3;j++) {
                    newX[i+offset][j] = fr.x[i][j];
                }
            }
            frout = fr;
            frout.x = newX;
            frout.natoms++;
            write_trxframe(trxout,&frout,gc);
            delete[] newX;
        }
        frameN++;
    } while(read_next_frame(oenv, status, &fr));
    
    if (trxout) {
        close_trx(trxout);
    }
    if (out_gro) {
        ffclose(out_gro);
    }
    
    return 0;
}

//.........这里部分代码省略.........
            remove(buf);
            xf[i].fp2 = gmx_ffopen(buf, "w");
        }
    }

    /* Read reference frame from tpx file to compute helix length */
    snew(xref, top->atoms.nr);
    read_tpx(ftp2fn(efTPR, NFILE, fnm),
             nullptr, nullptr, &natoms, xref, nullptr, nullptr);
    calc_hxprops(nres, bb, xref);
    do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
    sfree(xref);
    if (bDBG)
    {
        fprintf(stderr, "nca=%d, nbb=%d\n", nca, nbb);
        pr_bb(stdout, nres, bb);
    }

    gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);

    teller = 0;
    do
    {
        if ((teller++ % 10) == 0)
        {
            fprintf(stderr, "\rt=%.2f", t);
            fflush(stderr);
        }
        gmx_rmpbc(gpbc, natoms, box, x);


        calc_hxprops(nres, bb, x);
        if (bCheck)
        {
            do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, FALSE, 0, 0);
        }

        if (nca >= 5)
        {
            rms = fit_ahx(nres, bb, natoms, nall, allindex, x, nca, caindex, bFit);

            if (teller == 1)
            {
                write_sto_conf(opt2fn("-cz", NFILE, fnm), "Helix fitted to Z-Axis",
                               &(top->atoms), x, nullptr, ePBC, box);
            }

            xf[efhRAD].val   = radius(xf[efhRAD].fp2, nca, caindex, x);
            xf[efhTWIST].val = twist(nca, caindex, x);
            xf[efhRISE].val  = rise(nca, caindex, x);
            xf[efhLEN].val   = ahx_len(nca, caindex, x);
            xf[efhCD222].val = ellipticity(nres, bb);
            xf[efhDIP].val   = dip(nbb, bbindex, x, top->atoms.atom);
            xf[efhRMS].val   = rms;
            xf[efhCPHI].val  = ca_phi(nca, caindex, x);
            xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2, nres, bb);

            for (j = 0; (j <= efhCPHI); j++)
            {
                fprintf(xf[j].fp,   "%10g  %10g\n", t, xf[j].val);
            }

            av_phipsi(xf[efhPHI].fp, xf[efhPSI].fp, xf[efhPHI].fp2, xf[efhPSI].fp2,
                      t, nres, bb);
            av_hblen(xf[efhHB3].fp, xf[efhHB3].fp2,
                     xf[efhHB4].fp, xf[efhHB4].fp2,
                     xf[efhHB5].fp, xf[efhHB5].fp2,
                     t, nres, bb);
        }
    }
    while (read_next_x(oenv, status, &t, x, box));
    fprintf(stderr, "\n");

    gmx_rmpbc_done(gpbc);

    close_trx(status);

    for (i = 0; (i < nres); i++)
    {
        if (bb[i].nrms > 0)
        {
            fprintf(xf[efhRMSA].fp, "%10d  %10g\n", r0+i, bb[i].rmsa/bb[i].nrms);
        }
        fprintf(xf[efhAHX].fp, "%10d  %10g\n", r0+i, (bb[i].nhx*100.0)/static_cast<real>(teller));
        fprintf(xf[efhJCA].fp, "%10d  %10g\n",
                r0+i, 140.3+(bb[i].jcaha/static_cast<double>(teller)));
    }

    for (i = 0; (i < efhNR); i++)
    {
        xvgrclose(xf[i].fp);
        if (xf[i].bfp2)
        {
            xvgrclose(xf[i].fp2);
        }
        do_view(oenv, xf[i].filenm, "-nxy");
    }

    return 0;
}

//.........这里部分代码省略.........
    }
    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(efTPR, NFILE, fnm), &ePBC);

    snew(c1, nvec);
    for (i = 0; (i < nvec); i++)
    {
        c1[i] = nullptr;
    }
    n_alloc = 0;

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

    gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms);

    /* Start the loop over frames */
    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 */
        gmx_rmpbc_copy(gpbc, 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(oenv, status, &t, x, box));
    close_trx(status);
    fprintf(stderr, "\nDone with trajectory\n");

    gmx_rmpbc_done(gpbc);


    /* 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), oenv, "Rotational Correlation Function",
                    teller, nvec, c1, dt, mode, bAver);
    }

    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);

    return 0;
}

//.........这里部分代码省略.........
        {
            if (!tpr)
            {
                if (bTPRwarn)
                {
                    printf("You need a [REF].tpr[ref] file to analyse temperatures\n");
                    bTPRwarn = FALSE;
                }
            }
            else
            {
                v = fr.v;
                /* Loop over clusters and for each cluster compute 1/2 m v^2 */
                if (max_clust_ind >= 0)
                {
                    ekin = 0;
                    for (i = 0; (i < nindex); i++)
                    {
                        if (clust_index[i] == max_clust_ind)
                        {
                            ai    = index[i];
                            gmx_mtop_atomnr_to_atom(alook, ai, &atom);
                            ekin += 0.5*atom->m*iprod(v[ai], v[ai]);
                        }
                    }
                    temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
                    fprintf(tp, "%10.3f  %10.3f\n", fr.time, temp);
                }
            }
        }
        nframe++;
    }
    while (read_next_frame(oenv, status, &fr));
    close_trx(status);
    xvgrclose(fp);
    xvgrclose(gp);
    xvgrclose(hp);
    xvgrclose(tp);

    gmx_mtop_atomlookup_destroy(alook);

    if (max_clust_ind >= 0)
    {
        fp = gmx_ffopen(mcn, "w");
        fprintf(fp, "[ max_clust ]\n");
        for (i = 0; (i < nindex); i++)
        {
            if (clust_index[i] == max_clust_ind)
            {
                if (bMol)
                {
                    GMX_RELEASE_ASSERT(mols != NULL, "Cannot access index[] from NULL mols pointer");
                    for (j = mols->index[i]; (j < mols->index[i+1]); j++)
                    {
                        fprintf(fp, "%d\n", j+1);
                    }
                }
                else
                {
                    fprintf(fp, "%d\n", index[i]+1);
                }
            }
        }
        gmx_ffclose(fp);
    }

static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
                     real **value, real *time, real dt_remd, int isize,
                     int index[], real dt, const gmx_output_env_t *oenv)
{
    int           i, j, k, natoms, nnn;
    t_trxstatus **fp_in, **fp_out;
    gmx_bool      bCont, *bSet;
    real          t, first_time = 0;
    t_trxframe   *trx;

    snew(fp_in, nset);
    snew(trx, nset);
    snew(bSet, nset);
    natoms = -1;
    t      = -1;
    for (i = 0; (i < nset); i++)
    {
        nnn = read_first_frame(oenv, &(fp_in[i]), fnms[i], &(trx[i]),
                               TRX_NEED_X);
        if (natoms == -1)
        {
            natoms     = nnn;
            first_time = trx[i].time;
        }
        else if (natoms != nnn)
        {
            gmx_fatal(FARGS, "Trajectory file %s has %d atoms while previous trajs had %d atoms", fnms[i], nnn, natoms);
        }
        if (t == -1)
        {
            t = trx[i].time;
        }
        else if (t != trx[i].time)
        {
            gmx_fatal(FARGS, "Trajectory file %s has time %f while previous trajs had time %f", fnms[i], trx[i].time, t);
        }
    }

    snew(fp_out, nset);
    for (i = 0; (i < nset); i++)
    {
        fp_out[i] = open_trx(fnms_out[i], "w");
    }
    k = 0;
    if (std::round(time[k] - t) != 0)
    {
        gmx_fatal(FARGS, "First time in demuxing table does not match trajectories");
    }
    do
    {
        while ((k+1 < nval) && ((trx[0].time - time[k+1]) > dt_remd*0.1))
        {
            k++;
        }
        if (debug)
        {
            fprintf(debug, "trx[0].time = %g, time[k] = %g\n", trx[0].time, time[k]);
        }
        for (i = 0; (i < nset); i++)
        {
            bSet[i] = FALSE;
        }
        for (i = 0; (i < nset); i++)
        {
            j = std::round(value[i][k]);
            range_check(j, 0, nset);
            if (bSet[j])
            {
                gmx_fatal(FARGS, "Demuxing the same replica %d twice at time %f",
                          j, trx[0].time);
            }
            bSet[j] = TRUE;

            if (dt == 0 || bRmod(trx[i].time, first_time, dt))
            {
                if (index)
                {
                    write_trxframe_indexed(fp_out[j], &trx[i], isize, index, NULL);
                }
                else
                {
                    write_trxframe(fp_out[j], &trx[i], NULL);
                }
            }
        }

        bCont = (k < nval);
        for (i = 0; (i < nset); i++)
        {
            bCont = bCont && read_next_frame(oenv, fp_in[i], &trx[i]);
        }
    }
    while (bCont);

    for (i = 0; (i < nset); i++)
    {
        close_trx(fp_in[i]);
        close_trx(fp_out[i]);
    }
}

//.........这里部分代码省略.........
        real last_t = t;

        if(*dt < 0) {
            // User didn't provide dt, so use default:
            // Use trajectory timestep as dt and just store data for all trajectory frames

            real last_dt = *dt;

            do {
                cur_dt = t - last_t;
                last_t = t;

                // Check for consistency of trajectory timestep
                // (this is necessary for proper autocorrelation with default trajectory timestep,
                // otherwise specify a timestep in dt that is at least as large as the largest timestep
                // in the trajectory and is a multiple of all other timesteps in the trajectory)
                if(!GC_TIME_EQ(cur_dt, *dt) && nframes > 1) {
                    gk_log_fatal(FARGS, "Inconsistent time step in frame %d of %s: dt is %f vs. previous dt of %f.\n",
                        nframes, traj_fname, cur_dt, *dt);
                }
                else if(nframes == 1) {
                    // Use the dt between the first two frames of the trajectory as the expected dt.
                    *dt = cur_dt;
                }

                // Expand memory if needed
                if(nframes + 1 > cur_alloc) {
                    cur_alloc *= 2;
                    srenew(*unit_vecs, cur_alloc);
                }

                // Get unit vectors for the atom pairs in this frame
                snew((*unit_vecs)[nframes], npairs);
                for(int p = 0; p < npairs; ++p) {
                    gc_get_unit_vec(x, atompairs[2 * p], atompairs[2 * p + 1], (*unit_vecs)[nframes][p]);
                }

                ++nframes;

            } while(read_next_x(*oenv, status, &t,
#ifndef GRO_V5
                                natoms,
#endif
                                x, box));
        } // if dt < 0
        else { // if provided dt >= 0
            // Only store the data in the trajectory intervals matching the given dt

            last_t = t - *dt; // so that the first trajectory frame will pass dt boundary check
            int trajframe = 0;

            do {
                cur_dt = t - last_t;

                // if this trajectory frame is not on a dt boundary, skip it
                if(GC_TIME_EQ(cur_dt, *dt)) {
                    // Expand memory if needed
                    if(nframes + 1 > cur_alloc) {
                        cur_alloc *= 2;
                        srenew(*unit_vecs, cur_alloc);
                    }

                    // Get unit vectors for the atom pairs in this frame
                    snew((*unit_vecs)[nframes], npairs);
                    for(int p = 0; p < npairs; ++p) {
                        gc_get_unit_vec(x, atompairs[2 * p], atompairs[2 * p + 1], (*unit_vecs)[nframes][p]);
                    }

                    // Reset timestep counter
                    last_t = t;

                    ++nframes;
                }
                else if(cur_dt > *dt) {
                    gk_log_fatal(FARGS, "Error at frame %d of %s: "
                        "given dt = %f is not a multiple of the trajectory timestep, "
                        "or the trajectory has inconsistent timesteps.\n",
                        trajframe, traj_fname, *dt);
                }

                ++trajframe;
            } while(read_next_x(*oenv, status, &t,
#ifndef GRO_V5
                                natoms,
#endif
                                x, box));
        } // if provided dt >= 0
    } // if natoms > 0
    else {
        gk_log_fatal(FARGS, "No atoms found in %s!\n", traj_fname);
    }

    // Done with trajectory
    close_trx(status);
    sfree(x);

    srenew(*unit_vecs, nframes); // free excess memory

    return nframes;
}

//.........这里部分代码省略.........
             * we let the time start at zero, see the edit_files routine.
             */

            bNewFile = TRUE;

            if (!lastTimeSet)
            {
                lasttime    = 0;
                lastTimeSet = true;
            }
            printf("\n");
            printf("lasttime %g\n", lasttime);

            do
            {
                /* copy the input frame to the output frame */
                frout = fr;
                /* set the new time by adding the correct calculated above */
                frout.time += t_corr;
                if (ftpout == efTNG)
                {
                    frout.step += prevEndStep;
                }
                /* quit if we have reached the end of what should be written */
                if ((end > 0) && (frout.time > end+GMX_REAL_EPS))
                {
                    i = nfile_in;
                    break;
                }

                /* determine if we should write this frame (dt is handled elsewhere) */
                if (bCat) /* write all frames of all files */
                {
                    bWrite = TRUE;