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_;
}
//.........这里部分代码省略.........
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;
}
//.........这里部分代码省略.........
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
* I tried the code seems to work properly. B. Hess 2008-4-2.
*/
int main(int argn, char* args[])
{
static char *desc[] = {
"This script will read the xtc trajectory from an MD simulation",
"performed by GROMACS and output all data related to the analysis",
"of preferential interactions.",
"[PAR]",
"Generated: Thu Aug 20 16:49:21 EDT 2015.[BR]",
"Last updated: Thu Aug 27 17:26:01 EDT 2015.",
"[PAR]",
"Future Developments:[BR]",
"[BB]1.[bb] Partition calculation into groups (constituent amino acids).[BR]",
"[BB]2.[bb] Calculate residence times.[BR]",
"[BB]3.[bb] Include a measure for geometric orientation.[BR]",
"[BB]4.[bb] Calculate molecular distributions.[BR]",
"[BB]5.[bb] Parallelize it.",
"[PAR]",
"Remember to check the dependencie of the results on the number of block",
" used by using [TT]-block[tt].",
"Also, check that your simulation has been properly equilibrated by using [TT]-f0[tt].",
"In terms of reducing the computational cost of the analysis, use the option",
"[TT]-skip[tt] to only use every nr[TT]-th[tt] frame from the trajectory."
};
unsigned int protein_sequence=20; /* parse tpr.itp for info */
unsigned int N_co=4; /* look at topology for number of co-solvents */
unsigned int num_sol=7906; /* look at topology as well */
unsigned int skip_nr=1;
unsigned int granularity=2*100;
unsigned int f0=0;
bool bIndex=false;
t_pargs pa[] = {
{"-start", FALSE, etINT, {&f0}, "Start after n-th frame."},
{"-skip", FALSE, etINT, {&skip_nr}, "Only write every nr-th frame."},
{"-grid", FALSE, etINT, {&granularity}, "Number points to calculate."},
{"-seq", FALSE, etINT,{&protein_sequence}, "Number of amino acids in protein sequence."},
{"-Nco", FALSE, etINT,{&N_co}, "Number of cosolvent molecules."},
{"-Sol", FALSE, etINT,{&num_sol}, "Number of water molecules."}
};
char title[STRLEN];
t_topology top;
int ePBC;
rvec *xtop;
matrix box;
t_filenm fnm[] = { /* See filenm.h */
{ efTOP, NULL, NULL, ffREAD },
{ efTPS, NULL, NULL, ffREAD }, /* topology */
{ efTRX, "-f", NULL, ffREAD }, /* trajcetory */
{ efNDX, "-n", NULL, ffOPTRD } /* index file */
};
#define NFILE asize(fnm)
/* Interface. Adds default options. */
CopyRight(stderr,args[0]);
parse_common_args(&argn,args,PCA_CAN_TIME | PCA_CAN_VIEW, NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
read_tps_conf( ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box, TRUE);
sfree(xtop);
/*******************************************************************************************************************/
/*
int iter;
for (iter=0;iter<top.atoms.nr/100;++iter)
{
printf("Atom name: %s\tAtom charge: %f\n", *(top.atoms.atomname[iter]), top.atoms.atom[iter].q );
printf("Atom type: %d\tAtomic Residue NUmber: %d\n", top.atoms.atom[iter].type, top.atoms.atom[iter].resnr );
printf("Chain Identifier: %u\tNr of residues names: %d\n", top.atoms.atom[iter].ptype, top.atoms.nres );
printf("Residue name: %s\n\n", *(top.atoms.resname[iter]) );
}
*/
/*******************************************************************************************************************/
int status;
t_trxframe fr;
int flags = TRX_READ_X;
/* Count Number of Frames */
unsigned int bframes=0, frames=0;
int bwrite;
read_first_frame(&status, ftp2fn(efTRX,NFILE,fnm), &fr, flags);
do {
bwrite = bframes % skip_nr;
++bframes;
if ((bframes>f0) && (bwrite==0)){ ++frames; }
} while ( read_next_frame(status,&fr) );
/* Preparing Block Analysis */
bIndex = ftp2bSet(efNDX,NFILE,fnm);
std::vector<int> blocks_list;
if (bIndex)
{
FILE *fp = ffopen(opt2fn("-n",NFILE,fnm),"r");
int numblock;
while(fscanf(fp,"%d",&numblock)!=EOF)
{
blocks_list.push_back(numblock);
}
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论