本文整理汇总了C++中close_trj函数的典型用法代码示例。如果您正苦于以下问题:C++ close_trj函数的具体用法?C++ close_trj怎么用?C++ close_trj使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了close_trj函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: comp_trx
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");
}
开发者ID:TTarenzi,项目名称:MMCG-HAdResS,代码行数:32,代码来源:tpbcmp.c
示例2: write_bfactors
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);
}
开发者ID:alwanderer,项目名称:gromacs,代码行数:58,代码来源:gmx_order.c
示例3: scan_trj_files
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);
}
开发者ID:BioinformaticsArchive,项目名称:GromPy,代码行数:48,代码来源:gmx_trjcat.c
示例4: calc_tetra_order_parm
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);
}
开发者ID:tanigawa,项目名称:gromacs,代码行数:94,代码来源:gmx_order.cpp
示例5: main
//.........这里部分代码省略.........
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;
}
开发者ID:MelroLeandro,项目名称:gromacs,代码行数:101,代码来源:do_multiprot.c
示例6: gmx_velacc
//.........这里部分代码省略.........
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;
}
开发者ID:MrTheodor,项目名称:gromacs,代码行数:101,代码来源:gmx_velacc.cpp
示例7: gmx_trjcat
//.........这里部分代码省略.........
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 */
开发者ID:BioinformaticsArchive,项目名称:GromPy,代码行数:67,代码来源:gmx_trjcat.c
示例8: read_ang_dih
//.........这里部分代码省略.........
{
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;
}
开发者ID:exianshine,项目名称:gromacs,代码行数:101,代码来源:anadih.c
示例9: gmx_densmap
//.........这里部分代码省略.........
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)
开发者ID:BradleyDickson,项目名称:ABPenabledGROMACS,代码行数:67,代码来源:gmx_densmap.c
示例10: sgangle_plot_single
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);
}
开发者ID:alejandrox1,项目名称:gromacs_flatbottom,代码行数:83,代码来源:gmx_sgangle.c
示例11: gmx_rms
//.........这里部分代码省略.........
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);
}
开发者ID:rbharath,项目名称:gromacs,代码行数:67,代码来源:gmx_rms.c
示例12: main
//.........这里部分代码省略.........
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;
}
开发者ID:orbeckst,项目名称:g_count,代码行数:101,代码来源:g_zcoord.c
示例13: gmx_trjorder
//.........这里部分代码省略.........
}
}
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;
}
开发者ID:aalhossary,项目名称:gromacs-HREMD,代码行数:101,代码来源:gmx_trjorder.cpp
示例14: gmx_gyrate
//.........这里部分代码省略.........
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;
}
开发者ID:tanigawa,项目名称:gromacs,代码行数:101,代码来源:gmx_gyrate.cpp
示例15: gmx_disre
//.........这里部分代码省略.........
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;
}
开发者ID:andersx,项目名称:gmx-debug,代码行数:101,代码来源:gmx_disre.c
示例16: gmx_dos
//.........这里部分代码省略.........
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;
nframes++;
}
while (read_next_frame(oenv, status, &fr));
close_trj(status);
dt = (t1-t0)/(nframes-1);
if (nV > 0)
{
V = Vsum/nV;
}
if (bVerbose)
{
printf("Going to do %d fourier transforms of length %d. Hang on.\n",
gnx, nframes);
}
low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
FALSE, FALSE, -1, -1, 0, 0);
snew(dos, DOS_NR);
for (j = 0; (j < DOS_NR); j++)
{
snew(dos[j], nframes+4);
}
if (bVerbose)
{
printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
}
for (i = 0; (i < gnx); i += DIM)
{
mi = top.atoms.atom[i/DIM].m;
for (j = 0; (j < nframes/2); j++)
{
c1j = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
dos[VACF][j] += c1j/Natom;
dos[MVACF][j] += mi*c1j;
}
开发者ID:alwanderer,项目名称:gromacs,代码行数:67,代码来源:gmx_dos.c
示例17: gmx_saltbr
int gmx_saltbr(int argc, char *argv[])
{
const char *desc[] = {
"[TT]g_saltbr[tt] plots the distance between all combination of charged groups",
"as a function of time. The groups are combined in different ways.",
"A minimum distance can be given (i.e. a cut-off), such that groups",
"that are never closer than that distance will not be plotted.[PAR]",
"Output will be in a number of fixed filenames, [TT]min-min.xvg[tt], [TT]plus-min.xvg[tt]",
"and [TT]plus-plus.xvg[tt], or files for every individual ion pair if the [TT]
|
请发表评论