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