}
-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; (ftype<F_NRE); ftype++) {
- if (interaction_function[ftype].flags & IF_VSITE) {
- nra = interaction_function[ftype].nratoms;
- nr = ilist[ftype].nr;
- ia = ilist[ftype].iatoms;
-
- if (pbc_null) {
- vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
- } else {
- vsite_pbc = NULL;
- }
-
- for(i=0; (i<nr); ) {
- tp = ia[0];
- /*
- if (ftype != idef->functype[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; mb<mtop->nmolblock; mb++) {
- molb = &mtop->molblock[mb];
- molt = &mtop->moltype[molb->type];
- for(mol=0; mol<molb->nmol; 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; i<DIM; i++)
+ if (cr)
+ {
+ if (bDomDec)
+ {
+ dd_move_x_vsites(cr->dd, box, x);
+ }
+ else if (vsite->bPDvsitecomm)
{
- for(j=0; j<DIM; j++)
+ /* 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 != NULL)
{
- /* As xix is a linear combination of j and k, use that here */
- dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
+ unshift_self(graph, box, x);
+ }
+
+ move_construct_x(vsite->vsitecomm, x, cr);
+
+ if (graph != NULL)
+ {
+ shift_self(graph, box, x);
}
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff