NAMD: ComputeBondedCUDAKernel.cu Source File

00001 #ifdef WIN32
00002 #define _USE_MATH_DEFINES
00003 #define __thread __declspec(thread)
00004 #endif
00005 #include <math.h>
00006 #include <cuda.h>
00007 #if __CUDACC_VER_MAJOR__ >= 11
00008 #include <cub/cub.cuh>
00009 #else
00010 #include <namd_cub/cub.cuh>
00011 #endif
00012 #include "ComputeBondedCUDAKernel.h"
00013 
00014 #ifdef FUTURE_CUDADEVICE
00015 #include "CudaDevice.h"
00016 #else
00017 #include "DeviceCUDA.h"
00018 extern __thread DeviceCUDA *deviceCUDA;
00019 #endif
00020 
00021 __device__ inline long long int lliroundf(float f)
00022 {
00023     long long int l;
00024     asm("cvt.rni.s64.f32  %0, %1;" : "=l"(l) : "f"(f));
00025     return l;
00026 }
00027 
00028 __device__ inline unsigned long long int llitoulli(long long int l)
00029 {
00030     unsigned long long int u;
00031     asm("mov.b64    %0, %1;" : "=l"(u) : "l"(l));
00032     return u;
00033 }
00034 
00035 template <typename T>
00036 __forceinline__ __device__
00037 void convertForces(const float x, const float y, const float z,
00038   T &Tx, T &Ty, T &Tz) {
00039 
00040   Tx = (T)(x);
00041   Ty = (T)(y);
00042   Tz = (T)(z);
00043 }
00044 
00045 template <>
00046 __forceinline__ __device__
00047 void convertForces<long long int>(const float x, const float y, const float z,
00048   long long int &Tx, long long int &Ty, long long int &Tz) {
00049 
00050   Tx = lliroundf(x*float_to_force);
00051   Ty = lliroundf(y*float_to_force);
00052   Tz = lliroundf(z*float_to_force);
00053 }
00054 
00055 template <typename T>
00056 __forceinline__ __device__
00057 void calcComponentForces(
00058   float fij,
00059   const float dx, const float dy, const float dz,
00060   T &fxij, T &fyij, T &fzij) {
00061 
00062   fxij = (T)(fij*dx);
00063   fyij = (T)(fij*dy);
00064   fzij = (T)(fij*dz);
00065 
00066 }
00067 
00068 template <>
00069 __forceinline__ __device__
00070 void calcComponentForces<long long int>(
00071   float fij,
00072   const float dx, const float dy, const float dz,
00073   long long int &fxij,
00074   long long int &fyij,
00075   long long int &fzij) {
00076 
00077   fij *= float_to_force;
00078   fxij = lliroundf(fij*dx);
00079   fyij = lliroundf(fij*dy);
00080   fzij = lliroundf(fij*dz);
00081 
00082 }
00083 
00084 template <typename T>
00085 __forceinline__ __device__
00086 void calcComponentForces(
00087   float fij1,
00088   const float dx1, const float dy1, const float dz1,
00089   float fij2,
00090   const float dx2, const float dy2, const float dz2,
00091   T &fxij, T &fyij, T &fzij) {
00092 
00093   fxij = (T)(fij1*dx1 + fij2*dx2);
00094   fyij = (T)(fij1*dy1 + fij2*dy2);
00095   fzij = (T)(fij1*dz1 + fij2*dz2);
00096 }
00097 
00098 template <>
00099 __forceinline__ __device__
00100 void calcComponentForces<long long int>(
00101   float fij1,
00102   const float dx1,
00103   const float dy1,
00104   const float dz1,
00105   float fij2,
00106   const float dx2,
00107   const float dy2,
00108   const float dz2,
00109   long long int &fxij,
00110   long long int &fyij,
00111   long long int &fzij) {
00112 
00113   fij1 *= float_to_force;
00114   fij2 *= float_to_force;
00115   fxij = lliroundf(fij1*dx1 + fij2*dx2);
00116   fyij = lliroundf(fij1*dy1 + fij2*dy2);
00117   fzij = lliroundf(fij1*dz1 + fij2*dz2);
00118 }
00119 
00120 template <typename T>
00121 __forceinline__ __device__
00122 void storeForces(const T fx, const T fy, const T fz,
00123      const int ind, const int stride,
00124      T* force) {
00125   // The generic version can not be used
00126 }
00127 
00128 // Template specialization for 64bit integer = "long long int"
00129 template <>
00130 __forceinline__ __device__ 
00131 void storeForces <long long int> (const long long int fx,
00132           const long long int fy,
00133           const long long int fz,
00134           const int ind, const int stride,
00135           long long int* force) {
00136   atomicAdd((unsigned long long int *)&force[ind           ], llitoulli(fx));
00137   atomicAdd((unsigned long long int *)&force[ind + stride  ], llitoulli(fy));
00138   atomicAdd((unsigned long long int *)&force[ind + stride*2], llitoulli(fz));
00139 }
00140 
00141 //
00142 // Calculates bonds
00143 //
00144 template <typename T, bool doEnergy, bool doVirial>
00145 __device__ void bondForce(
00146   const int index,
00147   const CudaBond* __restrict__ bondList,
00148   const CudaBondValue* __restrict__ bondValues,
00149   const float4* __restrict__ xyzq,
00150   const int stride,
00151   const float3 lata, const float3 latb, const float3 latc,
00152   T* __restrict__ force, double &energy,
00153 #ifdef WRITE_FULL_VIRIALS
00154   ComputeBondedCUDAKernel::BondedVirial<double>& virial
00155 #else
00156   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virial
00157 #endif
00158   ) {
00159 
00160   CudaBond bl = bondList[index];
00161   CudaBondValue bondValue = bondValues[bl.itype];
00162   if (bondValue.x0 == 0.0f) return;
00163 
00164   float shx = bl.ioffsetXYZ.x*lata.x + bl.ioffsetXYZ.y*latb.x + bl.ioffsetXYZ.z*latc.x;
00165   float shy = bl.ioffsetXYZ.x*lata.y + bl.ioffsetXYZ.y*latb.y + bl.ioffsetXYZ.z*latc.y;
00166   float shz = bl.ioffsetXYZ.x*lata.z + bl.ioffsetXYZ.y*latb.z + bl.ioffsetXYZ.z*latc.z;
00167 
00168   float4 xyzqi = xyzq[bl.i];
00169   float4 xyzqj = xyzq[bl.j];
00170 
00171   float xij = xyzqi.x + shx - xyzqj.x;
00172   float yij = xyzqi.y + shy - xyzqj.y;
00173   float zij = xyzqi.z + shz - xyzqj.z;
00174 
00175   float r = sqrtf(xij*xij + yij*yij + zij*zij);
00176 
00177   float db = r - bondValue.x0;
00178   if (bondValue.x1) {
00179     // in this case, the bond represents a harmonic wall potential
00180     // where x0 is the lower wall and x1 is the upper
00181     db = (r > bondValue.x1 ? r - bondValue.x1 : (r > bondValue.x0 ? 0 : db));
00182   }
00183   float fij = db * bondValue.k * bl.scale;
00184  
00185   if (doEnergy) {
00186     energy += (double)(fij*db);
00187   }
00188   fij *= -2.0f/r;
00189   
00190   T T_fxij, T_fyij, T_fzij;
00191   calcComponentForces<T>(fij, xij, yij, zij, T_fxij, T_fyij, T_fzij);
00192   
00193   // Store forces
00194   storeForces<T>(T_fxij, T_fyij, T_fzij, bl.i, stride, force);
00195   storeForces<T>(-T_fxij, -T_fyij, -T_fzij, bl.j, stride, force);
00196 
00197   // Store virial
00198   if (doVirial) {
00199 #ifdef WRITE_FULL_VIRIALS
00200     float fxij = fij*xij;
00201     float fyij = fij*yij;
00202     float fzij = fij*zij;
00203     virial.xx = (fxij*xij);
00204     virial.xy = (fxij*yij);
00205     virial.xz = (fxij*zij);
00206     virial.yx = (fyij*xij);
00207     virial.yy = (fyij*yij);
00208     virial.yz = (fyij*zij);
00209     virial.zx = (fzij*xij);
00210     virial.zy = (fzij*yij);
00211     virial.zz = (fzij*zij);
00212 #endif
00213   }
00214 }
00215 
00216 //
00217 // Calculates modified exclusions
00218 //
00219 template <typename T, bool doEnergy, bool doVirial, bool doElect>
00220 __device__ void modifiedExclusionForce(
00221   const int index,
00222   const CudaExclusion* __restrict__ exclusionList,
00223   const bool doSlow,
00224   const float one_scale14,                // 1 - scale14
00225   const int vdwCoefTableWidth,
00226 #if __CUDA_ARCH__ >= 350
00227   const float2* __restrict__ vdwCoefTable,
00228 #else
00229   cudaTextureObject_t vdwCoefTableTex, 
00230 #endif
00231   cudaTextureObject_t forceTableTex, cudaTextureObject_t energyTableTex,
00232   const float4* __restrict__ xyzq,
00233   const int stride,
00234   const float3 lata, const float3 latb, const float3 latc,
00235   const float cutoff2,
00236   double &energyVdw,
00237   T* __restrict__ forceNbond, double &energyNbond,
00238   T* __restrict__ forceSlow, double &energySlow,
00239 #ifdef WRITE_FULL_VIRIALS
00240   ComputeBondedCUDAKernel::BondedVirial<double>& virialNbond,
00241   ComputeBondedCUDAKernel::BondedVirial<double>& virialSlow
00242 #else
00243   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virialNbond,
00244   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virialSlow
00245 #endif
00246   ) {
00247 
00248   CudaExclusion bl = exclusionList[index];
00249 
00250   float shx = bl.ioffsetXYZ.x*lata.x + bl.ioffsetXYZ.y*latb.x + bl.ioffsetXYZ.z*latc.x;
00251   float shy = bl.ioffsetXYZ.x*lata.y + bl.ioffsetXYZ.y*latb.y + bl.ioffsetXYZ.z*latc.y;
00252   float shz = bl.ioffsetXYZ.x*lata.z + bl.ioffsetXYZ.y*latb.z + bl.ioffsetXYZ.z*latc.z;
00253 
00254   float4 xyzqi = xyzq[bl.i];
00255   float4 xyzqj = xyzq[bl.j];
00256 
00257   float xij = xyzqi.x + shx - xyzqj.x;
00258   float yij = xyzqi.y + shy - xyzqj.y;
00259   float zij = xyzqi.z + shz - xyzqj.z;
00260 
00261   float r2 = xij*xij + yij*yij + zij*zij;
00262   if (r2 < cutoff2) {
00263 
00264     float rinv = rsqrtf(r2);
00265 
00266     float qq;
00267     if (doElect) qq = one_scale14 * xyzqi.w * xyzqj.w;
00268 
00269     int vdwIndex = bl.vdwtypei + bl.vdwtypej*vdwCoefTableWidth;
00270 #if __CUDA_ARCH__ >= 350
00271     float2 ljab = __ldg(&vdwCoefTable[vdwIndex]);
00272 #else
00273     float2 ljab = tex1Dfetch<float2>(vdwCoefTableTex, vdwIndex);
00274 #endif
00275 
00276     float4 fi = tex1D<float4>(forceTableTex, rinv);
00277     float4 ei;
00278     if (doEnergy) {
00279       ei = tex1D<float4>(energyTableTex, rinv);
00280       energyVdw   += (double)(ljab.x * ei.z + ljab.y * ei.y);
00281       if (doElect) {
00282         energyNbond += qq * ei.x;
00283         if (doSlow) energySlow  += qq * ei.w;
00284       }
00285     }
00286 
00287     float fNbond;
00288     if (doElect) {
00289       fNbond = -(ljab.x * fi.z + ljab.y * fi.y + qq * fi.x);
00290     } else {
00291       fNbond = -(ljab.x * fi.z + ljab.y * fi.y);
00292     }
00293     T T_fxij, T_fyij, T_fzij;
00294     calcComponentForces<T>(fNbond, xij, yij, zij, T_fxij, T_fyij, T_fzij);
00295     storeForces<T>(T_fxij, T_fyij, T_fzij, bl.i, stride, forceNbond);
00296     storeForces<T>(-T_fxij, -T_fyij, -T_fzij, bl.j, stride, forceNbond);
00297 
00298     float fSlow;
00299     if (doSlow && doElect) {
00300       fSlow = -qq * fi.w;
00301       T T_fxij, T_fyij, T_fzij;
00302       calcComponentForces<T>(fSlow, xij, yij, zij, T_fxij, T_fyij, T_fzij);
00303       storeForces<T>(T_fxij, T_fyij, T_fzij, bl.i, stride, forceSlow);
00304       storeForces<T>(-T_fxij, -T_fyij, -T_fzij, bl.j, stride, forceSlow);
00305     }
00306 
00307     // Store virial
00308     if (doVirial) {
00309 #ifdef WRITE_FULL_VIRIALS
00310       float fxij = fNbond*xij;
00311       float fyij = fNbond*yij;
00312       float fzij = fNbond*zij;
00313       virialNbond.xx = (fxij*xij);
00314       virialNbond.xy = (fxij*yij);
00315       virialNbond.xz = (fxij*zij);
00316       virialNbond.yx = (fyij*xij);
00317       virialNbond.yy = (fyij*yij);
00318       virialNbond.yz = (fyij*zij);
00319       virialNbond.zx = (fzij*xij);
00320       virialNbond.zy = (fzij*yij);
00321       virialNbond.zz = (fzij*zij);
00322 #endif
00323     }
00324 
00325     // Store virial
00326     if (doVirial && doSlow && doElect) {
00327 #ifdef WRITE_FULL_VIRIALS
00328       float fxij = fSlow*xij;
00329       float fyij = fSlow*yij;
00330       float fzij = fSlow*zij;
00331       virialSlow.xx = (fxij*xij);
00332       virialSlow.xy = (fxij*yij);
00333       virialSlow.xz = (fxij*zij);
00334       virialSlow.yx = (fyij*xij);
00335       virialSlow.yy = (fyij*yij);
00336       virialSlow.yz = (fyij*zij);
00337       virialSlow.zx = (fzij*xij);
00338       virialSlow.zy = (fzij*yij);
00339       virialSlow.zz = (fzij*zij);
00340 #endif
00341     }
00342 
00343   }
00344 }
00345 
00346 //
00347 // Calculates exclusions. Here doSlow = true
00348 //
00349 template <typename T, bool doEnergy, bool doVirial>
00350 __device__ void exclusionForce(
00351   const int index,
00352   const CudaExclusion* __restrict__ exclusionList,
00353   const float r2_delta, const int r2_delta_expc,
00354 #if __CUDA_ARCH__ >= 350
00355   const float* __restrict__ r2_table,
00356   const float4* __restrict__ exclusionTable,
00357 #else
00358   cudaTextureObject_t r2_table_tex,
00359   cudaTextureObject_t exclusionTableTex,
00360 #endif
00361   const float4* __restrict__ xyzq,
00362   const int stride,
00363   const float3 lata, const float3 latb, const float3 latc,
00364   const float cutoff2,
00365   T* __restrict__ forceSlow, double &energySlow,
00366 #ifdef WRITE_FULL_VIRIALS
00367   ComputeBondedCUDAKernel::BondedVirial<double>& virialSlow
00368 #else
00369   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virialSlow
00370 #endif
00371   ) {
00372 
00373   CudaExclusion bl = exclusionList[index];
00374 
00375   float shx = bl.ioffsetXYZ.x*lata.x + bl.ioffsetXYZ.y*latb.x + bl.ioffsetXYZ.z*latc.x;
00376   float shy = bl.ioffsetXYZ.x*lata.y + bl.ioffsetXYZ.y*latb.y + bl.ioffsetXYZ.z*latc.y;
00377   float shz = bl.ioffsetXYZ.x*lata.z + bl.ioffsetXYZ.y*latb.z + bl.ioffsetXYZ.z*latc.z;
00378 
00379   float4 xyzqi = xyzq[bl.i];
00380   float4 xyzqj = xyzq[bl.j];
00381 
00382   float xij = xyzqi.x + shx - xyzqj.x;
00383   float yij = xyzqi.y + shy - xyzqj.y;
00384   float zij = xyzqi.z + shz - xyzqj.z;
00385 
00386   float r2 = xij*xij + yij*yij + zij*zij;
00387   if (r2 < cutoff2) {
00388     r2 += r2_delta;
00389 
00390     union { float f; int i; } r2i;
00391     r2i.f = r2;
00392     int table_i = (r2i.i >> 17) + r2_delta_expc;  // table_i >= 0
00393 
00394 #if __CUDA_ARCH__ >= 350
00395     float r2_table_val = __ldg(&r2_table[table_i]);
00396 #else
00397     float r2_table_val = tex1Dfetch<float>(r2_table_tex, table_i);
00398 #endif
00399     float diffa = r2 - r2_table_val;
00400     float qq = xyzqi.w * xyzqj.w;
00401 
00402 #if __CUDA_ARCH__ >= 350
00403     float4 slow = __ldg(&exclusionTable[table_i]);
00404 #else
00405     float4 slow = tex1Dfetch<float4>(exclusionTableTex, table_i);
00406 #endif
00407 
00408     if (doEnergy) {
00409       energySlow += (double)(qq*(((diffa * (1.0f/6.0f)*slow.x + 0.25f*slow.y ) * diffa + 0.5f*slow.z ) * diffa + slow.w));
00410     }
00411 
00412     float fSlow = -qq*((diffa * slow.x + slow.y) * diffa + slow.z);
00413 
00414     T T_fxij, T_fyij, T_fzij;
00415     calcComponentForces<T>(fSlow, xij, yij, zij, T_fxij, T_fyij, T_fzij);
00416     storeForces<T>(T_fxij, T_fyij, T_fzij, bl.i, stride, forceSlow);
00417     storeForces<T>(-T_fxij, -T_fyij, -T_fzij, bl.j, stride, forceSlow);
00418 
00419     // Store virial
00420     if (doVirial) {
00421 #ifdef WRITE_FULL_VIRIALS
00422       float fxij = fSlow*xij;
00423       float fyij = fSlow*yij;
00424       float fzij = fSlow*zij;
00425       virialSlow.xx = (fxij*xij);
00426       virialSlow.xy = (fxij*yij);
00427       virialSlow.xz = (fxij*zij);
00428       virialSlow.yx = (fyij*xij);
00429       virialSlow.yy = (fyij*yij);
00430       virialSlow.yz = (fyij*zij);
00431       virialSlow.zx = (fzij*xij);
00432       virialSlow.zy = (fzij*yij);
00433       virialSlow.zz = (fzij*zij);
00434 #endif
00435     }
00436   }
00437 }
00438 
00439 template <typename T, bool doEnergy, bool doVirial>
00440 __device__ void angleForce(const int index,
00441   const CudaAngle* __restrict__ angleList,
00442   const CudaAngleValue* __restrict__ angleValues,
00443   const float4* __restrict__ xyzq,
00444   const int stride,
00445   const float3 lata, const float3 latb, const float3 latc,
00446   T* __restrict__ force, double &energy,
00447 #ifdef WRITE_FULL_VIRIALS
00448   ComputeBondedCUDAKernel::BondedVirial<double>& virial
00449 #else
00450   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virial
00451 #endif
00452   ) {
00453 
00454   CudaAngle al = angleList[index];
00455 
00456   float ishx = al.ioffsetXYZ.x*lata.x + al.ioffsetXYZ.y*latb.x + al.ioffsetXYZ.z*latc.x;
00457   float ishy = al.ioffsetXYZ.x*lata.y + al.ioffsetXYZ.y*latb.y + al.ioffsetXYZ.z*latc.y;
00458   float ishz = al.ioffsetXYZ.x*lata.z + al.ioffsetXYZ.y*latb.z + al.ioffsetXYZ.z*latc.z;
00459 
00460   float kshx = al.koffsetXYZ.x*lata.x + al.koffsetXYZ.y*latb.x + al.koffsetXYZ.z*latc.x;
00461   float kshy = al.koffsetXYZ.x*lata.y + al.koffsetXYZ.y*latb.y + al.koffsetXYZ.z*latc.y;
00462   float kshz = al.koffsetXYZ.x*lata.z + al.koffsetXYZ.y*latb.z + al.koffsetXYZ.z*latc.z;
00463 
00464   float xij = xyzq[al.i].x + ishx - xyzq[al.j].x;
00465   float yij = xyzq[al.i].y + ishy - xyzq[al.j].y;
00466   float zij = xyzq[al.i].z + ishz - xyzq[al.j].z;
00467 
00468   float xkj = xyzq[al.k].x + kshx - xyzq[al.j].x;
00469   float ykj = xyzq[al.k].y + kshy - xyzq[al.j].y;
00470   float zkj = xyzq[al.k].z + kshz - xyzq[al.j].z;
00471 
00472   float rij_inv = rsqrtf(xij*xij + yij*yij + zij*zij);
00473   float rkj_inv = rsqrtf(xkj*xkj + ykj*ykj + zkj*zkj);
00474 
00475   float xijr = xij*rij_inv;
00476   float yijr = yij*rij_inv;
00477   float zijr = zij*rij_inv;
00478   float xkjr = xkj*rkj_inv;
00479   float ykjr = ykj*rkj_inv;
00480   float zkjr = zkj*rkj_inv;
00481   float cos_theta = xijr*xkjr + yijr*ykjr + zijr*zkjr;
00482 
00483   CudaAngleValue angleValue = angleValues[al.itype];
00484   angleValue.k *= al.scale;
00485 
00486   float diff;
00487   if (angleValue.normal == 1) {
00488     // Restrict values of cos_theta to the interval [-0.999 ... 0.999]
00489     cos_theta = min(0.999f, max(-0.999f, cos_theta));
00490     float theta = acosf(cos_theta);
00491     diff = theta - angleValue.theta0;
00492   } else {
00493     diff = cos_theta - angleValue.theta0;
00494   }
00495 
00496   if (doEnergy) {
00497     energy += (double)(angleValue.k * diff * diff);
00498   }
00499   if (angleValue.normal == 1) {
00500     float inv_sin_theta = rsqrtf(1.0f - cos_theta*cos_theta);
00501     if (inv_sin_theta > 1.0e6) {
00502       diff = (diff < 0.0f) ? 1.0f : -1.0f;
00503     } else {
00504       diff *= -inv_sin_theta;
00505     }
00506   }
00507   diff *= -2.0f*angleValue.k;
00508 
00509   float dtxi = rij_inv*(xkjr - cos_theta*xijr);
00510   float dtxj = rkj_inv*(xijr - cos_theta*xkjr);
00511   float dtyi = rij_inv*(ykjr - cos_theta*yijr);
00512   float dtyj = rkj_inv*(yijr - cos_theta*ykjr);
00513   float dtzi = rij_inv*(zkjr - cos_theta*zijr);
00514   float dtzj = rkj_inv*(zijr - cos_theta*zkjr);
00515 
00516   T T_dtxi, T_dtyi, T_dtzi;
00517   T T_dtxj, T_dtyj, T_dtzj;
00518   calcComponentForces<T>(diff, dtxi, dtyi, dtzi, T_dtxi, T_dtyi, T_dtzi);
00519   calcComponentForces<T>(diff, dtxj, dtyj, dtzj, T_dtxj, T_dtyj, T_dtzj);
00520   T T_dtxk = -T_dtxi - T_dtxj;
00521   T T_dtyk = -T_dtyi - T_dtyj;
00522   T T_dtzk = -T_dtzi - T_dtzj;
00523   storeForces<T>(T_dtxk, T_dtyk, T_dtzk, al.j, stride, force);
00524 
00525   if (angleValue.k_ub) {
00526     float xik = xij - xkj;
00527     float yik = yij - ykj;
00528     float zik = zij - zkj;
00529     float rik_inv = rsqrtf(xik*xik + yik*yik + zik*zik);
00530     float rik = 1.0f/rik_inv;
00531     float diff_ub = rik - angleValue.r_ub;
00532     if (doEnergy) {
00533       energy += (double)(angleValue.k_ub * diff_ub * diff_ub);
00534     }
00535     diff_ub *= -2.0f*angleValue.k_ub*rik_inv;
00536     T T_dxik, T_dyik, T_dzik;
00537     calcComponentForces<T>(diff_ub, xik, yik, zik, T_dxik, T_dyik, T_dzik);
00538     T_dtxi += T_dxik;
00539     T_dtyi += T_dyik;
00540     T_dtzi += T_dzik;
00541     T_dtxj -= T_dxik;
00542     T_dtyj -= T_dyik;
00543     T_dtzj -= T_dzik;
00544   }
00545 
00546   storeForces<T>(T_dtxi, T_dtyi, T_dtzi, al.i, stride, force);
00547   storeForces<T>(T_dtxj, T_dtyj, T_dtzj, al.k, stride, force);
00548 
00549   // Store virial
00550   if (doVirial) {
00551 #ifdef WRITE_FULL_VIRIALS
00552     float fxi = ((float)T_dtxi)*force_to_float;
00553     float fyi = ((float)T_dtyi)*force_to_float;
00554     float fzi = ((float)T_dtzi)*force_to_float;
00555     float fxk = ((float)T_dtxj)*force_to_float;
00556     float fyk = ((float)T_dtyj)*force_to_float;
00557     float fzk = ((float)T_dtzj)*force_to_float;
00558     virial.xx = (fxi*xij) + (fxk*xkj);
00559     virial.xy = (fxi*yij) + (fxk*ykj);
00560     virial.xz = (fxi*zij) + (fxk*zkj);
00561     virial.yx = (fyi*xij) + (fyk*xkj);
00562     virial.yy = (fyi*yij) + (fyk*ykj);
00563     virial.yz = (fyi*zij) + (fyk*zkj);
00564     virial.zx = (fzi*xij) + (fzk*xkj);
00565     virial.zy = (fzi*yij) + (fzk*ykj);
00566     virial.zz = (fzi*zij) + (fzk*zkj);
00567 #endif
00568   }
00569 }
00570 
00571 //
00572 // Dihedral computation
00573 //
00574 // Out: df, e
00575 //
00576 template <bool doEnergy>
00577 __forceinline__ __device__
00578 void diheComp(const CudaDihedralValue* dihedralValues, int ic,
00579   const float sin_phi, const float cos_phi, const float scale, float& df, double& e) {
00580 
00581   df = 0.0f;
00582   if (doEnergy) e = 0.0;
00583 
00584   float phi = atan2f(sin_phi, cos_phi);
00585 
00586   bool lrep = true;
00587   while (lrep) {
00588     CudaDihedralValue dihedralValue = dihedralValues[ic];
00589     dihedralValue.k *= scale;
00590 
00591     // Last dihedral has n low bit set to 0
00592     lrep = (dihedralValue.n & 1);
00593     dihedralValue.n >>= 1;
00594     if (dihedralValue.n) {
00595       float nf = dihedralValue.n;
00596       float x = nf*phi - dihedralValue.delta;
00597       if (doEnergy) {
00598         float sin_x, cos_x;
00599         sincosf(x, &sin_x, &cos_x);
00600         e += (double)(dihedralValue.k*(1.0f + cos_x));
00601         df += (double)(nf*dihedralValue.k*sin_x);
00602       } else {
00603         float sin_x = sinf(x);
00604         df += (double)(nf*dihedralValue.k*sin_x);
00605       }
00606     } else {
00607       float diff = phi - dihedralValue.delta;
00608       if (diff < -(float)(M_PI)) diff += (float)(2.0*M_PI);
00609       if (diff >  (float)(M_PI)) diff -= (float)(2.0*M_PI);
00610       if (doEnergy) {
00611         e += (double)(dihedralValue.k*diff*diff);
00612       }
00613       df -= 2.0f*dihedralValue.k*diff;
00614     }
00615     ic++;
00616   }
00617 
00618 }
00619 
00620 
00621 template <typename T, bool doEnergy, bool doVirial>
00622 __device__ void diheForce(const int index,
00623   const CudaDihedral* __restrict__ diheList,
00624   const CudaDihedralValue* __restrict__ dihedralValues,
00625   const float4* __restrict__ xyzq,
00626   const int stride,
00627   const float3 lata, const float3 latb, const float3 latc,
00628   T* __restrict__ force, double &energy,
00629 #ifdef WRITE_FULL_VIRIALS
00630   ComputeBondedCUDAKernel::BondedVirial<double>& virial
00631 #else
00632   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virial
00633 #endif
00634   ) {
00635 
00636   CudaDihedral dl = diheList[index];
00637 
00638   float ishx = dl.ioffsetXYZ.x*lata.x + dl.ioffsetXYZ.y*latb.x + dl.ioffsetXYZ.z*latc.x;
00639   float ishy = dl.ioffsetXYZ.x*lata.y + dl.ioffsetXYZ.y*latb.y + dl.ioffsetXYZ.z*latc.y;
00640   float ishz = dl.ioffsetXYZ.x*lata.z + dl.ioffsetXYZ.y*latb.z + dl.ioffsetXYZ.z*latc.z;
00641 
00642   float jshx = dl.joffsetXYZ.x*lata.x + dl.joffsetXYZ.y*latb.x + dl.joffsetXYZ.z*latc.x;
00643   float jshy = dl.joffsetXYZ.x*lata.y + dl.joffsetXYZ.y*latb.y + dl.joffsetXYZ.z*latc.y;
00644   float jshz = dl.joffsetXYZ.x*lata.z + dl.joffsetXYZ.y*latb.z + dl.joffsetXYZ.z*latc.z;
00645 
00646   float lshx = dl.loffsetXYZ.x*lata.x + dl.loffsetXYZ.y*latb.x + dl.loffsetXYZ.z*latc.x;
00647   float lshy = dl.loffsetXYZ.x*lata.y + dl.loffsetXYZ.y*latb.y + dl.loffsetXYZ.z*latc.y;
00648   float lshz = dl.loffsetXYZ.x*lata.z + dl.loffsetXYZ.y*latb.z + dl.loffsetXYZ.z*latc.z;
00649 
00650   float xij = xyzq[dl.i].x + ishx - xyzq[dl.j].x;
00651   float yij = xyzq[dl.i].y + ishy - xyzq[dl.j].y;
00652   float zij = xyzq[dl.i].z + ishz - xyzq[dl.j].z;
00653 
00654   float xjk = xyzq[dl.j].x + jshx - xyzq[dl.k].x;
00655   float yjk = xyzq[dl.j].y + jshy - xyzq[dl.k].y;
00656   float zjk = xyzq[dl.j].z + jshz - xyzq[dl.k].z;
00657 
00658   float xlk = xyzq[dl.l].x + lshx - xyzq[dl.k].x;
00659   float ylk = xyzq[dl.l].y + lshy - xyzq[dl.k].y;
00660   float zlk = xyzq[dl.l].z + lshz - xyzq[dl.k].z;
00661 
00662   // A=F^G, B=H^G.
00663   float ax = yij*zjk - zij*yjk;
00664   float ay = zij*xjk - xij*zjk;
00665   float az = xij*yjk - yij*xjk;
00666   float bx = ylk*zjk - zlk*yjk;
00667   float by = zlk*xjk - xlk*zjk;
00668   float bz = xlk*yjk - ylk*xjk;
00669 
00670   float ra2 = ax*ax + ay*ay + az*az;
00671   float rb2 = bx*bx + by*by + bz*bz;
00672   float rg = sqrtf(xjk*xjk + yjk*yjk + zjk*zjk);
00673 
00674   //    if((ra2 <= rxmin2) .or. (rb2 <= rxmin2) .or. (rg <= rxmin)) then
00675   //          nlinear = nlinear + 1
00676   //       endif
00677 
00678   float rgr = 1.0f / rg;
00679   float ra2r = 1.0f / ra2;
00680   float rb2r = 1.0f / rb2;
00681   float rabr = sqrtf(ra2r*rb2r);
00682 
00683   float cos_phi = (ax*bx + ay*by + az*bz)*rabr;
00684   //
00685   // Note that sin(phi).G/|G|=B^A/(|A|.|B|)
00686   // which can be simplify to sin(phi)=|G|H.A/(|A|.|B|)
00687   float sin_phi = rg*rabr*(ax*xlk + ay*ylk + az*zlk);
00688   //
00689   //     Energy and derivative contributions.
00690 
00691   float df;
00692   double e;
00693   diheComp<doEnergy>(dihedralValues, dl.itype, sin_phi, cos_phi, dl.scale, df, e);
00694 
00695   if (doEnergy) energy += e;
00696 
00697   //
00698   //     Compute derivatives wrt catesian coordinates.
00699   //
00700   // GAA=dE/dphi.|G|/A^2, GBB=dE/dphi.|G|/B^2, FG=F.G, HG=H.G
00701   //  FGA=dE/dphi*F.G/(|G|A^2), HGB=dE/dphi*H.G/(|G|B^2)
00702 
00703   float fg = xij*xjk + yij*yjk + zij*zjk;
00704   float hg = xlk*xjk + ylk*yjk + zlk*zjk;
00705   ra2r *= df;
00706   rb2r *= df;
00707   float fga = fg*ra2r*rgr;
00708   float hgb = hg*rb2r*rgr;
00709   float gaa = ra2r*rg;
00710   float gbb = rb2r*rg;
00711   // DFi=dE/dFi, DGi=dE/dGi, DHi=dE/dHi.
00712 
00713   // Store forces
00714   T T_fix, T_fiy, T_fiz;
00715   calcComponentForces<T>(-gaa, ax, ay, az, T_fix, T_fiy, T_fiz);
00716   storeForces<T>(T_fix, T_fiy, T_fiz, dl.i, stride, force);
00717 
00718   T dgx, dgy, dgz;
00719   calcComponentForces<T>(fga, ax, ay, az, -hgb, bx, by, bz, dgx, dgy, dgz);
00720   T T_fjx = dgx - T_fix;
00721   T T_fjy = dgy - T_fiy;
00722   T T_fjz = dgz - T_fiz;
00723   storeForces<T>(T_fjx, T_fjy, T_fjz, dl.j, stride, force);
00724 
00725   T dhx, dhy, dhz;
00726   calcComponentForces<T>(gbb, bx, by, bz, dhx, dhy, dhz);
00727   T T_fkx = -dhx - dgx;
00728   T T_fky = -dhy - dgy;
00729   T T_fkz = -dhz - dgz;
00730   storeForces<T>(T_fkx, T_fky, T_fkz, dl.k, stride, force);
00731   storeForces<T>(dhx, dhy, dhz, dl.l, stride, force);
00732 
00733   // Store virial
00734   if (doVirial) {
00735 #ifdef WRITE_FULL_VIRIALS
00736     float fix = ((float)T_fix)*force_to_float;
00737     float fiy = ((float)T_fiy)*force_to_float;
00738     float fiz = ((float)T_fiz)*force_to_float;
00739     float fjx = ((float)dgx)*force_to_float;
00740     float fjy = ((float)dgy)*force_to_float;
00741     float fjz = ((float)dgz)*force_to_float;
00742     float flx = ((float)dhx)*force_to_float;
00743     float fly = ((float)dhy)*force_to_float;
00744     float flz = ((float)dhz)*force_to_float;
00745     virial.xx = (fix*xij) + (fjx*xjk) + (flx*xlk);
00746     virial.xy = (fix*yij) + (fjx*yjk) + (flx*ylk);
00747     virial.xz = (fix*zij) + (fjx*zjk) + (flx*zlk);
00748     virial.yx = (fiy*xij) + (fjy*xjk) + (fly*xlk);
00749     virial.yy = (fiy*yij) + (fjy*yjk) + (fly*ylk);
00750     virial.yz = (fiy*zij) + (fjy*zjk) + (fly*zlk);
00751     virial.zx = (fiz*xij) + (fjz*xjk) + (flz*xlk);
00752     virial.zy = (fiz*yij) + (fjz*yjk) + (flz*ylk);
00753     virial.zz = (fiz*zij) + (fjz*zjk) + (flz*zlk);
00754 #endif
00755   }
00756 
00757 }
00758 
00759 __device__ __forceinline__ float3 cross(const float3 v1, const float3 v2) {
00760  return make_float3(
00761   v1.y*v2.z-v2.y*v1.z,
00762   v2.x*v1.z-v1.x*v2.z,
00763   v1.x*v2.y-v2.x*v1.y
00764   );
00765 }
00766 
00767 __device__ __forceinline__ float dot(const float3 v1, const float3 v2) {
00768   return (v1.x * v2.x + v1.y * v2.y + v1.z * v2.z);
00769 }
00770 
00771 //
00772 // Calculates crossterms
00773 //
00774 template <typename T, bool doEnergy, bool doVirial>
00775 __device__ void crosstermForce(
00776   const int index,
00777   const CudaCrossterm* __restrict__ crosstermList,
00778   const CudaCrosstermValue* __restrict__ crosstermValues,
00779   const float4* __restrict__ xyzq,
00780   const int stride,
00781   const float3 lata, const float3 latb, const float3 latc,
00782   T* __restrict__ force, double &energy,
00783 #ifdef WRITE_FULL_VIRIALS
00784   ComputeBondedCUDAKernel::BondedVirial<double>& virial
00785 #else
00786   ComputeBondedCUDAKernel::BondedVirial* __restrict__ virial
00787 #endif
00788   ) {
00789 
00790   CudaCrossterm cl = crosstermList[index];
00791 
00792   // ----------------------------------------------------------------------------
00793   // Angle between 1 - 2 - 3 - 4
00794   //
00795   // 1 - 2
00796   float3 sh12 = make_float3(
00797     cl.offset12XYZ.x*lata.x + cl.offset12XYZ.y*latb.x + cl.offset12XYZ.z*latc.x,
00798     cl.offset12XYZ.x*lata.y + cl.offset12XYZ.y*latb.y + cl.offset12XYZ.z*latc.y,
00799     cl.offset12XYZ.x*lata.z + cl.offset12XYZ.y*latb.z + cl.offset12XYZ.z*latc.z);
00800 
00801   float4 xyzq1 = xyzq[cl.i1];
00802   float4 xyzq2 = xyzq[cl.i2];
00803 
00804   float3 r12 = make_float3(
00805     xyzq1.x + sh12.x - xyzq2.x,
00806     xyzq1.y + sh12.y - xyzq2.y,
00807     xyzq1.z + sh12.z - xyzq2.z);
00808 
00809   // 2 - 3
00810   float3 sh23 = make_float3(
00811     cl.offset23XYZ.x*lata.x + cl.offset23XYZ.y*latb.x + cl.offset23XYZ.z*latc.x,
00812     cl.offset23XYZ.x*lata.y + cl.offset23XYZ.y*latb.y + cl.offset23XYZ.z*latc.y,
00813     cl.offset23XYZ.x*lata.z + cl.offset23XYZ.y*latb.z + cl.offset23XYZ.z*latc.z);
00814 
00815   float4 xyzq3 = xyzq[cl.i3];
00816 
00817   float3 r23 = make_float3(
00818     xyzq2.x + sh23.x - xyzq3.x,
00819     xyzq2.y + sh23.y - xyzq3.y,
00820     xyzq2.z + sh23.z - xyzq3.z);
00821 
00822   // 3 - 4
00823   float3 sh34 = make_float3(
00824     cl.offset34XYZ.x*lata.x + cl.offset34XYZ.y*latb.x + cl.offset34XYZ.z*latc.x,
00825     cl.offset34XYZ.x*lata.y + cl.offset34XYZ.y*latb.y + cl.offset34XYZ.z*latc.y,
00826     cl.offset34XYZ.x*lata.z + cl.offset34XYZ.y*latb.z + cl.offset34XYZ.z*latc.z);
00827 
00828   float4 xyzq4 = xyzq[cl.i4];
00829 
00830   float3 r34 = make_float3(
00831     xyzq3.x + sh34.x - xyzq4.x,
00832     xyzq3.y + sh34.y - xyzq4.y,
00833     xyzq3.z + sh34.z - xyzq4.z);
00834 
00835   // Calculate the cross products
00836   float3 A = cross(r12, r23);
00837   float3 B = cross(r23, r34);
00838   float3 C = cross(r23, A);
00839 
00840   // Calculate the inverse distances
00841   float inv_rA = rsqrtf(dot(A, A));
00842   float inv_rB = rsqrtf(dot(B, B));
00843   float inv_rC = rsqrtf(dot(C, C));
00844 
00845   //  Calculate the sin and cos
00846   float cos_phi = dot(A, B)*(inv_rA*inv_rB);
00847   float sin_phi = dot(C, B)*(inv_rC*inv_rB);
00848 
00849   float phi = -atan2f(sin_phi,cos_phi);
00850 
00851   // ----------------------------------------------------------------------------
00852   // Angle between 5 - 6 - 7 - 8
00853   //
00854 
00855   // 5 - 6
00856   float3 sh56 = make_float3(
00857     cl.offset56XYZ.x*lata.x + cl.offset56XYZ.y*latb.x + cl.offset56XYZ.z*latc.x,
00858     cl.offset56XYZ.x*lata.y + cl.offset56XYZ.y*latb.y + cl.offset56XYZ.z*latc.y,
00859     cl.offset56XYZ.x*lata.z + cl.offset56XYZ.y*latb.z + cl.offset56XYZ.z*latc.z);
00860 
00861   float4 xyzq5 = xyzq[cl.i5];
00862   float4 xyzq6 = xyzq[cl.i6];
00863 
00864   float3 r56 = make_float3(
00865     xyzq5.x + sh56.x - xyzq6.x,
00866     xyzq5.y + sh56.y - xyzq6.y,
00867     xyzq5.z + sh56.z - xyzq6.z);
00868 
00869   // 6 - 7
00870   float3 sh67 = make_float3(
00871     cl.offset67XYZ.x*lata.x + cl.offset67XYZ.y*latb.x + cl.offset67XYZ.z*latc.x,
00872     cl.offset67XYZ.x*lata.y + cl.offset67XYZ.y*latb.y + cl.offset67XYZ.z*latc.y,
00873     cl.offset67XYZ.x*lata.z + cl.offset67XYZ.y*latb.z + cl.offset67XYZ.z*latc.z);
00874 
00875   float4 xyzq7 = xyzq[cl.i7];
00876 
00877   float3 r67 = make_float3(
00878     xyzq6.x + sh67.x - xyzq7.x,
00879     xyzq6.y + sh67.y - xyzq7.y,
00880     xyzq6.z + sh67.z - xyzq7.z);
00881 
00882   // 7 - 8
00883   float3 sh78 = make_float3(
00884     cl.offset78XYZ.x*lata.x + cl.offset78XYZ.y*latb.x + cl.offset78XYZ.z*latc.x,
00885     cl.offset78XYZ.x*lata.y + cl.offset78XYZ.y*latb.y + cl.offset78XYZ.z*latc.y,
00886     cl.offset78XYZ.x*lata.z + cl.offset78XYZ.y*latb.z + cl.offset78XYZ.z*latc.z);
00887 
00888   float4 xyzq8 = xyzq[cl.i8];
00889 
00890   float3 r78 = make_float3(
00891     xyzq7.x + sh78.x - xyzq8.x,
00892     xyzq7.y + sh78.y - xyzq8.y,
00893     xyzq7.z + sh78.z - xyzq8.z);
00894 
00895   // Calculate the cross products
00896   float3 D = cross(r56, r67);
00897   float3 E = cross(r67, r78);
00898   float3 F = cross(r67, D);
00899   
00900   // Calculate the inverse distances
00901   float inv_rD = rsqrtf(dot(D, D));
00902   float inv_rE = rsqrtf(dot(E, E));
00903   float inv_rF = rsqrtf(dot(F, F));
00904 
00905   //  Calculate the sin and cos
00906   float cos_psi = dot(D, E)*(inv_rD*inv_rE);
00907   float sin_psi = dot(F, E)*(inv_rF*inv_rE);
00908 
00909   float psi = -atan2f(sin_psi,cos_psi);
00910 
00911   // ----------------------------------------------------------------------------
00912   // Calculate the energy
00913 
00914   const float inv_h = (float)( (CudaCrosstermValue::dim) / (2.0*M_PI) );
00915 
00916   // Shift angles
00917   phi = fmod(phi + (float)M_PI, 2.0f*(float)M_PI);
00918   psi = fmod(psi + (float)M_PI, 2.0f*(float)M_PI);
00919 
00920   // distance measured in grid points between angle and smallest value
00921   float phi_h = phi * inv_h;
00922   float psi_h = psi * inv_h;
00923 
00924   // find smallest numbered grid point in stencil
00925   int iphi = (int)floor(phi_h);
00926   int ipsi = (int)floor(psi_h);
00927 
00928   const CudaCrosstermValue& cp = crosstermValues[cl.itype];
00929 
00930   // Load coefficients
00931   float4 c[4];
00932   c[0] = cp.c[iphi][ipsi][0];
00933   c[1] = cp.c[iphi][ipsi][1];
00934   c[2] = cp.c[iphi][ipsi][2];
00935   c[3] = cp.c[iphi][ipsi][3];
00936 
00937   float dphi = phi_h - (float)iphi;
00938   float dpsi = psi_h - (float)ipsi;
00939 
00940   if (doEnergy) {
00941     float energyf =          c[3].x + dphi*( c[3].y + dphi*( c[3].z + dphi*c[3].w ) );
00942     energyf = energyf*dpsi + c[2].x + dphi*( c[2].y + dphi*( c[2].z + dphi*c[2].w ) );
00943     energyf = energyf*dpsi + c[1].x + dphi*( c[1].y + dphi*( c[1].z + dphi*c[1].w ) );
00944     energyf = energyf*dpsi + c[0].x + dphi*( c[0].y + dphi*( c[0].z + dphi*c[0].w ) );
00945     energy += energyf*cl.scale;
00946   }
00947 
00948   float dEdphi =         3.0f*(c[0].w + dpsi*( c[1].w + dpsi*( c[2].w + dpsi*c[3].w ) ));
00949   dEdphi = dEdphi*dphi + 2.0f*(c[0].z + dpsi*( c[1].z + dpsi*( c[2].z + dpsi*c[3].z ) ));
00950   dEdphi = dEdphi*dphi +      (c[0].y + dpsi*( c[1].y + dpsi*( c[2].y + dpsi*c[3].y ) ));  // 13 muls
00951   dEdphi *= cl.scale*inv_h;
00952 
00953   float dEdpsi =         3.0f*(c[3].x + dphi*( c[3].y + dphi*( c[3].z + dphi*c[3].w ) ));
00954   dEdpsi = dEdpsi*dpsi + 2.0f*(c[2].x + dphi*( c[2].y + dphi*( c[2].z + dphi*c[2].w ) ));
00955   dEdpsi = dEdpsi*dpsi +      (c[1].x + dphi*( c[1].y + dphi*( c[1].z + dphi*c[1].w ) ));  // 13 muls
00956   dEdpsi *= cl.scale*inv_h;
00957 
00958   // float normCross1 = dot(A, A);
00959   float square_r23 = dot(r23, r23);
00960   float norm_r23 = sqrtf(square_r23);
00961   float inv_square_r23 = 1.0f/square_r23;
00962   float ff1 = dEdphi*norm_r23*inv_rA*inv_rA;
00963   float ff2 = -dot(r12, r23)*inv_square_r23;
00964   float ff3 = -dot(r34, r23)*inv_square_r23;
00965   float ff4 = -dEdphi*norm_r23*inv_rB*inv_rB;
00966 
00967   float3 f1 = make_float3(ff1*A.x, ff1*A.y, ff1*A.z);
00968   float3 f4 = make_float3(ff4*B.x, ff4*B.y, ff4*B.z);
00969   float3 t1 = make_float3( ff2*f1.x - ff3*f4.x, ff2*f1.y - ff3*f4.y, ff2*f1.z - ff3*f4.z );
00970   float3 f2 = make_float3(  t1.x - f1.x,  t1.y - f1.y,  t1.z - f1.z);
00971   float3 f3 = make_float3( -t1.x - f4.x, -t1.y - f4.y, -t1.z - f4.z);
00972 
00973   T T_f1x, T_f1y, T_f1z;
00974   T T_f2x, T_f2y, T_f2z;
00975   T T_f3x, T_f3y, T_f3z;
00976   T T_f4x, T_f4y, T_f4z;
00977   convertForces<T>(f1.x, f1.y, f1.z, T_f1x, T_f1y, T_f1z);
00978   convertForces<T>(f2.x, f2.y, f2.z, T_f2x, T_f2y, T_f2z);
00979   convertForces<T>(f3.x, f3.y, f3.z, T_f3x, T_f3y, T_f3z);
00980   convertForces<T>(f4.x, f4.y, f4.z, T_f4x, T_f4y, T_f4z);
00981   storeForces<T>(T_f1x, T_f1y, T_f1z, cl.i1, stride, force);
00982   storeForces<T>(T_f2x, T_f2y, T_f2z, cl.i2, stride, force);
00983   storeForces<T>(T_f3x, T_f3y, T_f3z, cl.i3, stride, force);
00984   storeForces<T>(T_f4x, T_f4y, T_f4z, cl.i4, stride, force);
00985 
00986   float square_r67 = dot(r67, r67);
00987   float norm_r67 = sqrtf(square_r67);
00988   float inv_square_r67 = 1.0f/(square_r67);
00989   ff1 = dEdpsi*norm_r67*inv_rD*inv_rD;
00990   ff2 = -dot(r56, r67)*inv_square_r67;
00991   ff3 = -dot(r78, r67)*inv_square_r67;
00992   ff4 = -dEdpsi*norm_r67*inv_rE*inv_rE;
00993   
00994   float3 f5 = make_float3( ff1*D.x, ff1*D.y, ff1*D.z );
00995   float3 f8 = make_float3( ff4*E.x, ff4*E.y, ff4*E.z );
00996   float3 t2 = make_float3( ff2*f5.x - ff3*f8.x, ff2*f5.y - ff3*f8.y, ff2*f5.z - ff3*f8.z );
00997   float3 f6 = make_float3( t2.x - f5.x,  t2.y - f5.y,  t2.z - f5.z );
00998   float3 f7 = make_float3(-t2.x - f8.x, -t2.y - f8.y, -t2.z - f8.z );
00999 
01000   T T_f5x, T_f5y, T_f5z;
01001   T T_f6x, T_f6y, T_f6z;
01002   T T_f7x, T_f7y, T_f7z;
01003   T T_f8x, T_f8y, T_f8z;
01004   convertForces<T>(f5.x, f5.y, f5.z, T_f5x, T_f5y, T_f5z);
01005   convertForces<T>(f6.x, f6.y, f6.z, T_f6x, T_f6y, T_f6z);
01006   convertForces<T>(f7.x, f7.y, f7.z, T_f7x, T_f7y, T_f7z);
01007   convertForces<T>(f8.x, f8.y, f8.z, T_f8x, T_f8y, T_f8z);
01008   storeForces<T>(T_f5x, T_f5y, T_f5z, cl.i5, stride, force);
01009   storeForces<T>(T_f6x, T_f6y, T_f6z, cl.i6, stride, force);
01010   storeForces<T>(T_f7x, T_f7y, T_f7z, cl.i7, stride, force);
01011   storeForces<T>(T_f8x, T_f8y, T_f8z, cl.i8, stride, force);
01012 
01013   // Store virial
01014   if (doVirial) {
01015 #ifdef WRITE_FULL_VIRIALS
01016     float3 s12 = make_float3( f1.x + f2.x, f1.y + f2.y, f1.z + f2.z );
01017     float3 s56 = make_float3( f5.x + f6.x, f5.y + f6.y, f5.z + f6.z );
01018     virial.xx = f1.x*r12.x + s12.x*r23.x - f4.x*r34.x + f5.x*r56.x + s56.x*r67.x - f8.x*r78.x;
01019     virial.xy = f1.x*r12.y + s12.x*r23.y - f4.x*r34.y + f5.x*r56.y + s56.x*r67.y - f8.x*r78.y;
01020     virial.xz = f1.x*r12.z + s12.x*r23.z - f4.x*r34.z + f5.x*r56.z + s56.x*r67.z - f8.x*r78.z;
01021     virial.yx = f1.y*r12.x + s12.y*r23.x - f4.y*r34.x + f5.y*r56.x + s56.y*r67.x - f8.y*r78.x;
01022     virial.yy = f1.y*r12.y + s12.y*r23.y - f4.y*r34.y + f5.y*r56.y + s56.y*r67.y - f8.y*r78.y;
01023     virial.yz = f1.y*r12.z + s12.y*r23.z - f4.y*r34.z + f5.y*r56.z + s56.y*r67.z - f8.y*r78.z;
01024     virial.zx = f1.z*r12.x + s12.z*r23.x - f4.z*r34.x + f5.z*r56.x + s56.z*r67.x - f8.z*r78.x;
01025     virial.zy = f1.z*r12.y + s12.z*r23.y - f4.z*r34.y + f5.z*r56.y + s56.z*r67.y - f8.z*r78.y;
01026     virial.zz = f1.z*r12.z + s12.z*r23.z - f4.z*r34.z + f5.z*r56.z + s56.z*r67.z - f8.z*r78.z;
01027   }
01028 #endif
01029 
01030 }
01031 
01032 #define BONDEDFORCESKERNEL_NUM_WARP 4
01033 //
01034 // Calculates all forces in a single kernel call
01035 //
01036 template <typename T, bool doEnergy, bool doVirial>
01037 __global__ void bondedForcesKernel(
01038   const int start,
01039 
01040   const int numBonds,
01041   const CudaBond* __restrict__ bonds,
01042   const CudaBondValue* __restrict__ bondValues,
01043 
01044   const int numAngles,
01045   const CudaAngle* __restrict__ angles,
01046   const CudaAngleValue* __restrict__ angleValues,
01047 
01048   const int numDihedrals,
01049   const CudaDihedral* __restrict__ dihedrals,
01050   const CudaDihedralValue* __restrict__ dihedralValues,
01051 
01052   const int numImpropers,
01053   const CudaDihedral* __restrict__ impropers,
01054   const CudaDihedralValue* __restrict__ improperValues,
01055 
01056   const int numExclusions,
01057   const CudaExclusion* __restrict__ exclusions,
01058 
01059   const int numCrossterms,
01060   const CudaCrossterm* __restrict__ crossterms,
01061   const CudaCrosstermValue* __restrict__ crosstermValues,
01062 
01063   const float cutoff2,
01064   const float r2_delta, const int r2_delta_expc,
01065 
01066   const float* __restrict__ r2_table,
01067   const float4* __restrict__ exclusionTable,
01068   cudaTextureObject_t r2_table_tex,
01069   cudaTextureObject_t exclusionTableTex,
01070 
01071   const float4* __restrict__ xyzq,
01072   const int stride,
01073   const float3 lata, const float3 latb, const float3 latc,
01074   T* __restrict__ force,
01075   T* __restrict__ forceSlow,
01076   double* __restrict__ energies_virials) {
01077 
01078   // Thread-block index
01079   int indexTB = start + blockIdx.x;
01080 
01081   const int numBondsTB     = (numBonds + blockDim.x - 1)/blockDim.x;
01082   const int numAnglesTB    = (numAngles + blockDim.x - 1)/blockDim.x;
01083   const int numDihedralsTB = (numDihedrals + blockDim.x - 1)/blockDim.x;
01084   const int numImpropersTB = (numImpropers + blockDim.x - 1)/blockDim.x;
01085   const int numExclusionsTB= (numExclusions + blockDim.x - 1)/blockDim.x;
01086   const int numCrosstermsTB= (numCrossterms + blockDim.x - 1)/blockDim.x;
01087 
01088   // Each thread computes single bonded interaction.
01089   // Each thread block computes single bonded type
01090   double energy;
01091   int energyIndex;
01092 
01093   if (doEnergy) {
01094     energy = 0.0;
01095     energyIndex = -1;
01096   }
01097 
01098 #ifdef WRITE_FULL_VIRIALS
01099   ComputeBondedCUDAKernel::BondedVirial<double> virial;
01100   int virialIndex;
01101   if (doVirial) {
01102     virial.xx = 0.0;
01103     virial.xy = 0.0;
01104     virial.xz = 0.0;
01105     virial.yx = 0.0;
01106     virial.yy = 0.0;
01107     virial.yz = 0.0;
01108     virial.zx = 0.0;
01109     virial.zy = 0.0;
01110     virial.zz = 0.0;
01111     virialIndex = ComputeBondedCUDAKernel::normalVirialIndex_XX;
01112   }
01113 #endif
01114 
01115   if (indexTB < numBondsTB) {
01116     int i = threadIdx.x + indexTB*blockDim.x;
01117     if (doEnergy) {
01118       energyIndex = ComputeBondedCUDAKernel::energyIndex_BOND;
01119     }
01120     if (i < numBonds) {
01121       bondForce<T, doEnergy, doVirial>
01122       (i, bonds, bondValues, xyzq,
01123         stride, lata, latb, latc,
01124         force, energy, virial);
01125     }
01126     goto reduce;
01127   }
01128   indexTB -= numBondsTB;
01129 
01130   if (indexTB < numAnglesTB) {
01131     int i = threadIdx.x + indexTB*blockDim.x;
01132     if (doEnergy) {
01133       energyIndex = ComputeBondedCUDAKernel::energyIndex_ANGLE;
01134     }
01135     if (i < numAngles) {
01136       angleForce<T, doEnergy, doVirial>
01137       (i, angles, angleValues, xyzq, stride,
01138         lata, latb, latc,
01139         force, energy, virial);
01140     }
01141     goto reduce;
01142   }
01143   indexTB -= numAnglesTB;
01144 
01145   if (indexTB < numDihedralsTB) {
01146     int i = threadIdx.x + indexTB*blockDim.x;
01147     if (doEnergy) {
01148       energyIndex = ComputeBondedCUDAKernel::energyIndex_DIHEDRAL;
01149     }
01150     if (doVirial) {
01151       virialIndex = ComputeBondedCUDAKernel::amdDiheVirialIndex_XX;
01152     }
01153     if (i < numDihedrals) {
01154       diheForce<T, doEnergy, doVirial>
01155       (i, dihedrals, dihedralValues, xyzq, stride,
01156         lata, latb, latc,
01157         force, energy, virial);
01158     }
01159     goto reduce;
01160   }
01161   indexTB -= numDihedralsTB;
01162 
01163   if (indexTB < numImpropersTB) {
01164     int i = threadIdx.x + indexTB*blockDim.x;
01165     if (doEnergy) {
01166       energyIndex = ComputeBondedCUDAKernel::energyIndex_IMPROPER;
01167     }
01168     if (i < numImpropers) {
01169       diheForce<T, doEnergy, doVirial>
01170       (i, impropers, improperValues, xyzq, stride,
01171         lata, latb, latc,
01172         force, energy, virial);
01173     }
01174     goto reduce;
01175   }
01176   indexTB -= numImpropersTB;
01177 
01178   if (indexTB < numExclusionsTB) {
01179     int i = threadIdx.x + indexTB*blockDim.x;
01180     if (doEnergy) {
01181       energyIndex = ComputeBondedCUDAKernel::energyIndex_ELECT_SLOW;
01182     }
01183     if (doVirial) {
01184       virialIndex = ComputeBondedCUDAKernel::slowVirialIndex_XX;
01185     }
01186     if (i < numExclusions) {
01187       exclusionForce<T, doEnergy, doVirial>
01188       (i, exclusions, r2_delta, r2_delta_expc,
01189 #if __CUDA_ARCH__ >= 350
01190         r2_table, exclusionTable,
01191 #else
01192         r2_table_tex, exclusionTableTex,
01193 #endif
01194         xyzq, stride, lata, latb, latc, cutoff2,
01195         forceSlow, energy, virial);
01196     }
01197     goto reduce;
01198   }
01199   indexTB -= numExclusionsTB;
01200 
01201   if (indexTB < numCrosstermsTB) {
01202     int i = threadIdx.x + indexTB*blockDim.x;
01203     if (doEnergy) {
01204       energyIndex = ComputeBondedCUDAKernel::energyIndex_CROSSTERM;
01205     }
01206     if (doVirial) {
01207       virialIndex = ComputeBondedCUDAKernel::amdDiheVirialIndex_XX;
01208     }
01209     if (i < numCrossterms) {
01210       crosstermForce<T, doEnergy, doVirial>
01211       (i, crossterms, crosstermValues,
01212         xyzq, stride, lata, latb, latc,
01213         force, energy, virial);
01214     }
01215     goto reduce;
01216   }
01217   // indexTB -= numCrosstermsTB;
01218 
01219   reduce:
01220 
01221   // Write energies to global memory
01222   if (doEnergy) {
01223     // energyIndex is constant within thread-block
01224     __shared__ double shEnergy[BONDEDFORCESKERNEL_NUM_WARP];
01225 #pragma unroll
01226     for (int i=16;i >= 1;i/=2) {
01227       energy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energy, i, 32);
01228     }
01229     int laneID = (threadIdx.x & (WARPSIZE - 1));
01230     int warpID = threadIdx.x / WARPSIZE;
01231     if (laneID == 0) {
01232       shEnergy[warpID] = energy;
01233     }
01234     BLOCK_SYNC;
01235     if (warpID == 0) {
01236       energy = (laneID < BONDEDFORCESKERNEL_NUM_WARP) ? shEnergy[laneID] : 0.0;
01237 #pragma unroll
01238       for (int i=16;i >= 1;i/=2) {
01239         energy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energy, i, 32);
01240       }
01241       if (laneID == 0) {
01242         atomicAdd(&energies_virials[energyIndex], energy);
01243       }
01244     }
01245   }
01246 
01247   // Write virials to global memory
01248 #ifdef WRITE_FULL_VIRIALS
01249   if (doVirial) {
01250 #pragma unroll
01251     for (int i=16;i >= 1;i/=2) {
01252       virial.xx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.xx, i, 32);
01253       virial.xy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.xy, i, 32);
01254       virial.xz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.xz, i, 32);
01255       virial.yx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.yx, i, 32);
01256       virial.yy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.yy, i, 32);
01257       virial.yz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.yz, i, 32);
01258       virial.zx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.zx, i, 32);
01259       virial.zy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.zy, i, 32);
01260       virial.zz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.zz, i, 32);
01261     }
01262     __shared__ ComputeBondedCUDAKernel::BondedVirial<double> shVirial[BONDEDFORCESKERNEL_NUM_WARP];
01263     int laneID = (threadIdx.x & (WARPSIZE - 1));
01264     int warpID = threadIdx.x / WARPSIZE;
01265     if (laneID == 0) {
01266       shVirial[warpID] = virial;
01267     }
01268     BLOCK_SYNC;
01269 
01270     if (warpID == 0) {
01271       virial.xx = 0.0;
01272       virial.xy = 0.0;
01273       virial.xz = 0.0;
01274       virial.yx = 0.0;
01275       virial.yy = 0.0;
01276       virial.yz = 0.0;
01277       virial.zx = 0.0;
01278       virial.zy = 0.0;
01279       virial.zz = 0.0;
01280       if (laneID < BONDEDFORCESKERNEL_NUM_WARP) virial = shVirial[laneID];
01281 #pragma unroll
01282       for (int i=16;i >= 1;i/=2) {
01283         virial.xx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.xx, i, 32);
01284         virial.xy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.xy, i, 32);
01285         virial.xz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.xz, i, 32);
01286         virial.yx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.yx, i, 32);
01287         virial.yy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.yy, i, 32);
01288         virial.yz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.yz, i, 32);
01289         virial.zx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.zx, i, 32);
01290         virial.zy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.zy, i, 32);
01291         virial.zz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virial.zz, i, 32);
01292       }   
01293 
01294       if (laneID == 0) {
01295 #ifdef USE_FP_VIRIAL
01296         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 0], llitoulli(virial.xx*double_to_virial));
01297         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 1], llitoulli(virial.xy*double_to_virial));
01298         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 2], llitoulli(virial.xz*double_to_virial));
01299         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 3], llitoulli(virial.yx*double_to_virial));
01300         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 4], llitoulli(virial.yy*double_to_virial));
01301         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 5], llitoulli(virial.yz*double_to_virial));
01302         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 6], llitoulli(virial.zx*double_to_virial));
01303         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 7], llitoulli(virial.zy*double_to_virial));
01304         atomicAdd((unsigned long long int *)&energies_virials[virialIndex + 8], llitoulli(virial.zz*double_to_virial));
01305 #else
01306         atomicAdd(&energies_virials[virialIndex + 0], virial.xx);
01307         atomicAdd(&energies_virials[virialIndex + 1], virial.xy);
01308         atomicAdd(&energies_virials[virialIndex + 2], virial.xz);
01309         atomicAdd(&energies_virials[virialIndex + 3], virial.yx);
01310         atomicAdd(&energies_virials[virialIndex + 4], virial.yy);
01311         atomicAdd(&energies_virials[virialIndex + 5], virial.yz);
01312         atomicAdd(&energies_virials[virialIndex + 6], virial.zx);
01313         atomicAdd(&energies_virials[virialIndex + 7], virial.zy);
01314         atomicAdd(&energies_virials[virialIndex + 8], virial.zz);
01315 #endif
01316       }
01317     }
01318   }
01319 #endif
01320 
01321 }
01322 
01323 template <typename T, bool doEnergy, bool doVirial, bool doElect>
01324 __global__ void modifiedExclusionForcesKernel(
01325   const int start,
01326 
01327   const int numModifiedExclusions,
01328   const CudaExclusion* __restrict__ modifiedExclusions,
01329 
01330   const bool doSlow,
01331   const float one_scale14,                // 1 - scale14
01332   const float cutoff2,
01333   const int vdwCoefTableWidth,
01334   const float2* __restrict__ vdwCoefTable,
01335   cudaTextureObject_t vdwCoefTableTex, 
01336   cudaTextureObject_t modifiedExclusionForceTableTex, cudaTextureObject_t modifiedExclusionEnergyTableTex,
01337 
01338   const float4* __restrict__ xyzq,
01339   const int stride,
01340   const float3 lata, const float3 latb, const float3 latc,
01341   T* __restrict__ forceNbond, T* __restrict__ forceSlow,
01342   double* __restrict__ energies_virials
01343   ) {
01344 
01345   // index
01346   int i = threadIdx.x + (start + blockIdx.x)*blockDim.x;
01347 
01348   double energyVdw, energyNbond, energySlow;
01349   if (doEnergy) {
01350     energyVdw = 0.0;
01351     if (doElect) {
01352       energyNbond = 0.0;
01353       energySlow = 0.0;
01354     }
01355   }
01356 
01357 #ifdef WRITE_FULL_VIRIALS
01358   ComputeBondedCUDAKernel::BondedVirial<double> virialNbond;
01359   ComputeBondedCUDAKernel::BondedVirial<double> virialSlow;
01360   if (doVirial) {
01361     virialNbond.xx = 0.0;
01362     virialNbond.xy = 0.0;
01363     virialNbond.xz = 0.0;
01364     virialNbond.yx = 0.0;
01365     virialNbond.yy = 0.0;
01366     virialNbond.yz = 0.0;
01367     virialNbond.zx = 0.0;
01368     virialNbond.zy = 0.0;
01369     virialNbond.zz = 0.0;
01370     if (doElect) {
01371       virialSlow.xx = 0.0;
01372       virialSlow.xy = 0.0;
01373       virialSlow.xz = 0.0;
01374       virialSlow.yx = 0.0;
01375       virialSlow.yy = 0.0;
01376       virialSlow.yz = 0.0;
01377       virialSlow.zx = 0.0;
01378       virialSlow.zy = 0.0;
01379       virialSlow.zz = 0.0;
01380     }
01381   }
01382 #endif
01383 
01384   if (i < numModifiedExclusions)
01385   {
01386     modifiedExclusionForce<T, doEnergy, doVirial, doElect>
01387     (i, modifiedExclusions, doSlow, one_scale14, vdwCoefTableWidth,
01388 #if __CUDA_ARCH__ >= 350
01389       vdwCoefTable,
01390 #else
01391       vdwCoefTableTex,
01392 #endif
01393       modifiedExclusionForceTableTex, modifiedExclusionEnergyTableTex,
01394       xyzq, stride, lata, latb, latc, cutoff2,
01395       energyVdw, forceNbond, energyNbond,
01396       forceSlow, energySlow, virialNbond, virialSlow);
01397   }
01398 
01399   // Write energies to global memory
01400   if (doEnergy) {
01401     __shared__ double shEnergyVdw[BONDEDFORCESKERNEL_NUM_WARP];
01402     __shared__ double shEnergyNbond[(doElect) ? BONDEDFORCESKERNEL_NUM_WARP : 1];
01403     __shared__ double shEnergySlow[(doElect) ? BONDEDFORCESKERNEL_NUM_WARP : 1];
01404 #pragma unroll
01405     for (int i=16;i >= 1;i/=2) {
01406       energyVdw   += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energyVdw, i, 32);
01407       if (doElect) {
01408         energyNbond += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energyNbond, i, 32);
01409         energySlow  += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energySlow, i, 32);
01410       }
01411     }
01412     int laneID = (threadIdx.x & (WARPSIZE - 1));
01413     int warpID = threadIdx.x / WARPSIZE;
01414     if (laneID == 0) {
01415       shEnergyVdw[warpID]   = energyVdw;
01416       if (doElect) {
01417         shEnergyNbond[warpID] = energyNbond;
01418         shEnergySlow[warpID]  = energySlow;
01419       }
01420     }
01421     BLOCK_SYNC;
01422     if (warpID == 0) {
01423       energyVdw   = (laneID < BONDEDFORCESKERNEL_NUM_WARP) ? shEnergyVdw[laneID] : 0.0;
01424       if (doElect) {
01425         energyNbond = (laneID < BONDEDFORCESKERNEL_NUM_WARP) ? shEnergyNbond[laneID] : 0.0;
01426         energySlow  = (laneID < BONDEDFORCESKERNEL_NUM_WARP) ? shEnergySlow[laneID] : 0.0;
01427       }
01428 #pragma unroll
01429       for (int i=16;i >= 1;i/=2) {
01430         energyVdw   += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energyVdw, i, 32);
01431         if (doElect) {
01432           energyNbond += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energyNbond, i, 32);
01433           energySlow  += WARP_SHUFFLE_XOR(WARP_FULL_MASK, energySlow, i, 32);
01434         }
01435       }
01436       if (laneID == 0) {
01437         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::energyIndex_LJ],         energyVdw);
01438         if (doElect) {
01439           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::energyIndex_ELECT],      energyNbond);
01440           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::energyIndex_ELECT_SLOW], energySlow);
01441         }
01442       }
01443     }
01444   }
01445 
01446   // Write virials to global memory
01447 #ifdef WRITE_FULL_VIRIALS
01448   if (doVirial) {
01449 #pragma unroll
01450     for (int i=16;i >= 1;i/=2) {
01451       virialNbond.xx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.xx, i, 32);
01452       virialNbond.xy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.xy, i, 32);
01453       virialNbond.xz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.xz, i, 32);
01454       virialNbond.yx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.yx, i, 32);
01455       virialNbond.yy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.yy, i, 32);
01456       virialNbond.yz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.yz, i, 32);
01457       virialNbond.zx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.zx, i, 32);
01458       virialNbond.zy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.zy, i, 32);
01459       virialNbond.zz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.zz, i, 32);
01460       if (doElect && doSlow) {
01461         virialSlow.xx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.xx, i, 32);
01462         virialSlow.xy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.xy, i, 32);
01463         virialSlow.xz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.xz, i, 32);
01464         virialSlow.yx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.yx, i, 32);
01465         virialSlow.yy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.yy, i, 32);
01466         virialSlow.yz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.yz, i, 32);
01467         virialSlow.zx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.zx, i, 32);
01468         virialSlow.zy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.zy, i, 32);
01469         virialSlow.zz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.zz, i, 32);
01470       }
01471     }
01472     __shared__ ComputeBondedCUDAKernel::BondedVirial<double> shVirialNbond[BONDEDFORCESKERNEL_NUM_WARP];
01473     __shared__ ComputeBondedCUDAKernel::BondedVirial<double> shVirialSlow[(doElect) ? BONDEDFORCESKERNEL_NUM_WARP : 1];
01474     int laneID = (threadIdx.x & (WARPSIZE - 1));
01475     int warpID = threadIdx.x / WARPSIZE;
01476     if (laneID == 0) {
01477       shVirialNbond[warpID] = virialNbond;
01478       shVirialSlow[warpID] = virialSlow;
01479     }
01480     BLOCK_SYNC;
01481 
01482     virialNbond.xx = 0.0;
01483     virialNbond.xy = 0.0;
01484     virialNbond.xz = 0.0;
01485     virialNbond.yx = 0.0;
01486     virialNbond.yy = 0.0;
01487     virialNbond.yz = 0.0;
01488     virialNbond.zx = 0.0;
01489     virialNbond.zy = 0.0;
01490     virialNbond.zz = 0.0;
01491     if (doElect && doSlow) {
01492       virialSlow.xx = 0.0;
01493       virialSlow.xy = 0.0;
01494       virialSlow.xz = 0.0;
01495       virialSlow.yx = 0.0;
01496       virialSlow.yy = 0.0;
01497       virialSlow.yz = 0.0;
01498       virialSlow.zx = 0.0;
01499       virialSlow.zy = 0.0;
01500       virialSlow.zz = 0.0;
01501     }
01502 
01503     if (warpID == 0) {
01504       if (laneID < BONDEDFORCESKERNEL_NUM_WARP) {
01505         virialNbond = shVirialNbond[laneID];
01506         virialSlow = shVirialSlow[laneID];
01507       }
01508 #pragma unroll
01509       for (int i=16;i >= 1;i/=2) {
01510         virialNbond.xx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.xx, i, 32);
01511         virialNbond.xy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.xy, i, 32);
01512         virialNbond.xz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.xz, i, 32);
01513         virialNbond.yx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.yx, i, 32);
01514         virialNbond.yy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.yy, i, 32);
01515         virialNbond.yz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.yz, i, 32);
01516         virialNbond.zx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.zx, i, 32);
01517         virialNbond.zy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.zy, i, 32);
01518         virialNbond.zz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialNbond.zz, i, 32);
01519         if (doElect && doSlow) {
01520           virialSlow.xx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.xx, i, 32);
01521           virialSlow.xy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.xy, i, 32);
01522           virialSlow.xz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.xz, i, 32);
01523           virialSlow.yx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.yx, i, 32);
01524           virialSlow.yy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.yy, i, 32);
01525           virialSlow.yz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.yz, i, 32);
01526           virialSlow.zx += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.zx, i, 32);
01527           virialSlow.zy += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.zy, i, 32);
01528           virialSlow.zz += WARP_SHUFFLE_XOR(WARP_FULL_MASK, virialSlow.zz, i, 32);
01529         }
01530       }
01531 
01532       if (laneID == 0)
01533       {
01534 #ifdef USE_FP_VIRIAL
01535         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_XX], llitoulli(virialNbond.xx*double_to_virial));
01536         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_XY], llitoulli(virialNbond.xy*double_to_virial));
01537         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_XZ], llitoulli(virialNbond.xz*double_to_virial));
01538         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_YX], llitoulli(virialNbond.yx*double_to_virial));
01539         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_YY], llitoulli(virialNbond.yy*double_to_virial));
01540         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_YZ], llitoulli(virialNbond.yz*double_to_virial));
01541         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_ZX], llitoulli(virialNbond.zx*double_to_virial));
01542         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_ZY], llitoulli(virialNbond.zy*double_to_virial));
01543         atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_ZZ], llitoulli(virialNbond.zz*double_to_virial));
01544         if (doElect && doSlow) {
01545           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_XX], llitoulli(virialSlow.xx*double_to_virial));
01546           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_XY], llitoulli(virialSlow.xy*double_to_virial));
01547           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_XZ], llitoulli(virialSlow.xz*double_to_virial));
01548           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_YX], llitoulli(virialSlow.yx*double_to_virial));
01549           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_YY], llitoulli(virialSlow.yy*double_to_virial));
01550           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_YZ], llitoulli(virialSlow.yz*double_to_virial));
01551           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_ZX], llitoulli(virialSlow.zx*double_to_virial));
01552           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_ZY], llitoulli(virialSlow.zy*double_to_virial));
01553           atomicAdd((unsigned long long int *)&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_ZZ], llitoulli(virialSlow.zz*double_to_virial));
01554         }
01555 #else
01556         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_XX], virialNbond.xx);
01557         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_XY], virialNbond.xy);
01558         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_XZ], virialNbond.xz);
01559         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_YX], virialNbond.yx);
01560         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_YY], virialNbond.yy);
01561         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_YZ], virialNbond.yz);
01562         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_ZX], virialNbond.zx);
01563         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_ZY], virialNbond.zy);
01564         atomicAdd(&energies_virials[ComputeBondedCUDAKernel::nbondVirialIndex_ZZ], virialNbond.zz);
01565         if (doElect && doSlow) {
01566           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_XX], virialSlow.xx);
01567           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_XY], virialSlow.xy);
01568           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_XZ], virialSlow.xz);
01569           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_YX], virialSlow.yx);
01570           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_YY], virialSlow.yy);
01571           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_YZ], virialSlow.yz);
01572           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_ZX], virialSlow.zx);
01573           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_ZY], virialSlow.zy);
01574           atomicAdd(&energies_virials[ComputeBondedCUDAKernel::slowVirialIndex_ZZ], virialSlow.zz);
01575         }
01576 #endif
01577       }
01578     }
01579   }
01580 #endif
01581 
01582 }
01583 
01584 // ##############################################################################################
01585 // ##############################################################################################
01586 // ##############################################################################################
01587 
01588 //
01589 // Class constructor
01590 //
01591 ComputeBondedCUDAKernel::ComputeBondedCUDAKernel(int deviceID, CudaNonbondedTables& cudaNonbondedTables) :
01592 deviceID(deviceID), cudaNonbondedTables(cudaNonbondedTables) {
01593 
01594   cudaCheck(cudaSetDevice(deviceID));
01595 
01596   tupleData = NULL;
01597   tupleDataSize = 0;
01598 
01599   numBonds = 0;
01600   numAngles = 0;
01601   numDihedrals = 0;
01602   numImpropers = 0;
01603   numModifiedExclusions = 0;
01604   numExclusions = 0;
01605   numCrossterms = 0;
01606 
01607   bondValues = NULL;
01608   angleValues = NULL;
01609   dihedralValues = NULL;
01610   improperValues = NULL;
01611   crosstermValues = NULL;
01612 
01613   xyzq = NULL;
01614   xyzqSize = 0;
01615 
01616   forces = NULL;
01617   forcesSize = 0;
01618 
01619   allocate_device<double>(&energies_virials, energies_virials_SIZE);
01620 }
01621 
01622 //
01623 // Class destructor
01624 //
01625 ComputeBondedCUDAKernel::~ComputeBondedCUDAKernel() {
01626   cudaCheck(cudaSetDevice(deviceID));
01627 
01628   deallocate_device<double>(&energies_virials);
01629   // deallocate_device<BondedVirial>(&virial);
01630 
01631   if (tupleData != NULL) deallocate_device<char>(&tupleData);
01632   if (xyzq != NULL) deallocate_device<float4>(&xyzq);
01633   if (forces != NULL) deallocate_device<FORCE_TYPE>(&forces);
01634 
01635   if (bondValues != NULL) deallocate_device<CudaBondValue>(&bondValues);
01636   if (angleValues != NULL) deallocate_device<CudaAngleValue>(&angleValues);
01637   if (dihedralValues != NULL) deallocate_device<CudaDihedralValue>(&dihedralValues);
01638   if (improperValues != NULL) deallocate_device<CudaDihedralValue>(&improperValues);
01639   if (crosstermValues != NULL) deallocate_device<CudaCrosstermValue>(&crosstermValues);
01640 }
01641 
01642 void ComputeBondedCUDAKernel::setupBondValues(int numBondValues, CudaBondValue* h_bondValues) {
01643   allocate_device<CudaBondValue>(&bondValues, numBondValues);
01644   copy_HtoD_sync<CudaBondValue>(h_bondValues, bondValues, numBondValues);
01645 }
01646 
01647 void ComputeBondedCUDAKernel::setupAngleValues(int numAngleValues, CudaAngleValue* h_angleValues) {
01648   allocate_device<CudaAngleValue>(&angleValues, numAngleValues);
01649   copy_HtoD_sync<CudaAngleValue>(h_angleValues, angleValues, numAngleValues);
01650 }
01651 
01652 void ComputeBondedCUDAKernel::setupDihedralValues(int numDihedralValues, CudaDihedralValue* h_dihedralValues) {
01653   allocate_device<CudaDihedralValue>(&dihedralValues, numDihedralValues);
01654   copy_HtoD_sync<CudaDihedralValue>(h_dihedralValues, dihedralValues, numDihedralValues);
01655 }
01656 
01657 void ComputeBondedCUDAKernel::setupImproperValues(int numImproperValues, CudaDihedralValue* h_improperValues) {
01658   allocate_device<CudaDihedralValue>(&improperValues, numImproperValues);
01659   copy_HtoD_sync<CudaDihedralValue>(h_improperValues, improperValues, numImproperValues);
01660 }
01661 
01662 void ComputeBondedCUDAKernel::setupCrosstermValues(int numCrosstermValues, CudaCrosstermValue* h_crosstermValues) {
01663   allocate_device<CudaCrosstermValue>(&crosstermValues, numCrosstermValues);
01664   copy_HtoD_sync<CudaCrosstermValue>(h_crosstermValues, crosstermValues, numCrosstermValues);
01665 }
01666 
01667 //
01668 // Update bonded lists
01669 //
01670 void ComputeBondedCUDAKernel::update(
01671   const int numBondsIn,
01672   const int numAnglesIn,
01673   const int numDihedralsIn,
01674   const int numImpropersIn,
01675   const int numModifiedExclusionsIn,
01676   const int numExclusionsIn,
01677   const int numCrosstermsIn,
01678   const char* h_tupleData,
01679   cudaStream_t stream) {
01680 
01681   numBonds              = numBondsIn;
01682   numAngles             = numAnglesIn;
01683   numDihedrals          = numDihedralsIn;
01684   numImpropers          = numImpropersIn;
01685   numModifiedExclusions = numModifiedExclusionsIn;
01686   numExclusions         = numExclusionsIn;
01687   numCrossterms         = numCrosstermsIn;
01688 
01689   int numBondsWA     = warpAlign(numBonds);
01690   int numAnglesWA    = warpAlign(numAngles);
01691   int numDihedralsWA = warpAlign(numDihedrals);
01692   int numImpropersWA = warpAlign(numImpropers);
01693   int numModifiedExclusionsWA = warpAlign(numModifiedExclusions);
01694   int numExclusionsWA         = warpAlign(numExclusions);
01695   int numCrosstermsWA         = warpAlign(numCrossterms);
01696 
01697   int sizeTot = numBondsWA*sizeof(CudaBond) + numAnglesWA*sizeof(CudaAngle) + 
01698   numDihedralsWA*sizeof(CudaDihedral) + numImpropersWA*sizeof(CudaDihedral) +
01699   numModifiedExclusionsWA*sizeof(CudaExclusion) + numExclusionsWA*sizeof(CudaExclusion) + 
01700   numCrosstermsWA*sizeof(CudaCrossterm);
01701 
01702   reallocate_device<char>(&tupleData, &tupleDataSize, sizeTot, 1.4f);
01703   copy_HtoD<char>(h_tupleData, tupleData, sizeTot, stream);
01704 
01705   // Setup pointers
01706   int pos = 0;
01707   bonds = (CudaBond *)&tupleData[pos];
01708   pos += numBondsWA*sizeof(CudaBond);
01709 
01710   angles = (CudaAngle* )&tupleData[pos];
01711   pos += numAnglesWA*sizeof(CudaAngle);
01712 
01713   dihedrals = (CudaDihedral* )&tupleData[pos];
01714   pos += numDihedralsWA*sizeof(CudaDihedral);
01715 
01716   impropers = (CudaDihedral* )&tupleData[pos];
01717   pos += numImpropersWA*sizeof(CudaDihedral);
01718 
01719   modifiedExclusions = (CudaExclusion* )&tupleData[pos];
01720   pos += numModifiedExclusionsWA*sizeof(CudaExclusion);
01721 
01722   exclusions = (CudaExclusion* )&tupleData[pos];
01723   pos += numExclusionsWA*sizeof(CudaExclusion);
01724 
01725   crossterms = (CudaCrossterm* )&tupleData[pos];
01726   pos += numCrosstermsWA*sizeof(CudaCrossterm);
01727 }
01728 
01729 //
01730 // Return stride for force array
01731 //
01732 int ComputeBondedCUDAKernel::getForceStride(const int atomStorageSize) {
01733 #ifdef USE_STRIDED_FORCE
01734   // Align stride to 256 bytes
01735   return ((atomStorageSize*sizeof(FORCE_TYPE) - 1)/256 + 1)*256/sizeof(FORCE_TYPE);
01736 #else
01737   return 1;
01738 #endif
01739 }
01740 
01741 //
01742 // Return size of single force array
01743 //
01744 int ComputeBondedCUDAKernel::getForceSize(const int atomStorageSize) {
01745 #ifdef USE_STRIDED_FORCE
01746   return (3*getForceStride(atomStorageSize));
01747 #else
01748   return (3*atomStorageSize);
01749 #endif
01750 }
01751 
01752 //
01753 // Return size of the all force arrays
01754 //
01755 int ComputeBondedCUDAKernel::getAllForceSize(const int atomStorageSize, const bool doSlow) {
01756 
01757   int forceSize = getForceSize(atomStorageSize);
01758 
01759   if (numModifiedExclusions > 0 || numExclusions > 0) {
01760     if (doSlow) {
01761       // All three force arrays [normal, nbond, slow]
01762       forceSize *= 3;
01763     } else {
01764       // Two force arrays [normal, nbond]
01765       forceSize *= 2;
01766     }
01767   }
01768 
01769   return forceSize;
01770 }
01771 
01772 //
01773 // Compute bonded forces
01774 //
01775 void ComputeBondedCUDAKernel::bondedForce(
01776   const double scale14, const int atomStorageSize,
01777   const bool doEnergy, const bool doVirial, const bool doSlow,
01778   const float3 lata, const float3 latb, const float3 latc,
01779   const float cutoff2, const float r2_delta, const int r2_delta_expc,
01780   const float4* h_xyzq, FORCE_TYPE* h_forces,
01781   double *h_energies_virials,
01782   cudaStream_t stream) {
01783 
01784   int forceStorageSize = getAllForceSize(atomStorageSize, true);
01785   int forceCopySize = getAllForceSize(atomStorageSize, doSlow);
01786   int forceStride = getForceStride(atomStorageSize);
01787 
01788   int forceSize = getForceSize(atomStorageSize);
01789   int startNbond = forceSize;
01790   int startSlow = 2*forceSize;
01791 
01792   // Re-allocate coordinate and force arrays if neccessary
01793   reallocate_device<float4>(&xyzq, &xyzqSize, atomStorageSize, 1.4f);
01794   reallocate_device<FORCE_TYPE>(&forces, &forcesSize, forceStorageSize, 1.4f);
01795 
01796   // Copy coordinates to device
01797   copy_HtoD<float4>(h_xyzq, xyzq, atomStorageSize, stream);
01798 
01799   // Clear force array
01800   clear_device_array<FORCE_TYPE>(forces, forceCopySize, stream);
01801   if (doEnergy || doVirial) {
01802     clear_device_array<double>(energies_virials, energies_virials_SIZE, stream);
01803   }
01804 
01805   float one_scale14 = (float)(1.0 - scale14);
01806 
01807   // If doSlow = false, these exclusions are not calculated
01808   int numExclusionsDoSlow = doSlow ? numExclusions : 0;
01809 
01810   int nthread = BONDEDFORCESKERNEL_NUM_WARP * WARPSIZE;
01811 
01812   int numBondsTB     = (numBonds + nthread - 1)/nthread;
01813   int numAnglesTB    = (numAngles + nthread - 1)/nthread;
01814   int numDihedralsTB = (numDihedrals + nthread - 1)/nthread;
01815   int numImpropersTB = (numImpropers + nthread - 1)/nthread;
01816   int numExclusionsTB= (numExclusionsDoSlow + nthread - 1)/nthread;
01817   int numCrosstermsTB= (numCrossterms + nthread - 1)/nthread;
01818 
01819   int nblock = numBondsTB + numAnglesTB + numDihedralsTB + numImpropersTB + 
01820   numExclusionsTB + numCrosstermsTB;
01821   int shmem_size = 0;
01822 
01823   // printf("%d %d %d %d %d %d nblock %d\n",
01824   //   numBonds, numAngles, numDihedrals, numImpropers, numModifiedExclusions, numExclusions, nblock);
01825 
01826   int max_nblock = deviceCUDA->getMaxNumBlocks();
01827 
01828   int start = 0;
01829   while (start < nblock)
01830   {
01831     int nleft = nblock - start;
01832     int nblock_use = min(max_nblock, nleft);
01833 
01834 #define CALL(DOENERGY, DOVIRIAL) \
01835   bondedForcesKernel<FORCE_TYPE, DOENERGY, DOVIRIAL> \
01836   <<< nblock_use, nthread, shmem_size, stream >>> \
01837   (start, numBonds, bonds, bondValues, \
01838     numAngles, angles, angleValues, \
01839     numDihedrals, dihedrals, dihedralValues, \
01840     numImpropers, impropers, improperValues, \
01841     numExclusionsDoSlow, exclusions, \
01842     numCrossterms, crossterms, crosstermValues, \
01843     cutoff2, \
01844     r2_delta, r2_delta_expc, \
01845     cudaNonbondedTables.get_r2_table(), cudaNonbondedTables.getExclusionTable(), \
01846     cudaNonbondedTables.get_r2_table_tex(), cudaNonbondedTables.getExclusionTableTex(), \
01847     xyzq, forceStride, \
01848     lata, latb, latc, \
01849     forces, &forces[startSlow], energies_virials);
01850 
01851     if (!doEnergy && !doVirial) CALL(0, 0);
01852     if (!doEnergy && doVirial)  CALL(0, 1);
01853     if (doEnergy && !doVirial)  CALL(1, 0);
01854     if (doEnergy && doVirial)   CALL(1, 1);
01855 
01856 #undef CALL
01857     cudaCheck(cudaGetLastError());
01858 
01859     start += nblock_use;
01860   }
01861 
01862   nthread = BONDEDFORCESKERNEL_NUM_WARP * WARPSIZE;
01863   nblock = (numModifiedExclusions + nthread - 1)/nthread;
01864 
01865   bool doElect = (one_scale14 == 0.0f) ? false : true;
01866 
01867   start = 0;
01868   while (start < nblock)
01869   {
01870     int nleft = nblock - start;
01871     int nblock_use = min(max_nblock, nleft);
01872 
01873 #define CALL(DOENERGY, DOVIRIAL, DOELECT) \
01874   modifiedExclusionForcesKernel<FORCE_TYPE, DOENERGY, DOVIRIAL, DOELECT> \
01875   <<< nblock_use, nthread, shmem_size, stream >>> \
01876   (start, numModifiedExclusions, modifiedExclusions, \
01877     doSlow, one_scale14, cutoff2, \
01878     cudaNonbondedTables.getVdwCoefTableWidth(), cudaNonbondedTables.getExclusionVdwCoefTable(), \
01879     cudaNonbondedTables.getExclusionVdwCoefTableTex(), \
01880     cudaNonbondedTables.getModifiedExclusionForceTableTex(), cudaNonbondedTables.getModifiedExclusionEnergyTableTex(), \
01881     xyzq, forceStride, lata, latb, latc, \
01882     &forces[startNbond], &forces[startSlow], energies_virials);
01883 
01884     if (!doEnergy && !doVirial && !doElect) CALL(0, 0, 0);
01885     if (!doEnergy && doVirial && !doElect)  CALL(0, 1, 0);
01886     if (doEnergy && !doVirial && !doElect)  CALL(1, 0, 0);
01887     if (doEnergy && doVirial && !doElect)   CALL(1, 1, 0);
01888 
01889     if (!doEnergy && !doVirial && doElect)  CALL(0, 0, 1);
01890     if (!doEnergy && doVirial && doElect)   CALL(0, 1, 1);
01891     if (doEnergy && !doVirial && doElect)   CALL(1, 0, 1);
01892     if (doEnergy && doVirial && doElect)    CALL(1, 1, 1);
01893 
01894 #undef CALL
01895     cudaCheck(cudaGetLastError());
01896 
01897     start += nblock_use;
01898   }
01899 
01900   copy_DtoH<FORCE_TYPE>(forces, h_forces, forceCopySize, stream);
01901   if (doEnergy || doVirial) {
01902     copy_DtoH<double>(energies_virials, h_energies_virials, energies_virials_SIZE, stream);
01903   }
01904 
01905 }