Homec4science
Diffusion lammps (runner) Browse

lammps/lib/gpu/re_squared_lj.cu
84b5b0bef886runner

Last Change
Raw File

re_squared_lj.cu
View Options

// **************************************************************************
// re_squared_lj.cu
// -------------------
// W. Michael Brown
//
// Device code for RE-Squared - Lennard-Jones potential acceleration
//
// __________________________________________________________________________
// This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
// __________________________________________________________________________
//
// begin : Fri May 06 2011
// email : brownw@ornl.gov
// ***************************************************************************/
#ifndef RE_SQUARED_LJ_CU
#define RE_SQUARED_LJ_CU
#ifdef NV_KERNEL
#include "ellipsoid_extra.h"
#endif
#define SBBITS 30
#define NEIGHMASK 0x3FFFFFFF
__inline int sbmask(int j) { return j >> SBBITS & 3; }
__kernel void kernel_ellipsoid_sphere(__global numtyp4* x_,__global numtyp4 *q,
__global numtyp4* shape, __global numtyp4* well,
__global numtyp *splj, __global numtyp2* sig_eps,
const int ntypes, __global int *dev_nbor, const int stride,
__global acctyp4 *ans, const int astride,
__global acctyp *engv, __global int *err_flag,
const int eflag, const int vflag, const int inum,
const int t_per_atom) {
int tid=THREAD_ID_X;
int ii=mul24((int)BLOCK_ID_X,(int)(BLOCK_SIZE_X)/t_per_atom);
ii+=tid/t_per_atom;
int offset=tid%t_per_atom;
__local numtyp sp_lj[4];
sp_lj[0]=splj[0];
sp_lj[1]=splj[1];
sp_lj[2]=splj[2];
sp_lj[3]=splj[3];
__local numtyp b_alpha, cr60, solv_f_a, solv_f_r;
b_alpha=(numtyp)45.0/(numtyp)56.0;
cr60=pow((numtyp)60.0,(numtyp)1.0/(numtyp)3.0);
solv_f_a = (numtyp)3.0/((numtyp)16.0*atan((numtyp)1.0)*-(numtyp)36.0);
solv_f_r = (numtyp)3.0/((numtyp)16.0*atan((numtyp)1.0)*(numtyp)2025.0);
acctyp energy=(acctyp)0;
acctyp4 f;
f.x=(acctyp)0;
f.y=(acctyp)0;
f.z=(acctyp)0;
acctyp4 tor;
tor.x=(acctyp)0;
tor.y=(acctyp)0;
tor.z=(acctyp)0;
acctyp virial[6];
for (int i=0; i<6; i++)
virial[i]=(acctyp)0;
if (ii<inum) {
__global int *nbor=dev_nbor+ii;
int i=*nbor;
nbor+=stride;
int numj=*nbor;
nbor+=stride;
__global int *nbor_end=nbor+mul24(stride,numj);
nbor+=mul24(offset,stride);
int n_stride=mul24(t_per_atom,stride);
numtyp4 ix=x_[i];
int itype=ix.w;
numtyp a[9]; // Rotation matrix (lab->body)
numtyp aTe[9]; // A'*E
numtyp lA_0[9], lA_1[9], lA_2[9]; // -A*rotation generator (x,y, or z)
numtyp4 ishape;
ishape=shape[itype];
numtyp ilshape=ishape.x*ishape.y*ishape.z;
{
gpu_quat_to_mat_trans(q,i,a);
gpu_transpose_times_diag3(a,well[itype],aTe);
gpu_rotation_generator_x(a,lA_0);
gpu_rotation_generator_y(a,lA_1);
gpu_rotation_generator_z(a,lA_2);
}
numtyp factor_lj;
for ( ; nbor<nbor_end; nbor+=n_stride) {
int j=*nbor;
factor_lj = sp_lj[sbmask(j)];
j &= NEIGHMASK;
numtyp4 jx=x_[j];
int jtype=jx.w;
// Compute r12
numtyp r[3], rhat[3];
numtyp rnorm;
r[0] = jx.x-ix.x;
r[1] = jx.y-ix.y;
r[2] = jx.z-ix.z;
rnorm = gpu_dot3(r,r);
rnorm = rsqrt(rnorm);
rhat[0] = r[0]*rnorm;
rhat[1] = r[1]*rnorm;
rhat[2] = r[2]*rnorm;
numtyp sigma, epsilon;
int mtype=mul24(ntypes,itype)+jtype;
sigma = sig_eps[mtype].x;
epsilon = sig_eps[mtype].y*factor_lj;
numtyp aTs[9];
numtyp4 scorrect;
numtyp half_sigma=sigma*(numtyp)0.5;
scorrect.x = ishape.x+half_sigma;
scorrect.y = ishape.y+half_sigma;
scorrect.z = ishape.z+half_sigma;
scorrect.x = scorrect.x * scorrect.x * (numtyp)0.5;
scorrect.y = scorrect.y * scorrect.y * (numtyp)0.5;
scorrect.z = scorrect.z * scorrect.z * (numtyp)0.5;
gpu_transpose_times_diag3(a,scorrect,aTs);
// energy
numtyp gamma[9], s[3];
gpu_times3(aTs,a,gamma);
gpu_mldivide3(gamma,rhat,s,err_flag);
numtyp sigma12 = rsqrt((numtyp)0.5*gpu_dot3(s,rhat));
numtyp temp[9], w[3];
gpu_times3(aTe,a,temp);
temp[0] += (numtyp)1.0;
temp[4] += (numtyp)1.0;
temp[8] += (numtyp)1.0;
gpu_mldivide3(temp,rhat,w,err_flag);
numtyp h12 = (numtyp)1.0/rnorm-sigma12;
numtyp chi = (numtyp)2.0*gpu_dot3(rhat,w);
numtyp sigh = sigma/h12;
numtyp tprod = chi*sigh;
numtyp Ua, Ur;
numtyp h12p3 = h12*h12*h12;
numtyp sigmap3 = sigma*sigma*sigma;
numtyp stemp = h12*(numtyp)0.5;
Ua = (ishape.x+stemp)*(ishape.y+stemp)*(ishape.z+stemp)*h12p3/(numtyp)8.0;
Ua = ((numtyp)1.0+(numtyp)3.0*tprod)*ilshape/Ua;
Ua = epsilon*Ua*sigmap3*solv_f_a;
stemp = h12/cr60;
Ur = (ishape.x+stemp)*(ishape.y+stemp)*(ishape.z+stemp)*h12p3/
(numtyp)60.0;
Ur = ((numtyp)1.0+b_alpha*tprod)*ilshape/Ur;
numtyp sigh6=sigh*sigh*sigh;
sigh6*=sigh6;
Ur = epsilon*Ur*sigmap3*sigh6*solv_f_r;
energy+=Ua+Ur;
// force
numtyp fourw[3], spr[3];
numtyp sec = sigma*chi;
numtyp sigma12p3 = sigma12*sigma12*sigma12;
fourw[0] = (numtyp)4.0*w[0];
fourw[1] = (numtyp)4.0*w[1];
fourw[2] = (numtyp)4.0*w[2];
spr[0] = (numtyp)0.5*sigma12p3*s[0];
spr[1] = (numtyp)0.5*sigma12p3*s[1];
spr[2] = (numtyp)0.5*sigma12p3*s[2];
stemp = (numtyp)1.0/(ishape.x*(numtyp)2.0+h12)+
(numtyp)1.0/(ishape.y*(numtyp)2.0+h12)+
(numtyp)1.0/(ishape.z*(numtyp)2.0+h12)+
(numtyp)3.0/h12;
numtyp hsec = (numtyp)1.0/(h12+(numtyp)3.0*sec);
numtyp dspu = (numtyp)1.0/h12-hsec+stemp;
numtyp pbsu = (numtyp)3.0*sigma*hsec;
stemp = (numtyp)1.0/(ishape.x*cr60+h12)+
(numtyp)1.0/(ishape.y*cr60+h12)+
(numtyp)1.0/(ishape.z*cr60+h12)+
(numtyp)3.0/h12;
hsec = (numtyp)1.0/(h12+b_alpha*sec);
numtyp dspr = (numtyp)7.0/h12-hsec+stemp;
numtyp pbsr = b_alpha*sigma*hsec;
#pragma unroll
for (int i=0; i<3; i++) {
numtyp u[3];
u[0] = -rhat[i]*rhat[0];
u[1] = -rhat[i]*rhat[1];
u[2] = -rhat[i]*rhat[2];
u[i] += (numtyp)1.0;
u[0] *= rnorm;
u[1] *= rnorm;
u[2] *= rnorm;
numtyp dchi = gpu_dot3(u,fourw);
numtyp dh12 = rhat[i]+gpu_dot3(u,spr);
numtyp dUa = pbsu*dchi-dh12*dspu;
numtyp dUr = pbsr*dchi-dh12*dspr;
numtyp force=dUr*Ur+dUa*Ua;
if (i==0) {
f.x+=force;
if (vflag>0)
virial[0]+=-r[0]*force;
} else if (i==1) {
f.y+=force;
if (vflag>0) {
virial[1]+=-r[1]*force;
virial[3]+=-r[0]*force;
}
} else {
f.z+=force;
if (vflag>0) {
virial[2]+=-r[2]*force;
virial[4]+=-r[0]*force;
virial[5]+=-r[1]*force;
}
}
}
// torque on i
numtyp fwae[3];
gpu_row_times3(fourw,aTe,fwae);
{
numtyp tempv[3], p[3], lAtwo[9];
gpu_times_column3(lA_0,rhat,p);
gpu_times_column3(lA_0,w,tempv);
numtyp dchi = -gpu_dot3(fwae,tempv);
gpu_times3(aTs,lA_0,lAtwo);
gpu_times_column3(lAtwo,spr,tempv);
numtyp dh12 = -gpu_dot3(s,tempv);
numtyp dUa = pbsu*dchi-dh12*dspu;
numtyp dUr = pbsr*dchi-dh12*dspr;
tor.x -= (dUa*Ua+dUr*Ur);
}
{
numtyp tempv[3], p[3], lAtwo[9];
gpu_times_column3(lA_1,rhat,p);
gpu_times_column3(lA_1,w,tempv);
numtyp dchi = -gpu_dot3(fwae,tempv);
gpu_times3(aTs,lA_1,lAtwo);
gpu_times_column3(lAtwo,spr,tempv);
numtyp dh12 = -gpu_dot3(s,tempv);
numtyp dUa = pbsu*dchi-dh12*dspu;
numtyp dUr = pbsr*dchi-dh12*dspr;
tor.y -= (dUa*Ua+dUr*Ur);
}
{
numtyp tempv[3], p[3], lAtwo[9];
gpu_times_column3(lA_2,rhat,p);
gpu_times_column3(lA_2,w,tempv);
numtyp dchi = -gpu_dot3(fwae,tempv);
gpu_times3(aTs,lA_2,lAtwo);
gpu_times_column3(lAtwo,spr,tempv);
numtyp dh12 = -gpu_dot3(s,tempv);
numtyp dUa = pbsu*dchi-dh12*dspu;
numtyp dUr = pbsr*dchi-dh12*dspr;
tor.z -= (dUa*Ua+dUr*Ur);
}
} // for nbor
} // if ii
// Reduce answers
if (t_per_atom>1) {
__local acctyp red_acc[7][BLOCK_PAIR];
red_acc[0][tid]=f.x;
red_acc[1][tid]=f.y;
red_acc[2][tid]=f.z;
red_acc[3][tid]=tor.x;
red_acc[4][tid]=tor.y;
red_acc[5][tid]=tor.z;
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<6; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
f.x=red_acc[0][tid];
f.y=red_acc[1][tid];
f.z=red_acc[2][tid];
tor.x=red_acc[3][tid];
tor.y=red_acc[4][tid];
tor.z=red_acc[5][tid];
if (eflag>0 || vflag>0) {
for (int r=0; r<6; r++)
red_acc[r][tid]=virial[r];
red_acc[6][tid]=energy;
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<7; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
for (int r=0; r<6; r++)
virial[r]=red_acc[r][tid];
energy=red_acc[6][tid];
}
}
// Store answers
if (ii<inum && offset==0) {
__global acctyp *ap1=engv+ii;
if (eflag>0) {
*ap1+=energy;
ap1+=astride;
}
if (vflag>0) {
for (int i=0; i<6; i++) {
*ap1+=virial[i];
ap1+=astride;
}
}
acctyp4 old=ans[ii];
old.x+=f.x;
old.y+=f.y;
old.z+=f.z;
ans[ii]=old;
old=ans[ii+astride];
old.x+=tor.x;
old.y+=tor.y;
old.z+=tor.z;
ans[ii+astride]=old;
} // if ii
}
__kernel void kernel_sphere_ellipsoid(__global numtyp4 *x_,__global numtyp4 *q,
__global numtyp4* shape,__global numtyp4* well,
__global numtyp *splj, __global numtyp2* sig_eps,
const int ntypes, __global int *dev_nbor,
const int stride, __global acctyp4 *ans,
__global acctyp *engv, __global int *err_flag,
const int eflag, const int vflag,const int start,
const int inum, const int t_per_atom) {
int tid=THREAD_ID_X;
int ii=mul24((int)BLOCK_ID_X,(int)(BLOCK_SIZE_X)/t_per_atom);
ii+=tid/t_per_atom+start;
int offset=tid%t_per_atom;
__local numtyp sp_lj[4];
sp_lj[0]=splj[0];
sp_lj[1]=splj[1];
sp_lj[2]=splj[2];
sp_lj[3]=splj[3];
__local numtyp b_alpha, cr60, solv_f_a, solv_f_r;
b_alpha=(numtyp)45.0/(numtyp)56.0;
cr60=pow((numtyp)60.0,(numtyp)1.0/(numtyp)3.0);
solv_f_a = (numtyp)3.0/((numtyp)16.0*atan((numtyp)1.0)*-(numtyp)36.0);
solv_f_r = (numtyp)3.0/((numtyp)16.0*atan((numtyp)1.0)*(numtyp)2025.0);
acctyp energy=(acctyp)0;
acctyp4 f;
f.x=(acctyp)0;
f.y=(acctyp)0;
f.z=(acctyp)0;
acctyp virial[6];
for (int i=0; i<6; i++)
virial[i]=(acctyp)0;
if (ii<inum) {
__global int *nbor=dev_nbor+ii;
int j=*nbor;
nbor+=stride;
int numj=*nbor;
nbor+=stride;
__global int *nbor_end=nbor+mul24(stride,numj);
nbor+=mul24(offset,stride);
int n_stride=mul24(t_per_atom,stride);
numtyp4 jx=x_[j];
int jtype=jx.w;
numtyp factor_lj;
for ( ; nbor<nbor_end; nbor+=n_stride) {
int i=*nbor;
factor_lj = sp_lj[sbmask(i)];
i &= NEIGHMASK;
numtyp4 ix=x_[i];
int itype=ix.w;
numtyp a[9]; // Rotation matrix (lab->body)
numtyp aTe[9]; // A'*E
numtyp4 ishape;
ishape=shape[itype];
gpu_quat_to_mat_trans(q,i,a);
gpu_transpose_times_diag3(a,well[itype],aTe);
// Compute r12
numtyp r[3], rhat[3];
numtyp rnorm;
r[0] = ix.x-jx.x;
r[1] = ix.y-jx.y;
r[2] = ix.z-jx.z;
rnorm = gpu_dot3(r,r);
rnorm = rsqrt(rnorm);
rhat[0] = r[0]*rnorm;
rhat[1] = r[1]*rnorm;
rhat[2] = r[2]*rnorm;
numtyp sigma, epsilon;
int mtype=mul24(ntypes,itype)+jtype;
sigma = sig_eps[mtype].x;
epsilon = sig_eps[mtype].y*factor_lj;
numtyp aTs[9];
numtyp4 scorrect;
numtyp half_sigma=sigma * (numtyp)0.5;
scorrect.x = ishape.x+half_sigma;
scorrect.y = ishape.y+half_sigma;
scorrect.z = ishape.z+half_sigma;
scorrect.x = scorrect.x * scorrect.x * (numtyp)0.5;
scorrect.y = scorrect.y * scorrect.y * (numtyp)0.5;
scorrect.z = scorrect.z * scorrect.z * (numtyp)0.5;
gpu_transpose_times_diag3(a,scorrect,aTs);
// energy
numtyp gamma[9], s[3];
gpu_times3(aTs,a,gamma);
gpu_mldivide3(gamma,rhat,s,err_flag);
numtyp sigma12 = rsqrt((numtyp)0.5*gpu_dot3(s,rhat));
numtyp temp[9], w[3];
gpu_times3(aTe,a,temp);
temp[0] += (numtyp)1.0;
temp[4] += (numtyp)1.0;
temp[8] += (numtyp)1.0;
gpu_mldivide3(temp,rhat,w,err_flag);
numtyp h12 = (numtyp)1.0/rnorm-sigma12;
numtyp chi = (numtyp)2.0*gpu_dot3(rhat,w);
numtyp sigh = sigma/h12;
numtyp tprod = chi*sigh;
numtyp Ua, Ur;
numtyp h12p3 = h12*h12*h12;
numtyp sigmap3 = sigma*sigma*sigma;
numtyp stemp = h12/(numtyp)2.0;
Ua = (ishape.x+stemp)*(ishape.y+stemp)*(ishape.z+stemp)*h12p3/(numtyp)8.0;
numtyp ilshape=ishape.x*ishape.y*ishape.z;
Ua = ((numtyp)1.0+(numtyp)3.0*tprod)*ilshape/Ua;
Ua = epsilon*Ua*sigmap3*solv_f_a;
stemp = h12/cr60;
Ur = (ishape.x+stemp)*(ishape.y+stemp)*(ishape.z+stemp)*h12p3/
(numtyp)60.0;
Ur = ((numtyp)1.0+b_alpha*tprod)*ilshape/Ur;
numtyp sigh6=sigh*sigh*sigh;
sigh6*=sigh6;
Ur = epsilon*Ur*sigmap3*sigh6*solv_f_r;
energy+=Ua+Ur;
// force
numtyp fourw[3], spr[3];
numtyp sec = sigma*chi;
numtyp sigma12p3 = sigma12*sigma12*sigma12;
fourw[0] = (numtyp)4.0*w[0];
fourw[1] = (numtyp)4.0*w[1];
fourw[2] = (numtyp)4.0*w[2];
spr[0] = (numtyp)0.5*sigma12p3*s[0];
spr[1] = (numtyp)0.5*sigma12p3*s[1];
spr[2] = (numtyp)0.5*sigma12p3*s[2];
stemp = (numtyp)1.0/(ishape.x*(numtyp)2.0+h12)+
(numtyp)1.0/(ishape.y*(numtyp)2.0+h12)+
(numtyp)1.0/(ishape.z*(numtyp)2.0+h12)+
(numtyp)3.0/h12;
numtyp hsec = (numtyp)1.0/(h12+(numtyp)3.0*sec);
numtyp dspu = (numtyp)1.0/h12-hsec+stemp;
numtyp pbsu = (numtyp)3.0*sigma*hsec;
stemp = (numtyp)1.0/(ishape.x*cr60+h12)+
(numtyp)1.0/(ishape.y*cr60+h12)+
(numtyp)1.0/(ishape.z*cr60+h12)+
(numtyp)3.0/h12;
hsec = (numtyp)1.0/(h12+b_alpha*sec);
numtyp dspr = (numtyp)7.0/h12-hsec+stemp;
numtyp pbsr = b_alpha*sigma*hsec;
#pragma unroll
for (int i=0; i<3; i++) {
numtyp u[3];
u[0] = -rhat[i]*rhat[0];
u[1] = -rhat[i]*rhat[1];
u[2] = -rhat[i]*rhat[2];
u[i] += (numtyp)1.0;
u[0] *= rnorm;
u[1] *= rnorm;
u[2] *= rnorm;
numtyp dchi = gpu_dot3(u,fourw);
numtyp dh12 = rhat[i]+gpu_dot3(u,spr);
numtyp dUa = pbsu*dchi-dh12*dspu;
numtyp dUr = pbsr*dchi-dh12*dspr;
numtyp force=dUr*Ur+dUa*Ua;
if (i==0) {
f.x+=force;
if (vflag>0)
virial[0]+=-r[0]*force;
} else if (i==1) {
f.y+=force;
if (vflag>0) {
virial[1]+=-r[1]*force;
virial[3]+=-r[0]*force;
}
} else {
f.z+=force;
if (vflag>0) {
virial[2]+=-r[2]*force;
virial[4]+=-r[0]*force;
virial[5]+=-r[1]*force;
}
}
}
} // for nbor
} // if ii
// Reduce answers
if (t_per_atom>1) {
__local acctyp red_acc[7][BLOCK_PAIR];
red_acc[0][tid]=f.x;
red_acc[1][tid]=f.y;
red_acc[2][tid]=f.z;
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<3; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
f.x=red_acc[0][tid];
f.y=red_acc[1][tid];
f.z=red_acc[2][tid];
if (eflag>0 || vflag>0) {
for (int r=0; r<6; r++)
red_acc[r][tid]=virial[r];
red_acc[6][tid]=energy;
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<7; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
for (int r=0; r<6; r++)
virial[r]=red_acc[r][tid];
energy=red_acc[6][tid];
}
}
// Store answers
if (ii<inum && offset==0) {
__global acctyp *ap1=engv+ii;
if (eflag>0) {
*ap1=energy;
ap1+=inum;
}
if (vflag>0) {
for (int i=0; i<6; i++) {
*ap1=virial[i];
ap1+=inum;
}
}
ans[ii]=f;
} // if ii
}
__kernel void kernel_lj(__global numtyp4 *x_, __global numtyp4 *lj1,
__global numtyp4* lj3, const int lj_types,
__global numtyp *gum,
const int stride, __global int *dev_ij,
__global acctyp4 *ans, __global acctyp *engv,
__global int *err_flag, const int eflag,
const int vflag, const int start, const int inum,
const int t_per_atom) {
int tid=THREAD_ID_X;
int ii=mul24((int)BLOCK_ID_X,(int)(BLOCK_SIZE_X)/t_per_atom);
ii+=tid/t_per_atom+start;
int offset=tid%t_per_atom;
__local numtyp sp_lj[4];
sp_lj[0]=gum[0];
sp_lj[1]=gum[1];
sp_lj[2]=gum[2];
sp_lj[3]=gum[3];
acctyp energy=(acctyp)0;
acctyp4 f;
f.x=(acctyp)0;
f.y=(acctyp)0;
f.z=(acctyp)0;
acctyp virial[6];
for (int i=0; i<6; i++)
virial[i]=(acctyp)0;
if (ii<inum) {
__global int *nbor=dev_ij+ii;
int i=*nbor;
nbor+=stride;
int numj=*nbor;
nbor+=stride;
__global int *list_end=nbor+mul24(stride,numj);
nbor+=mul24(offset,stride);
int n_stride=mul24(t_per_atom,stride);
numtyp4 ix=x_[i];
int itype=ix.w;
numtyp factor_lj;
for ( ; nbor<list_end; nbor+=n_stride) {
int j=*nbor;
factor_lj = sp_lj[sbmask(j)];
j &= NEIGHMASK;
numtyp4 jx=x_[j];
int jtype=jx.w;
// Compute r12
numtyp delx = ix.x-jx.x;
numtyp dely = ix.y-jx.y;
numtyp delz = ix.z-jx.z;
numtyp r2inv = delx*delx+dely*dely+delz*delz;
int ii=itype*lj_types+jtype;
if (r2inv<lj1[ii].z && lj1[ii].w==SPHERE_SPHERE) {
r2inv=(numtyp)1.0/r2inv;
numtyp r6inv = r2inv*r2inv*r2inv;
numtyp force = r2inv*r6inv*(lj1[ii].x*r6inv-lj1[ii].y);
force*=factor_lj;
f.x+=delx*force;
f.y+=dely*force;
f.z+=delz*force;
if (eflag>0) {
numtyp e=r6inv*(lj3[ii].x*r6inv-lj3[ii].y);
energy+=factor_lj*(e-lj3[ii].z);
}
if (vflag>0) {
virial[0] += delx*delx*force;
virial[1] += dely*dely*force;
virial[2] += delz*delz*force;
virial[3] += delx*dely*force;
virial[4] += delx*delz*force;
virial[5] += dely*delz*force;
}
}
} // for nbor
} // if ii
// Reduce answers
if (t_per_atom>1) {
__local acctyp red_acc[6][BLOCK_PAIR];
red_acc[0][tid]=f.x;
red_acc[1][tid]=f.y;
red_acc[2][tid]=f.z;
red_acc[3][tid]=energy;
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<4; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
f.x=red_acc[0][tid];
f.y=red_acc[1][tid];
f.z=red_acc[2][tid];
energy=red_acc[3][tid];
if (vflag>0) {
for (int r=0; r<6; r++)
red_acc[r][tid]=virial[r];
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<6; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
for (int r=0; r<6; r++)
virial[r]=red_acc[r][tid];
}
}
// Store answers
if (ii<inum && offset==0) {
__global acctyp *ap1=engv+ii;
if (eflag>0) {
*ap1+=energy;
ap1+=inum;
}
if (vflag>0) {
for (int i=0; i<6; i++) {
*ap1+=virial[i];
ap1+=inum;
}
}
acctyp4 old=ans[ii];
old.x+=f.x;
old.y+=f.y;
old.z+=f.z;
ans[ii]=old;
} // if ii
}
__kernel void kernel_lj_fast(__global numtyp4 *x_, __global numtyp4 *lj1_in,
__global numtyp4* lj3_in, __global numtyp *gum,
const int stride, __global int *dev_ij,
__global acctyp4 *ans, __global acctyp *engv,
__global int *err_flag, const int eflag,
const int vflag, const int start, const int inum,
const int t_per_atom) {
int tid=THREAD_ID_X;
int ii=mul24((int)BLOCK_ID_X,(int)(BLOCK_SIZE_X)/t_per_atom);
ii+=tid/t_per_atom+start;
int offset=tid%t_per_atom;
__local numtyp sp_lj[4];
__local numtyp4 lj1[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
__local numtyp4 lj3[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
if (tid<4)
sp_lj[tid]=gum[tid];
if (tid<MAX_SHARED_TYPES*MAX_SHARED_TYPES) {
lj1[tid]=lj1_in[tid];
if (eflag>0)
lj3[tid]=lj3_in[tid];
}
acctyp energy=(acctyp)0;
acctyp4 f;
f.x=(acctyp)0;
f.y=(acctyp)0;
f.z=(acctyp)0;
acctyp virial[6];
for (int i=0; i<6; i++)
virial[i]=(acctyp)0;
__syncthreads();
if (ii<inum) {
__global int *nbor=dev_ij+ii;
int i=*nbor;
nbor+=stride;
int numj=*nbor;
nbor+=stride;
__global int *list_end=nbor+mul24(stride,numj);
nbor+=mul24(offset,stride);
int n_stride=mul24(t_per_atom,stride);
numtyp4 ix=x_[i];
int iw=ix.w;
int itype=mul24((int)MAX_SHARED_TYPES,iw);
numtyp factor_lj;
for ( ; nbor<list_end; nbor+=n_stride) {
int j=*nbor;
factor_lj = sp_lj[sbmask(j)];
j &= NEIGHMASK;
numtyp4 jx=x_[j];
int mtype=itype+jx.w;
// Compute r12
numtyp delx = ix.x-jx.x;
numtyp dely = ix.y-jx.y;
numtyp delz = ix.z-jx.z;
numtyp r2inv = delx*delx+dely*dely+delz*delz;
if (r2inv<lj1[mtype].z && lj1[mtype].w==SPHERE_SPHERE) {
r2inv=(numtyp)1.0/r2inv;
numtyp r6inv = r2inv*r2inv*r2inv;
numtyp force = factor_lj*r2inv*r6inv*(lj1[mtype].x*r6inv-lj1[mtype].y);
f.x+=delx*force;
f.y+=dely*force;
f.z+=delz*force;
if (eflag>0) {
numtyp e=r6inv*(lj3[mtype].x*r6inv-lj3[mtype].y);
energy+=factor_lj*(e-lj3[mtype].z);
}
if (vflag>0) {
virial[0] += delx*delx*force;
virial[1] += dely*dely*force;
virial[2] += delz*delz*force;
virial[3] += delx*dely*force;
virial[4] += delx*delz*force;
virial[5] += dely*delz*force;
}
}
} // for nbor
} // if ii
// Reduce answers
if (t_per_atom>1) {
__local acctyp red_acc[6][BLOCK_PAIR];
red_acc[0][tid]=f.x;
red_acc[1][tid]=f.y;
red_acc[2][tid]=f.z;
red_acc[3][tid]=energy;
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<4; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
f.x=red_acc[0][tid];
f.y=red_acc[1][tid];
f.z=red_acc[2][tid];
energy=red_acc[3][tid];
if (vflag>0) {
for (int r=0; r<6; r++)
red_acc[r][tid]=virial[r];
for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
if (offset < s) {
for (int r=0; r<6; r++)
red_acc[r][tid] += red_acc[r][tid+s];
}
}
for (int r=0; r<6; r++)
virial[r]=red_acc[r][tid];
}
}
// Store answers
if (ii<inum && offset==0) {
__global acctyp *ap1=engv+ii;
if (eflag>0) {
*ap1+=energy;
ap1+=inum;
}
if (vflag>0) {
for (int i=0; i<6; i++) {
*ap1+=virial[i];
ap1+=inum;
}
}
acctyp4 old=ans[ii];
old.x+=f.x;
old.y+=f.y;
old.z+=f.z;
ans[ii]=old;
} // if ii
}
#endif
c4science · Help