KERNEL void hybridForce(int numParticles, int paddedNumParticles, GLOBAL mm_long* RESTRICT force, GLOBAL mm_long* RESTRICT force0, GLOBAL mm_long* RESTRICT force1, GLOBAL mm_long* RESTRICT dforce0, GLOBAL mm_long* RESTRICT dforce1, GLOBAL int* RESTRICT invAtomOrder, GLOBAL int* RESTRICT inner0InvAtomOrder, GLOBAL int* RESTRICT inner1InvAtomOrder, real dEdu0, real dEdu1) { for (int i = GLOBAL_ID; i < numParticles; i += GLOBAL_SIZE) { int index = invAtomOrder[i]; int index0 = inner0InvAtomOrder[i]; int index1 = inner1InvAtomOrder[i]; mm_long fx0 = force0[index0]+dforce0[index]; mm_long fy0 = force0[index0+paddedNumParticles]+dforce0[index+paddedNumParticles]; mm_long fz0 = force0[index0+paddedNumParticles*2]+dforce0[index+paddedNumParticles*2]; mm_long fx1 = force1[index1]+dforce1[index]; mm_long fy1 = force1[index1+paddedNumParticles]+dforce1[index+paddedNumParticles]; mm_long fz1 = force1[index1+paddedNumParticles*2]+dforce1[index+paddedNumParticles*2]; force[index] += (mm_long) (dEdu0*fx0 + dEdu1*fx1); force[index+paddedNumParticles] += (mm_long) (dEdu0*fy0 + dEdu1*fy1); force[index+paddedNumParticles*2] += (mm_long) (dEdu0*fz0 + dEdu1*fz1); } } //reset variable displacement forces KERNEL void resetDisplForce(int numParticles, int paddedNumParticles, GLOBAL mm_ulong* RESTRICT dforce0, GLOBAL mm_ulong* RESTRICT dforce1) { mm_ulong zero = 0; for (int index = GLOBAL_ID; index < numParticles; index += GLOBAL_SIZE) { dforce0[index] = zero; dforce0[index+paddedNumParticles] = zero; dforce0[index+paddedNumParticles*2] = zero; dforce1[index] = zero; dforce1[index+paddedNumParticles] = zero; dforce1[index+paddedNumParticles*2] = zero; } } //add forces due to variable displacements KERNEL void displForce(int numParticles, int paddedNumParticles, GLOBAL mm_long* RESTRICT force0, GLOBAL mm_long* RESTRICT force1, GLOBAL mm_long* RESTRICT dforce0, GLOBAL mm_long* RESTRICT dforce1, GLOBAL int4* displParticles, GLOBAL int* RESTRICT atomOrder, GLOBAL int* RESTRICT invAtomOrder, GLOBAL int* RESTRICT inner0InvAtomOrder, GLOBAL int* RESTRICT inner1InvAtomOrder) { GLOBAL mm_ulong* df0 = (GLOBAL mm_ulong*) dforce0; GLOBAL mm_ulong* df1 = (GLOBAL mm_ulong*) dforce1; for (int index = GLOBAL_ID; index < numParticles; index += GLOBAL_SIZE) { int atom = atomOrder[index]; int pj1 = displParticles[atom].x; int pi1 = displParticles[atom].y; int pj0 = displParticles[atom].z; int pi0 = displParticles[atom].w; int index0 = inner0InvAtomOrder[atom]; int index1 = inner1InvAtomOrder[atom]; if (pj1 >= 0 && pi1 >= 0) { int j1 = invAtomOrder[pj1]; int i1 = invAtomOrder[pi1]; ATOMIC_ADD(&df1[j1], (mm_ulong) force1[index1]); ATOMIC_ADD(&df1[j1+paddedNumParticles], (mm_ulong) force1[index1+paddedNumParticles]); ATOMIC_ADD(&df1[j1+paddedNumParticles*2], (mm_ulong) force1[index1+paddedNumParticles*2]); ATOMIC_ADD(&df1[i1], (mm_ulong) -force1[index1]); ATOMIC_ADD(&df1[i1+paddedNumParticles], (mm_ulong) -force1[index1+paddedNumParticles]); ATOMIC_ADD(&df1[i1+paddedNumParticles*2],(mm_ulong) -force1[index1+paddedNumParticles*2]); } if (pj0 >= 0 && pi0 >= 0) { int j0 = invAtomOrder[pj0]; int i0 = invAtomOrder[pi0]; ATOMIC_ADD(&df0[j0], (mm_ulong) force0[index0]); ATOMIC_ADD(&df0[j0+paddedNumParticles], (mm_ulong) force0[index0+paddedNumParticles]); ATOMIC_ADD(&df0[j0+paddedNumParticles*2],(mm_ulong) force0[index0+paddedNumParticles*2]); ATOMIC_ADD(&df0[i0], (mm_ulong) -force0[index0]); ATOMIC_ADD(&df0[i0+paddedNumParticles], (mm_ulong) -force0[index0+paddedNumParticles]); ATOMIC_ADD(&df0[i0+paddedNumParticles*2], (mm_ulong) -force0[index0+paddedNumParticles*2]); } } } KERNEL void copyState(int numParticles, GLOBAL real4* RESTRICT posq, GLOBAL real4* RESTRICT posq0, GLOBAL real4* RESTRICT posq1, GLOBAL real4* RESTRICT displacement0, GLOBAL real4* RESTRICT displacement1, GLOBAL int4* displParticles, GLOBAL int* RESTRICT atomOrder, GLOBAL int* RESTRICT invAtomOrder, GLOBAL int* RESTRICT inner0InvAtomOrder, GLOBAL int* RESTRICT inner1InvAtomOrder #ifdef USE_MIXED_PRECISION , GLOBAL real4* RESTRICT posqCorrection, GLOBAL real4* RESTRICT posq0Correction, GLOBAL real4* RESTRICT posq1Correction #endif ) { for (int i = GLOBAL_ID; i < numParticles; i += GLOBAL_SIZE) { int atom = atomOrder[i]; //default fixed lab frame displacement real4 displ0 = displacement0[atom]; real4 displ1 = displacement1[atom]; //override with variable displacements if set int pj1 = displParticles[atom].x; int pi1 = displParticles[atom].y; int pj0 = displParticles[atom].z; int pi0 = displParticles[atom].w; if (pj1 >= 0 && pi1 >= 0) { // variable system coordinate displacements int indexj1 = invAtomOrder[pj1]; int indexi1 = invAtomOrder[pi1]; displ1 = make_real4((real) posq[indexj1].x- posq[indexi1].x, (real) posq[indexj1].y- posq[indexi1].y, (real) posq[indexj1].z- posq[indexi1].z, (real) 0); if (pj0 >= 0 && pi0 >= 0) { int indexj0 = invAtomOrder[pj0]; int indexi0 = invAtomOrder[pi0]; displ0 = make_real4((real) posq[indexj0].x - posq[indexi0].x, (real) posq[indexj0].y - posq[indexi0].y, (real) posq[indexj0].z - posq[indexi0].z, (real) 0); } else { displ0 = make_real4((real) 0, (real) 0, (real) 0, (real) 0); } } int index0 = inner0InvAtomOrder[atom]; int index1 = inner1InvAtomOrder[atom]; real4 p0 = posq[i] + make_real4((real) displ0.x, (real) displ0.y, (real) displ0.z, 0); real4 p1 = posq[i] + make_real4((real) displ1.x, (real) displ1.y, (real) displ1.z, 0); p0.w = posq0[i].w; p1.w = posq1[i].w; posq0[index0] = p0; posq1[index1] = p1; #ifdef USE_MIXED_PRECISION posq0Correction[index0] = posqCorrection[i]; posq1Correction[index1] = posqCorrection[i]; #endif } }