#include <math.h>
/*#define HAVE_NN_LOOPS*/
-/* Set environment variable CFLAGS = "-fopenmp" when running
- * configure and define DOUSEOPENMP to make use of parallelized
- * calculation of autocorrelation function.
- * It also adds a new option -nthreads which sets the number of threads.
- * */
-/*#define DOUSEOPENMP*/
-
-#ifdef DOUSEOPENMP
-#define HAVE_OPENMP
-#include "omp.h"
-#endif
+
+#include "gmx_omp.h"
#include "statutil.h"
#include "copyrite.h"
{"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
#define MAXHH 4
-#ifdef HAVE_OPENMP
+#ifdef GMX_OPENMP
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==0)
#else
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==(threadNr))
/* Not found apparently. Add it to the list! */
/* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
-/* Unfortunately this needs to be critical it seems. */
-#ifdef HAVE_OPENMP
#pragma omp critical
-#endif
{
if (!per->p2i) {
fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
E = 0;
hb->hbE.E[d][a][h][frame] = E;
-#ifdef HAVE_OPENMP
+
#pragma omp critical
-#endif
{
hb->hbE.Etot[frame] += E;
}
}
}
+
}
static void inc_nhbonds(t_donors *ddd,int d, int h)
k = 0;
if (hb->bHBmap) {
- if (hb->hbmap[id][ia] == NULL) {
- snew(hb->hbmap[id][ia],1);
- snew(hb->hbmap[id][ia]->h,hb->maxhydro);
- snew(hb->hbmap[id][ia]->g,hb->maxhydro);
+
+#pragma omp critical
+ {
+ if (hb->hbmap[id][ia] == NULL) {
+ snew(hb->hbmap[id][ia],1);
+ snew(hb->hbmap[id][ia]->h,hb->maxhydro);
+ snew(hb->hbmap[id][ia]->g,hb->maxhydro);
+ }
+ add_ff(hb,id,k,ia,frame,ihb,p);
}
- add_ff(hb,id,k,ia,frame,ihb,p);
}
/* Strange construction with frame >=0 is a relic from old code
inc_nhbonds(&(hb->d),d,h);
}
-/* Now a redundant function. It might find use at some point though. */
-static gmx_bool in_list(atom_id selection,int isize,atom_id *index)
-{
- int i;
- gmx_bool bFound;
-
- bFound=FALSE;
- for(i=0; (i<isize) && !bFound; i++)
- if(selection == index[i])
- bFound=TRUE;
-
- return bFound;
-}
-
static char *mkatomname(t_atoms *atoms,int i)
{
static char buf[32];
int i,j,nra,n;
t_functype func_type;
t_ilist *interaction;
- atom_id nr1,nr2;
+ atom_id nr1,nr2,nr3;
gmx_bool stop;
if (!ddd->dptr) {
/* check out this functype */
if (func_type == F_SETTLE) {
- nr1=interaction->iatoms[i+1];
+ nr1 = interaction->iatoms[i+1];
+ nr2 = interaction->iatoms[i+2];
+ nr3 = interaction->iatoms[i+3];
if (ISINGRP(datable[nr1])) {
- if (ISINGRP(datable[nr1+1])) {
+ if (ISINGRP(datable[nr2])) {
datable[nr1] |= DON;
add_dh(ddd,nr1,nr1+1,grp,datable);
}
- if (ISINGRP(datable[nr1+2])) {
+ if (ISINGRP(datable[nr3])) {
datable[nr1] |= DON;
add_dh(ddd,nr1,nr1+2,grp,datable);
}
return grid;
}
-static void control_pHist(t_hbdata *hb, int nframes)
-{
- int i,j,k;
- PSTYPE p;
- for (i=0;i<hb->d.nrd;i++)
- for (j=0;j<hb->a.nra;j++)
- if (hb->per->pHist[i][j].len != 0)
- for (k=hb->hbmap[i][j][0].n0; k<nframes; k++) {
- p = getPshift(hb->per->pHist[i][j], k);
- if (p>hb->per->nper)
- fprintf(stderr, "Weird stuff in pHist[%i][%i].p at frame %i: p=%i\n",
- i,j,k,p);
- }
-}
-
static void reset_nhbonds(t_donors *ddd)
{
int i,j;
* This could be implemented slightly more efficient, but the code
* would get much more complicated.
*/
-#define B(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? 0 : (x-1))
-#define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n-1 : (x+1))
-#define GRIDMOD(j,n) (j+n)%(n)
-#define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n,bTric) \
- for(zz=B(n[ZZ],zo,bTric,FALSE); zz<=E(n[ZZ],zo,bTric,FALSE); zz++) { \
- z=GRIDMOD(zz,n[ZZ]); \
- for(yy=B(n[YY],yo,bTric,z==0||z==n[ZZ]-1); \
- yy<=E(n[YY],yo,bTric,z==0||z==n[ZZ]-1); yy++) { \
- y=GRIDMOD(yy,n[YY]); \
- for(xx=B(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); \
- xx<=E(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); xx++) { \
- x=GRIDMOD(xx,n[XX]);
-#define ENDLOOPGRIDINNER \
- } \
- } \
- } \
- \
+static inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
+{
+ return ((n==1) ? x : bTric && bEdge ? 0 : (x-1));
+}
+static inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
+{
+ return ((n==1) ? x : bTric && bEdge ? (n-1) : (x+1));
+}
+static inline int grid_mod(int j, int n)
+{
+ return (j+n) % (n);
+}
static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
{
g=NULL;
}
-static void pbc_correct(rvec dx,matrix box,rvec hbox)
-{
- int m;
- for(m=DIM-1; m>=0; m--) {
- if ( dx[m] < -hbox[m] )
- rvec_inc(dx,box[m]);
- else if ( dx[m] >= hbox[m] )
- rvec_dec(dx,box[m]);
- }
-}
-
void pbc_correct_gem(rvec dx,matrix box,rvec hbox)
{
int m;
if (bDA || (!bDA && (rha2 <= rc2))) {
rvec_sub(x[d],x[hh],r_dh);
if (bBox) {
- if (hb->bGem)
- pbc_correct_gem(r_dh,box,hbox);
- else
- pbc_correct_gem(r_dh,box,hbox);
+ pbc_correct_gem(r_dh,box,hbox);
}
if (!bDA)
"Ac(t)",
"Cc\\scontact,hb\\v{}\\z{}(t)",
"-dAc\\sfs\\v{}\\z{}/dt" };
- gmx_bool bNorm=FALSE;
+ gmx_bool bNorm=FALSE, bOMP=FALSE;
double nhb = 0;
int nhbi=0;
real *rhbex=NULL,*ht,*gt,*ght,*dght,*kt;
real *ct,*p_ct,tail,tail2,dtail,ct_fac,ght_fac,*cct;
const real tol = 1e-3;
int nframes = hb->nframes,nf;
- unsigned int **h,**g;
+ unsigned int **h=NULL,**g=NULL;
int nh,nhbonds,nhydro,ngh;
t_hbond *hbh;
PSTYPE p, *pfound = NULL, np;
t_E *E;
double *ctdouble, *timedouble, *fittedct;
double fittolerance=0.1;
+ int *dondata=NULL, thisThread;
enum {AC_NONE, AC_NN, AC_GEM, AC_LUZAR};
-
-#ifdef HAVE_OPENMP
- int *dondata=NULL, thisThread;
+#ifdef GMX_OPENMP
+ bOMP = TRUE;
+#else
+ bOMP = FALSE;
#endif
-
printf("Doing autocorrelation ");
/* Decide what kind of ACF calculations to do. */
if (bGemFit)
sprintf(legGem[(bBallistic ? 3:2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
- } else {
+ }
+ else
+ {
acType = AC_LUZAR;
printf("according to the theory of Luzar and Chandler.\n");
}
nn = nframes/2;
- if (acType != AC_NN ||
-#ifndef HAVE_OPENMP
- TRUE
-#else
- FALSE
-#endif
- ) {
+ if (acType != AC_NN || bOMP) {
snew(h,hb->maxhydro);
snew(g,hb->maxhydro);
}
ngh = 0;
anhb = 0;
- /* ------------------------------------------------
- * I got tired of waiting for the acf calculations
- * and parallelized it with openMP
- * set environment variable CFLAGS = "-fopenmp" when running
- * configure and define DOUSEOPENMP to make use of it.
- */
-
-#ifdef HAVE_OPENMP /* ================================================= \
- * Set up the OpenMP stuff, |
- * like the number of threads and such |
- */
- if (acType != AC_LUZAR)
+ if (acType != AC_LUZAR && bOMP)
{
-/* #if (_OPENMP >= 200805) /\* =====================\ *\/ */
-/* nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit()); */
-/* #else */
- nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
-/* #endif /\* _OPENMP >= 200805 ====================/ *\/ */
+ nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
- omp_set_num_threads(nThreads);
+ gmx_omp_set_num_threads(nThreads);
snew(dondata, nThreads);
for (i=0; i<nThreads; i++)
dondata[i] = -1;
fprintf(stderr, "%-7s", tmpstr);
}
}
- fprintf(stderr, "\n"); /* | */
- } /* | */
-#endif /* HAVE_OPENMP ===================================================/ */
+ fprintf(stderr, "\n");
+ }
/* Build the ACF according to acType */
#ifdef HAVE_NN_LOOPS
/* Here we're using the estimated energy for the hydrogen bonds. */
snew(ct,nn);
-#ifdef HAVE_OPENMP /* ==================================\ */
+
#pragma omp parallel \
- private(i, j, k, nh, E, rhbex, thisThread), \
+ private(i, j, k, nh, E, rhbex, thisThread) \
default(shared)
{
#pragma omp barrier
- thisThread = omp_get_thread_num();
+ thisThread = gmx_omp_get_thread_num();
rhbex = NULL;
-#endif /* ==============================================/ */
snew(rhbex, n2);
memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
-#ifdef HAVE_OPENMP /* ################################################## \
- * #
- * #
- */
#pragma omp barrier
#pragma omp for schedule (dynamic)
-#endif
for (i=0; i<hb->d.nrd; i++) /* loop over donors */
{
-#ifdef HAVE_OPENMP /* ====== Write some output ======\ */
+ if (bOMP)
+ {
#pragma omp critical
+ {
+ dondata[thisThread] = i;
+ parallel_print(dondata, nThreads);
+ }
+ }
+ else
{
- dondata[thisThread] = i;
- parallel_print(dondata, nThreads);
+ fprintf(stderr, "\r %i", i);
}
-#else
- fprintf(stderr, "\r %i", i);
-#endif /* ===========================================/ */
for (j=0; j<hb->a.nra; j++) /* loop over acceptors */
{
low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rhbex),hb->time[1]-hb->time[0],
eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
-#ifdef HAVE_OPENMP
#pragma omp critical
-#endif
{
for(k=0; (k<nn); k++)
ct[k] += rhbex[k];
} /* i loop */
sfree(rhbex);
#pragma omp barrier
-#ifdef HAVE_OPENMP
- /* # */
- } /* End of parallel block # */
- /* ##############################################################/ */
- sfree(dondata);
-#endif
+ }
+
+ if (bOMP)
+ {
+ sfree(dondata);
+ }
normalizeACF(ct, NULL, 0, nn);
snew(ctdouble, nn);
snew(timedouble, nn);
hb->time[j]-hb->time[0],
ct[j],
ctdouble[j]);
- fclose(fp);
+ xvgrclose(fp);
sfree(ct);
sfree(ctdouble);
sfree(timedouble);
case AC_GEM:
snew(ct,2*n2);
memset(ct,0,2*n2*sizeof(real));
-#ifndef HAVE_OPENMP
+#ifndef GMX_OPENMP
fprintf(stderr, "Donor:\n");
#define __ACDATA ct
#else
#define __ACDATA p_ct
#endif
-#ifdef HAVE_OPENMP /* =========================================\
- * */
-#pragma omp parallel default(none) \
+#pragma omp parallel \
private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
pfound, poff, rHbExGem, p, ihb, mMax, \
thisThread, p_ct) \
- shared(hb, dondata, ct, nn, nThreads, n2, stderr, bNorm, \
- nframes, bMerge, bContact)
+ default(shared)
{ /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
h = NULL;
g = NULL;
- thisThread = omp_get_thread_num();
+ thisThread = gmx_omp_get_thread_num();
snew(h,hb->maxhydro);
snew(g,hb->maxhydro);
mMax = INT_MIN;
/* I'm using a chunk size of 1, since I expect \
* the overhead to be really small compared \
* to the actual calculations \ */
-#pragma omp for schedule(dynamic,1) nowait /* \ */
-#endif /* HAVE_OPENMP =========================================/ */
-
+#pragma omp for schedule(dynamic,1) nowait
for (i=0; i<hb->d.nrd; i++) {
-#ifdef HAVE_OPENMP
+
+ if (bOMP)
+ {
#pragma omp critical
+ {
+ dondata[thisThread] = i;
+ parallel_print(dondata, nThreads);
+ }
+ }
+ else
{
- dondata[thisThread] = i;
- parallel_print(dondata, nThreads);
+ fprintf(stderr, "\r %i", i);
}
-#else
- fprintf(stderr, "\r %i", i);
-#endif
-
for (k=0; k<hb->a.nra; k++) {
for (nh=0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) {
hbh = hb->hbmap[i][k];
pHist = &(hb->per->pHist[i][k]);
if (ISHB(hbh->history[nh]) && pHist->len != 0) {
-/* No need for a critical section */
-/* #ifdef HAVE_OPENMP */
-/* #pragma omp critical */
-/* #endif */
{
h[nh] = hbh->h[nh];
g[nh] = hb->per->gemtype==gemAD ? hbh->g[nh] : NULL;
srenew(poff,np);
}
-/* This shouldn't have to be critical, right? */
-/* #ifdef HAVE_OPENMP */
-/* #pragma omp critical */
-/* #endif */
{
if (rHbExGem != NULL && rHbExGem[m] != NULL) {
/* This must be done, as this array was most likey
sfree(h);
sfree(g);
-#ifdef HAVE_OPENMP /* =======================================\ */
-#pragma omp critical
+
+ if (bOMP)
{
- for (i=0; i<nn; i++)
- ct[i] += p_ct[i];
+#pragma omp critical
+ {
+ for (i=0; i<nn; i++)
+ ct[i] += p_ct[i];
+ }
+ sfree(p_ct);
}
- sfree(p_ct);
} /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
- sfree(dondata);
-#endif /* HAVE_OPENMP =======================================/ */
+ if (bOMP)
+ {
+ sfree(dondata);
+ }
normalizeACF(ct, NULL, 0, nn);
fprintf(fp," %10g", fittedct[j]);
fprintf(fp,"\n");
}
- fclose(fp);
+ xvgrclose(fp);
sfree(ctdouble);
sfree(timedouble);
fp = opt2FILE("-hbn",nfile,fnm,"w");
if (opt2bSet("-g",nfile,fnm)) {
fplog = ffopen(opt2fn("-g",nfile,fnm),"w");
- if (bContact)
- fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
- else
- fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
+ fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
}
else
fplog = NULL;
ffclose(fplog);
}
-#ifdef HAVE_OPENMP
/* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
* It mimics add_frames() and init_frame() to some extent. */
static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
* even though the data its members point to will change,
* hence no need for re-syncing. */
}
-#endif
int gmx_hbond(int argc,char *argv[])
{
"which should contain exactly one atom. In this case, only hydrogen",
"bonds between atoms within the shell distance from the one atom are",
"considered.[PAR]",
+
+ "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
+ "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
+ "If contact kinetics are analyzed by using the -contact option, then",
+ "n(t) can be defined as either all pairs that are not within contact distance r at time t",
+ "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
+ "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
+ "See mentioned literature for more details and definitions."
+ "[PAR]",
/* "It is also possible to analyse specific hydrogen bonds with",
"[TT]-sel[tt]. This index file must contain a group of atom triplets",
"Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
{ "-diff", FALSE, etREAL, {&D},
"Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
-#ifdef HAVE_OPENMP
+#ifdef GMX_OPENMP
{ "-nthreads", FALSE, etINT, {&nThreads},
"Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
#endif
matrix box;
real t,ccut,dist=0.0,ang=0.0;
double max_nhb,aver_nhb,aver_dist;
- int h=0,i,j,k=0,l,start,end,id,ja,ogrp,nsel;
+ int h=0,i=0,j,k=0,l,start,end,id,ja,ogrp,nsel;
int xi,yi,zi,ai;
int xj,yj,zj,aj,xjj,yjj,zjj;
int xk,yk,zk,ak,xkk,ykk,zkk;
gmx_bool bSelected,bHBmap,bStop,bTwo,was,bBox,bTric;
- int *adist,*rdist;
+ int *adist,*rdist,*aptr,*rprt;
int grp,nabin,nrbin,bin,resdist,ihb;
char **leg;
- t_hbdata *hb;
+ t_hbdata *hb,*hbptr;
FILE *fp,*fpins=NULL,*fpnhb=NULL;
t_gridcell ***grid;
t_ncell *icell,*jcell,*kcell;
int threadNr=0;
gmx_bool bGem, bNN, bParallel;
t_gemParams *params=NULL;
+ gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
- CopyRight(stdout,argv[0]);
+ t_hbdata **p_hb=NULL; /* one per thread, then merge after the frame loop */
+ int **p_adist=NULL, **p_rdist=NULL; /* a histogram for each thread. */
+
+#ifdef GMX_OPENMP
+ bOMP = TRUE;
+#else
+ bOMP = FALSE;
+#endif
+
+ CopyRight(stderr,argv[0]);
npargs = asize(pa);
ppa = add_acf_pargs(&npargs,pa);
gmx_fatal(FARGS,"Can not analyze contact between H and A: turn off -noda");
}
}
-
-#ifndef HAVE_LIBGSL
- /* Don't pollute stdout with information about external libraries.
- *
- * printf("NO GSL! Can't find and take away ballistic term in ACF without GSL\n.");
- */
-#endif
/* Initiate main data structure! */
bHBmap = (opt2bSet("-ac",NFILE,fnm) ||
opt2bSet("-hbm",NFILE,fnm) ||
bGem);
-#ifdef HAVE_OPENMP
- /* Same thing here. There is no reason whatsoever to write the specific version of
- * OpenMP used for compilation to stdout for normal usage.
- *
- * printf("Compiled with OpenMP (%i)\n", _OPENMP);
- */
-#endif
-
- /* if (bContact && bGem) */
- /* gmx_fatal(FARGS, "Can't do reversible geminate recombination with -contact yet."); */
-
if (opt2bSet("-nhbdist",NFILE,fnm)) {
const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
fpnhb = xvgropen(opt2fn("-nhbdist",NFILE,fnm),
bParallel = FALSE;
-#ifndef HAVE_OPENMP
+#ifndef GMX_OPENMP
#define __ADIST adist
#define __RDIST rdist
#define __HBDATA hb
-#else /* HAVE_OPENMP ================================================== \
+#else /* GMX_OPENMP ================================================== \
* Set up the OpenMP stuff, |
* like the number of threads and such |
* Also start the parallel loop. |
#define __ADIST p_adist[threadNr]
#define __RDIST p_rdist[threadNr]
#define __HBDATA p_hb[threadNr]
+#endif
+ if (bOMP)
+ {
+ bParallel = !bSelected;
- bParallel = !bSelected;
+ if (bParallel)
+ {
+ actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
- if (bParallel)
- {
-/* #if (_OPENMP > 200805) */
-/* actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit()); */
-/* #else */
- actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
-/* #endif */
- omp_set_num_threads(actual_nThreads);
- printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
- fflush(stdout);
- }
- else
- {
- actual_nThreads = 1;
- }
+ gmx_omp_set_num_threads(actual_nThreads);
+ printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
+ fflush(stdout);
+ }
+ else
+ {
+ actual_nThreads = 1;
+ }
- t_hbdata **p_hb; /* one per thread, then merge after the frame loop */
- int **p_adist, **p_rdist; /* a histogram for each thread. */
- snew(p_hb, actual_nThreads);
- snew(p_adist, actual_nThreads);
- snew(p_rdist, actual_nThreads);
- for (i=0; i<actual_nThreads; i++)
- {
- snew(p_hb[i], 1);
- snew(p_adist[i], nabin+1);
- snew(p_rdist[i], nrbin+1);
-
- p_hb[i]->max_frames = 0;
- p_hb[i]->nhb = NULL;
- p_hb[i]->ndist = NULL;
- p_hb[i]->n_bound = NULL;
- p_hb[i]->time = NULL;
- p_hb[i]->nhx = NULL;
-
- p_hb[i]->bHBmap = hb->bHBmap;
- p_hb[i]->bDAnr = hb->bDAnr;
- p_hb[i]->bGem = hb->bGem;
- p_hb[i]->wordlen = hb->wordlen;
- p_hb[i]->nframes = hb->nframes;
- p_hb[i]->maxhydro = hb->maxhydro;
- p_hb[i]->danr = hb->danr;
- p_hb[i]->d = hb->d;
- p_hb[i]->a = hb->a;
- p_hb[i]->hbmap = hb->hbmap;
- p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
- p_hb[i]->per = hb->per;
+ snew(p_hb, actual_nThreads);
+ snew(p_adist, actual_nThreads);
+ snew(p_rdist, actual_nThreads);
+ for (i=0; i<actual_nThreads; i++)
+ {
+ snew(p_hb[i], 1);
+ snew(p_adist[i], nabin+1);
+ snew(p_rdist[i], nrbin+1);
+
+ p_hb[i]->max_frames = 0;
+ p_hb[i]->nhb = NULL;
+ p_hb[i]->ndist = NULL;
+ p_hb[i]->n_bound = NULL;
+ p_hb[i]->time = NULL;
+ p_hb[i]->nhx = NULL;
+
+ p_hb[i]->bHBmap = hb->bHBmap;
+ p_hb[i]->bDAnr = hb->bDAnr;
+ p_hb[i]->bGem = hb->bGem;
+ p_hb[i]->wordlen = hb->wordlen;
+ p_hb[i]->nframes = hb->nframes;
+ p_hb[i]->maxhydro = hb->maxhydro;
+ p_hb[i]->danr = hb->danr;
+ p_hb[i]->d = hb->d;
+ p_hb[i]->a = hb->a;
+ p_hb[i]->hbmap = hb->hbmap;
+ p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
+ p_hb[i]->per = hb->per;
#ifdef HAVE_NN_LOOPS
- p_hb[i]->hbE = hb->hbE;
+ p_hb[i]->hbE = hb->hbE;
#endif
- p_hb[i]->nrhb = 0;
- p_hb[i]->nrdist = 0;
+ p_hb[i]->nrhb = 0;
+ p_hb[i]->nrdist = 0;
+ }
}
/* Make a thread pool here,
* instead of forking anew at every frame. */
#pragma omp parallel \
- private(i, j, h, ii, jj, hh, E, \
+ firstprivate(i) \
+ private(j, h, ii, jj, hh, E, \
xi, yi, zi, xj, yj, zj, threadNr, \
dist, ang, peri, icell, jcell, \
grp, ogrp, ai, aj, xjj, yjj, zjj, \
xk, yk, zk, ihb, id, resdist, \
- xkk, ykk, zkk, kcell, ak, k, bTric) \
- default(none) \
- shared(hb, p_hb, p_adist, p_rdist, actual_nThreads, \
- x, bBox, box, hbox, rcut, r2cut, rshell, \
- shatom, ngrid, grid, nframes, t, \
- bParallel, bNN, index, bMerge, bContact, \
- bTwo, bDA,ccut, abin, rbin, top, \
- bSelected, bDebug, stderr, nsel, \
- bGem, oenv, fnm, fpnhb, trrStatus, natoms, \
- status, nabin, nrbin, adist, rdist, debug)
+ xkk, ykk, zkk, kcell, ak, k, bTric, \
+ bEdge_xjj, bEdge_yjj) \
+ default(shared)
{ /* Start of parallel region */
- threadNr = omp_get_thread_num();
-#endif /* HAVE_OPENMP ================================================= */
+ threadNr = gmx_omp_get_thread_num();
+
do
{
+
bTric = bBox && TRICLINIC(box);
-#ifdef HAVE_OPENMP
- sync_hbdata(hb, p_hb[threadNr], nframes, t);
+ if (bOMP)
+ {
+ sync_hbdata(hb, p_hb[threadNr], nframes, t);
+ }
#pragma omp single
-#endif
{
build_grid(hb,x,x[shatom], bBox,box,hbox, (rcut>r2cut)?rcut:r2cut,
rshell, ngrid,grid);
if (hb->bDAnr)
count_da_grid(ngrid, grid, hb->danr[nframes]);
} /* omp single */
-#ifdef HAVE_OPENMP
- p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
-#endif
+
+ if (bOMP)
+ {
+ p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
+ }
+
if (bNN)
{
#ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
/* Loop over all atom pairs and estimate interaction energy */
-#ifdef HAVE_OPENMP /* ------- */
+
#pragma omp single
-#endif /* HAVE_OPENMP ------- */
{
addFramesNN(hb, nframes);
}
-#ifdef HAVE_OPENMP /* ---------------- */
+
#pragma omp barrier
#pragma omp for schedule(dynamic)
-#endif /* HAVE_OPENMP ---------------- */
for (i=0; i<hb->d.nrd; i++)
{
for(j=0;j<hb->a.nra; j++)
{
if (bSelected)
{
-#ifdef HAVE_OPENMP
+
#pragma omp single
-#endif
{
/* Do not parallelize this just yet. */
/* int ii; */
} /* if (bSelected) */
else
{
-#ifdef HAVE_OPENMP
+
#pragma omp single
{
-#endif
- if (bGem)
- calcBoxProjection(box, hb->per->P);
+ if (bGem)
+ calcBoxProjection(box, hb->per->P);
+
+ /* loop over all gridcells (xi,yi,zi) */
+ /* Removed confusing macro, DvdS 27/12/98 */
- /* loop over all gridcells (xi,yi,zi) */
- /* Removed confusing macro, DvdS 27/12/98 */
-#ifdef HAVE_OPENMP
}
/* The outer grid loop will have to do for now. */
#pragma omp for schedule(dynamic)
-#endif
for(xi=0; xi<ngrid[XX]; xi++)
for(yi=0; (yi<ngrid[YY]); yi++)
for(zi=0; (zi<ngrid[ZZ]); zi++) {
i = icell->atoms[ai];
/* loop over all adjacent gridcells (xj,yj,zj) */
- /* This is a macro!!! */
- LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
- jcell=&(grid[zj][yj][xj].a[ogrp]);
- /* loop over acceptor atoms from other group (ogrp)
- * in this adjacent gridcell (jcell)
- */
- for (aj=0; (aj<jcell->nr); aj++) {
- j = jcell->atoms[aj];
-
- /* check if this once was a h-bond */
- peri = -1;
- ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
- hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
-
- if (ihb) {
- /* add to index if not already there */
- /* Add a hbond */
- add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
-
- /* make angle and distance distributions */
- if (ihb == hbHB && !bContact) {
- if (dist>rcut)
- gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
- ang*=RAD2DEG;
- __ADIST[(int)( ang/abin)]++;
- __RDIST[(int)(dist/rbin)]++;
- if (!bTwo) {
- int id,ia;
- if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
- gmx_fatal(FARGS,"Invalid donor %d",i);
- if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
- gmx_fatal(FARGS,"Invalid acceptor %d",j);
- resdist=abs(top.atoms.atom[i].resind-
- top.atoms.atom[j].resind);
- if (resdist >= max_hx)
- resdist = max_hx-1;
- __HBDATA->nhx[nframes][resdist]++;
+ for(zjj = grid_loop_begin(ngrid[ZZ],zi,bTric,FALSE);
+ zjj <= grid_loop_end(ngrid[ZZ],zi,bTric,FALSE);
+ zjj++)
+ {
+ zj = grid_mod(zjj,ngrid[ZZ]);
+ bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
+ for(yjj = grid_loop_begin(ngrid[YY],yi,bTric,bEdge_yjj);
+ yjj <= grid_loop_end(ngrid[YY],yi,bTric,bEdge_yjj);
+ yjj++)
+ {
+ yj = grid_mod(yjj,ngrid[YY]);
+ bEdge_xjj =
+ (yj == 0) || (yj == ngrid[YY] - 1) ||
+ (zj == 0) || (zj == ngrid[ZZ] - 1);
+ for(xjj = grid_loop_begin(ngrid[XX],xi,bTric,bEdge_xjj);
+ xjj <= grid_loop_end(ngrid[XX],xi,bTric,bEdge_xjj);
+ xjj++)
+ {
+ xj = grid_mod(xjj,ngrid[XX]);
+ jcell=&(grid[zj][yj][xj].a[ogrp]);
+ /* loop over acceptor atoms from other group (ogrp)
+ * in this adjacent gridcell (jcell)
+ */
+ for (aj=0; (aj<jcell->nr); aj++) {
+ j = jcell->atoms[aj];
+
+ /* check if this once was a h-bond */
+ peri = -1;
+ ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
+ hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
+
+ if (ihb) {
+ /* add to index if not already there */
+ /* Add a hbond */
+ add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
+
+ /* make angle and distance distributions */
+ if (ihb == hbHB && !bContact) {
+ if (dist>rcut)
+ gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
+ ang*=RAD2DEG;
+ __ADIST[(int)( ang/abin)]++;
+ __RDIST[(int)(dist/rbin)]++;
+ if (!bTwo) {
+ int id,ia;
+ if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
+ gmx_fatal(FARGS,"Invalid donor %d",i);
+ if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
+ gmx_fatal(FARGS,"Invalid acceptor %d",j);
+ resdist=abs(top.atoms.atom[i].resind-
+ top.atoms.atom[j].resind);
+ if (resdist >= max_hx)
+ resdist = max_hx-1;
+ __HBDATA->nhx[nframes][resdist]++;
+ }
+ }
+
}
- }
-
- }
- } /* for aj */
- }
- ENDLOOPGRIDINNER;
+ } /* for aj */
+ } /* for xjj */
+ } /* for yjj */
+ } /* for zjj */
} /* for ai */
} /* for grp */
} /* for xi,yi,zi */
} /* if (bSelected) {...} else */
-#ifdef HAVE_OPENMP /* ---------------------------- */
+
/* Better wait for all threads to finnish using x[] before updating it. */
- k = nframes; /* */
-#pragma omp barrier /* */
-#pragma omp critical /* */
- { /* */
+ k = nframes;
+#pragma omp barrier
+#pragma omp critical
+ {
/* Sum up histograms and counts from p_hb[] into hb */
- { /* */
+ if (bOMP)
+ {
hb->nhb[k] += p_hb[threadNr]->nhb[k];
hb->ndist[k] += p_hb[threadNr]->ndist[k];
- for (j=0; j<max_hx; j++) /**/
+ for (j=0; j<max_hx; j++)
hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
- } /* */
- } /* */
- /* */
+ }
+ }
+
/* Here are a handful of single constructs
* to share the workload a bit. The most
* important one is of course the last one,
* where there's a potential bottleneck in form
* of slow I/O. */
-#pragma omp single /* ++++++++++++++++, */
-#endif /* HAVE_OPENMP ----------------+------------*/
- { /* + */
- if (hb != NULL) /* */
- { /* + */
+#pragma omp barrier
+#pragma omp single
+ {
+ if (hb != NULL)
+ {
analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,k,t,oenv);
- } /* + */
- } /* + */
-#ifdef HAVE_OPENMP /* + */
-#pragma omp single /* +++ +++ */
-#endif /* + */
- { /* + */
- if (fpnhb) /* + */
+ }
+ }
+
+#pragma omp single
+ {
+ if (fpnhb)
do_nhb_dist(fpnhb,hb,t);
- } /* + */
+ }
} /* if (bNN) {...} else + */
-#ifdef HAVE_OPENMP /* + */
-#pragma omp single /* +++ +++ */
-#endif /* + */
- { /* + */
+
+#pragma omp single
+ {
trrStatus = (read_next_x(oenv,status,&t,natoms,x,box));
- nframes++; /* + */
- } /* + */
-#ifdef HAVE_OPENMP /* +++++++++++++++++ */
+ nframes++;
+ }
+
#pragma omp barrier
-#endif
} while (trrStatus);
-#ifdef HAVE_OPENMP
-#pragma omp critical
+ if (bOMP)
{
- hb->nrhb += p_hb[threadNr]->nrhb;
- hb->nrdist += p_hb[threadNr]->nrdist;
- }
- /* Free parallel datastructures */
- sfree(p_hb[threadNr]->nhb);
- sfree(p_hb[threadNr]->ndist);
- sfree(p_hb[threadNr]->nhx);
+#pragma omp critical
+ {
+ hb->nrhb += p_hb[threadNr]->nrhb;
+ hb->nrdist += p_hb[threadNr]->nrdist;
+ }
+ /* Free parallel datastructures */
+ sfree(p_hb[threadNr]->nhb);
+ sfree(p_hb[threadNr]->ndist);
+ sfree(p_hb[threadNr]->nhx);
#pragma omp for
- for (i=0; i<nabin; i++)
- for (j=0; j<actual_nThreads; j++)
+ for (i=0; i<nabin; i++)
+
for (j=0; j<actual_nThreads; j++)
- adist[i] += p_adist[j][i];
+ adist[i] += p_adist[j][i];
#pragma omp for
- for (i=0; i<=nrbin; i++)
- for (j=0; j<actual_nThreads; j++)
- rdist[i] += p_rdist[j][i];
+ for (i=0; i<=nrbin; i++)
+
for (j=0; j<actual_nThreads; j++)
+ rdist[i] += p_rdist[j][i];
- sfree(p_adist[threadNr]);
- sfree(p_rdist[threadNr]);
+ sfree(p_adist[threadNr]);
+ sfree(p_rdist[threadNr]);
+ }
} /* End of parallel region */
- sfree(p_adist);
- sfree(p_rdist);
-#endif
+ if (bOMP)
+ {
+ sfree(p_adist);
+ sfree(p_rdist);
+ }
if(nframes <2 && (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm)))
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff