Add shared memory optimization (d024b3d3) · Commits · 郑智淋 / lammps

lib/gpu/lal_lj_tip4p_long.cpp

+18 −7

Original line number	Diff line number	Diff line
		@@ -74,11 +74,11 @@ int LJTIP4PLongT::init(const int ntypes,
		// If atom type constants fit in shared memory use fast kernel
		int lj_types=ntypes;
		shared_types=false;
		// int max_shared_types=this->device->max_shared_types();
		// if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
		// lj_types=max_shared_types;
		// shared_types=true;
		// }
		int max_shared_types=this->device->max_shared_types();
		if (lj_types<=max_shared_types && this->_block_size>=max_shared_types) {
		lj_types=max_shared_types;
		shared_types=true;
		}
		_lj_types=lj_types;

		// Allocate a host write buffer for data initialization
		@@ -202,7 +202,6 @@ void LJTIP4PLongT::loop(const bool _eflag, const bool _vflag) {

		GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/
		(BX/this->_threads_per_atom)));
		this->k_pair.set_size(GX,BX);
		if (vflag) {
		this->ansO.resize_ib(ainum*3);
		} else {
		@@ -210,6 +209,17 @@ void LJTIP4PLongT::loop(const bool _eflag, const bool _vflag) {
		}
		this->ansO.zero();
		this->device->gpu->sync();
		if(shared_types) {
		this->k_pair_fast.set_size(GX,BX);
		this->k_pair_fast.run(&this->atom->x, &lj1, &lj3, &_lj_types, &sp_lj,
		&this->nbor->dev_nbor, &this->_nbor_data->begin(),
		&this->ans->force, &this->ans->engv, &eflag, &vflag,
		&ainum, &nbor_pitch, &this->_threads_per_atom,
		&hneight, &m, &TypeO, &TypeH, &alpha,
		&this->atom->q, &cutsq, &_qqrd2e, &_g_ewald,
		&cut_coulsq, &cut_coulsqplus, &this->ansO);
		} else {
		this->k_pair.set_size(GX,BX);
		this->k_pair.run(&this->atom->x, &lj1, &lj3, &_lj_types, &sp_lj,
		&this->nbor->dev_nbor, &this->_nbor_data->begin(),
		&this->ans->force, &this->ans->engv, &eflag, &vflag,
		@@ -217,6 +227,7 @@ void LJTIP4PLongT::loop(const bool _eflag, const bool _vflag) {
		&hneight, &m, &TypeO, &TypeH, &alpha,
		&this->atom->q, &cutsq, &_qqrd2e, &_g_ewald,
		&cut_coulsq, &cut_coulsqplus, &this->ansO);
		}
		GX=static_cast<int>(ceil(static_cast<double>(this->ans->inum())/BX));
		this->k_pair_distrib.set_size(GX,BX);
		this->k_pair_distrib.run(&this->atom->x, &this->ans->force, &this->ans->engv,

lib/gpu/lal_lj_tip4p_long.cu

+325 −1

Original line number	Diff line number	Diff line
		@@ -572,4 +572,328 @@ __kernel void k_lj_tip4p_long(const __global numtyp4 *restrict x_,
		} // if ii
		}

		__kernel void k_lj_tip4p_long_fast() {}
		__kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
		const __global numtyp4 *restrict lj1_in,
		const __global numtyp4 *restrict lj3_in,
		const int lj_types,
		const __global numtyp *restrict sp_lj_in,
		const __global int * dev_nbor,
		const __global int * dev_packed,
		__global acctyp4 *restrict ans,
		__global acctyp *restrict engv,
		const int eflag, const int vflag, const int inum,
		const int nbor_pitch, const int t_per_atom,
		__global int *restrict hneigh,
		__global numtyp4 *restrict m,
		const int typeO, const int typeH,
		const numtyp alpha,
		const __global numtyp *restrict q_,
		const __global numtyp *restrict cutsq,
		const numtyp qqrd2e, const numtyp g_ewald,
		const numtyp cut_coulsq, const numtyp cut_coulsqplus,
		__global acctyp4 *restrict ansO) {
		int tid, ii, offset;
		atom_info(t_per_atom,ii,tid,offset);

		const numtyp eq_zero = 1e-6;

		__local numtyp4 lj1[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
		__local numtyp4 lj3[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
		__local numtyp sp_lj[8];
		if (tid<8)
		sp_lj[tid]=sp_lj_in[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;
		acctyp e_coul = (acctyp)0;
		acctyp4 f, fO;
		f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
		fO.x=(acctyp)0; fO.y=(acctyp)0; fO.z=(acctyp)0;
		acctyp virial[6],vO[6];
		for (int i=0; i<6; i++) {
		virial[i]=(acctyp)0;
		vO[i]=(acctyp)0;
		}

		__syncthreads();
		if (ii<inum) {
		int i, numj, nbor, nbor_end;
		__local int n_stride;
		nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
		n_stride,nbor_end,nbor);

		numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
		numtyp qtmp; fetch(qtmp,i,q_tex);
		int itype = ix.w;
		numtyp4 x1 = ix;

		int non_local_oxy = 0;
		int iH1, iH2, iO;

		if(itype == typeO) {
		iO = i;
		iH1 = hneigh[i*4 ];
		iH2 = hneigh[i*4+1];
		x1 = m[iO];
		} else {
		iO = hneigh[i *4 ];
		iH1 = hneigh[iO*4 ];
		iH2 = hneigh[iO*4+1];
		if (iO >= inum) {
		non_local_oxy = 1;
		}
		}

		for ( ; nbor<nbor_end; nbor+=n_stride) {
		int j=dev_packed[nbor];

		numtyp factor_lj,factor_coul;
		factor_lj = sp_lj[sbmask(j)];
		factor_coul = (numtyp)1.0-sp_lj[sbmask(j)+4];
		j &= NEIGHMASK;

		numtyp4 jx; fetch4(jx,j,pos_tex); //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 rsq = delxdelx+delydely+delz*delz;

		int mtype = itype*lj_types+jtype;
		if (rsq < lj1[mtype].z) { // cut_ljsq
		numtyp r2inv = ucl_recip(rsq);
		numtyp r6inv = r2invr2invr2inv;
		numtyp forcelj = r6inv(lj1[mtype].xr6inv-lj1[mtype].y);
		forcelj = r2invfactor_lj;

		f.x += delx*forcelj;
		f.y += dely*forcelj;
		f.z += delz*forcelj;

		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] += delxdelxforcelj;
		virial[1] += delydelyforcelj;
		virial[2] += delzdelzforcelj;
		virial[3] += delxdelyforcelj;
		virial[4] += delxdelzforcelj;
		virial[5] += delydelzforcelj;
		}
		} // if LJ

		if (rsq < cut_coulsqplus) { //cut_coulsqplus
		int jH1, jH2, jO;
		numtyp qj; fetch(qj,j,q_tex);
		numtyp4 x2 = jx;
		if(itype == typeO \|\| jtype == typeO) {
		if (jtype == typeO) {
		jO = j;
		jH1 = hneigh[j*4 ];
		jH2 = hneigh[j*4+1];
		x2 = m[j];
		}
		delx = x1.x-x2.x;
		dely = x1.y-x2.y;
		delz = x1.z-x2.z;
		rsq = delxdelx+delydely+delz*delz;
		}
		if (rsq < cut_coulsq) {
		numtyp r2inv = ucl_recip(rsq);
		numtyp r = ucl_rsqrt(r2inv);
		numtyp grij = g_ewald * r;
		numtyp expm2 = ucl_exp(-grij*grij);
		numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
		numtyp _erfc = t * (A1+t(A2+t(A3+t(A4+tA5)))) * expm2;

		numtyp prefactor = qj;
		prefactor = qqrd2eqtmp/r;
		numtyp force_coul = r2invprefactor (_erfc + EWALD_Fgrijexpm2 - factor_coul);

		if (itype == typeH) {
		f.x += delx * force_coul;
		f.y += dely * force_coul;
		f.z += delz * force_coul;
		f.w += 0;
		} else {
		fO.x += delx * force_coul;
		fO.y += dely * force_coul;
		fO.z += delz * force_coul;
		fO.w += 0;
		}
		if (eflag>0) {
		e_coul += prefactor*(_erfc-factor_coul);
		}
		if (vflag>0) {
		acctyp4 fd;
		fd.x = delx*force_coul;
		fd.y = dely*force_coul;
		fd.z = delz*force_coul;
		if (itype == typeH) {
		if (jtype == typeH) {
		virial[0] += delx*fd.x;
		virial[1] += dely*fd.y;
		virial[2] += delz*fd.z;
		virial[3] += delx*fd.y;
		virial[4] += delx*fd.z;
		virial[5] += dely*fd.z;
		} else {
		numtyp cO = 1 - alpha, cH = 0.5*alpha;
		numtyp4 vdj;
		numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
		numtyp4 xjH2; fetch4(xjH2,jH2,pos_tex);
		numtyp4 xjO; fetch4(xjO,jO,pos_tex);
		vdj.x = xjO.xcO + xjH1.xcH + xjH2.x*cH;
		vdj.y = xjO.ycO + xjH1.ycH + xjH2.y*cH;
		vdj.z = xjO.zcO + xjH1.zcH + xjH2.z*cH;
		//vdj.w = vdj.w;
		virial[0] += (ix.x - vdj.x)*fd.x;
		virial[1] += (ix.y - vdj.y)*fd.y;
		virial[2] += (ix.z - vdj.z)*fd.z;
		virial[3] += (ix.x - vdj.x)*fd.y;
		virial[4] += (ix.x - vdj.x)*fd.z;
		virial[5] += (ix.y - vdj.y)*fd.z;
		}
		} else {
		numtyp cO = 1 - alpha, cH = 0.5*alpha;
		numtyp4 vdi, vdj;
		numtyp4 xH1; fetch4(xH1,iH1,pos_tex);
		numtyp4 xH2; fetch4(xH2,iH2,pos_tex);
		numtyp4 xO; fetch4(xO,iO,pos_tex);
		vdi.x = xO.xcO + xH1.xcH + xH2.x*cH;
		vdi.y = xO.ycO + xH1.ycH + xH2.y*cH;
		vdi.z = xO.zcO + xH1.zcH + xH2.z*cH;
		//vdi.w = vdi.w;
		if (jtype != typeH) {
		numtyp4 xjH1; fetch4(xjH1,jH1,pos_tex);
		numtyp4 xjH2; fetch4(xjH2,jH2,pos_tex);
		numtyp4 xjO; fetch4(xjO,jO,pos_tex);
		vdj.x = xjO.xcO + xjH1.xcH + xjH2.x*cH;
		vdj.y = xjO.ycO + xjH1.ycH + xjH2.y*cH;
		vdj.z = xjO.zcO + xjH1.zcH + xjH2.z*cH;
		//vdj.w = vdj.w;
		} else vdj = jx;
		vO[0] += 0.5(vdi.x - vdj.x)fd.x;
		vO[1] += 0.5(vdi.y - vdj.y)fd.y;
		vO[2] += 0.5(vdi.z - vdj.z)fd.z;
		vO[3] += 0.5(vdi.x - vdj.x)fd.y;
		vO[4] += 0.5(vdi.x - vdj.x)fd.z;
		vO[5] += 0.5(vdi.y - vdj.y)fd.z;
		}
		}
		}
		if (non_local_oxy == 1) {
		if (iO == j) {
		x2 = ix;
		qj = qtmp;
		}
		numtyp4 x1m = m[iO];
		delx = x1m.x-x2.x;
		dely = x1m.y-x2.y;
		delz = x1m.z-x2.z;
		rsq = delxdelx+delydely+delz*delz;
		if (rsq < cut_coulsq) {
		numtyp r2inv = ucl_recip(rsq);
		numtyp r = ucl_rsqrt(r2inv);
		numtyp grij = g_ewald * r;
		numtyp expm2 = ucl_exp(-grij*grij);
		numtyp t = ucl_recip((numtyp)1.0 + EWALD_P*grij);
		numtyp _erfc = t * (A1+t(A2+t(A3+t(A4+tA5)))) * expm2;

		numtyp prefactor = qj;
		prefactor = qqrd2ex1m.w/r;
		numtyp force_coul = r2invprefactor (_erfc + EWALD_Fgrijexpm2 - factor_coul);

		numtyp cO = 1 - alpha, cH = 0.5*alpha;
		numtyp4 fd;
		fd.x = delx * force_coul * cH;
		fd.y = dely * force_coul * cH;
		fd.z = delz * force_coul * cH;

		f.x += fd.x;
		f.y += fd.y;
		f.z += fd.z;

		if (eflag>0) {
		e_coul += prefactor(_erfc-factor_coul) (acctyp)0.5 * alpha;
		}
		if (vflag>0) {
		numtyp4 xH1; fetch4(xH1,iH1,pos_tex);
		numtyp4 xH2; fetch4(xH2,iH2,pos_tex);
		numtyp4 xO; fetch4(xO,iO,pos_tex);

		virial[0] += ((xO.xcO + xH1.xcH + xH2.xcH) - x2.x) fd.x;
		virial[1] += ((xO.ycO + xH1.ycH + xH2.ycH) - x2.y) fd.y;
		virial[2] += ((xO.zcO + xH1.zcH + xH2.zcH) - x2.z) fd.z;
		virial[3] += ((xO.xcO + xH1.xcH + xH2.xcH) - x2.x) fd.y;
		virial[4] += ((xO.xcO + xH1.xcH + xH2.xcH) - x2.x) fd.z;
		virial[5] += ((xO.ycO + xH1.ycH + xH2.ycH) - x2.y) fd.z;
		}
		}
		}
		} // if cut_coulsqplus
		} // for nbor
		if (t_per_atom>1) {
		#if (ARCH < 300)
		__local acctyp red_acc[6][BLOCK_PAIR];
		red_acc[0][tid]=fO.x;
		red_acc[1][tid]=fO.y;
		red_acc[2][tid]=fO.z;
		red_acc[3][tid]=fO.w;
		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];
		}
		}
		fO.x=red_acc[0][tid];
		fO.y=red_acc[1][tid];
		fO.z=red_acc[2][tid];
		fO.w=red_acc[3][tid];
		if (vflag>0) {
		for (int r=0; r<6; r++) red_acc[r][tid]=vO[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++) vO[r]=red_acc[r][tid];
		}
		#else
		for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
		fO.x += shfl_xor(fO.x, s, t_per_atom);
		fO.y += shfl_xor(fO.y, s, t_per_atom);
		fO.z += shfl_xor(fO.z, s, t_per_atom);
		fO.w += shfl_xor(fO.w, s, t_per_atom);
		}
		if (vflag>0) {
		for (unsigned int s=t_per_atom/2; s>0; s>>=1) {
		for (int r=0; r<6; r++)
		vO[r] += shfl_xor(vO[r], s, t_per_atom);
		}
		}
		#endif
		}
		if(offset == 0) {
		ansO[i] = fO;
		if (vflag>0) {
		ansO[inum + i].x = vO[0];
		ansO[inum + i].y = vO[1];
		ansO[inum + i].z = vO[2];
		ansO[inum*2 + i].x = vO[3];
		ansO[inum*2 + i].y = vO[4];
		ansO[inum*2 + i].z = vO[5];
		}
		}
		store_answers_q(f,energy,e_coul,virial,ii,inum,tid,t_per_atom,offset,eflag,
		vflag,ans,engv);
		} // if ii
		}

Admin message