本文整理汇总了C++中read_next_x函数的典型用法代码示例。如果您正苦于以下问题:C++ read_next_x函数的具体用法?C++ read_next_x怎么用?C++ read_next_x使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了read_next_x函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: new_data
bool new_data(t_xrama *xr)
{
if (!read_next_x(xr->traj,&xr->t,xr->natoms,xr->x,xr->box))
return FALSE;
calc_dihs(xr);
return TRUE;
}
开发者ID:Chadi-akel,项目名称:cere,代码行数:9,代码来源:nrama.c
示例2: new_data
gmx_bool new_data(t_xrama *xr)
{
if (!read_next_x(xr->oenv, xr->traj, &xr->t, xr->x, xr->box))
{
return FALSE;
}
calc_dihs(xr);
return TRUE;
}
开发者ID:FoldingAtHome,项目名称:gromacs,代码行数:11,代码来源:nrama.c
示例3: step_man
static bool step_man(t_manager *man, int *nat)
{
static int ncount = 0;
static bool bWarn = false;
bool bEof;
const char *warn;
if (!man->natom)
{
fprintf(stderr, "Not initiated yet!");
exit(1);
}
bEof = read_next_x(man->oenv, man->status, &man->time, man->x, man->box);
*nat = man->natom;
if (ncount == man->nSkip)
{
switch (man->molw->boxtype)
{
case esbTri:
put_atoms_in_triclinic_unitcell(ecenterDEF, man->box, man->natom, man->x);
break;
case esbTrunc:
warn = put_atoms_in_compact_unitcell(man->molw->ePBC, ecenterDEF, man->box,
man->natom, man->x);
if (warn && !bWarn)
{
fprintf(stderr, "\n%s\n", warn);
bWarn = true;
}
break;
case esbRect:
case esbNone:
default:
break;
}
if (man->bPbc)
{
gmx_rmpbc(man->gpbc, man->natom, man->box, man->x);
reset_mols(&(man->top.mols), man->box, man->x);
}
ncount = 0;
}
else
{
if (man->nSkip > 0)
{
ncount++;
return step_man(man, nat);
}
}
return bEof;
}
开发者ID:ElsevierSoftwareX,项目名称:SOFTX-D-15-00003,代码行数:53,代码来源:manager.cpp
示例4: step_man
static bool step_man(t_manager *man, int *nat)
{
static int ncount = 0;
bool bEof;
if (!man->natom)
{
std::fprintf(stderr, "Not initiated yet!");
std::exit(1);
}
bEof = read_next_x(man->oenv, man->status, &man->time, man->x, man->box);
*nat = man->natom;
if (ncount == man->nSkip)
{
auto atomsArrayRef = gmx::arrayRefFromArray(reinterpret_cast<gmx::RVec *>(man->x), man->natom);
switch (man->molw->boxtype)
{
case esbTri:
put_atoms_in_triclinic_unitcell(ecenterDEF, man->box, atomsArrayRef);
break;
case esbTrunc:
put_atoms_in_compact_unitcell(man->molw->ePBC, ecenterDEF, man->box,
atomsArrayRef);
break;
case esbRect:
case esbNone:
default:
break;
}
if (man->bPbc)
{
gmx_rmpbc(man->gpbc, man->natom, man->box, man->x);
reset_mols(&(man->top.mols), man->box, man->x);
}
ncount = 0;
}
else
{
if (man->nSkip > 0)
{
ncount++;
return step_man(man, nat);
}
}
return bEof;
}
开发者ID:friforever,项目名称:gromacs,代码行数:47,代码来源:manager.cpp
示例5: gmx_rotacf
//.........这里部分代码省略.........
{ efNDX, NULL, NULL, ffREAD },
{ efXVG, "-o", "rotacf", ffWRITE }
};
#define NFILE asize(fnm)
int npargs;
t_pargs *ppa;
CopyRight(stderr,argv[0]);
npargs = asize(pa);
ppa = add_acf_pargs(&npargs,pa);
parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL);
rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
if (bVec)
nvec = isize/2;
else
nvec = isize/3;
if (((isize % 3) != 0) && !bVec)
gmx_fatal(FARGS,"number of index elements not multiple of 3, "
"these can not be atom triplets\n");
if (((isize % 2) != 0) && bVec)
gmx_fatal(FARGS,"number of index elements not multiple of 2, "
"these can not be atom doublets\n");
top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
snew(c1,nvec);
for (i=0; (i<nvec); i++)
c1[i]=NULL;
n_alloc=0;
natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
snew(x_s,natoms);
/* Start the loop over frames */
t1 = t0 = t;
teller = 0;
do {
if (teller >= n_alloc) {
n_alloc+=100;
for (i=0; (i<nvec); i++)
srenew(c1[i],DIM*n_alloc);
}
t1 = t;
/* Remove periodicity */
rm_pbc(&(top->idef),ePBC,natoms,box,x,x_s);
/* Compute crossproducts for all vectors, if triplets.
* else, just get the vectors in case of doublets.
*/
if (bVec == FALSE) {
for (i=0; (i<nvec); i++) {
ai=index[3*i];
aj=index[3*i+1];
ak=index[3*i+2];
rvec_sub(x_s[ai],x_s[aj],xij);
rvec_sub(x_s[aj],x_s[ak],xjk);
cprod(xij,xjk,n);
for(m=0; (m<DIM); m++)
c1[i][DIM*teller+m]=n[m];
}
}
else {
for (i=0; (i<nvec); i++) {
ai=index[2*i];
aj=index[2*i+1];
rvec_sub(x_s[ai],x_s[aj],n);
for(m=0; (m<DIM); m++)
c1[i][DIM*teller+m]=n[m];
}
}
/* Increment loop counter */
teller++;
} while (read_next_x(status,&t,natoms,x,box));
close_trj(status);
fprintf(stderr,"\nDone with trajectory\n");
/* Autocorrelation function */
if (teller < 2)
fprintf(stderr,"Not enough frames for correlation function\n");
else {
dt=(t1 - t0)/(teller-1);
mode = eacVector;
do_autocorr(ftp2fn(efXVG,NFILE,fnm),"Rotational Correlation Function",
teller,nvec,c1,dt,mode,bAver);
}
do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
thanx(stderr);
return 0;
}
开发者ID:BioinformaticsArchive,项目名称:GromPy,代码行数:101,代码来源:gmx_rotacf.c
示例6: 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
示例7: gmx_mdmat
//.........这里部分代码省略.........
trxnat = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
nframes = 0;
rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;
rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;
gpbc = gmx_rmpbc_init(&top.idef, ePBC, trxnat, box);
if (bFrames)
{
out = opt2FILE("-frames", NFILE, fnm, "w");
}
do
{
gmx_rmpbc(gpbc, trxnat, box, x);
nframes++;
calc_mat(nres, natoms, rndx, x, index, truncate, mdmat, nmat, ePBC, box);
for (i = 0; (i < nres); i++)
{
for (j = 0; (j < natoms); j++)
{
if (nmat[i][j])
{
totnmat[i][j]++;
}
}
}
for (i = 0; (i < nres); i++)
{
for (j = 0; (j < nres); j++)
{
totmdmat[i][j] += mdmat[i][j];
}
}
if (bFrames)
{
sprintf(label, "t=%.0f ps", t);
write_xpm(out, 0, label, "Distance (nm)", "Residue Index", "Residue Index",
nres, nres, resnr, resnr, mdmat, 0, truncate, rlo, rhi, &nlevels);
}
}
while (read_next_x(oenv, status, &t, trxnat, x, box));
fprintf(stderr, "\n");
close_trj(status);
gmx_rmpbc_done(gpbc);
if (bFrames)
{
ffclose(out);
}
fprintf(stderr, "Processed %d frames\n", nframes);
for (i = 0; (i < nres); i++)
{
for (j = 0; (j < nres); j++)
{
totmdmat[i][j] /= nframes;
}
}
write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean smallest distance",
"Distance (nm)", "Residue Index", "Residue Index",
nres, nres, resnr, resnr, totmdmat, 0, truncate, rlo, rhi, &nlevels);
if (bCalcN)
{
tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);
fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
"Increase in number of contacts", "Residue", "Ratio", oenv);
fprintf(fp, "@ legend on\n");
fprintf(fp, "@ legend box on\n");
fprintf(fp, "@ legend loctype view\n");
fprintf(fp, "@ legend 0.75, 0.8\n");
fprintf(fp, "@ legend string 0 \"Total/mean\"\n");
fprintf(fp, "@ legend string 1 \"Total\"\n");
fprintf(fp, "@ legend string 2 \"Mean\"\n");
fprintf(fp, "@ legend string 3 \"# atoms\"\n");
fprintf(fp, "@ legend string 4 \"Mean/# atoms\"\n");
fprintf(fp, "#%3s %8s %3s %8s %3s %8s\n",
"res", "ratio", "tot", "mean", "natm", "mean/atm");
for (i = 0; (i < nres); i++)
{
if (mean_n[i] == 0)
{
ratio = 1;
}
else
{
ratio = tot_n[i]/mean_n[i];
}
fprintf(fp, "%3d %8.3f %3d %8.3f %3d %8.3f\n",
i+1, ratio, tot_n[i], mean_n[i], natm[i], mean_n[i]/natm[i]);
}
ffclose(fp);
}
thanx(stderr);
return 0;
}
开发者ID:yhalcyon,项目名称:Gromacs,代码行数:101,代码来源:gmx_mdmat.c
示例8: dotrj
/* scan trajectory and accumulate data */
static int dotrj(hist_t *hsrms, bb_t *bb, int rescnt, FILE *fplog,
trj_t *trj, const char *fxtc, rvec *xrefbb,
int nat, matrix box, double *wtot)
{
int nfr, match, natoms;
int fpxtc = 0; /* gromacs file handle for trj */
rvec *x = NULL, *xbb;
double t1, w, rmsd;
real t = 0.f;
xnew(xbb, rescnt * 3);
/* loop over xtc frames */
for (nfr = 0; ; nfr++) {
if (nfr == 0) {
/* first frame, get the file handle fpxtc, allocate x
* and check the number of atoms */
natoms = read_first_x(&fpxtc, fxtc, &t, &x, box);
if (natoms != nat) {
fprintf(stderr, "natoms mismatch! top %d, xtc %d\n", nat, natoms);
return -1;
}
} else {
if (!read_next_x(fpxtc, &t, natoms, x, box))
break;
}
if (ctmd) { /* each frame has a weight */
w = 1.;
} else {
/* for each xtc frame, we search a corresponding
* frame in trj at the same time */
for (match = 0; trj->pos < trj->n; trj->pos++) {
t1 = trj->fr[trj->pos].t;
if (fabs(t1 - t) < 0.5) {
match = 1;
break;
} else if (t1 > t + 0.5) {
fprintf(stderr, "time difference %g (trace), %g (xtc)\n", t1, t);
return -1;
}
}
if (!match) {
fprintf(stderr, "cannot find a frame for t = %g\n", t);
break;
}
if (!trj->fr[trj->pos].in) {
w = 0.0;
} else {
w = exp(trj->fr[trj->pos].lnw);
}
}
/* skip trivial frames */
if (fplog == NULL && w < 1e-6) continue;
/* extract backbone indices */
xgetbb(xbb, x, bb, rescnt);
/* fit x to xref */
rmsd = rv3_rmsd(xbb, NULL, xrefbb, vmass, rescnt*3, 0, NULL, NULL);
if (fplog) {
fprintf(fplog, "%10.0f %8.4f %g %d\n", t, rmsd, w, trj->fr[trj->pos].dataid);
}
hs_add(hsrms, &rmsd, w, HIST_VERBOSE);
*wtot += w;
}
free(xbb);
return 0;
}
开发者ID:3ki5tj,项目名称:gmxuser,代码行数:71,代码来源:rmsdis.c
示例9: gmx_covar
//.........这里部分代码省略.........
/* Prepare reference frame */
if (bPBC) {
gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,box);
gmx_rmpbc(gpbc,atoms->nr,box,xref);
}
if (bFit)
reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);
snew(x,natoms);
snew(xav,natoms);
ndim=natoms*DIM;
if (sqrt(GMX_LARGE_INT_MAX)<ndim) {
gmx_fatal(FARGS,"Number of degrees of freedoms to large for matrix.\n");
}
snew(mat,ndim*ndim);
fprintf(stderr,"Calculating the average structure ...\n");
nframes0 = 0;
nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
if (nat != atoms->nr)
fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat);
do {
nframes0++;
/* calculate x: a fitted struture of the selected atoms */
if (bPBC)
gmx_rmpbc(gpbc,nat,box,xread);
if (bFit) {
reset_x(nfit,ifit,nat,NULL,xread,w_rls);
do_fit(nat,w_rls,xref,xread);
}
for (i=0; i<natoms; i++)
rvec_inc(xav[i],xread[index[i]]);
} while (read_next_x(oenv,status,&t,nat,xread,box));
close_trj(status);
inv_nframes = 1.0/nframes0;
for(i=0; i<natoms; i++)
for(d=0; d<DIM; d++) {
xav[i][d] *= inv_nframes;
xread[index[i]][d] = xav[i][d];
}
write_sto_conf_indexed(opt2fn("-av",NFILE,fnm),"Average structure",
atoms,xread,NULL,epbcNONE,zerobox,natoms,index);
sfree(xread);
fprintf(stderr,"Constructing covariance matrix (%dx%d) ...\n",(int)ndim,(int)ndim);
nframes=0;
nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
tstart = t;
do {
nframes++;
tend = t;
/* calculate x: a (fitted) structure of the selected atoms */
if (bPBC)
gmx_rmpbc(gpbc,nat,box,xread);
if (bFit) {
reset_x(nfit,ifit,nat,NULL,xread,w_rls);
do_fit(nat,w_rls,xref,xread);
}
if (bRef)
for (i=0; i<natoms; i++)
rvec_sub(xread[index[i]],xref[index[i]],x[i]);
else
for (i=0; i<natoms; i++)
rvec_sub(xread[index[i]],xav[i],x[i]);
开发者ID:BradleyDickson,项目名称:ABPenabledGROMACS,代码行数:67,代码来源:gmx_covar.c
示例10: gmx_dist
//.........这里部分代码省略.........
pbc = NULL;
gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
do {
/* initialisation for correct distance calculations */
if (pbc) {
set_pbc(pbc,ePBC,box);
/* make molecules whole again */
gmx_rmpbc(gpbc,natoms,box,x);
}
/* calculate center of masses */
for(g=0;(g<ngrps);g++) {
if (isize[g] == 1) {
copy_rvec(x[index[g][0]],com[g]);
} else {
for(d=0;(d<DIM);d++) {
com[g][d]=0;
for(i=0;(i<isize[g]);i++) {
com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
}
com[g][d] /= mass[g];
}
}
}
if (!bCutoff) {
/* write to output */
fprintf(fp,"%12.7f ",t);
for(g=0;(g<ngrps/2);g++) {
if (pbc)
pbc_dx(pbc,com[2*g],com[2*g+1],dx);
else
rvec_sub(com[2*g],com[2*g+1],dx);
fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
norm(dx),dx[XX],dx[YY],dx[ZZ]);
}
fprintf(fp,"\n");
} else {
for(i=0;(i<isize[1]);i++) {
j=index[1][i];
if (pbc)
pbc_dx(pbc,x[j],com[0],dx);
else
rvec_sub(x[j],com[0],dx);
dist2 = norm2(dx);
if (dist2<cut2) {
if (bPrintDist) {
res=top->atoms.atom[j].resind;
fprintf(stdout,"\rt: %g %d %s %d %s %g (nm)\n",
t,top->atoms.resinfo[res].nr,*top->atoms.resinfo[res].name,
j+1,*top->atoms.atomname[j],sqrt(dist2));
}
if (bLifeTime)
contact_time[i]++;
} else {
if (bLifeTime) {
if (contact_time[i]) {
add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
contact_time[i] = 0;
}
}
}
}
}
teller++;
} while (read_next_x(oenv,status,&t,natoms,x,box));
gmx_rmpbc_done(gpbc);
if (!bCutoff)
ffclose(fp);
close_trj(status);
if (bCutoff && bLifeTime) {
/* Add the contacts still present in the last frame */
for(i=0; i<isize[1]; i++)
if (contact_time[i])
add_contact_time(&ccount,&ccount_nalloc,contact_time[i]-1);
sprintf(buf,"%s - %s within %g nm",
grpname[0],grpname[1],cut);
fp = xvgropen(opt2fn("-lt",NFILE,fnm),
buf,"Time (ps)","Number of contacts",oenv);
for(i=0; i<min(ccount_nalloc,teller-1); i++) {
/* Account for all subintervals of longer intervals */
sum = 0;
for(j=i; j<ccount_nalloc; j++)
sum += (j-i+1)*ccount[j];
fprintf(fp,"%10.3f %10.3f\n",i*(t-t0)/(teller-1),sum/(double)(teller-i));
}
ffclose(fp);
}
thanx(stderr);
return 0;
}
开发者ID:BradleyDickson,项目名称:ABPenabledGROMACS,代码行数:101,代码来源:gmx_dist.c
示例11: 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
示例12: 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
示例13: calc_order
//.........这里部分代码省略.........
{
/* store per-molecule order parameter
* To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
* following is for Scd order: */
(*slOrder)[j][i] += -1* (1.0/3.0 * (3 * cossum[XX] - 1) + 1.0/3.0 * 0.5 * (3.0 * cossum[YY] - 1));
}
if (distcalc)
{
if (radial)
{
/* bin order parameter by arc distance from reference group*/
arcdist = gmx_angle(dvec, direction);
(*distvals)[j][i] += arcdist;
}
else if (i == 1)
{
/* Want minimum lateral distance to first group calculated */
tmpdist = trace(box); /* should be max value */
for (k = 0; k < distsize; k++)
{
pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
/* at the moment, just remove dvec[axis] */
dvec[axis] = 0;
tmpdist = std::min(tmpdist, norm2(dvec));
}
//fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
(*distvals)[j][i] += std::sqrt(tmpdist);
}
}
} /* end loop j, over all atoms in group */
for (m = 0; m < DIM; m++)
{
(*order)[i][m] += (frameorder[m]/size);
}
if (!permolecule)
{ /*Skip following if doing per-molecule*/
for (k = 0; k < nslices; k++)
{
if (slCount[k]) /* if no elements, nothing has to be added */
{
(*slOrder)[k][i] += slFrameorder[k]/slCount[k];
slFrameorder[k] = 0; slCount[k] = 0;
}
}
} /* end loop i, over all groups in indexfile */
}
nr_frames++;
}
while (read_next_x(oenv, status, &t, x0, box));
/*********** done with status file **********/
fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
gmx_rmpbc_done(gpbc);
/* average over frames */
for (i = 1; i < ngrps - 1; i++)
{
svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX],
(*order)[i][YY], (*order)[i][ZZ]);
if (bSliced || permolecule)
{
for (k = 0; k < nslices; k++)
{
(*slOrder)[k][i] /= nr_frames;
}
}
if (distcalc)
{
for (k = 0; k < nslices; k++)
{
(*distvals)[k][i] /= nr_frames;
}
}
}
if (bUnsat)
{
fprintf(stderr, "Average angle between double bond and normal: %f\n",
180*sdbangle/(nr_frames * size*M_PI));
}
sfree(x0); /* free memory used by coordinate arrays */
sfree(x1);
if (comidx != NULL)
{
sfree(comidx);
}
if (distidx != NULL)
{
sfree(distidx);
}
if (grpname != NULL)
{
sfree(grpname);
}
}
开发者ID:tanigawa,项目名称:gromacs,代码行数:101,代码来源:gmx_order.cpp
示例14: 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
示例15: gmx_vanhove
//.........这里部分代码省略.........
}
read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,boxtop,
FALSE);
get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
nalloc = 0;
time = NULL;
sbox = NULL;
sx = NULL;
clear_mat(avbox);
natom=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
nfr = 0;
do {
if (nfr >= nalloc) {
nalloc += 100;
srenew(time,nalloc);
srenew(sbox,nalloc);
srenew(sx,nalloc);
}
time[nfr] = t;
copy_mat(box,sbox[nfr]);
/* This assumes that the off-diagonal box elements
* are not affected by jumps across the periodic boundaries.
*/
m_add(avbox,box,avbox);
snew(sx[nfr],isize);
for(i=0; i<isize; i++)
copy_rvec(x[index[i]],sx[nfr][i]);
nfr++;
} while (read_next_x(oenv,status,&t,natom,x,box));
/* clean up */
sfree(x);
close_trj(status);
fprintf(stderr,"Read %d frames\n",nfr);
dt = (time[nfr-1] - time[0])/(nfr - 1);
/* Some ugly rounding to get nice nice times in the output */
dt = (int)(10000.0*dt + 0.5)/10000.0;
invbin = 1.0/rbin;
if (matfile) {
if (fmmax <= 0 || fmmax >= nfr)
fmmax = nfr - 1;
snew(mcount,fmmax);
nbin = (int)(rmax*invbin + 0.5);
if (sbin == 0) {
mat_nx = fmmax + 1;
} else {
invsbin = 1.0/sbin;
mat_nx = sqrt(fmmax*dt)*invsbin + 1;
}
snew(mat,mat_nx);
for(f=0; f<mat_nx; f++)
snew(mat[f],nbin);
rmax2 = sqr(nbin*rbin);
/* Initialize time zero */
mat[0][0] = nfr*isize;
mcount[0] += nfr;
} else {
开发者ID:andersx,项目名称:gmx-debug,代码行数:67,代码来源:gmx_vanhove.c
示例16: gmx_sans
//.........这里部分代码省略.........
}
pr->grn = prframecurrent->grn;
pr->binwidth = prframecurrent->binwidth;
/* summ up gr and fill r */
for (i = 0; i < prframecurrent->grn; i++)
{
pr->gr[i] += prframecurrent->gr[i];
pr->r[i] = prframecurrent->r[i];
}
/* normalize histo */
normalize_probability(prframecurrent->grn, prframecurrent->gr);
/* convert p(r) to sq */
sqframecurrent = convert_histogram_to_intensity_curve(prframecurrent, start_q, end_q, q_step);
/* print frame data if needed */
if (opt2fn_null("-prframe", NFILE, fnm))
{
snew(hdr, 25);
snew(suffix, GMX_PATH_MAX);
/* prepare header */
sprintf(hdr, "g(r), t = %f", t);
/* prepare output filename */
fnmdup = dup_tfn(NFILE, fnm);
sprintf(suffix, "-t%.2f", t);
add_suffix_to_output_names(fnmdup, NFILE, suffix);
fp = xvgropen(opt2fn_null("-prframe", NFILE, fnmdup), hdr, "Distance (nm)", "Probability", oenv);
for (i = 0; i < prframecurrent->grn; i++)
{
fprintf(fp, "%10.6f%10.6f\n", prframecurrent->r[i], prframecurrent->gr[i]);
}
done_filenms(NFILE, fnmdup);
fclose(fp);
sfree(hdr);
sfree(suffix);
sfree(fnmdup);
}
if (opt2fn_null("-sqframe", NFILE, fnm))
{
snew(hdr, 25);
snew(suffix, GMX_PATH_MAX);
/* prepare header */
sprintf(hdr, "I(q), t = %f", t);
/* prepare output filename */
fnmdup = dup_tfn(NFILE, fnm);
sprintf(suffix, "-t%.2f", t);
add_suffix_to_output_names(fnmdup, NFILE, suffix);
fp = xvgropen(opt2fn_null("-sqframe", NFILE, fnmdup), hdr, "q (nm^-1)", "s(q)/s(0)", oenv);
for (i = 0; i < sqframecurrent->qn; i++)
{
fprintf(fp, "%10.6f%10.6f\n", sqframecurrent->q[i], sqframecurrent->s[i]);
}
done_filenms(NFILE, fnmdup);
fclose(fp);
sfree(hdr);
sfree(suffix);
sfree(fnmdup);
}
/* free pr structure */
sfree(prframecurrent->gr);
sfree(prframecurrent->r);
sfree(prframecurrent);
/* free sq structure */
sfree(sqframecurrent->q);
sfree(sqframecurrent->s);
sfree(sqframecurrent);
}
while (read_next_x(oenv, status, &t, x, box));
close_trj(status);
/* normalize histo */
normalize_probability(pr->grn, pr->gr);
sq = convert_histogram_to_intensity_curve(pr, start_q, end_q, q_step);
/* prepare pr.xvg */
fp = xvgropen(opt2fn_null("-pr", NFILE, fnm), "G(r)", "Distance (nm)", "Probability", oenv);
for (i = 0; i < pr->grn; i++)
{
fprintf(fp, "%10.6f%10.6f\n", pr->r[i], pr->gr[i]);
}
xvgrclose(fp);
/* prepare sq.xvg */
fp = xvgropen(opt2fn_null("-sq", NFILE, fnm), "I(q)", "q (nm^-1)", "s(q)/s(0)", oenv);
for (i = 0; i < sq->qn; i++)
{
fprintf(fp, "%10.6f%10.6f\n", sq->q[i], sq->s[i]);
}
xvgrclose(fp);
/*
* Clean up memory
*/
sfree(pr->gr);
sfree(pr->r);
sfree(pr);
sfree(sq->q);
sfree(sq->s);
sfree(sq);
please_cite(stdout, "Garmay2012");
return 0;
}
开发者ID:ElsevierSoftwareX,项目名称:SOFTX-D-15-00003,代码行数:101,代码来源:gmx_sans.c
示例17: gmx_genconf
//.........这里部分代码省略.........
if (bRandom)
{
x[ndx+l][m] = xrot[l][m];
v[ndx+l][m] = vrot[l][m];
}
else
{
x[ndx+l][m] = xx[l][m];
v[ndx+l][m] = v[l][m];
}
}
if (ePBC == epbcSCREW && i % 2 == 1)
{
/* Rotate around x axis */
for (m = YY; m <= ZZ; m++)
{
x[ndx+l][m] = box[YY][m] + box[ZZ][m] - x[ndx+l][m];
v[ndx+l][m] = -v[ndx+l][m];
}
}
for (m = 0; (m < DIM); m++)
{
x[ndx+l][m] += shift[m];
}
atoms->atom[ndx+l].resind = nrdx + atoms->atom[l].resind;
atoms->atomname[ndx+l] = atoms->atomname[l];
}
for (l = 0; (l < nres); l++)
{
atoms->resinfo[nrdx+l] = atoms->resinfo[l];
if (bRenum)
{
atoms->resinfo[nrdx+l].nr += nrdx;
}
}
if (bTRX)
{
if (!read_next_x(oenv, status, &t, natoms, xx, boxx) &&
((i+1)*(j+1)*(k+1) < vol))
{
gmx_fatal(FARGS, "Not enough frames in trajectory");
}
}
}
}
}
if (bTRX)
{
close_trj(status);
}
/* make box bigger */
for (m = 0; (m < DIM); m++)
{
box[m][m] += dist[m];
}
svmul(nx, box[XX], box[XX]);
svmul(ny, box[YY], box[YY]);
svmul(nz, box[ZZ], box[ZZ]);
if (ePBC == epbcSCREW && nx % 2 == 0)
{
/* With an even number of boxes in x we can forgot about the screw */
ePBC = epbcXYZ;
}
/* move_x(natoms*vol,x,box); */ /* put atoms in box? */
atoms->nr *= vol;
atoms->nres *= vol;
/*depending on how you look at it, this is either a nasty hack or the way it should work*/
if (bRenum)
{
for (i = 0; i < atoms->nres; i++)
{
atoms->resinfo[i].nr = i+1;
}
}
if (bShuffle)
{
randwater(0, atoms->nr/nmolat, nmolat, x, v, &seed);
}
else if (bSort)
{
sortwater(0, atoms->nr/nmolat, nmolat, x, v);
}
else if (opt2parg_bSet("-block", asize(pa), pa))
{
mkcompact(0, atoms->nr/nmolat, nmolat, x, v, nblock, box);
}
write_sto_conf(opt2fn("-o", NFILE, fnm), title, atoms, x, v, ePBC, box);
thanx(stderr);
|
请发表评论