本文整理汇总了C++中close_trx函数的典型用法代码示例。如果您正苦于以下问题:C++ close_trx函数的具体用法?C++ close_trx怎么用?C++ close_trx使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了close_trx函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: close_trx
void
TrajectoryAnalysisRunnerCommon::Impl::finishTrajectory()
{
if (bTrajOpen_)
{
close_trx(status_);
bTrajOpen_ = false;
}
if (gpbc_ != NULL)
{
gmx_rmpbc_done(gpbc_);
gpbc_ = NULL;
}
}
开发者ID:smendozabarrera,项目名称:gromacs,代码行数:14,代码来源:runnercommon.cpp
示例2: init_gmx
static void init_gmx(t_x11 *x11, char *program, int nfile, t_filenm fnm[],
const output_env_t oenv)
{
Pixmap pm;
t_gmx *gmx;
XSizeHints hints;
int w0, h0;
int natom, natom_trx;
t_topology top;
int ePBC;
matrix box;
t_trxframe fr;
t_trxstatus *status;
char quote[256];
snew(gmx, 1);
snew(gmx->wd, 1);
ePBC = read_tpx_top(ftp2fn(efTPR, nfile, fnm),
NULL, box, &natom, NULL, NULL, &top);
read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_DONT_SKIP);
close_trx(status);
natom_trx = fr.natoms;
/* Creates a simple window */
w0 = DisplayWidth(x11->disp, x11->screen)-132;
h0 = DisplayHeight(x11->disp, x11->screen)-140;
bromacs(quote, 255);
InitWin(gmx->wd, 0, 0, w0, h0, 3, quote);
gmx->wd->self = XCreateSimpleWindow(x11->disp, x11->root,
gmx->wd->x, gmx->wd->y,
gmx->wd->width, gmx->wd->height,
gmx->wd->bwidth, WHITE, BLACK);
pm = XCreatePixmapFromBitmapData(x11->disp, x11->root,
(char *)gromacs_bits, gromacs_width,
gromacs_height,
WHITE, BLACK, 1);
hints.flags = PMinSize;
hints.min_width = 2*EWIDTH+40;
hints.min_height = EHEIGHT+LDHEIGHT+LEGHEIGHT+40;
XSetStandardProperties(x11->disp, gmx->wd->self, gmx->wd->text, program,
pm, NULL, 0, &hints);
x11->RegisterCallback(x11, gmx->wd->self, x11->root, MainCallBack, gmx);
x11->SetInputMask(x11, gmx->wd->self,
ButtonPressMask | ButtonReleaseMask |
OwnerGrabButtonMask | ExposureMask |
StructureNotifyMask);
/* The order of creating windows is important here! */
/* Manager */
gmx->man = init_man(x11, gmx->wd->self, 0, 0, 1, 1, WHITE, BLACK, ePBC, box, oenv);
gmx->logo = init_logo(x11, gmx->wd->self, false);
/* Now put all windows in the proper place */
move_gmx(x11, gmx, w0, h0, false);
XMapWindow(x11->disp, gmx->wd->self);
map_man(x11, gmx->man);
/* Pull Down menu */
gmx->pd = init_pd(x11, gmx->wd->self, gmx->wd->width,
x11->fg, x11->bg,
MSIZE, gmx_pd_size, gmx_pd, MenuTitle);
/* Dialogs & Filters */
gmx->filter = init_filter(&(top.atoms), ftp2fn_null(efNDX, nfile, fnm),
natom_trx);
init_dlgs(x11, gmx);
/* Now do file operations */
set_file(x11, gmx->man, ftp2fn(efTRX, nfile, fnm), ftp2fn(efTPR, nfile, fnm));
ShowDlg(gmx->dlgs[edFilter]);
}
开发者ID:aalhossary,项目名称:gromacs-HREMD,代码行数:78,代码来源:view.cpp
示例3: gmx_trjcat
//.........这里部分代码省略.........
* I tried the code seems to work properly. B. Hess 2008-4-2.
*/
}
/* Or, if time is set explicitly, we check for overlap/gap */
if(cont_type[i]==TIME_EXPLICIT)
if( ( i < nfile_in ) &&
( frout.time < settime[i]-1.5*timestep ) )
fprintf(stderr, "WARNING: Frames around t=%f %s have a different "
"spacing than the rest,\n"
"might be a gap or overlap that couldn't be corrected "
"automatically.\n",convert_time(frout.time),time_unit());
}
/* if we don't have a timestep in the current file, use the old one */
if ( timest[i] != 0 )
timestep = timest[i];
read_first_frame(&status,fnms[i],&fr,FLAGS);
if(!fr.bTime) {
fr.time=0;
fprintf(stderr,"\nWARNING: Couldn't find a time in the frame.\n");
}
if(cont_type[i]==TIME_EXPLICIT)
t_corr=settime[i]-fr.time;
/* t_corr is the amount we want to change the time.
* If the user has chosen not to change the time for
* this part of the trajectory t_corr remains at
* the value it had in the last part, changing this
* by the same amount.
* If no value was given for the first trajectory part
* we let the time start at zero, see the edit_files routine.
*/
bNewFile=TRUE;
printf("\n");
if (lasttime != NOTSET)
printf("lasttime %g\n", lasttime);
do {
/* copy the input frame to the output frame */
frout=fr;
/* set the new time by adding the correct calculated above */
frout.time += t_corr;
/* quit if we have reached the end of what should be written */
if((end > 0) && (frout.time > end+GMX_REAL_EPS)) {
i=nfile_in;
break;
}
/* determine if we should write this frame (dt is handled elsewhere) */
if (bCat) /* write all frames of all files */
bWrite = TRUE;
else if ( bKeepLast ) /* write till last frame of this traj
and skip first frame(s) of next traj */
bWrite = ( frout.time > lasttime+0.5*timestep );
else /* write till first frame of next traj */
bWrite = ( frout.time < settime[i+1]-0.5*timestep );
if( bWrite && (frout.time >= begin) ) {
frame++;
if (frame_out == -1)
first_time = frout.time;
lasttime = frout.time;
if (dt==0 || bRmod(frout.time,first_time,dt)) {
frame_out++;
last_ok_t=frout.time;
if(bNewFile) {
fprintf(stderr,"\nContinue writing frames from %s t=%g %s, "
"frame=%d \n",
fnms[i],convert_time(frout.time),time_unit(),frame);
bNewFile=FALSE;
}
if (bIndex)
write_trxframe_indexed(trxout,&frout,isize,index);
else
write_trxframe(trxout,&frout);
if ( ((frame % 10) == 0) || (frame < 10) )
fprintf(stderr," -> frame %6d time %8.3f %s \r",
frame_out,convert_time(frout.time),time_unit());
}
}
} while( read_next_frame(status,&fr));
close_trj(status);
earliersteps+=step;
}
if (trxout >= 0)
close_trx(trxout);
fprintf(stderr,"\nLast frame written was %d, time %f %s\n",
frame,convert_time(last_ok_t),time_unit());
}
thanx(stderr);
return 0;
}
开发者ID:BioinformaticsArchive,项目名称:GromPy,代码行数:101,代码来源:gmx_trjcat.c
示例4: gmx_polystat
//.........这里部分代码省略.........
sum_gyro_tot += sum_gyro;
if (outp)
{
i = -1;
for (j = 0; j < nat_min/2; j += 2)
{
sum_inp[j] /= ninp[j];
if (i == -1 && sum_inp[j] <= std::exp(-1.0))
{
i = j;
}
}
if (i == -1)
{
pers = j;
}
else
{
/* Do linear interpolation on a log scale */
pers = i - 2.0
+ 2.0*(std::log(sum_inp[i-2]) + 1.0)/(std::log(sum_inp[i-2]) - std::log(sum_inp[i]));
}
fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
sum_pers_tot += pers;
}
frame++;
}
while (read_next_x(oenv, status, &t, x, box));
gmx_rmpbc_done(gpbc);
close_trx(status);
xvgrclose(out);
if (outv)
{
xvgrclose(outv);
}
if (outp)
{
xvgrclose(outp);
}
sum_eed2_tot /= frame;
sum_gyro_tot /= frame;
sum_pers_tot /= frame;
fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
std::sqrt(sum_eed2_tot));
fprintf(stdout, "\nAverage radius of gyration: %.3f (nm)\n",
std::sqrt(sum_gyro_tot));
if (opt2bSet("-p", NFILE, fnm))
{
fprintf(stdout, "\nAverage persistence length: %.2f bonds\n",
sum_pers_tot);
}
/* Handle printing of internal distances. */
if (outi)
{
if (output_env_get_print_xvgr_codes(oenv))
{
fprintf(outi, "@ xaxes scale Logarithmic\n");
}
ymax = -1;
开发者ID:HITS-MBM,项目名称:gromacs-fda,代码行数:67,代码来源:gmx_polystat.cpp
示例5: main
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)",
//.........这里部分代码省略.........
开发者ID:cudabigdata,项目名称:gromacs,代码行数:101,代码来源:g_anavel.c
示例6: gmx_nmens
//.........这里部分代码省略.........
printf("Select eigenvectors for output, end your selection with 0\n");
nout = -1;
iout = NULL;
do
{
nout++;
srenew(iout, nout+1);
if (1 != scanf("%d", &iout[nout]))
{
gmx_fatal(FARGS, "Error reading user input");
}
iout[nout]--;
}
while (iout[nout] >= 0);
printf("\n");
}
/* make an index of the eigenvectors which are present */
snew(outvec, nout);
noutvec = 0;
for (i = 0; i < nout; i++)
{
j = 0;
while ((j < nvec) && (eignr[j] != iout[i]))
{
j++;
}
if ((j < nvec) && (eignr[j] == iout[i]))
{
outvec[noutvec] = j;
iout[noutvec] = iout[i];
noutvec++;
}
}
fprintf(stderr, "%d eigenvectors selected for output\n", noutvec);
if (seed == -1)
{
seed = (int)gmx_rng_make_seed();
rng = gmx_rng_init(seed);
}
else
{
rng = gmx_rng_init(seed);
}
fprintf(stderr, "Using seed %d and a temperature of %g K\n", seed, temp);
snew(xout1, natoms);
snew(xout2, atoms->nr);
out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
jran = (unsigned long)((real)im*gmx_rng_uniform_real(rng));
gmx_rng_destroy(rng);
for (s = 0; s < nstruct; s++)
{
for (i = 0; i < natoms; i++)
{
copy_rvec(xav[i], xout1[i]);
}
for (j = 0; j < noutvec; j++)
{
v = outvec[j];
/* (r-0.5) n times: var_n = n * var_1 = n/12
n=4: var_n = 1/3, so multiply with 3 */
rfac = sqrt(3.0 * BOLTZ*temp/eigval[iout[j]]);
rhalf = 2.0*rfac;
rfac = rfac/(real)im;
jran = (jran*ia+ic) & im;
jr = (real)jran;
jran = (jran*ia+ic) & im;
jr += (real)jran;
jran = (jran*ia+ic) & im;
jr += (real)jran;
jran = (jran*ia+ic) & im;
jr += (real)jran;
disp = rfac * jr - rhalf;
for (i = 0; i < natoms; i++)
{
for (d = 0; d < DIM; d++)
{
xout1[i][d] += disp*eigvec[v][i][d]*invsqrtm[i];
}
}
}
for (i = 0; i < natoms; i++)
{
copy_rvec(xout1[i], xout2[index[i]]);
}
t = s+1;
write_trx(out, natoms, index, atoms, 0, t, box, xout2, NULL, NULL);
fprintf(stderr, "\rGenerated %d structures", s+1);
}
fprintf(stderr, "\n");
close_trx(out);
return 0;
}
开发者ID:drmaruyama,项目名称:gromacs,代码行数:101,代码来源:gmx_nmens.c
示例7: 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
示例8: calc_potential
//.........这里部分代码省略.........
the sphere
if (slice > (*nslices))
{
fprintf(stderr,"Warning: slice = %d\n",slice);
}
*/
(*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
}
else
{
z = x0[index[n][i]][axis];
z = z + fudge_z;
if (z < 0)
{
z += box[axis][axis];
}
if (z > box[axis][axis])
{
z -= box[axis][axis];
}
/* determine which slice atom is in */
slice = static_cast<int>((z / (*slWidth)));
(*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
}
}
}
nr_frames++;
}
while (read_next_x(oenv, status, &t, x0, box));
gmx_rmpbc_done(gpbc);
/*********** done with status file **********/
close_trx(status);
/* slCharge now contains the total charge per slice, summed over all
frames. Now divide by nr_frames and integrate twice
*/
if (bSpherical)
{
fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential"
"in spherical coordinates\n", nr_frames);
}
else
{
fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential\n",
nr_frames);
}
for (n = 0; n < nr_grps; n++)
{
for (i = 0; i < *nslices; i++)
{
if (bSpherical)
{
/* charge per volume is now the summed charge, divided by the nr
of frames and by the volume of the slice it's in, 4pi r^2 dr
*/
slVolume = 4*M_PI * gmx::square(i) * gmx::square(*slWidth) * *slWidth;
if (slVolume == 0)
{
(*slCharge)[n][i] = 0;
}
else
开发者ID:friforever,项目名称:gromacs,代码行数:67,代码来源:gmx_potential.cpp
示例9: gmx_nmtraj
//.........这里部分代码省略.........
}
else
{
for (i = 0; (i < natoms); i++)
{
invsqrtm[i] = 1.0;
}
}
snew(xout, natoms);
snew(amplitude, nmodes);
printf("mode phases: %g %g\n", phases[0], phases[1]);
for (i = 0; i < nmodes; i++)
{
kmode = out_eigidx[i];
this_eigvec = eigvec[kmode];
if ( (kmode >= 6) && (eigval[kmode] > 0))
{
/* Derive amplitude from temperature and eigenvalue if we can */
/* Convert eigenvalue to angular frequency, in units s^(-1) */
omega = std::sqrt(eigval[kmode]*1.0E21/(AVOGADRO*AMU));
/* Harmonic motion will be x=x0 + A*sin(omega*t)*eigenvec.
* The velocity is thus:
*
* v = A*omega*cos(omega*t)*eigenvec.
*
* And the average kinetic energy the integral of mass*v*v/2 over a
* period:
*
* (1/4)*mass*A*omega*eigenvec
*
* For t =2*pi*n, all energy will be kinetic, and v=A*omega*eigenvec.
* The kinetic energy will be sum(0.5*mass*v*v) if we temporarily set A to 1,
* and the average over a period half of this.
*/
Ekin = 0;
for (k = 0; k < natoms; k++)
{
m = atoms->atom[k].m;
for (d = 0; d < DIM; d++)
{
vel = omega*this_eigvec[k][d];
Ekin += 0.5*0.5*m*vel*vel;
}
}
/* Convert Ekin from amu*(nm/s)^2 to J, i.e., kg*(m/s)^2
* This will also be proportional to A^2
*/
Ekin *= AMU*1E-18;
/* Set the amplitude so the energy is kT/2 */
amplitude[i] = std::sqrt(0.5*BOLTZMANN*temp/Ekin);
}
else
{
amplitude[i] = refamplitude;
}
}
out = open_trx(ftp2fn(efTRO, NFILE, fnm), "w");
/* Write a sine oscillation around the average structure,
* modulated by the eigenvector with selected amplitude.
*/
for (i = 0; i < nframes; i++)
{
fraction = static_cast<real>(i)/nframes;
for (j = 0; j < natoms; j++)
{
copy_rvec(xav[j], xout[j]);
}
for (k = 0; k < nmodes; k++)
{
kmode = out_eigidx[k];
this_eigvec = eigvec[kmode];
for (j = 0; j < natoms; j++)
{
for (d = 0; d < DIM; d++)
{
xout[j][d] += amplitude[k]*std::sin(2*M_PI*(fraction+phases[k]/360.0))*this_eigvec[j][d];
}
}
}
write_trx(out, natoms, dummy, atoms, i, static_cast<real>(i)/nframes, box, xout, NULL, NULL);
}
fprintf(stderr, "\n");
close_trx(out);
return 0;
}
开发者ID:mpharrigan,项目名称:gromacs,代码行数:101,代码来源:gmx_nmtraj.cpp
示例10: gmx_insert_dummy_atom
//.........这里部分代码省略.........
out_gro = ffopen(out_file,"w");
}
else if (ftp == efXTC) {
trxout = open_trx(out_file,"w");
}
/* read file and loop through frames */
int frameN = 0;
do {
if (ftp == efGRO) {
fprintf(out_gro,"Dummy atom inserted into %s, FRAME %i\n",traj_file,frameN);
fprintf(out_gro,"%i\n",top.atoms.nr+1);
float CD[3] = { fr.x[a1][0], fr.x[a1][1], fr.x[a1][2] };
float NE[3] = { fr.x[a2][0], fr.x[a2][1], fr.x[a2][2] };
float MP[3] ;
for (int i=0;i<3;i++) {
MP[i] = (CD[i]+NE[i])*0.5;
}
// GRO format:
// RESID, RESNAME, ATOM, INDEX, X, Y, Z, vX, vY, vZ
//"%5d%-5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f"
int index = 1;
if (insert_at <= 0) {
fprintf(out_gro,"%5d%-5s%5s%5d%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n",
0,"TCHG","TCHG",0,MP[0],MP[1],MP[2],0.0f,0.0f,0.0f);
index++;
}
/* Loop over atoms */
int i;
int resid_offset = 0;
for (i=0;i<top.atoms.nr;i++){
if (insert_at == index) {
fprintf(out_gro,"%5d%-5s%5s%5d%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n",
top.atoms.atom[i-1].resind+2,"TCHG","TCHG",index,MP[0],MP[1],MP[2],0.0f,0.0f,0.0f);
index++;
resid_offset++;
}
// Ignoring velocities since I'm using this with mdrun -rerun for
// force calculations only, which don't care about velocities
fprintf(out_gro,"%5d%-5s%5s%5d%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n",
top.atoms.atom[i].resind+1 + resid_offset,
*top.atoms.resinfo[top.atoms.atom[i].resind].name,
*top.atoms.atomname[i],
index,
fr.x[i][0], fr.x[i][1], fr.x[i][2],
0.0f,0.0f,0.0f);
index++;
if (index > 99999) {
index = 0;
}
}
/* Get box information */
write_hconf_box(out_gro,1,box);
}
else if (ftp == efXTC) {
float CD[3] = { fr.x[a1][0], fr.x[a1][1], fr.x[a1][2] };
float NE[3] = { fr.x[a2][0], fr.x[a2][1], fr.x[a2][2] };
rvec MP ;
for (int i=0;i<3;i++) {
MP[i] = (CD[i]+NE[i])*0.5;
}
rvec * newX = new rvec [top.atoms.nr+1];
int i = 0;
int offset = 0;
if (insert_at <= 0) {
for (int j=0;j<3;j++) {
newX[i][j] = MP[j];
}
offset++;
}
for (i=0; i<top.atoms.nr; i++) {
if (insert_at == i) {
for (int j=0;j<3;j++) {
newX[i][j] = MP[j];
}
offset++;
}
for (int j=0;j<3;j++) {
newX[i+offset][j] = fr.x[i][j];
}
}
frout = fr;
frout.x = newX;
frout.natoms++;
write_trxframe(trxout,&frout,gc);
delete[] newX;
}
frameN++;
} while(read_next_frame(oenv, status, &fr));
if (trxout) {
close_trx(trxout);
}
if (out_gro) {
ffclose(out_gro);
}
return 0;
}
开发者ID:awritchie,项目名称:my_gmx,代码行数:101,代码来源:gmx_insert_dummy_atom.cpp
示例11: gmx_helix
//.........这里部分代码省略.........
remove(buf);
xf[i].fp2 = gmx_ffopen(buf, "w");
}
}
/* Read reference frame from tpx file to compute helix length */
snew(xref, top->atoms.nr);
read_tpx(ftp2fn(efTPR, NFILE, fnm),
nullptr, nullptr, &natoms, xref, nullptr, nullptr);
calc_hxprops(nres, bb, xref);
do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
sfree(xref);
if (bDBG)
{
fprintf(stderr, "nca=%d, nbb=%d\n", nca, nbb);
pr_bb(stdout, nres, bb);
}
gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
teller = 0;
do
{
if ((teller++ % 10) == 0)
{
fprintf(stderr, "\rt=%.2f", t);
fflush(stderr);
}
gmx_rmpbc(gpbc, natoms, box, x);
calc_hxprops(nres, bb, x);
if (bCheck)
{
do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, FALSE, 0, 0);
}
if (nca >= 5)
{
rms = fit_ahx(nres, bb, natoms, nall, allindex, x, nca, caindex, bFit);
if (teller == 1)
{
write_sto_conf(opt2fn("-cz", NFILE, fnm), "Helix fitted to Z-Axis",
&(top->atoms), x, nullptr, ePBC, box);
}
xf[efhRAD].val = radius(xf[efhRAD].fp2, nca, caindex, x);
xf[efhTWIST].val = twist(nca, caindex, x);
xf[efhRISE].val = rise(nca, caindex, x);
xf[efhLEN].val = ahx_len(nca, caindex, x);
xf[efhCD222].val = ellipticity(nres, bb);
xf[efhDIP].val = dip(nbb, bbindex, x, top->atoms.atom);
xf[efhRMS].val = rms;
xf[efhCPHI].val = ca_phi(nca, caindex, x);
xf[efhPPRMS].val = pprms(xf[efhPPRMS].fp2, nres, bb);
for (j = 0; (j <= efhCPHI); j++)
{
fprintf(xf[j].fp, "%10g %10g\n", t, xf[j].val);
}
av_phipsi(xf[efhPHI].fp, xf[efhPSI].fp, xf[efhPHI].fp2, xf[efhPSI].fp2,
t, nres, bb);
av_hblen(xf[efhHB3].fp, xf[efhHB3].fp2,
xf[efhHB4].fp, xf[efhHB4].fp2,
xf[efhHB5].fp, xf[efhHB5].fp2,
t, nres, bb);
}
}
while (read_next_x(oenv, status, &t, x, box));
fprintf(stderr, "\n");
gmx_rmpbc_done(gpbc);
close_trx(status);
for (i = 0; (i < nres); i++)
{
if (bb[i].nrms > 0)
{
fprintf(xf[efhRMSA].fp, "%10d %10g\n", r0+i, bb[i].rmsa/bb[i].nrms);
}
fprintf(xf[efhAHX].fp, "%10d %10g\n", r0+i, (bb[i].nhx*100.0)/static_cast<real>(teller));
fprintf(xf[efhJCA].fp, "%10d %10g\n",
r0+i, 140.3+(bb[i].jcaha/static_cast<double>(teller)));
}
for (i = 0; (i < efhNR); i++)
{
xvgrclose(xf[i].fp);
if (xf[i].bfp2)
{
xvgrclose(xf[i].fp2);
}
do_view(oenv, xf[i].filenm, "-nxy");
}
return 0;
}
开发者ID:HITS-MBM,项目名称:gromacs-fda,代码行数:101,代码来源:gmx_helix.cpp
示例12: gmx_rotacf
//.........这里部分代码省略.........
}
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(efTPR, NFILE, fnm), &ePBC);
snew(c1, nvec);
for (i = 0; (i < nvec); i++)
{
c1[i] = nullptr;
}
n_alloc = 0;
natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
snew(x_s, natoms);
gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms);
/* Start the loop over frames */
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 */
gmx_rmpbc_copy(gpbc, 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(oenv, status, &t, x, box));
close_trx(status);
fprintf(stderr, "\nDone with trajectory\n");
gmx_rmpbc_done(gpbc);
/* 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), oenv, "Rotational Correlation Function",
teller, nvec, c1, dt, mode, bAver);
}
do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);
return 0;
}
开发者ID:HITS-MBM,项目名称:gromacs-fda,代码行数:101,代码来源:gmx_rotacf.cpp
示例13: clust_size
//.........这里部分代码省略.........
{
if (!tpr)
{
if (bTPRwarn)
{
printf("You need a [REF].tpr[ref] file to analyse temperatures\n");
bTPRwarn = FALSE;
}
}
else
{
v = fr.v;
/* Loop over clusters and for each cluster compute 1/2 m v^2 */
if (max_clust_ind >= 0)
{
ekin = 0;
for (i = 0; (i < nindex); i++)
{
if (clust_index[i] == max_clust_ind)
{
ai = index[i];
gmx_mtop_atomnr_to_atom(alook, ai, &atom);
ekin += 0.5*atom->m*iprod(v[ai], v[ai]);
}
}
temp = (ekin*2.0)/(3.0*tfac*max_clust_size*BOLTZ);
fprintf(tp, "%10.3f %10.3f\n", fr.time, temp);
}
}
}
nframe++;
}
while (read_next_frame(oenv, status, &fr));
close_trx(status);
xvgrclose(fp);
xvgrclose(gp);
xvgrclose(hp);
xvgrclose(tp);
gmx_mtop_atomlookup_destroy(alook);
if (max_clust_ind >= 0)
{
fp = gmx_ffopen(mcn, "w");
fprintf(fp, "[ max_clust ]\n");
for (i = 0; (i < nindex); i++)
{
if (clust_index[i] == max_clust_ind)
{
if (bMol)
{
GMX_RELEASE_ASSERT(mols != NULL, "Cannot access index[] from NULL mols pointer");
for (j = mols->index[i]; (j < mols->index[i+1]); j++)
{
fprintf(fp, "%d\n", j+1);
}
}
else
{
fprintf(fp, "%d\n", index[i]+1);
}
}
}
gmx_ffclose(fp);
}
开发者ID:mpharrigan,项目名称:gromacs,代码行数:66,代码来源:gmx_clustsize.cpp
示例14: do_demux
static void do_demux(int nset, char *fnms[], char *fnms_out[], int nval,
real **value, real *time, real dt_remd, int isize,
int index[], real dt, const gmx_output_env_t *oenv)
{
int i, j, k, natoms, nnn;
t_trxstatus **fp_in, **fp_out;
gmx_bool bCont, *bSet;
real t, first_time = 0;
t_trxframe *trx;
snew(fp_in, nset);
snew(trx, nset);
snew(bSet, nset);
natoms = -1;
t = -1;
for (i = 0; (i < nset); i++)
{
nnn = read_first_frame(oenv, &(fp_in[i]), fnms[i], &(trx[i]),
TRX_NEED_X);
if (natoms == -1)
{
natoms = nnn;
first_time = trx[i].time;
}
else if (natoms != nnn)
{
gmx_fatal(FARGS, "Trajectory file %s has %d atoms while previous trajs had %d atoms", fnms[i], nnn, natoms);
}
if (t == -1)
{
t = trx[i].time;
}
else if (t != trx[i].time)
{
gmx_fatal(FARGS, "Trajectory file %s has time %f while previous trajs had time %f", fnms[i], trx[i].time, t);
}
}
snew(fp_out, nset);
for (i = 0; (i < nset); i++)
{
fp_out[i] = open_trx(fnms_out[i], "w");
}
k = 0;
if (std::round(time[k] - t) != 0)
{
gmx_fatal(FARGS, "First time in demuxing table does not match trajectories");
}
do
{
while ((k+1 < nval) && ((trx[0].time - time[k+1]) > dt_remd*0.1))
{
k++;
}
if (debug)
{
fprintf(debug, "trx[0].time = %g, time[k] = %g\n", trx[0].time, time[k]);
}
for (i = 0; (i < nset); i++)
{
bSet[i] = FALSE;
}
for (i = 0; (i < nset); i++)
{
j = std::round(value[i][k]);
range_check(j, 0, nset);
if (bSet[j])
{
gmx_fatal(FARGS, "Demuxing the same replica %d twice at time %f",
j, trx[0].time);
}
bSet[j] = TRUE;
if (dt == 0 || bRmod(trx[i].time, first_time, dt))
{
if (index)
{
write_trxframe_indexed(fp_out[j], &trx[i], isize, index, NULL);
}
else
{
write_trxframe(fp_out[j], &trx[i], NULL);
}
}
}
bCont = (k < nval);
for (i = 0; (i < nset); i++)
{
bCont = bCont && read_next_frame(oenv, fp_in[i], &trx[i]);
}
}
while (bCont);
for (i = 0; (i < nset); i++)
{
close_trx(fp_in[i]);
close_trx(fp_out[i]);
}
}
开发者ID:kmtu,项目名称:gromacs,代码行数:100,代码来源:gmx_trjcat.cpp
示例15: gc_traj2uvecs
//.........这里部分代码省略.........
real last_t = t;
if(*dt < 0) {
// User didn't provide dt, so use default:
// Use trajectory timestep as dt and just store data for all trajectory frames
real last_dt = *dt;
do {
cur_dt = t - last_t;
last_t = t;
// Check for consistency of trajectory timestep
// (this is necessary for proper autocorrelation with default trajectory timestep,
// otherwise specify a timestep in dt that is at least as large as the largest timestep
// in the trajectory and is a multiple of all other timesteps in the trajectory)
if(!GC_TIME_EQ(cur_dt, *dt) && nframes > 1) {
gk_log_fatal(FARGS, "Inconsistent time step in frame %d of %s: dt is %f vs. previous dt of %f.\n",
nframes, traj_fname, cur_dt, *dt);
}
else if(nframes == 1) {
// Use the dt between the first two frames of the trajectory as the expected dt.
*dt = cur_dt;
}
// Expand memory if needed
if(nframes + 1 > cur_alloc) {
cur_alloc *= 2;
srenew(*unit_vecs, cur_alloc);
}
// Get unit vectors for the atom pairs in this frame
snew((*unit_vecs)[nframes], npairs);
for(int p = 0; p < npairs; ++p) {
gc_get_unit_vec(x, atompairs[2 * p], atompairs[2 * p + 1], (*unit_vecs)[nframes][p]);
}
++nframes;
} while(read_next_x(*oenv, status, &t,
#ifndef GRO_V5
natoms,
#endif
x, box));
} // if dt < 0
else { // if provided dt >= 0
// Only store the data in the trajectory intervals matching the given dt
last_t = t - *dt; // so that the first trajectory frame will pass dt boundary check
int trajframe = 0;
do {
cur_dt = t - last_t;
// if this trajectory frame is not on a dt boundary, skip it
if(GC_TIME_EQ(cur_dt, *dt)) {
// Expand memory if needed
if(nframes + 1 > cur_alloc) {
cur_alloc *= 2;
srenew(*unit_vecs, cur_alloc);
}
// Get unit vectors for the atom pairs in this frame
snew((*unit_vecs)[nframes], npairs);
for(int p = 0; p < npairs; ++p) {
gc_get_unit_vec(x, atompairs[2 * p], atompairs[2 * p + 1], (*unit_vecs)[nframes][p]);
}
// Reset timestep counter
last_t = t;
++nframes;
}
else if(cur_dt > *dt) {
gk_log_fatal(FARGS, "Error at frame %d of %s: "
"given dt = %f is not a multiple of the trajectory timestep, "
"or the trajectory has inconsistent timesteps.\n",
trajframe, traj_fname, *dt);
}
++trajframe;
} while(read_next_x(*oenv, status, &t,
#ifndef GRO_V5
natoms,
#endif
x, box));
} // if provided dt >= 0
} // if natoms > 0
else {
gk_log_fatal(FARGS, "No atoms found in %s!\n", traj_fname);
}
// Done with trajectory
close_trx(status);
sfree(x);
srenew(*unit_vecs, nframes); // free excess memory
return nframes;
}
开发者ID:luky1971,项目名称:g_correlate,代码行数:101,代码来源:correlate.c
示例16: gmx_trjcat
//.........这里部分代码省略.........
* we let the time start at zero, see the edit_files routine.
*/
bNewFile = TRUE;
if (!lastTimeSet)
{
lasttime = 0;
lastTimeSet = true;
}
printf("\n");
printf("lasttime %g\n", lasttime);
do
{
/* copy the input frame to the output frame */
frout = fr;
/* set the new time by adding the correct calculated above */
frout.time += t_corr;
if (ftpout == efTNG)
{
frout.step += prevEndStep;
}
/* quit if we have reached the end of what should be written */
if ((end > 0) && (frout.time > end+GMX_REAL_EPS))
{
i = nfile_in;
break;
}
/* determine if we should write this frame (dt is handled elsewhere) */
if (bCat) /* write all frames of all files */
{
bWrite = TRUE;
|
请发表评论