C++ close_trj函数代码示例

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

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

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

//.........这里部分代码省略.........
    if (gnx > natoms) {
	gmx_fatal(FARGS,"\nTrajectory does not match selected group!");
    }

    fo = xvgropen(opt2fn("-o",NFILE,fnm),"RMSD","Time (ps)","RMSD (nm)",oenv);
    frc = xvgropen(opt2fn("-rc",NFILE,fnm),"Number of Residues in the alignment","Time (ps)","Residues",oenv);
    
    do {
	t = output_env_conv_time(oenv,fr.time);
	gmx_rmpbc(gpbc,natoms,fr.box,fr.x);
	tapein=gmx_ffopen(pdbfile,"w");
	write_pdbfile_indexed(tapein,NULL,atoms,fr.x,ePBC,fr.box,' ',-1,gnx,index,NULL,TRUE); 
	gmx_ffclose(tapein);
	system(multiprot);
	remove(pdbfile);
	process_multiprot_output(fn, &rmsd, &nres2,rotangles,translation,bCountres,countres);
	fprintf(fo,"%12.7f",t);
	fprintf(fo," %12.7f\n",rmsd);
	fprintf(frc,"%12.7f",t);
	fprintf(frc,"%12d\n",nres2);
	if (bTrjout) {
	    rotate_conf(natoms,fr.x,NULL,rotangles[XX],rotangles[YY],rotangles[ZZ]);
	    for(i=0; i<natoms; i++) {
		rvec_inc(fr.x[i],translation);
	    }
	    switch(outftp) {
		case efTRJ:
		case efTRR:
		case efG87:
		case efXTC:
		    write_trxframe(trxout,&fr,NULL);
		    break;
		case efGRO:
		case efG96:
		case efPDB:
		    sprintf(out_title,"Generated by do_multiprot : %s t= %g %s",
			    title,output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
		    switch(outftp) {
			case efGRO: 
			    write_hconf_p(out,out_title,&useatoms,prec2ndec(fr.prec),
					  fr.x,NULL,fr.box);
			    break;
			case efPDB:
			    fprintf(out,"REMARK    GENERATED BY DO_MULTIPROT\n");
			    sprintf(out_title,"%s t= %g %s",title,output_env_conv_time(oenv,fr.time),output_env_get_time_unit(oenv));
			    /* if reading from pdb, we want to keep the original 
			       model numbering else we write the output frame
			       number plus one, because model 0 is not allowed in pdb */
			    if (ftp==efPDB && fr.step > model_nr) {
				model_nr = fr.step;
			    }
			    else {
				model_nr++;
			    }
			    write_pdbfile(out,out_title,&useatoms,fr.x,ePBC,fr.box,' ',model_nr,NULL,TRUE);
			    break;
			case efG96:
			    fr.title = out_title;
			    fr.bTitle = (nframe == 0);
			    fr.bAtoms = FALSE;
			    fr.bStep = TRUE;
			    fr.bTime = TRUE;
			    write_g96_conf(out,&fr,-1,NULL);
		    }
		    break;
	    }
	}
	nframe++;
    } while(read_next_frame(oenv,status,&fr));
    if (bCountres) {
	fres=  xvgropen(opt2fn("-cr",NFILE,fnm),"Number of frames in which the residues are aligned to","Residue","Number",oenv);
	for (i=0;i<ratoms.nres;i++) {
	    fprintf(fres,"%10d  %12d\n",countres[i].resnr,countres[i].count);
	}
	gmx_ffclose(fres);
    }
    gmx_ffclose(fo);
    gmx_ffclose(frc);
    fprintf(stderr,"\n");
    close_trj(status);
    if (trxout != NULL) {
	close_trx(trxout);
    }
    else if (out != NULL) {
	gmx_ffclose(out);
    }
    view_all(oenv,NFILE, fnm);
    sfree(xr);
    if (bCountres) {
	sfree(countres);
    }
    free_t_atoms(&ratoms,TRUE);
    if (bTrjout) {
	if (outftp==efPDB || outftp==efGRO || outftp==efG96) {
	    free_t_atoms(&useatoms,TRUE);
	}
    }
    gmx_thanx(stderr);
    return 0;
}

//.........这里部分代码省略.........
    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];
                    }
                    mv_mol[XX] += mass*fr.v[j][XX];
                    mv_mol[YY] += mass*fr.v[j][YY];
                    mv_mol[ZZ] += mass*fr.v[j][ZZ];
                }
                c1[i][counter_dim+XX] = mv_mol[XX];
                c1[i][counter_dim+YY] = mv_mol[YY];
                c1[i][counter_dim+ZZ] = mv_mol[ZZ];
            }
        }
        else
        {
            for (i = 0; i < gnx; i++)
            {
                if (bMass)
                {
                    mass = top.atoms.atom[index[i]].m;
                }
                else
                {
                    mass = 1;
                }
                c1[i][counter_dim+XX] = mass*fr.v[index[i]][XX];
                c1[i][counter_dim+YY] = mass*fr.v[index[i]][YY];
                c1[i][counter_dim+ZZ] = mass*fr.v[index[i]][ZZ];
            }
        }

        t1 = fr.time;

        counter++;
    }
    while (read_next_frame(oenv, status, &fr));

    close_trj(status);

    if (counter >= 4)
    {
        /* Compute time step between frames */
        dt = (t1-t0)/(counter-1);
        do_autocorr(opt2fn("-o", NFILE, fnm), oenv,
                    bMass ?
                    "Momentum Autocorrelation Function" :
                    "Velocity Autocorrelation Function",
                    counter, gnx, c1, dt, eacVector, TRUE);

        do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");

        if (opt2bSet("-os", NFILE, fnm))
        {
            calc_spectrum(counter/2, (real *) (c1[0]), (t1-t0)/2, opt2fn("-os", NFILE, fnm),
                          oenv, bRecip);
            do_view(oenv, opt2fn("-os", NFILE, fnm), "-nxy");
        }
    }
    else
    {
        fprintf(stderr, "Not enough frames in trajectory - no output generated.\n");
    }

    return 0;
}

//.........这里部分代码省略.........
      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
	   * I tried the code seems to work properly. B. Hess 2008-4-2.
	   */
	}
	/* Or, if time is set explicitly, we check for overlap/gap */

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

//.........这里部分代码省略.........
                    copy_rvec(x[ind[i][0]],xcom[i]);
                } else {
                    /* Calculate the center of mass */
                    clear_rvec(xcom[i]);
                    mtot = 0;
                    for(j=0; j<gnx[i]; j++) {
                        k = ind[i][j];
                        m = top.atoms.atom[k].m;
                        for(l=0; l<DIM; l++)
                            xcom[i][l] += m*x[k][l];
                        mtot += m;
                    }
                    svmul(1/mtot,xcom[i],xcom[i]);
                }
            }
            pbc_dx(&pbc,xcom[1],xcom[0],direction);
            for(i=0; i<DIM; i++)
                center[i] = xcom[0][i] + 0.5*direction[i];
            unitv(direction,direction);
            for(i=0; i<nindex; i++) {
                j = index[i];
                pbc_dx(&pbc,x[j],center,dx);
                axial = iprod(dx,direction);
                r = sqrt(norm2(dx) - axial*axial);
                if (axial>=-amax && axial<amax && r<rmax) {
                    if (bMirror)
                        r += rmax;
                    grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
                }
            }
        }
        nfr++;
    } while(read_next_x(oenv,status,&t,natoms,x,box));
    close_trj(status);

    /* normalize gridpoints */
    maxgrid = 0;
    if (!bRadial) {
        for (i=0; i<n1; i++) {
            for (j=0; j<n2; j++) {
                grid[i][j] /= nfr;
                if (grid[i][j] > maxgrid)
                    maxgrid = grid[i][j];
            }
        }
    } else {
        for (i=0; i<n1; i++) {
            vol_old = 0;
            for (j=0; j<nradial; j++) {
                switch (nmpower) {
                case -3:
                    vol = M_PI*(j+1)*(j+1)/(invspz*invspz*invspa);
                    break;
                case -2:
                    vol =            (j+1)/(invspz*invspa);
                    break;
                default:
                    vol =             j+1;
                    break;
                }
                if (bMirror)
                    k = j + nradial;
                else
                    k = j;
                grid[i][k] /= nfr*(vol - vol_old);
                if (bMirror)

void sgangle_plot_single(char *fn,char *afile,char *dfile, 
			 char *d1file, char *d2file,
			 atom_id index1[], int gnx1, char *grpn1,
			 atom_id index2[], int gnx2, char *grpn2,
			 t_topology *top,int ePBC)
{
  FILE         
    *sg_angle,           /* xvgr file with angles */
    *sg_distance = NULL, /* xvgr file with distances */
    *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
    *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
  real         
    t,                   /* time */
    angle,               /* cosine of angle between two groups */
    distance,            /* distance between two groups. */
    distance1,           /* distance between plane and one of two atoms */
    distance2;           /* same for second of two atoms */
  int        status,natoms,teller=0;
  int        i;
  rvec       *x0;   /* coordinates, and coordinates corrected for pb */
  rvec       *xzero;
  matrix     box;        
  char       buf[256];   /* for xvgr title */
  
  if ((natoms = read_first_x(&status,fn,&t,&x0,box)) == 0)
    gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
  
  sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
  sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ");
  
  if (dfile) {
    sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
    sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)");
  }
  
  if (d1file) {
    sprintf(buf,"Distance between plane and first atom of vector");
    sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)");
  }
  
  if (d2file) {
    sprintf(buf,"Distance between plane and second atom of vector");
    sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)");
  }
  
  snew(xzero,natoms);

  do {
    teller++;
    
    rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0);
    if (teller==1) {
      for(i=0;i<natoms;i++)
	copy_rvec(x0[i],xzero[i]);
    }
    
    
    calc_angle_single(ePBC,box,
		      xzero,x0,index1,index2,gnx1,gnx2,&angle,
		      &distance,&distance1,&distance2);
    
    fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
    if (dfile)
      fprintf(sg_distance,"%12g  %12g\n",t,distance);
    if (d1file)
      fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
    if (d2file)
      fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
    
  } while (read_next_x(status,&t,natoms,x0,box));
  
  fprintf(stderr,"\n");
  close_trj(status);
  fclose(sg_angle);
  if (dfile)
    fclose(sg_distance);
  if (d1file)
    fclose(sg_distance1);
  if (d2file)
    fclose(sg_distance2);
  
  sfree(x0);
}

//.........这里部分代码省略.........
            if (bFit)
            {
                /*do the least squares fit to mirror of original structure*/
                do_fit(natoms, w_rls, xm, x);
            }

            for (j = 0; j < nrms; j++)
            {
                rlsm[j][teller] =
                    calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
            }
        }
        time[teller] = output_env_conv_time(oenv, t);

        teller++;
        if (teller >= maxframe)
        {
            maxframe += NFRAME;
            srenew(time, maxframe);
            for (j = 0; (j < nrms); j++)
            {
                srenew(rls[j], maxframe);
            }
            if (bMirror)
            {
                for (j = 0; (j < nrms); j++)
                {
                    srenew(rlsm[j], maxframe);
                }
            }
        }
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);

    if (bFile2)
    {
        snew(time2, maxframe2);

        fprintf(stderr, "\nWill read second trajectory file\n");
        snew(mat_x2, NFRAME);
        natoms_trx2 =
            read_first_x(oenv, &status, opt2fn("-f2", NFILE, fnm), &t, &x, box);
        if (natoms_trx2 != natoms_trx)
        {
            gmx_fatal(FARGS,
                      "Second trajectory (%d atoms) does not match the first one"
                      " (%d atoms)", natoms_trx2, natoms_trx);
        }
        tel_mat2 = 0;
        teller2  = 0;
        do
        {
            if (bPBC)
            {
                gmx_rmpbc(gpbc, natoms, box, x);
            }

            if (bReset)
            {
                reset_x(ifit, ind_fit, natoms, NULL, x, w_rls);
            }
            if (ewhat == ewRhoSc)
            {
                norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
            }

//.........这里部分代码省略.........
      if (k_frame >= max_frames) {
        max_frames += XPM_INITIAL_FRAMES;
        if (!(srenew(zbins,max_frames)) ||
            !(srenew(x_axis,max_frames)) ) {
          msg("Out of memory for xpm graph: I abort and try to process what "
              "I have got so far (frame = %g).\n", k_frame);
          break;
        }
      }
      snew(zbins[k_frame],nbins);
      x_axis[k_frame]=t/1000;  // in ns

      // initialize bins to zero
      for (i=0; i<nbins; i++) 
        zbins[k_frame][i] = 0.0;
    }

    /* We are not interested in molecules, just in the z-coordinates 
       of atoms */
    for(i=0; i<gnx; i++) {
      if ( bInCavity (x[gindex[i]], &cavity) ) {  // only look at cylinder
	if (bXPM) {
	  if ( (x[gindex[i]][ZZ] >= minzval) && (x[gindex[i]][ZZ] <= maxzval))
	    zbins[k_frame][(int)floorl((x[gindex[i]][ZZ] - minzval) / spacing)]=1.0;
	}
      } 
      else {
	/* particle OUTSIDE cylinder: mark coordinates as 'invalid'
	   but keep them in the output file */
	for(j=0; j<3; j++) {
	  x[gindex[i]][j] = invalid;
	}
      }
      dump_data (dfiles,"\t%10g", x[gindex[i]][ZZ]);
    }

    /* some printout for debugging/development - MB 
    for (i=0; i<nbins; i++) 
    for (i=0; i<10; i++) 
      printf("%1d",(int)zbins[k_frame][i]);
    printf("\n");
    */

    dump_data (dfiles, "\n");
    k_frame++;
  } while(read_next_x(oenv,status,&t,natoms,x,box));
  close_trj(status);

  if (bXPM) {
    /* colors for matrix */
    rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0; // no water color
    rhi.r = 0.4, rhi.g = 0.4, rhi.b = 1.0; // water color
  
    /* create labels */
    snew(y_axis,nbins+1);
    for(i=0; i<=nbins; i++) {                   // n+1 labels with MAT_SPATIAL_Y
      y_axis[i]=10 * (minzval + (i * spacing)); // Yes, I WANT Angstroms !!!
    }
    /* write xpm matrix */
    xpmOut = fopen(opt2fn("-mat", NFILE, fnm),"w");
    write_xpm(xpmOut,
	      MAT_SPATIAL_Y,
	      "g_zcoord matrix",     // title
	      "existence",           // label for the continuous legend
	      "timeframe / ns",      // label for the x-axis
	      "z coordinate / A",    // label for the y-axis - where the F**k is the
                          	      // bloody Angstrom
	      k_frame, nbins,        // size of the matrix
	      x_axis,            // the x-ticklabels
	      // would be nice to have the times here
	      y_axis,            // the y-ticklabels
	      // would be nice to have the zcoordinate label in nm or A
	      zbins,             // element x,y is matrix[x][y]
	      0.0,               // output lower than lo is set to lo
	      1.0,               // output higher than hi is set to hi
	      rlo,               // rgb value for level lo
	      rhi,               // rgb value for level hi
	      &nlevels);         // number of color levels for the output
    fclose(xpmOut);
  };

  /* clean up a bit */
  fprintf(stderr,"\n");

  if (bDAT) fclose(fData);
  if (bXVG) fclose(out);
  
  if (bXPM) {
    // free memory from xpm matrix data
    for (i=0; i<k_frame; i++)
      sfree(zbins[i]);
    sfree(zbins);
    sfree(x_axis);
  }


  thanx(stdout);
  
  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;
}

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