From: Mark Abraham Date: Tue, 19 Feb 2013 13:01:57 +0000 (+0100) Subject: Merge release-4-5-patches into release-4-6 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=ceb646365debf0b2c90e4f5a25e1b01b48a314d5;p=alexxy%2Fgromacs.git Merge release-4-5-patches into release-4-6 Avoided the 4.5.6 release patch by using git merge -s ours 10b109beb3ffa95e37386d38140d2c0fd7769f20 git merge release-4-5-patches tree=$(git log -1 HEAD --pretty=%T) git reset --hard $(git cat-file commit HEAD | sed '1,/^$/d' | \ git commit-tree $tree -p 2b05689 -p 8d6cc14) Conflict resolution was straightforward, except that the vsite code was somewhat reorganized in the meantime. Conflicts: Makefile.am src/contrib/Makefile.am src/contrib/addquote.c src/gmxlib/copyrite.c src/mdlib/vsite.c Change-Id: Ib90ac4e55fac308c433b15c344caf5ee41f0146c --- ceb646365debf0b2c90e4f5a25e1b01b48a314d5 diff --cc src/gmxlib/copyrite.c index fe3a2df9fe,f4d65ab56b..ae32be48e8 --- a/src/gmxlib/copyrite.c +++ b/src/gmxlib/copyrite.c @@@ -152,72 -119,53 +152,60 @@@ gmx_bool be_cool(void void space(FILE *out, int n) { - fprintf(out,"%*s",n,""); + fprintf(out, "%*s", n, ""); } - void f(char *a) - { - int i; - int len = strlen(a); - - for (i = 0; i < len; i++) - { - a[i] = ~a[i]; - } - } - -static void sp_print(FILE *out,const char *s) +static void sp_print(FILE *out, const char *s) { - int slen; - - slen=strlen(s); - space(out,(80-slen)/2); - fprintf(out,"%s\n",s); + int slen; + + slen = strlen(s); + space(out, (80-slen)/2); + fprintf(out, "%s\n", s); } -static void ster_print(FILE *out,const char *s) +static void ster_print(FILE *out, const char *s) { - int slen; - char buf[128]; - - snprintf(buf,128,":-) %s (-:",s); - slen=strlen(buf); - space(out,(80-slen)/2); - fprintf(out,"%s\n",buf); + int slen; + char buf[128]; + + snprintf(buf, 128, ":-) %s (-:", s); + slen = strlen(buf); + space(out, (80-slen)/2); + fprintf(out, "%s\n", buf); } -static void pukeit(const char *db,const char *defstring, char *retstring, - int retsize, int *cqnum) +static void pukeit(const char *db, const char *defstring, char *retstring, + int retsize, int *cqnum) { - FILE *fp; - char **help; - int i,nhlp; - int seed; - - if (be_cool() && ((fp = low_libopen(db,FALSE)) != NULL)) { - nhlp=fget_lines(fp,&help); - /* for libraries we can use the low-level close routines */ - ffclose(fp); - seed=time(NULL); - *cqnum=nhlp*rando(&seed); - if (strlen(help[*cqnum]) >= STRLEN) - help[*cqnum][STRLEN-1] = '\0'; - strncpy(retstring,help[*cqnum],retsize); - for(i=0; (i= STRLEN) + { + help[*cqnum][STRLEN-1] = '\0'; + } + strncpy(retstring, help[*cqnum], retsize); - f(retstring); + for (i = 0; (i < nhlp); i++) + { + sfree(help[i]); + } + sfree(help); + } + else + { + strncpy(retstring, defstring, retsize); + } } void bromacs(char *retstring, int retsize) diff --cc src/mdlib/vsite.c index ab250de74f,2eadddecb1..a2d541ea25 --- a/src/mdlib/vsite.c +++ b/src/mdlib/vsite.c @@@ -422,224 -406,400 +422,230 @@@ static int constr_vsiten(t_iatom *ia, t } -void construct_vsites(FILE *log,gmx_vsite_t *vsite, - rvec x[],t_nrnb *nrnb, - real dt,rvec *v, - t_iparams ip[],t_ilist ilist[], - int ePBC,gmx_bool bMolPBC,t_graph *graph, - t_commrec *cr,matrix box) +void construct_vsites_thread(gmx_vsite_t *vsite, + rvec x[], t_nrnb *nrnb, + real dt, rvec *v, + t_iparams ip[], t_ilist ilist[], + t_pbc *pbc_null) { - rvec xpbc,xv,vv,dx; - real a1,b1,c1,inv_dt; - int i,inc,ii,nra,nr,tp,ftype; - t_iatom avsite,ai,aj,ak,al,pbc_atom; - t_iatom *ia; - t_pbc pbc,*pbc_null,*pbc_null2; - gmx_bool bDomDec; - int *vsite_pbc,ishift; - rvec reftmp,vtmp,rtmp; - - bDomDec = cr && DOMAINDECOMP(cr); - - /* We only need to do pbc when we have inter-cg vsites */ - if (ePBC != epbcNONE && (bDomDec || bMolPBC) && vsite->n_intercg_vsite) { - /* This is wasting some CPU time as we now do this multiple times - * per MD step. But how often do we have vsites with full pbc? - */ - pbc_null = set_pbc_dd(&pbc,ePBC,cr!=NULL ? cr->dd : NULL,FALSE,box); - } else { - pbc_null = NULL; - } - - if (cr) { - if (bDomDec) { - dd_move_x_vsites(cr->dd,box,x); - } else if (vsite->bPDvsitecomm) { - /* I'm not sure whether the periodicity and shift are guaranteed - * to be consistent between different nodes when running e.g. polymers - * in parallel. In this special case we thus unshift/shift, - * but only when necessary. This is to make sure the coordinates - * we move don't end up a box away... - */ - if (graph) - unshift_self(graph,box,x); - - move_construct_x(vsite->vsitecomm,x,cr); - - if (graph) - shift_self(graph,box,x); - } - } - - if (v) { - inv_dt = 1.0/dt; - } else { - inv_dt = 1.0; - } - - pbc_null2 = NULL; - for(ftype=0; (ftypevsite_pbc_loc[ftype-F_VSITE2]; - } else { - vsite_pbc = NULL; - } - - for(i=0; (ifunctype[tp]) - gmx_incons("Function types for vsites wrong"); - */ - - /* The vsite and constructing atoms */ - avsite = ia[1]; - ai = ia[2]; - - /* Constants for constructing vsites */ - a1 = ip[tp].vsite.a; - /* Check what kind of pbc we need to use */ - if (vsite_pbc) { - pbc_atom = vsite_pbc[i/(1+nra)]; - if (pbc_atom > -2) { - if (pbc_atom >= 0) { - /* We need to copy the coordinates here, - * single for single atom cg's pbc_atom is the vsite itself. - */ - copy_rvec(x[pbc_atom],xpbc); - } - pbc_null2 = pbc_null; - } else { - pbc_null2 = NULL; - } - } else { - pbc_atom = -2; - } - /* Copy the old position */ - copy_rvec(x[avsite],xv); - - /* Construct the vsite depending on type */ - inc = nra+1; - switch (ftype) { - case F_VSITE2: - aj = ia[3]; - constr_vsite2(x[ai],x[aj],x[avsite],a1,pbc_null2); - break; - case F_VSITE3: - aj = ia[3]; - ak = ia[4]; - b1 = ip[tp].vsite.b; - constr_vsite3(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2); - break; - case F_VSITE3FD: - aj = ia[3]; - ak = ia[4]; - b1 = ip[tp].vsite.b; - constr_vsite3FD(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2); - break; - case F_VSITE3FAD: - aj = ia[3]; - ak = ia[4]; - b1 = ip[tp].vsite.b; - constr_vsite3FAD(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2); - break; - case F_VSITE3OUT: - aj = ia[3]; - ak = ia[4]; - b1 = ip[tp].vsite.b; - c1 = ip[tp].vsite.c; - constr_vsite3OUT(x[ai],x[aj],x[ak],x[avsite],a1,b1,c1,pbc_null2); - break; - case F_VSITE4FD: - aj = ia[3]; - ak = ia[4]; - al = ia[5]; - b1 = ip[tp].vsite.b; - c1 = ip[tp].vsite.c; - constr_vsite4FD(x[ai],x[aj],x[ak],x[al],x[avsite],a1,b1,c1, - pbc_null2); - break; - case F_VSITE4FDN: - aj = ia[3]; - ak = ia[4]; - al = ia[5]; - b1 = ip[tp].vsite.b; - c1 = ip[tp].vsite.c; - constr_vsite4FDN(x[ai],x[aj],x[ak],x[al],x[avsite],a1,b1,c1, - pbc_null2); - break; - case F_VSITEN: - inc = constr_vsiten(ia,ip,x,pbc_null2); - break; - default: - gmx_fatal(FARGS,"No such vsite type %d in %s, line %d", - ftype,__FILE__,__LINE__); - } - - if (pbc_atom >= 0) { - /* Match the pbc of this vsite to the rest of its charge group */ - ishift = pbc_dx_aiuc(pbc_null,x[avsite],xpbc,dx); - if (ishift != CENTRAL) - rvec_add(xpbc,dx,x[avsite]); - } - if (v) { - /* Calculate velocity of vsite... */ - rvec_sub(x[avsite],xv,vv); - svmul(inv_dt,vv,v[avsite]); - } - - /* Increment loop variables */ - i += inc; - ia += inc; - } - } - } -} + gmx_bool bPBCAll; + rvec xpbc, xv, vv, dx; + real a1, b1, c1, inv_dt; + int i, inc, ii, nra, nr, tp, ftype; + t_iatom avsite, ai, aj, ak, al, pbc_atom; + t_iatom *ia; + t_pbc *pbc_null2; + int *vsite_pbc, ishift; + rvec reftmp, vtmp, rtmp; + + if (v != NULL) + { + inv_dt = 1.0/dt; + } + else + { + inv_dt = 1.0; + } -static void spread_vsite2(t_iatom ia[],real a, - rvec x[],rvec f[],rvec fshift[], - t_pbc *pbc,t_graph *g) -{ - rvec fi,fj,dx; - t_iatom av,ai,aj; - ivec di; - real b; - int siv,sij; - - av = ia[1]; - ai = ia[2]; - aj = ia[3]; - - svmul(1-a,f[av],fi); - svmul( a,f[av],fj); - /* 7 flop */ - - rvec_inc(f[ai],fi); - rvec_inc(f[aj],fj); - /* 6 Flops */ - - if (g) { - ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,av),di); - siv = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),di); - sij = IVEC2IS(di); - } else if (pbc) { - siv = pbc_dx_aiuc(pbc,x[ai],x[av],dx); - sij = pbc_dx_aiuc(pbc,x[ai],x[aj],dx); - } else { - siv = CENTRAL; - sij = CENTRAL; - } - - if (fshift && (siv != CENTRAL || sij != CENTRAL)) { - rvec_inc(fshift[siv],f[av]); - rvec_dec(fshift[CENTRAL],fi); - rvec_dec(fshift[sij],fj); - } - - /* TOTAL: 13 flops */ -} + bPBCAll = (pbc_null != NULL && !vsite->bHaveChargeGroups); -void construct_vsites_mtop(FILE *log,gmx_vsite_t *vsite, - gmx_mtop_t *mtop,rvec x[]) -{ - int as,mb,mol; - gmx_molblock_t *molb; - gmx_moltype_t *molt; - - as = 0; - for(mb=0; mbnmolblock; mb++) { - molb = &mtop->molblock[mb]; - molt = &mtop->moltype[molb->type]; - for(mol=0; molnmol; mol++) { - construct_vsites(log,vsite,x+as,NULL,0.0,NULL, - mtop->ffparams.iparams,molt->ilist, - epbcNONE,TRUE,NULL,NULL,NULL); - as += molt->atoms.nr; - } - } -} + pbc_null2 = NULL; + vsite_pbc = NULL; + for (ftype = 0; (ftype < F_NRE); ftype++) + { + if ((interaction_function[ftype].flags & IF_VSITE) && + ilist[ftype].nr > 0) + { + nra = interaction_function[ftype].nratoms; + inc = 1 + nra; + nr = ilist[ftype].nr; + ia = ilist[ftype].iatoms; -static void spread_vsite3(t_iatom ia[],real a,real b, - rvec x[],rvec f[],rvec fshift[], - t_pbc *pbc,t_graph *g) -{ - rvec fi,fj,fk,dx; - atom_id av,ai,aj,ak; - ivec di; - int siv,sij,sik; - - av = ia[1]; - ai = ia[2]; - aj = ia[3]; - ak = ia[4]; - - svmul(1-a-b,f[av],fi); - svmul( a,f[av],fj); - svmul( b,f[av],fk); - /* 11 flops */ - - rvec_inc(f[ai],fi); - rvec_inc(f[aj],fj); - rvec_inc(f[ak],fk); - /* 9 Flops */ - - if (g) { - ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ia[1]),di); - siv = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),di); - sij = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ak),di); - sik = IVEC2IS(di); - } else if (pbc) { - siv = pbc_dx_aiuc(pbc,x[ai],x[av],dx); - sij = pbc_dx_aiuc(pbc,x[ai],x[aj],dx); - sik = pbc_dx_aiuc(pbc,x[ai],x[ak],dx); - } else { - siv = CENTRAL; - sij = CENTRAL; - sik = CENTRAL; - } - - if (fshift && (siv!=CENTRAL || sij!=CENTRAL || sik!=CENTRAL)) { - rvec_inc(fshift[siv],f[av]); - rvec_dec(fshift[CENTRAL],fi); - rvec_dec(fshift[sij],fj); - rvec_dec(fshift[sik],fk); - } - - /* TOTAL: 20 flops */ + if (bPBCAll) + { + pbc_null2 = pbc_null; + } + else if (pbc_null != NULL) + { + vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2]; + } + + for (i = 0; i < nr; ) + { + tp = ia[0]; + + /* The vsite and constructing atoms */ + avsite = ia[1]; + ai = ia[2]; - aj = ia[3]; + + /* Constants for constructing vsites */ + a1 = ip[tp].vsite.a; + /* Check what kind of pbc we need to use */ + if (bPBCAll) + { + /* No charge groups, vsite follows its own pbc */ + pbc_atom = avsite; + copy_rvec(x[avsite], xpbc); + } + else if (vsite_pbc != NULL) + { + pbc_atom = vsite_pbc[i/(1+nra)]; + if (pbc_atom > -2) + { + if (pbc_atom >= 0) + { + /* We need to copy the coordinates here, + * single for single atom cg's pbc_atom + * is the vsite itself. + */ + copy_rvec(x[pbc_atom], xpbc); + } + pbc_null2 = pbc_null; + } + else + { + pbc_null2 = NULL; + } + } + else + { + pbc_atom = -2; + } + /* Copy the old position */ + copy_rvec(x[avsite], xv); + + /* Construct the vsite depending on type */ + switch (ftype) + { + case F_VSITE2: ++ aj = ia[3]; + constr_vsite2(x[ai], x[aj], x[avsite], a1, pbc_null2); + break; + case F_VSITE3: ++ aj = ia[3]; + ak = ia[4]; + b1 = ip[tp].vsite.b; + constr_vsite3(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2); + break; + case F_VSITE3FD: ++ aj = ia[3]; + ak = ia[4]; + b1 = ip[tp].vsite.b; + constr_vsite3FD(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2); + break; + case F_VSITE3FAD: ++ aj = ia[3]; + ak = ia[4]; + b1 = ip[tp].vsite.b; + constr_vsite3FAD(x[ai], x[aj], x[ak], x[avsite], a1, b1, pbc_null2); + break; + case F_VSITE3OUT: ++ aj = ia[3]; + ak = ia[4]; + b1 = ip[tp].vsite.b; + c1 = ip[tp].vsite.c; + constr_vsite3OUT(x[ai], x[aj], x[ak], x[avsite], a1, b1, c1, pbc_null2); + break; + case F_VSITE4FD: ++ aj = ia[3]; + ak = ia[4]; + al = ia[5]; + b1 = ip[tp].vsite.b; + c1 = ip[tp].vsite.c; + constr_vsite4FD(x[ai], x[aj], x[ak], x[al], x[avsite], a1, b1, c1, + pbc_null2); + break; + case F_VSITE4FDN: ++ aj = ia[3]; + ak = ia[4]; + al = ia[5]; + b1 = ip[tp].vsite.b; + c1 = ip[tp].vsite.c; + constr_vsite4FDN(x[ai], x[aj], x[ak], x[al], x[avsite], a1, b1, c1, + pbc_null2); + break; + case F_VSITEN: + inc = constr_vsiten(ia, ip, x, pbc_null2); + break; + default: + gmx_fatal(FARGS, "No such vsite type %d in %s, line %d", + ftype, __FILE__, __LINE__); + } + + if (pbc_atom >= 0) + { + /* Match the pbc of this vsite to the rest of its charge group */ + ishift = pbc_dx_aiuc(pbc_null, x[avsite], xpbc, dx); + if (ishift != CENTRAL) + { + rvec_add(xpbc, dx, x[avsite]); + } + } + if (v != NULL) + { + /* Calculate velocity of vsite... */ + rvec_sub(x[avsite], xv, vv); + svmul(inv_dt, vv, v[avsite]); + } + + /* Increment loop variables */ + i += inc; + ia += inc; + } + } + } } -static void spread_vsite3FD(t_iatom ia[],real a,real b, - rvec x[],rvec f[],rvec fshift[], - gmx_bool VirCorr,matrix dxdf, - t_pbc *pbc,t_graph *g) +void construct_vsites(FILE *log, gmx_vsite_t *vsite, + rvec x[], t_nrnb *nrnb, + real dt, rvec *v, + t_iparams ip[], t_ilist ilist[], + int ePBC, gmx_bool bMolPBC, t_graph *graph, + t_commrec *cr, matrix box) { - real fx,fy,fz,c,invl,fproj,a1; - rvec xvi,xij,xjk,xix,fv,temp; - t_iatom av,ai,aj,ak; - int svi,sji,skj,d; - ivec di; - - av = ia[1]; - ai = ia[2]; - aj = ia[3]; - ak = ia[4]; - copy_rvec(f[av],fv); - - sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij); - skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk); - /* 6 flops */ - - /* xix goes from i to point x on the line jk */ - xix[XX]=xij[XX]+a*xjk[XX]; - xix[YY]=xij[YY]+a*xjk[YY]; - xix[ZZ]=xij[ZZ]+a*xjk[ZZ]; - /* 6 flops */ - - invl=gmx_invsqrt(iprod(xix,xix)); - c=b*invl; - /* 4 + ?10? flops */ - - fproj=iprod(xix,fv)*invl*invl; /* = (xix . f)/(xix . xix) */ - - temp[XX]=c*(fv[XX]-fproj*xix[XX]); - temp[YY]=c*(fv[YY]-fproj*xix[YY]); - temp[ZZ]=c*(fv[ZZ]-fproj*xix[ZZ]); - /* 16 */ - - /* c is already calculated in constr_vsite3FD - storing c somewhere will save 26 flops! */ - - a1=1-a; - f[ai][XX] += fv[XX] - temp[XX]; - f[ai][YY] += fv[YY] - temp[YY]; - f[ai][ZZ] += fv[ZZ] - temp[ZZ]; - f[aj][XX] += a1*temp[XX]; - f[aj][YY] += a1*temp[YY]; - f[aj][ZZ] += a1*temp[ZZ]; - f[ak][XX] += a*temp[XX]; - f[ak][YY] += a*temp[YY]; - f[ak][ZZ] += a*temp[ZZ]; - /* 19 Flops */ - - if (g) { - ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di); - sji = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di); - skj = IVEC2IS(di); - } else if (pbc) { - svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi); - } else { - svi = CENTRAL; - } - - if (fshift && (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL)) { - rvec_dec(fshift[svi],fv); - fshift[CENTRAL][XX] += fv[XX] - (1 + a)*temp[XX]; - fshift[CENTRAL][YY] += fv[YY] - (1 + a)*temp[YY]; - fshift[CENTRAL][ZZ] += fv[ZZ] - (1 + a)*temp[ZZ]; - fshift[ sji][XX] += temp[XX]; - fshift[ sji][YY] += temp[YY]; - fshift[ sji][ZZ] += temp[ZZ]; - fshift[ skj][XX] += a*temp[XX]; - fshift[ skj][YY] += a*temp[YY]; - fshift[ skj][ZZ] += a*temp[ZZ]; - } + t_pbc pbc, *pbc_null; + gmx_bool bDomDec; + int nthreads; - if (VirCorr) + bDomDec = cr && DOMAINDECOMP(cr); + + /* We only need to do pbc when we have inter-cg vsites */ + if (ePBC != epbcNONE && (bDomDec || bMolPBC) && vsite->n_intercg_vsite) { - /* When VirCorr=TRUE, the virial for the current forces is not - * calculated from the redistributed forces. This means that - * the effect of non-linear virtual site constructions on the virial - * needs to be added separately. This contribution can be calculated - * in many ways, but the simplest and cheapest way is to use - * the first constructing atom ai as a reference position in space: - * subtract (xv-xi)*fv and add (xj-xi)*fj + (xk-xi)*fk. + /* This is wasting some CPU time as we now do this multiple times + * per MD step. But how often do we have vsites with full pbc? */ - rvec xiv; - int i,j; - - pbc_rvec_sub(pbc,x[av],x[ai],xiv); + pbc_null = set_pbc_dd(&pbc, ePBC, cr != NULL ? cr->dd : NULL, FALSE, box); + } + else + { + pbc_null = NULL; + } - for(i=0; idd, box, x); + } + else if (vsite->bPDvsitecomm) { - for(j=0; jvsitecomm, x, cr); + + if (graph != NULL) + { + shift_self(graph, box, x); } } }