C++ read_first_frame函数代码示例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    do
    {

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

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

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

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

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

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

    /* Correlation stuff */
    snew(c1, gnx);
    for (i = 0; (i < gnx); i++)
    {
        c1[i] = NULL;
    }

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

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

        t1 = fr.time;

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

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

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

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

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

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

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

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

        if (indexOK)
        {

            if (bDatout)
            {

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

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


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

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

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

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

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

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

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

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

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

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

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

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

    if (bMol)
    {
        if (!bTop)
        {
            gmx_fatal(FARGS, "Need a topology to determine the molecules");
        }
        snew(normm, top.atoms.nr);
        precalc(top, normm);
        index_atom2mol(&gnx, index, &top.mols);
    }

    /* Correlation stuff */
    snew(c1, gnx);
    for (i = 0; (i < gnx); i++)
    {
        c1[i] = NULL;
    }

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

    n_alloc = 0;
    counter = 0;
    do
    {
        if (counter >= n_alloc)
        {
            n_alloc += 100;
            for (i = 0; i < gnx; i++)
            {
                srenew(c1[i], DIM*n_alloc);
            }
        }
        counter_dim = DIM*counter;
        if (bMol)
        {
            for (i = 0; i < gnx; i++)
            {
                clear_rvec(mv_mol);
                k = top.mols.index[index[i]];
                l = top.mols.index[index[i]+1];
                for (j = k; j < l; j++)
                {
                    if (bMass)
                    {
                        mass = top.atoms.atom[j].m;
                    }
                    else
                    {
                        mass = normm[j];
                    }

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

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

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

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

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

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

}

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

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

    snew(impl_->fr, 1);

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

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

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