From d657191b2ffeba74ed2ed06e4289188014bf4e23 Mon Sep 17 00:00:00 2001 From: David van der Spoel Date: Fri, 2 Jul 2010 14:13:01 +0200 Subject: [PATCH] Made removing pbc non-static. This affects many analysis tools and small parts of the trjana library. --- include/rmpbc.h | 4 ++-- src/gmxlib/mshift.c | 51 ++++++++++++++++++++++++++------------------ src/gmxlib/nrama.c | 7 ++++-- src/gmxlib/rmpbc.c | 10 --------- src/gmxlib/trajana/trajana.c | 16 +++++++++----- src/ngmx/manager.c | 4 ++-- src/ngmx/manager.h | 4 +++- src/tools/do_dssp.c | 9 +++++--- src/tools/gmx_anaeig.c | 22 +++++++++++++++---- src/tools/gmx_bundle.c | 8 +++++-- src/tools/gmx_covar.c | 14 +++++++----- src/tools/gmx_current.c | 12 +++++++---- src/tools/gmx_density.c | 14 ++++++++---- src/tools/gmx_dipoles.c | 6 +++++- src/tools/gmx_disre.c | 10 +++++++-- src/tools/gmx_dist.c | 5 ++++- src/tools/gmx_filter.c | 10 ++++++--- src/tools/gmx_gyrate.c | 7 +++++- src/tools/gmx_h2order.c | 7 ++++-- src/tools/gmx_helix.c | 8 ++++++- src/tools/gmx_helixorient.c | 9 ++++++-- src/tools/gmx_mdmat.c | 7 +++++- src/tools/gmx_mindist.c | 15 +++++++++---- src/tools/gmx_msd.c | 18 +++++++++++----- src/tools/gmx_multipoles.c | 11 +++++++--- src/tools/gmx_order.c | 27 ++++++++++++++--------- src/tools/gmx_polystat.c | 8 ++++++- src/tools/gmx_potential.c | 11 +++++++--- src/tools/gmx_principal.c | 10 ++++++++- src/tools/gmx_rdf.c | 10 +++++++-- src/tools/gmx_relax.c | 7 +++++- src/tools/gmx_rms.c | 13 +++++++---- src/tools/gmx_rmsf.c | 15 +++++++++---- src/tools/gmx_rotacf.c | 10 +++++++-- src/tools/gmx_rotmat.c | 17 ++++++++++++--- src/tools/gmx_sas.c | 11 ++++++++-- src/tools/gmx_sdf.c | 12 +++++++---- src/tools/gmx_sgangle.c | 16 +++++++++++--- src/tools/gmx_spatial.c | 9 +++++++- src/tools/gmx_spol.c | 10 +++++++-- src/tools/gmx_traj.c | 6 +++++- src/tools/gmx_trjconv.c | 14 ++++++++---- src/tools/gmx_trjorder.c | 7 ++++-- 43 files changed, 360 insertions(+), 141 deletions(-) diff --git a/include/rmpbc.h b/include/rmpbc.h index 74bb24ac9d..006f2c01b6 100644 --- a/include/rmpbc.h +++ b/include/rmpbc.h @@ -64,8 +64,8 @@ extern "C" { * When ePBC=-1, the type of pbc is guessed from the box matrix. */ - extern void rm_pbc(t_idef *idef,int ePBC,int natoms, - matrix box,rvec x[],rvec x_s[]); + /*extern void rm_pbc(t_idef *idef,int ePBC,int natoms, + matrix box,rvec x[],rvec x_s[]);*/ /* Convenience function that still holds a static variable. */ extern void rm_gropbc(t_atoms *atoms,rvec x[],matrix box); diff --git a/src/gmxlib/mshift.c b/src/gmxlib/mshift.c index b92d72880c..b650767b75 100644 --- a/src/gmxlib/mshift.c +++ b/src/gmxlib/mshift.c @@ -378,6 +378,9 @@ void mk_graph_ilist(FILE *fplog, g->nbound++; } + g->negc = 0; + g->egc = NULL; + sfree(nbond); if (gmx_debug_at) @@ -506,7 +509,7 @@ static void mk_1shift_screw(matrix box,rvec hbox, } } -static int mk_grey(FILE *log,t_graph *g,int *AtomI, +static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI, int npbcdim,matrix box,rvec x[],int *nerror) { int m,j,ng,ai,aj,g0; @@ -521,13 +524,11 @@ static int mk_grey(FILE *log,t_graph *g,int *AtomI, ng=0; ai=*AtomI; - range_check(ai,0,g->negc); g0=g->start; /* Loop over all the bonds */ for(j=0; (jnedge[ai]); j++) { aj=g->edge[ai][j]-g0; - range_check(aj,0,g->negc); /* If there is a white one, make it grey and set pbc */ if (g->bScrewPBC) mk_1shift_screw(box,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj); @@ -536,10 +537,10 @@ static int mk_grey(FILE *log,t_graph *g,int *AtomI, else mk_1shift(npbcdim,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj); - if (g->egc[aj] == egcolWhite) { + if (egc[aj] == egcolWhite) { if (aj < *AtomI) *AtomI = aj; - g->egc[aj] = egcolGrey; + egc[aj] = egcolGrey; copy_ivec(is_aj,g->ishift[aj]); @@ -564,7 +565,7 @@ static int mk_grey(FILE *log,t_graph *g,int *AtomI, return ng; } -static int first_colour(int fC,egCol Col,t_graph *g) +static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[]) /* Return the first node with colour Col starting at fC. * return -1 if none found. */ @@ -572,7 +573,7 @@ static int first_colour(int fC,egCol Col,t_graph *g) int i; for(i=fC; (innodes); i++) - if ((g->nedge[i] > 0) && (g->egc[i]==Col)) + if ((g->nedge[i] > 0) && (egc[i]==Col)) return i; return -1; @@ -580,13 +581,13 @@ static int first_colour(int fC,egCol Col,t_graph *g) void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[]) { + static int nerror_tot = 0; int npbcdim; - int ng,i; + int ng,nnodes,i; int nW,nG,nB; /* Number of Grey, Black, White */ int fW,fG; /* First of each category */ int nerror=0; - printf("Start. g->nnodes = %d\n",g->nnodes); g->bScrewPBC = (ePBC == epbcSCREW); if (ePBC == epbcXY) @@ -606,12 +607,12 @@ void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[]) if (!g->nbound) return; - if (g->nnodes > g->negc) { - g->negc = g->nnodes; + nnodes=g->nnodes; + if (nnodes > g->negc) { + g->negc = nnodes; srenew(g->egc,g->negc); } - for(i=0; (innodes); i++) - g->egc[i] = egcolWhite; + memset(g->egc,0,(size_t)(nnodes*sizeof(g->egc[0]))); nW=g->nbound; nG=0; @@ -630,7 +631,7 @@ void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[]) * number than before, because no nodes are made white * in the loop */ - if ((fW=first_colour(fW,egcolWhite,g)) == -1) + if ((fW=first_colour(fW,egcolWhite,g,g->egc)) == -1) gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW); /* Make the first white node grey */ @@ -645,7 +646,7 @@ void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[]) nW,nG,nB,nW+nG+nB); #endif while (nG > 0) { - if ((fG=first_colour(fG,egcolGrey,g)) == -1) + if ((fG=first_colour(fG,egcolGrey,g,g->egc)) == -1) gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG); /* Make the first grey node black */ @@ -656,21 +657,29 @@ void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[]) /* Make all the neighbours of this black node grey * and set their periodicity */ - ng=mk_grey(log,g,&fG,npbcdim,box,x,&nerror); + ng=mk_grey(log,nnodes,g->egc,g,&fG,npbcdim,box,x,&nerror); /* ng is the number of white nodes made grey */ nG+=ng; nW-=ng; } } if (nerror > 0) { - fprintf(stderr,"There were %d inconsistent shifts. Check your topology\n", - nerror); - if (log) { - fprintf(log,"There were %d inconsistent shifts. Check your topology\n", + nerror_tot++; + if (nerror_tot <= 100) { + fprintf(stderr,"There were %d inconsistent shifts. Check your topology\n", nerror); + if (log) { + fprintf(log,"There were %d inconsistent shifts. Check your topology\n", + nerror); + } + } + if (nerror_tot == 100) { + fprintf(stderr,"Will stop reporting inconsistent shifts\n"); + if (log) { + fprintf(log,"Will stop reporting inconsistent shifts\n"); + } } } - printf("End. g->nnodes = %d\n",g->nnodes); } /************************************************************ diff --git a/src/gmxlib/nrama.c b/src/gmxlib/nrama.c index 50e2a9266f..05143f0e56 100644 --- a/src/gmxlib/nrama.c +++ b/src/gmxlib/nrama.c @@ -73,8 +73,11 @@ static void calc_dihs(t_xrama *xr) rvec r_ij,r_kj,r_kl,m,n; real sign; t_dih *dd; - - rm_pbc(xr->idef,xr->ePBC,xr->natoms,xr->box,xr->x,xr->x); + gmx_rmpbc_t gpbc=NULL; + + gpbc = gmx_rmpbc_init(xr->idef,xr->ePBC,xr->natoms,xr->box); + gmx_rmpbc(gpbc,xr->box,xr->x,xr->x); + gmx_rmpbc_done(gpbc); for(i=0; (indih); i++) { dd=&(xr->dih[i]); diff --git a/src/gmxlib/rmpbc.c b/src/gmxlib/rmpbc.c index 6e9fde19e7..33e543b410 100644 --- a/src/gmxlib/rmpbc.c +++ b/src/gmxlib/rmpbc.c @@ -92,16 +92,6 @@ void gmx_rmpbc(gmx_rmpbc_t gpbc,matrix box,rvec x[],rvec x_s[]) copy_rvec(x[i],x_s[i]); } -void rm_pbc(t_idef *idef,int ePBC,int natoms, - matrix box,rvec x[],rvec x_s[]) -{ - static gmx_rmpbc_t gpbc=NULL; - - if (NULL == gpbc) - gpbc = gmx_rmpbc_init(idef,ePBC,natoms,box); - gmx_rmpbc(gpbc,box,x,x_s); -} - void rm_gropbc(t_atoms *atoms,rvec x[],matrix box) { real dist; diff --git a/src/gmxlib/trajana/trajana.c b/src/gmxlib/trajana/trajana.c index d7bfc6e9f7..6f870821d9 100644 --- a/src/gmxlib/trajana/trajana.c +++ b/src/gmxlib/trajana/trajana.c @@ -1614,7 +1614,8 @@ int gmx_ana_do(gmx_ana_traj_t *d, int flags, gmx_analysisfunc analyze, void *dat t_pbc pbc; t_pbc *ppbc; int rc; - + gmx_rmpbc_t gpbc=NULL; + rc = init_first_frame(d); if (rc != 0) { @@ -1626,14 +1627,16 @@ int gmx_ana_do(gmx_ana_traj_t *d, int flags, gmx_analysisfunc analyze, void *dat { d->bRmPBC = FALSE; } - + if (d->bRmPBC) + { + gpbc = gmx_rmpbc_init(&d->top->idef,d->ePBC,d->fr->natoms,d->fr->box); + } d->nframes = 0; do { if (d->bRmPBC) { - rm_pbc(&d->top->idef, d->ePBC, d->fr->natoms, d->fr->box, - d->fr->x, d->fr->x); + gmx_rmpbc(gpbc,d->fr->box,d->fr->x, d->fr->x); } if (ppbc) { @@ -1658,7 +1661,10 @@ int gmx_ana_do(gmx_ana_traj_t *d, int flags, gmx_analysisfunc analyze, void *dat d->nframes++; } while (d->trjfile && read_next_frame(d->oenv, d->status, d->fr)); - + if (d->bRmPBC) + { + gmx_rmpbc_done(gpbc); + } if (d->trjfile) { close_trj(d->status); diff --git a/src/ngmx/manager.c b/src/ngmx/manager.c index 92839b8d25..29fd19f07a 100644 --- a/src/ngmx/manager.c +++ b/src/ngmx/manager.c @@ -198,6 +198,7 @@ void set_file(t_x11 *x11,t_manager *man,const char *trajectory, snew(man->bHydro,sh.natoms); snew(bB,sh.natoms); read_tpx_top(status,NULL,man->box,&man->natom,NULL,NULL,NULL,&man->top); + man->gpbc = gmx_rmpbc_init(&man->top.idef,-1,man->natom,man->box); man->natom= read_first_x(man->oenv,&man->status,trajectory,&(man->time),&(man->x), @@ -337,8 +338,7 @@ static bool step_man(t_manager *man,int *nat) break; } if (man->bPbc) { - rm_pbc(&(man->top.idef),man->molw->ePBC, - man->natom,man->box,man->x,man->x); + gmx_rmpbc(man->gpbc,man->box,man->x,man->x); reset_mols(&(man->top.mols),man->box,man->x); } ncount=0; diff --git a/src/ngmx/manager.h b/src/ngmx/manager.h index 4f961e192a..0fc964c4b1 100644 --- a/src/ngmx/manager.h +++ b/src/ngmx/manager.h @@ -44,6 +44,7 @@ #include "nleg.h" #include "buttons.h" #include "statutil.h" +#include "rmpbc.h" /* Some window sizes */ #define EWIDTH 200 @@ -119,7 +120,8 @@ typedef struct { bool bPlus; /* Draw plus for single atom */ int nSkip; /* Skip n steps after each frame */ int nWait; /* Wait n ms after each frame */ - + gmx_rmpbc_t gpbc; /* For removing peridiocity */ + t_windata wd; /* The manager subwindow */ t_windata title; /* Title window */ t_3dview *view; /* The 3d struct */ diff --git a/src/tools/do_dssp.c b/src/tools/do_dssp.c index a20fdee4b5..77c3d4802c 100644 --- a/src/tools/do_dssp.c +++ b/src/tools/do_dssp.c @@ -399,7 +399,7 @@ int main(int argc,char *argv[]) char dssp[256]; const char *dptr; output_env_t oenv; - + gmx_rmpbc_t gpbc=NULL; t_filenm fnm[] = { { efTRX, "-f", NULL, ffREAD }, @@ -495,6 +495,8 @@ int main(int argc,char *argv[]) snew(norm_av_area, atoms->nres); accr=NULL; naccr=0; + + gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box); do { t = output_env_conv_time(oenv,t); if (bDoAccSurf && nframe>=naccr) { @@ -503,7 +505,7 @@ int main(int argc,char *argv[]) for(i=naccr-10; inres-1); } - rm_pbc(&(top.idef),ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); tapein=ffopen(pdbfile,"w"); write_pdbfile_indexed(tapein,NULL,atoms,x,ePBC,box,0,-1,gnx,index,NULL); ffclose(tapein); @@ -535,7 +537,8 @@ int main(int argc,char *argv[]) close_trj(status); if (fTArea) ffclose(fTArea); - + gmx_rmpbc_done(gpbc); + prune_ss_legend(&mat); ss=opt2FILE("-o",NFILE,fnm,"w"); diff --git a/src/tools/gmx_anaeig.c b/src/tools/gmx_anaeig.c index 33022bd14f..f506b1bba1 100644 --- a/src/tools/gmx_anaeig.c +++ b/src/tools/gmx_anaeig.c @@ -415,6 +415,7 @@ static void project(const char *trajfile,t_topology *top,int ePBC,matrix topbox, real t,inp,**inprod=NULL,min=0,max=0; char str[STRLEN],str2[STRLEN],**ylabel,*c; real fact; + gmx_rmpbc_t gpbc=NULL; snew(x,natoms); @@ -442,12 +443,18 @@ static void project(const char *trajfile,t_topology *top,int ePBC,matrix topbox, if (nat>atoms->nr) gmx_fatal(FARGS,"the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)",nat,atoms->nr); snew(all_at,nat); + + if (top) + gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box); + if (top) + gmx_rmpbc_done(gpbc); + for(i=0; iidef),ePBC,nat,box,xread,xread); + gmx_rmpbc(gpbc,box,xread,xread); if (nframes>=snew_size) { snew_size+=100; for(i=0; inr); + if (top) + gmx_rmpbc_done(gpbc); + if (projfile) { snew(ylabel,noutvec); @@ -866,7 +876,8 @@ int gmx_anaeig(int argc,char *argv[]) int neig1,neig2; double **xvgdata; output_env_t oenv; - + gmx_rmpbc_t gpbc=NULL; + t_filenm fnm[] = { { efTRN, "-v", "eigenvec", ffREAD }, { efTRN, "-v2", "eigenvec2", ffOPTRD }, @@ -988,13 +999,15 @@ int gmx_anaeig(int argc,char *argv[]) nfit=0; ifit=NULL; w_rls=NULL; + gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,topbox); + if (!bTPS) bTop=FALSE; else { bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm), title,&top,&ePBC,&xtop,NULL,topbox,bM); atoms=&top.atoms; - rm_pbc(&(top.idef),ePBC,atoms->nr,topbox,xtop,xtop); + gmx_rmpbc(gpbc,topbox,xtop,xtop); /* Fitting is only required for the projection */ if (bProj && bFit1) { if (xref1 == NULL) { @@ -1027,7 +1040,8 @@ int gmx_anaeig(int argc,char *argv[]) } } } - + gmx_rmpbc_done(gpbc); + if (bIndex) { printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms); get_index(atoms,indexfile,1,&i,&index,&grpname); diff --git a/src/tools/gmx_bundle.c b/src/tools/gmx_bundle.c index 5ea38694f7..0772cec7bd 100644 --- a/src/tools/gmx_bundle.c +++ b/src/tools/gmx_bundle.c @@ -218,6 +218,8 @@ int gmx_bundle(int argc,char *argv[]) bool bKink; rvec va,vb,vc,vr,vl; output_env_t oenv; + gmx_rmpbc_t gpbc=NULL; + #define NLEG asize(leg) t_filenm fnm[] = { { efTRX, "-f", NULL, ffREAD }, @@ -307,9 +309,10 @@ int gmx_bundle(int argc,char *argv[]) fpdb = NULL; read_first_frame(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_NEED_X); - + gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box); + do { - rm_pbc(&top.idef,ePBC,fr.natoms,fr.box,fr.x,fr.x); + gmx_rmpbc(gpbc,fr.box,fr.x,fr.x); calc_axes(fr.x,top.atoms.atom,gnx,index,!bZ,&bun); t = output_env_conv_time(oenv,fr.time); fprintf(flen," %10g",t); @@ -366,6 +369,7 @@ int gmx_bundle(int argc,char *argv[]) if (fpdb ) dump_axes(fpdb,&fr,&outatoms,&bun); } while(read_next_frame(oenv,status,&fr)); + gmx_rmpbc_done(gpbc); close_trx(status); diff --git a/src/tools/gmx_covar.c b/src/tools/gmx_covar.c index f0cbccd533..78da9e032b 100644 --- a/src/tools/gmx_covar.c +++ b/src/tools/gmx_covar.c @@ -138,6 +138,7 @@ int gmx_covar(int argc,char *argv[]) t_rgb rlo,rmi,rhi; real *tmp; output_env_t oenv; + gmx_rmpbc_t gpbc=NULL; t_filenm fnm[] = { { efTRX, "-f", NULL, ffREAD }, @@ -217,10 +218,12 @@ int gmx_covar(int argc,char *argv[]) w_rls[ifit[i]]=1.0; } } - + /* Prepare reference frame */ - if (bPBC) - rm_pbc(&(top.idef),ePBC,atoms->nr,box,xref,xref); + if (bPBC) { + gpbc = gmx_rmpbc_init(&top.idef,ePBC,atoms->nr,box); + gmx_rmpbc(gpbc,box,xref,xref); + } if (bFit) reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls); @@ -241,7 +244,7 @@ int gmx_covar(int argc,char *argv[]) nframes0++; /* calculate x: a fitted struture of the selected atoms */ if (bPBC) - rm_pbc(&(top.idef),ePBC,nat,box,xread,xread); + gmx_rmpbc(gpbc,box,xread,xread); if (bFit) { reset_x(nfit,ifit,nat,NULL,xread,w_rls); do_fit(nat,w_rls,xref,xread); @@ -270,7 +273,7 @@ int gmx_covar(int argc,char *argv[]) tend = t; /* calculate x: a (fitted) structure of the selected atoms */ if (bPBC) - rm_pbc(&(top.idef),ePBC,nat,box,xread,xread); + gmx_rmpbc(gpbc,box,xread,xread); if (bFit) { reset_x(nfit,ifit,nat,NULL,xread,w_rls); do_fit(nat,w_rls,xref,xread); @@ -296,6 +299,7 @@ int gmx_covar(int argc,char *argv[]) } while (read_next_x(oenv,status,&t,nat,xread,box) && (bRef || nframes < nframes0)); close_trj(status); + gmx_rmpbc_done(gpbc); fprintf(stderr,"Read %d frames\n",nframes); diff --git a/src/tools/gmx_current.c b/src/tools/gmx_current.c index a76dfbad29..8f67f7f550 100644 --- a/src/tools/gmx_current.c +++ b/src/tools/gmx_current.c @@ -344,7 +344,8 @@ static void dielectric(FILE *fmj,FILE *fmd,FILE *outf,FILE *fcur,FILE *mcor, real err=0.0; real *xfit; real *yfit; - + gmx_rmpbc_t gpbc=NULL; + /* * indices for EH */ @@ -377,7 +378,8 @@ static void dielectric(FILE *fmj,FILE *fmd,FILE *outf,FILE *fcur,FILE *mcor, clear_rvec(mjd_tmp); clear_rvec(mdvec); clear_rvec(tmp); - + gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box); + do{ refr=(real)(nfr+1); @@ -426,8 +428,8 @@ static void dielectric(FILE *fmj,FILE *fmd,FILE *outf,FILE *fcur,FILE *mcor, } - rm_pbc(&top.idef,ePBC,fr.natoms,fr.box,fr.x,fr.x); - + gmx_rmpbc(gpbc,fr.box,fr.x,fr.x); + calc_mj(top,ePBC,fr.box,bNoJump,nmols,indexm,fr.x,mtrans[nfr],mass2,qmol); for(i=0;iidef,ePBC,top->atoms.nr,box); /*********** Start processing trajectory ***********/ do { - rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x0,x0); + gmx_rmpbc(gpbc,box,x0,x0); if (bCenter) center_coords(&top->atoms,box,x0,axis); @@ -224,7 +226,8 @@ void calc_electron_density(const char *fn, atom_id **index, int gnx[], } nr_frames++; } while (read_next_x(oenv,status,&t,natoms,x0,box)); - + gmx_rmpbc_done(gpbc); + /*********** done with status file **********/ close_trj(status); @@ -262,6 +265,7 @@ void calc_density(const char *fn, atom_id **index, int gnx[], real t, z; char *buf; /* for tmp. keeping atomname */ + gmx_rmpbc_t gpbc=NULL; switch(axis) { case 0: @@ -289,9 +293,10 @@ void calc_density(const char *fn, atom_id **index, int gnx[], for (i = 0; i < nr_grps; i++) snew((*slDensity)[i], *nslices); + gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box); /*********** Start processing trajectory ***********/ do { - rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x0,x0); + gmx_rmpbc(gpbc,box,x0,x0); if (bCenter) center_coords(&top->atoms,box,x0,axis); @@ -315,7 +320,8 @@ void calc_density(const char *fn, atom_id **index, int gnx[], nr_frames++; } while (read_next_x(oenv,status,&t,natoms,x0,box)); - + gmx_rmpbc_done(gpbc); + /*********** done with status file **********/ close_trj(status); diff --git a/src/tools/gmx_dipoles.c b/src/tools/gmx_dipoles.c index 6e73fdf775..2224c00c7f 100644 --- a/src/tools/gmx_dipoles.c +++ b/src/tools/gmx_dipoles.c @@ -691,6 +691,7 @@ static void do_dip(t_topology *top,int ePBC,real volume, rvec *slab_dipoles=NULL; t_atom *atom=NULL; t_block *mols=NULL; + gmx_rmpbc_t gpbc=NULL; gnx_tot = gnx[0]; if (ncos > 1) { @@ -848,6 +849,7 @@ static void do_dip(t_topology *top,int ePBC,real volume, gkrbin = mk_gkrbin(rcut,rcmax,bPhi,ndegrees); } + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natom,box); /* Start while loop over frames */ t1 = t0 = t; @@ -876,7 +878,7 @@ static void do_dip(t_topology *top,int ePBC,real volume, M_av[m] = 0; M_av2[m] = 0; } - rm_pbc(&(top->idef),ePBC,natom,box,x,x); + gmx_rmpbc(gpbc,box,x,x); muframelsq = gmx_stats_init(); /* Begin loop of all molecules in frame */ @@ -1070,6 +1072,8 @@ static void do_dip(t_topology *top,int ePBC,real volume, bCont = read_next_x(oenv,status,&t,natom,x,box); } while (bCont); + gmx_rmpbc_done(gpbc); + if (!bMU) close_trj(status); diff --git a/src/tools/gmx_disre.c b/src/tools/gmx_disre.c index 62c013a753..e67d3c1450 100644 --- a/src/tools/gmx_disre.c +++ b/src/tools/gmx_disre.c @@ -606,7 +606,8 @@ int gmx_disre(int argc,char *argv[]) int my_clust; FILE *fplog; output_env_t oenv; - + gmx_rmpbc_t gpbc=NULL; + t_filenm fnm[] = { { efTPX, NULL, NULL, ffREAD }, { efTRX, "-f", NULL, ffREAD }, @@ -713,13 +714,16 @@ int gmx_disre(int argc,char *argv[]) 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 - rm_pbc(&top->idef,ir.ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); } if (clust) { @@ -761,6 +765,8 @@ int gmx_disre(int argc,char *argv[]) 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]), diff --git a/src/tools/gmx_dist.c b/src/tools/gmx_dist.c index 12726592ae..ca23546dbc 100644 --- a/src/tools/gmx_dist.c +++ b/src/tools/gmx_dist.c @@ -105,6 +105,7 @@ int gmx_dist(int argc,char *argv[]) int *contact_time=NULL,*ccount=NULL,ccount_nalloc=0,sum; char buf[STRLEN]; output_env_t oenv; + gmx_rmpbc_t gpbc=NULL; char *leg[4] = { "|d|","d\\sx\\N","d\\sy\\N","d\\sz\\N" }; @@ -181,12 +182,13 @@ int gmx_dist(int argc,char *argv[]) else 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 */ - rm_pbc(&top->idef,ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); } /* calculate center of masses */ for(g=0;(gidef,ePBC,top->atoms.nr,box); /*********** Start processing trajectory ***********/ do { *slWidth = box[axis][axis]/(*nslices); teller++; - rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x0,x0); + gmx_rmpbc(gpbc,box,x0,x0); if (bMicel) calc_xcm(x0, nmic, micel, top->atoms.atom, com, FALSE); @@ -197,6 +199,7 @@ void calc_h2order(const char *fn, atom_id index[], int ngx, rvec **slDipole, /*********** done with status file **********/ fprintf(stderr,"\nRead trajectory. Printing parameters to file\n"); + gmx_rmpbc_done(gpbc); for (i = 0; i < *nslices; i++) /* average over frames */ { diff --git a/src/tools/gmx_helix.c b/src/tools/gmx_helix.c index 1683c57dbd..92e207906e 100644 --- a/src/tools/gmx_helix.c +++ b/src/tools/gmx_helix.c @@ -222,6 +222,7 @@ int gmx_helix(int argc,char *argv[]) real t; real rms,fac; matrix box; + gmx_rmpbc_t gpbc=NULL; bool bRange; t_filenm fnm[] = { { efTPX, NULL, NULL, ffREAD }, @@ -289,12 +290,15 @@ int gmx_helix(int argc,char *argv[]) pr_bb(stdout,nres,bb); } + gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box); + snew(xav,natoms); teller=0; do { if ((teller++ % 10) == 0) fprintf(stderr,"\rt=%.2f",t); - rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x,x); + gmx_rmpbc(gpbc,box,x,x); + calc_hxprops(nres,bb,x,box); if (bCheck) @@ -334,6 +338,8 @@ int gmx_helix(int argc,char *argv[]) } while (read_next_x(oenv,status,&t,natoms,x,box)); fprintf(stderr,"\n"); + gmx_rmpbc_done(gpbc); + close_trj(status); if (otrj) { diff --git a/src/tools/gmx_helixorient.c b/src/tools/gmx_helixorient.c index 4768f8a780..09280de96d 100644 --- a/src/tools/gmx_helixorient.c +++ b/src/tools/gmx_helixorient.c @@ -135,6 +135,7 @@ int gmx_helixorient(int argc,char *argv[]) int ePBC; output_env_t oenv; + gmx_rmpbc_t gpbc=NULL; static bool bSC=FALSE; static bool bIncremental = FALSE; @@ -249,12 +250,14 @@ int gmx_helixorient(int argc,char *argv[]) unitaxes[1][1]=1; unitaxes[2][2]=1; + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); + do { /* initialisation for correct distance calculations */ set_pbc(&pbc,ePBC,box); /* make molecules whole again */ - rm_pbc(&top->idef,ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); /* copy coords to our smaller arrays */ for(i=0;iidef,ePBC,natoms,box); + bFirst=TRUE; do { - if (top) { - rm_pbc(&(top->idef),ePBC,natoms,box,x,x); - } + if (NULL != top) + gmx_rmpbc(gpbc,box,x,x); + periodic_dist(box,x,n,index,&rmin,&rmax,ind_min); if (rmin < rmint) { rmint = rmin; @@ -147,6 +151,9 @@ static void periodic_mindist_plot(const char *trxfn,const char *outfn, output_env_conv_time(oenv,t),rmin,rmax,norm(box[0]),norm(box[1]),norm(box[2])); bFirst=FALSE; } while(read_next_x(oenv,status,&t,natoms,x,box)); + + if (NULL != top) + gmx_rmpbc_done(gpbc); ffclose(out); diff --git a/src/tools/gmx_msd.c b/src/tools/gmx_msd.c index 666ebfcead..357b5cd31c 100644 --- a/src/tools/gmx_msd.c +++ b/src/tools/gmx_msd.c @@ -569,6 +569,7 @@ int corr_loop(t_corr *curr,const char *fn,t_topology *top,int ePBC, #define prev (1-cur) matrix box; bool bFirst; + gmx_rmpbc_t gpbc=NULL; natoms = read_first_x(oenv,&status,fn,&curr->t0,&(x[cur]),box); #ifdef DEBUG @@ -595,6 +596,9 @@ int corr_loop(t_corr *curr,const char *fn,t_topology *top,int ePBC, if (x_pdb) *x_pdb = NULL; + if (bMol) + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); + /* the loop over all frames */ do { @@ -666,7 +670,7 @@ int corr_loop(t_corr *curr,const char *fn,t_topology *top,int ePBC, /* make the molecules whole */ if (bMol) - rm_pbc(&top->idef,ePBC,natoms,box,x[cur],x[cur]); + gmx_rmpbc(gpbc,box,x[cur],x[cur]); /* first remove the periodic boundary condition crossings */ for(i=0;ingrp;i++) @@ -701,8 +705,12 @@ int corr_loop(t_corr *curr,const char *fn,t_topology *top,int ePBC, fprintf(stderr,"\nUsed %d restart points spaced %g %s over %g %s\n\n", curr->nrestart, output_env_conv_time(oenv,dt), output_env_get_time_unit(oenv), - output_env_conv_time(oenv,curr->time[curr->nframes-1]), output_env_get_time_unit(oenv) ); + output_env_conv_time(oenv,curr->time[curr->nframes-1]), + output_env_get_time_unit(oenv) ); + if (bMol) + gmx_rmpbc_done(gpbc); + close_trj(status); return natoms; @@ -753,7 +761,7 @@ void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file, int *gnx_com=NULL; /* the COM removal group size */ atom_id **index_com=NULL; /* the COM removal group atom indices */ char **grpname_com=NULL; /* the COM removal group name */ - + snew(gnx,nrgrp); snew(index,nrgrp); snew(grpname,nrgrp); @@ -770,7 +778,7 @@ void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file, fprintf(stderr, "\nNow select a group for center of mass removal:\n"); get_index(&top->atoms, ndx_file, 1, gnx_com, index_com, grpname_com); } - + if (mol_file) index_atom2mol(&gnx[0],index[0],&top->mols); @@ -783,7 +791,7 @@ void do_corr(const char *trx_file, const char *ndx_file, const char *msd_file, (mol_file!=NULL) ? calc1_mol : (bMW ? calc1_mw : calc1_norm), bTen, gnx_com, index_com, dt,t_pdb, pdb_file ? &x : NULL,box,oenv); - + /* Correct for the number of points */ for(j=0; (jngrp); j++) { for(i=0; (inframes); i++) { diff --git a/src/tools/gmx_multipoles.c b/src/tools/gmx_multipoles.c index 9b00e15279..4daa01b3ee 100644 --- a/src/tools/gmx_multipoles.c +++ b/src/tools/gmx_multipoles.c @@ -628,6 +628,7 @@ void do_multipoles(char *trjfn,char *topfn,char *molndxfn,bool bFull) tensor3 *m3; tensor4 *m4; matrix trans; + gmx_rmpbc_t gpbc=NULL; top = read_top(topfn); rd_index(molndxfn,1,&gnx,&grpindex,&grpname); @@ -639,12 +640,13 @@ void do_multipoles(char *trjfn,char *topfn,char *molndxfn,bool bFull) snew(m2,gnx); snew(m3,gnx); snew(m4,gnx); - + + gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box); + /* Start while loop over frames */ do { /* PvM, bug in rm_pbc??? Does not work for virtual sites ... - rm_pbc(&(top->idef),top->atoms.nr,box,x,x_s); */ - + gmx_rmpbc(gpbc,box,x,x_s); */ /* Begin loop of all molecules in index file */ for(i=0; (iidef,ePBC,top->atoms.nr,box); /*********** Start processing trajectory ***********/ do { if (bSliced) *slWidth = box[axis][axis]/nslices; teller++; - set_pbc(&pbc,ePBC,box); - rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x0,x1); + set_pbc(&pbc,ePBC,box); + gmx_rmpbc(gpbc,box,x0,x1); /* Now loop over all groups. There are ngrps groups, the order parameter can be calculated for grp 1 to grp ngrps - 1. For each group, loop over all @@ -608,7 +614,8 @@ void calc_order(const char *fn, atom_id *index, atom_id *a, rvec **order, /*********** 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]); diff --git a/src/tools/gmx_polystat.c b/src/tools/gmx_polystat.c index 5403c532a5..a310644c8a 100644 --- a/src/tools/gmx_polystat.c +++ b/src/tools/gmx_polystat.c @@ -162,6 +162,7 @@ int gmx_polystat(int argc,char *argv[]) " eig1", " eig2", " eig3", "", "", "" }; char **legp,buf[STRLEN]; + gmx_rmpbc_t gpbc=NULL; CopyRight(stderr,argv[0]); parse_common_args(&argc,argv, @@ -253,8 +254,11 @@ int gmx_polystat(int argc,char *argv[]) sum_eed2_tot = 0; sum_gyro_tot = 0; sum_pers_tot = 0; + + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); + do { - rm_pbc(&(top->idef),ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); sum_eed2 = 0; for(d=0; didef,ePBC,top->atoms.nr,box); + /*********** Start processing trajectory ***********/ do { *slWidth = box[axis][axis]/(*nslices); teller++; - - rm_pbc(&(top->idef),ePBC,top->atoms.nr,box,x0,x0); + gmx_rmpbc(gpbc,box,x0,x0); /* calculate position of center of mass based on group 1 */ calc_xcm(x0, gnx[0], index[0], top->atoms.atom, xcm, FALSE); @@ -201,6 +204,8 @@ void calc_potential(const char *fn, atom_id **index, int gnx[], nr_frames++; } while (read_next_x(oenv,status,&t,natoms,x0,box)); + gmx_rmpbc_done(gpbc); + /*********** done with status file **********/ close_trj(status); diff --git a/src/tools/gmx_principal.c b/src/tools/gmx_principal.c index 4746daf110..6721307a24 100644 --- a/src/tools/gmx_principal.c +++ b/src/tools/gmx_principal.c @@ -102,6 +102,8 @@ int gmx_principal(int argc,char *argv[]) FILE * fmoi; matrix axes,box; output_env_t oenv; + gmx_rmpbc_t gpbc=NULL; + t_filenm fnm[] = { { efTRX, "-f", NULL, ffREAD }, @@ -130,9 +132,12 @@ int gmx_principal(int argc,char *argv[]) natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box); + gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box); + do { - rm_pbc(&(top.idef),ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); + calc_principal_axes(&top,x,index,gnx,axes,moi); fprintf(axis1,"%15.10f %15.10f %15.10f %15.10f\n",t,axes[XX][XX],axes[YY][XX],axes[ZZ][XX]); @@ -142,6 +147,9 @@ int gmx_principal(int argc,char *argv[]) } while(read_next_x(oenv,status,&t,natoms,x,box)); + gmx_rmpbc_done(gpbc); + + close_trj(status); ffclose(axis1); ffclose(axis2); diff --git a/src/tools/gmx_rdf.c b/src/tools/gmx_rdf.c index 95360c4d1b..82818a857d 100644 --- a/src/tools/gmx_rdf.c +++ b/src/tools/gmx_rdf.c @@ -191,7 +191,7 @@ static void do_rdf(const char *fnNDX,const char *fnTPS,const char *fnTRX, t_blocka *excl; t_atom *atom=NULL; t_pbc pbc; - + gmx_rmpbc_t gpbc=NULL; int *is=NULL,**coi=NULL,cur,mol,i1,res,a; excl=NULL; @@ -353,12 +353,15 @@ static void do_rdf(const char *fnNDX,const char *fnTPS,const char *fnTRX, snew(x_i1,max_i); nframes = 0; invvol_sum = 0; + if (bPBC && (NULL != top)) + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); + do { /* Must init pbc every step because of pressure coupling */ copy_mat(box,box_pbc); if (bPBC) { if (top != NULL) - rm_pbc(&top->idef,ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); if (bXY) { check_box_c(box); clear_rvec(box_pbc[ZZ]); @@ -465,6 +468,9 @@ static void do_rdf(const char *fnNDX,const char *fnTPS,const char *fnTRX, } while (read_next_x(oenv,status,&t,natoms,x,box)); fprintf(stderr,"\n"); + if (bPBC && (NULL != top)) + gmx_rmpbc_done(gpbc); + close_trj(status); sfree(x); diff --git a/src/tools/gmx_relax.c b/src/tools/gmx_relax.c index 6b410b556c..8baf414b1a 100644 --- a/src/tools/gmx_relax.c +++ b/src/tools/gmx_relax.c @@ -389,6 +389,7 @@ void spectrum(bool bVerbose, rvec **corr; real **Corr; t_sij *spec; + gmx_rmpbc_t gpbc=NULL; snew(spec,npair); @@ -401,6 +402,8 @@ void spectrum(bool bVerbose, natoms = read_first_x(&status,trj,&t0,&x,box); if (natoms > nat) gmx_fatal(FARGS,"Not enough atoms in trajectory"); + gpbc = gmx_rmpbc_init(idef,ePBC,natoms,box); + do { if (nframes >= maxframes) { fprintf(stderr,"\nThere are more than the %d frames you told me!", @@ -410,7 +413,7 @@ void spectrum(bool bVerbose, t1 = t; if (bVerbose) fprintf(stderr,"\rframe: %d",nframes); - rm_pbc(idef,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); if (bFit) do_fit(natoms,w_rls,xp,x); @@ -437,6 +440,8 @@ void spectrum(bool bVerbose, if (bVerbose) fprintf(stderr,"\n"); + gmx_rmpbc_done(gpbc); + fp=ffopen("ylm.out","w"); calc_aver(fp,nframes,npair,pair,spec,maxdist); ffclose(fp); diff --git a/src/tools/gmx_rms.c b/src/tools/gmx_rms.c index a34472baa6..034983c0da 100644 --- a/src/tools/gmx_rms.c +++ b/src/tools/gmx_rms.c @@ -245,6 +245,8 @@ int gmx_rms(int argc, char *argv[]) char *gn_fit, **gn_rms; t_rgb rlo, rhi; output_env_t oenv; + gmx_rmpbc_t gpbc=NULL; + t_filenm fnm[] = { { efTPS, NULL, NULL, ffREAD }, @@ -429,8 +431,10 @@ int gmx_rms(int argc, char *argv[]) } } /* Prepare reference frame */ - if (bPBC) - rm_pbc(&(top.idef),ePBC,top.atoms.nr,box,xp,xp); + if (bPBC) { + gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box); + gmx_rmpbc(gpbc,box,xp,xp); + } if (bReset) reset_x(ifit,ind_fit,top.atoms.nr,NULL,xp,w_rls); if (bMirror) { @@ -529,7 +533,7 @@ int gmx_rms(int argc, char *argv[]) teller = 0; do { if (bPBC) - rm_pbc(&(top.idef),ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); if (bReset) reset_x(ifit,ind_fit,natoms,NULL,x,w_rls); @@ -611,7 +615,7 @@ int gmx_rms(int argc, char *argv[]) teller2 = 0; do { if (bPBC) - rm_pbc(&(top.idef),ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); if (bReset) reset_x(ifit,ind_fit,natoms,NULL,x,w_rls); @@ -649,6 +653,7 @@ int gmx_rms(int argc, char *argv[]) tel_mat2=tel_mat; freq2=freq; } + gmx_rmpbc_done(gpbc); if (bMat || bBond) { /* calculate RMS matrix */ diff --git a/src/tools/gmx_rmsf.c b/src/tools/gmx_rmsf.c index 09d36632d9..959ef4a311 100644 --- a/src/tools/gmx_rmsf.c +++ b/src/tools/gmx_rmsf.c @@ -234,6 +234,7 @@ int gmx_rmsf(int argc,char *argv[]) int d; real count=0; rvec xcm; + gmx_rmpbc_t gpbc=NULL; output_env_t oenv; @@ -300,9 +301,11 @@ int gmx_rmsf(int argc,char *argv[]) copy_mat(box,pdbbox); } - if (bFit) + if (bFit) { sub_xcm(xref,isize,index,top.atoms.atom,xcm,FALSE); - + gpbc = gmx_rmpbc_init(&top.idef,ePBC,natom,box); + } + natom = read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box); /* Now read the trj again to compute fluctuations */ @@ -310,8 +313,8 @@ int gmx_rmsf(int argc,char *argv[]) do { if (bFit) { /* Remove periodic boundary */ - rm_pbc(&(top.idef),ePBC,natom,box,x,x); - + gmx_rmpbc(gpbc,box,x,x); + /* Set center of mass to zero */ sub_xcm(x,isize,index,top.atoms.atom,xcm,FALSE); @@ -343,6 +346,10 @@ int gmx_rmsf(int argc,char *argv[]) } while(read_next_x(oenv,status,&t,natom,x,box)); close_trj(status); + if (bFit) + gmx_rmpbc_done(gpbc); + + invcount = 1.0/count; snew(Uaver,DIM*DIM); totmass = 0; diff --git a/src/tools/gmx_rotacf.c b/src/tools/gmx_rotacf.c index 3285ae41c8..c8a53d08e7 100644 --- a/src/tools/gmx_rotacf.c +++ b/src/tools/gmx_rotacf.c @@ -98,6 +98,7 @@ int gmx_rotacf(int argc,char *argv[]) int i,m,teller,n_alloc,natoms,nvec,ai,aj,ak; unsigned long mode; real t,t0,t1,dt; + gmx_rmpbc_t gpbc=NULL; t_topology *top; int ePBC; t_filenm fnm[] = { @@ -143,6 +144,8 @@ int gmx_rotacf(int argc,char *argv[]) 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,box); + /* Start the loop over frames */ t1 = t0 = t; teller = 0; @@ -155,8 +158,8 @@ int gmx_rotacf(int argc,char *argv[]) t1 = t; /* Remove periodicity */ - rm_pbc(&(top->idef),ePBC,natoms,box,x,x_s); - + gmx_rmpbc(gpbc,box,x,x_s); + /* Compute crossproducts for all vectors, if triplets. * else, just get the vectors in case of doublets. */ @@ -187,6 +190,9 @@ int gmx_rotacf(int argc,char *argv[]) close_trj(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"); diff --git a/src/tools/gmx_rotmat.c b/src/tools/gmx_rotmat.c index 67161ffc91..e9fed2f718 100644 --- a/src/tools/gmx_rotmat.c +++ b/src/tools/gmx_rotmat.c @@ -71,6 +71,8 @@ static void get_refx(output_env_t oenv,const char *trxfn,int nfitdim,int skip, real xf; matrix box,R; real *w_rls; + gmx_rmpbc_t gpbc=NULL; + nfr_all = 0; nfr = 0; @@ -89,12 +91,13 @@ static void get_refx(output_env_t oenv,const char *trxfn,int nfitdim,int skip, w_rls[a] = (bMW ? top->atoms.atom[index[a]].m : 1.0); tot_mass += w_rls[a]; } + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); do { if (nfr_all % skip == 0) { - rm_pbc(&top->idef,ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); snew(xi[nfr],gnx); for(i=0; iidef),ePBC,natoms,box); + do { teller++; - rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0); + gmx_rmpbc(gpbc,box,x0,x0); calc_angle(ePBC,box,x0,index1,index2,gnx1,gnx2,&angle, &distance,&distance1,&distance2); @@ -241,7 +245,9 @@ void sgangle_plot(const char *fn,const char *afile,const char *dfile, fprintf(sg_distance2,"%12g %12g\n",t,distance1); } while (read_next_x(oenv,status,&t,natoms,x0,box)); - + + gmx_rmpbc_done(gpbc); + fprintf(stderr,"\n"); close_trj(status); ffclose(sg_angle); @@ -363,7 +369,9 @@ void sgangle_plot_single(const char *fn,const char *afile,const char *dfile, rvec *xzero; matrix box; char buf[256]; /* for xvgr title */ + gmx_rmpbc_t gpbc=NULL; + if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0) gmx_fatal(FARGS,"Could not read coordinates from statusfile\n"); @@ -386,11 +394,12 @@ void sgangle_plot_single(const char *fn,const char *afile,const char *dfile, } snew(xzero,natoms); + gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box); do { teller++; - rm_pbc(&(top->idef),ePBC,natoms,box,x0,x0); + gmx_rmpbc(gpbc,box,x0,x0); if (teller==1) { for(i=0;imols.index; atom = top->atoms.atom; + + gpbc = gmx_rmpbc_init(&top->idef,ir->ePBC,natoms,box); /* start analysis of trajectory */ do { /* make molecules whole again */ - rm_pbc(&top->idef,ir->ePBC,natoms,box,x,x); - + gmx_rmpbc(gpbc,box,x,x); + set_pbc(&pbc,ir->ePBC,box); if (bCom) calc_com_pbc(nrefat,top,x,&pbc,index[0],xref,ir->ePBC,box); @@ -296,6 +300,8 @@ int gmx_spol(int argc,char *argv[]) nf++; } while (read_next_x(oenv,status,&t,natoms,x,box)); + + gmx_rmpbc_done(gpbc); /* clean up */ sfree(x); diff --git a/src/tools/gmx_traj.c b/src/tools/gmx_traj.c index f9d98aa678..2598d2cd24 100644 --- a/src/tools/gmx_traj.c +++ b/src/tools/gmx_traj.c @@ -530,6 +530,7 @@ int gmx_traj(int argc,char *argv[]) matrix topbox; t_trxstatus *status; t_trxstatus *status_out=NULL; + gmx_rmpbc_t gpbc=NULL; int i,j,n; int nr_xfr,nr_vfr,nr_ffr; char **grpname; @@ -721,6 +722,8 @@ int gmx_traj(int argc,char *argv[]) nr_xfr = 0; nr_vfr = 0; nr_ffr = 0; + + gpbc = gmx_rmpbc_init(&top.idef,ePBC,fr.natoms,fr.box); do { time = output_env_conv_time(oenv,fr.time); @@ -735,7 +738,7 @@ int gmx_traj(int argc,char *argv[]) } if (fr.bX && bCom) - rm_pbc(&(top.idef),ePBC,fr.natoms,fr.box,fr.x,fr.x); + gmx_rmpbc(gpbc,fr.box,fr.x,fr.x); if (bVD && fr.bV) update_histo(isize[0],index[0],fr.v,&nvhisto,&vhisto,binwidth); @@ -803,6 +806,7 @@ int gmx_traj(int argc,char *argv[]) } while(read_next_frame(oenv,status,&fr)); + gmx_rmpbc_done(gpbc); /* clean up a bit */ close_trj(status); diff --git a/src/tools/gmx_trjconv.c b/src/tools/gmx_trjconv.c index dcdb1f127e..c962361484 100644 --- a/src/tools/gmx_trjconv.c +++ b/src/tools/gmx_trjconv.c @@ -763,6 +763,7 @@ int gmx_trjconv(int argc,char *argv[]) real tshift=0,t0=-1,dt=0.001,prec; bool bFit,bFitXY,bPFit,bReset; int nfitdim; + gmx_rmpbc_t gpbc=NULL; bool bRmPBC,bPBCWhole,bPBCcomRes,bPBCcomMol,bPBCcomAtom,bPBC,bNoJump,bCluster; bool bCopy,bDoIt,bIndex,bTDump,bSetTime,bTPS=FALSE,bDTset=FALSE; bool bExec,bTimeStep=FALSE,bDumpFrame=FALSE,bSetPrec,bNeedPrec; @@ -846,6 +847,7 @@ int gmx_trjconv(int argc,char *argv[]) if (bFit || bReset) nfitdim = (fit_enum==efFitXY || fit_enum==efResetXY) ? 2 : 3; bRmPBC = bFit || bPBCWhole || bPBCcomRes || bPBCcomMol; + if (bSetUR) { if (!(bPBCcomRes || bPBCcomMol || bPBCcomAtom)) { fprintf(stderr, @@ -945,6 +947,8 @@ int gmx_trjconv(int argc,char *argv[]) if (bCONECT) gc = gmx_conect_generate(&top); + if (bRmPBC) + gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,top_box); } /* get frame number index */ @@ -1011,7 +1015,7 @@ int gmx_trjconv(int argc,char *argv[]) /* Restore reference structure and set to origin, store original location (to put structure back) */ if (bRmPBC) - rm_pbc(&(top.idef),ePBC,atoms->nr,top_box,xp,xp); + gmx_rmpbc(gpbc,top_box,xp,xp); copy_rvec(xp[index[0]],x_shift); reset_x_ndim(nfitdim,ifit,ind_fit,atoms->nr,NULL,xp,w_rls); rvec_dec(x_shift,xp[index[0]]); @@ -1202,7 +1206,7 @@ int gmx_trjconv(int argc,char *argv[]) /* Now modify the coords according to the flags, for normal fit, this is only done for output frames */ if (bRmPBC) - rm_pbc(&(top.idef),ePBC,natoms,fr.box,fr.x,fr.x); + gmx_rmpbc(gpbc,fr.box,fr.x,fr.x); reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls); do_fit(natoms,w_rls,xp,fr.x); @@ -1280,7 +1284,7 @@ int gmx_trjconv(int argc,char *argv[]) for PFit we did this already! */ if (bRmPBC) - rm_pbc(&(top.idef),ePBC,natoms,fr.box,fr.x,fr.x); + gmx_rmpbc(gpbc,fr.box,fr.x,fr.x); if (bReset) { reset_x_ndim(nfitdim,ifit,ind_fit,natoms,NULL,fr.x,w_rls); @@ -1486,7 +1490,9 @@ int gmx_trjconv(int argc,char *argv[]) fprintf(stderr,"\n"); close_trj(status); - + if (bRmPBC) + gmx_rmpbc_done(gpbc); + if (trxout) close_trx(trxout); else if (out != NULL) diff --git a/src/tools/gmx_trjorder.c b/src/tools/gmx_trjorder.c index 4e87837ff5..dbbd57d332 100644 --- a/src/tools/gmx_trjorder.c +++ b/src/tools/gmx_trjorder.c @@ -128,6 +128,7 @@ int gmx_trjorder(int argc,char *argv[]) rvec *x,*xsol,xcom,dx; matrix box; t_pbc pbc; + gmx_rmpbc_t gpbc; real t,totmass,mass,rcut2=0,n2; int natoms,nwat,ncut; char **grpname,title[256]; @@ -211,8 +212,9 @@ int gmx_trjorder(int argc,char *argv[]) } out = open_trx(opt2fn("-o",NFILE,fnm),"w"); } + gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box); do { - rm_pbc(&top.idef,ePBC,natoms,box,x,x); + gmx_rmpbc(gpbc,box,x,x); set_pbc(&pbc,ePBC,box); if (ref_a == -1) { @@ -308,7 +310,8 @@ int gmx_trjorder(int argc,char *argv[]) close_trx(out); if (fp) ffclose(fp); - + gmx_rmpbc_done(gpbc); + thanx(stderr); return 0; -- 2.11.4.GIT