major refactor for restart, data file handling. removal of dead code.

This commit is contained in:
Axel Kohlmeyer
2024-06-19 11:08:31 -04:00
parent 6ada6b7bf2
commit 45508baee5
7 changed files with 90 additions and 118 deletions

View File

@ -47,8 +47,7 @@ int DPDCoulSlaterLongT::init(const int ntypes,
double **host_cutsq, double **host_a0,
double **host_gamma, double **host_sigma,
double **host_cut_dpd, double **host_cut_dpdsq,
double **host_cut_slatersq, double **host_scale,
double *host_special_lj,
double **host_cut_slatersq, double *host_special_lj,
const bool tstat_only,
const int nlocal, const int nall,
const int max_nbors, const int maxspecial,
@ -115,7 +114,7 @@ int DPDCoulSlaterLongT::init(const int ntypes,
cutsq.alloc(lj_types*lj_types,*(this->ucl_device),UCL_READ_ONLY);
this->atom->type_pack4(ntypes,lj_types,cutsq,host_write,host_cutsq,
host_cut_dpdsq, host_scale, host_cut_slatersq);
host_cut_dpdsq,host_cut_slatersq);
double special_sqrt[4];
special_sqrt[0] = sqrt(host_special_lj[0]);

View File

@ -289,8 +289,8 @@ __kernel void k_dpd_coul_slater_long(const __global numtyp4 *restrict x_,
// apply Slater electrostatic force if distance below Slater cutoff
// and the two species have a slater coeff
// cutsq[mtype].w -> Coulombic squared cutoff
if ( cutsq[mtype].w != 0.0 && rsq < cutsq[mtype].w){
// cutsq[mtype].z -> Coulombic squared cutoff
if ( cutsq[mtype].z != 0.0 && rsq < cutsq[mtype].z){
numtyp r2inv=ucl_recip(rsq);
numtyp _erfc;
numtyp grij = g_ewald * r;
@ -426,7 +426,7 @@ __kernel void k_dpd_coul_slater_long_fast(const __global numtyp4 *restrict x_,
int mtype=itype+jx.w;
/// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.z = scale, cutsq.w = cut_slatersq
/// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.z = cut_slatersq
if (rsq<cutsq[mtype].x) {
numtyp r=ucl_sqrt(rsq);
numtyp force_dpd = (numtyp)0.0;
@ -474,8 +474,8 @@ __kernel void k_dpd_coul_slater_long_fast(const __global numtyp4 *restrict x_,
// apply Slater electrostatic force if distance below Slater cutoff
// and the two species have a slater coeff
// cutsq[mtype].w -> Coulombic squared cutoff
if ( cutsq[mtype].w != 0.0 && rsq < cutsq[mtype].w){
// cutsq[mtype].z -> Coulombic squared cutoff
if ( cutsq[mtype].z != 0.0 && rsq < cutsq[mtype].z){
numtyp r2inv=ucl_recip(rsq);
numtyp _erfc;
numtyp grij = g_ewald * r;

View File

@ -37,13 +37,12 @@ class DPDCoulSlaterLong : public BaseDPD<numtyp, acctyp> {
* - -3 if there is an out of memory error
* - -4 if the GPU library was not compiled for GPU
* - -5 Double precision is not supported on card **/
int init(const int ntypes, double **host_cutsq, double **host_a0,
double **host_gamma, double **host_sigma,
double **host_cut_dpd, double **host_cut_dpdsq, double **host_cut_slatersq,
double **host_scale, double *host_special_lj, bool tstat_only, const int nlocal,
const int nall, const int max_nbors, const int maxspecial,
const double cell_size, const double gpu_split, FILE *screen, double *host_special_coul,
const double qqrd2e, const double g_ewald, const double lamda);
int init(const int ntypes, double **host_cutsq, double **host_a0, double **host_gamma,
double **host_sigma, double **host_cut_dpd, double **host_cut_dpdsq,
double **host_cut_slatersq, double *host_special_lj, bool tstat_only, const int nlocal,
const int nall, const int max_nbors, const int maxspecial, const double cell_size,
const double gpu_split, FILE *screen, double *host_special_coul, const double qqrd2e,
const double g_ewald, const double lamda);
/// Clear all host and device data
/** \note This is called at the beginning of the init() routine **/
@ -66,7 +65,7 @@ class DPDCoulSlaterLong : public BaseDPD<numtyp, acctyp> {
/// coeff.x = a0, coeff.y = gamma, coeff.z = sigma, coeff.w = cut_dpd
UCL_D_Vec<numtyp4> coeff;
/// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.z = scale, cutsq.w = cut_slatersq
/// cutsq.x = cutsq, cutsq.y = cut_dpdsq, cutsq.w = cut_slatersq
UCL_D_Vec<numtyp4> cutsq;
/// Special LJ values

View File

@ -27,12 +27,14 @@ static DPDCoulSlaterLong<PRECISION,ACC_PRECISION> DPDCMF;
// ---------------------------------------------------------------------------
// Allocate memory on host and device and copy constants to device
// ---------------------------------------------------------------------------
int dpd_coul_slater_long_gpu_init(const int ntypes, double **host_cutsq, double **host_a0, double **host_gamma,
double **host_sigma, double **host_cut_dpd, double **host_cut_dpdsq, double **host_cut_slatersq,
double **host_scale, double *special_lj, const int inum,
const int nall, const int max_nbors, const int maxspecial, const double cell_size,
int &gpu_mode, FILE *screen, double *host_special_coul,
const double qqrd2e, const double g_ewald, const double lamda) {
int dpd_coul_slater_long_gpu_init(const int ntypes, double **host_cutsq, double **host_a0,
double **host_gamma, double **host_sigma, double **host_cut_dpd,
double **host_cut_dpdsq, double **host_cut_slatersq,
double *special_lj, const int inum, const int nall,
const int max_nbors, const int maxspecial,
const double cell_size, int &gpu_mode, FILE *screen,
double *host_special_coul, const double qqrd2e,
const double g_ewald, const double lamda) {
DPDCMF.clear();
gpu_mode=DPDCMF.device->gpu_mode();
double gpu_split=DPDCMF.device->particle_split();
@ -55,11 +57,10 @@ int dpd_coul_slater_long_gpu_init(const int ntypes, double **host_cutsq, double
int init_ok=0;
if (world_me==0)
init_ok=DPDCMF.init(ntypes, host_cutsq, host_a0, host_gamma, host_sigma,
host_cut_dpd, host_cut_dpdsq, host_cut_slatersq,
host_scale, special_lj, false, inum, nall, max_nbors,
maxspecial, cell_size, gpu_split, screen,
host_special_coul,qqrd2e, g_ewald, lamda);
init_ok=DPDCMF.init(ntypes, host_cutsq, host_a0, host_gamma, host_sigma, host_cut_dpd,
host_cut_dpdsq, host_cut_slatersq, special_lj, false, inum, nall,
max_nbors, maxspecial, cell_size, gpu_split, screen, host_special_coul,
qqrd2e, g_ewald, lamda);
DPDCMF.device->world_barrier();
if (message)
@ -75,11 +76,10 @@ int dpd_coul_slater_long_gpu_init(const int ntypes, double **host_cutsq, double
fflush(screen);
}
if (gpu_rank==i && world_me!=0)
init_ok=DPDCMF.init(ntypes, host_cutsq, host_a0, host_gamma, host_sigma,
host_cut_dpd, host_cut_dpdsq, host_cut_slatersq,
host_scale, special_lj, false, inum, nall, max_nbors,
maxspecial, cell_size, gpu_split, screen,
host_special_coul,qqrd2e, g_ewald, lamda);
init_ok=DPDCMF.init(ntypes, host_cutsq, host_a0, host_gamma, host_sigma, host_cut_dpd,
host_cut_dpdsq, host_cut_slatersq, special_lj, false, inum, nall,
max_nbors, maxspecial, cell_size, gpu_split, screen, host_special_coul,
qqrd2e, g_ewald, lamda);
DPDCMF.device->serialize_init();
if (message)

View File

@ -41,12 +41,13 @@ static constexpr double EPSILON = 1.0e-10;
/* ---------------------------------------------------------------------- */
PairDPDCoulSlaterLong::PairDPDCoulSlaterLong(LAMMPS *lmp) : Pair(lmp)
PairDPDCoulSlaterLong::PairDPDCoulSlaterLong(LAMMPS *lmp) :
Pair(lmp), cut_dpd(nullptr), cut_dpdsq(nullptr), cut_slatersq(nullptr),
a0(nullptr), gamma(nullptr), sigma(nullptr), random(nullptr)
{
writedata = 1;
ewaldflag = pppmflag = 1;
qdist = 0.0;
random = nullptr;
respa_enable = 0;
}
/* ---------------------------------------------------------------------- */
@ -60,14 +61,12 @@ PairDPDCoulSlaterLong::~PairDPDCoulSlaterLong()
memory->destroy(cutsq);
memory->destroy(cut_dpd);
memory->destroy(cut_dpdsq);
memory->destroy(cut_slater);
memory->destroy(cut_slatersq);
memory->destroy(cut);
memory->destroy(a0);
memory->destroy(gamma);
memory->destroy(sigma);
memory->destroy(scale);
}
if (random) delete random;
@ -178,7 +177,7 @@ void PairDPDCoulSlaterLong::compute(int eflag, int vflag)
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*(1-slater_term);
forcecoul *= r2inv;
@ -227,12 +226,9 @@ void PairDPDCoulSlaterLong::allocate()
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(scale,n+1,n+1,"pair:scale");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(cut_dpd,n+1,n+1,"pair:cut_dpd");
memory->create(cut_dpdsq,n+1,n+1,"pair:cut_dpdsq");
memory->create(cut_slater,n+1,n+1,"pair:cut_slater");
memory->create(cut_slatersq,n+1,n+1,"pair:cut_slatersq");
memory->create(a0,n+1,n+1,"pair:a0");
memory->create(gamma,n+1,n+1,"pair:gamma");
@ -256,9 +252,11 @@ void PairDPDCoulSlaterLong::settings(int narg, char **arg)
seed = utils::inumeric(FLERR,arg[2],false,lmp);
lamda = utils::numeric(FLERR,arg[3],false,lmp);
cut_coul = utils::numeric(FLERR,arg[4],false,lmp);
// initialize Marsaglia RNG with processor-unique seed
if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
if (seed <= 0)
error->all(FLERR,"Invalid random seed {} for pair_style dpd/coul/slater/long command", seed);
delete random;
random = new RanMars(lmp,seed + comm->me);
@ -294,7 +292,7 @@ void PairDPDCoulSlaterLong::coeff(int narg, char **arg)
if (narg > 4) {
bool do_slater = utils::logical(FLERR,arg[4],false,lmp);
if (do_slater) cut_two = cut_coul+2.0*qdist;
if (do_slater) cut_two = cut_coul;
}
if (narg > 5) cut_one = utils::numeric(FLERR,arg[5],false,lmp);
@ -305,10 +303,9 @@ void PairDPDCoulSlaterLong::coeff(int narg, char **arg)
a0[i][j] = a0_one;
gamma[i][j] = gamma_one;
cut_dpd[i][j] = cut_one;
cut_slater[i][j] = cut_two;
cut_slatersq[i][j] = cut_two * cut_two;
cut[i][j] = MAX(cut_one, cut_two);
setflag[i][j] = 1;
scale[i][j] = 1.0;
count++;
}
}
@ -360,20 +357,16 @@ double PairDPDCoulSlaterLong::init_one(int i, int j)
sigma[i][j] = sqrt(2.0*force->boltz*temperature*gamma[i][j]);
cut_dpdsq[i][j] = cut_dpd[i][j] * cut_dpd[i][j];
cut_dpdsq[j][i] = cut_dpdsq[i][j];
cut_slatersq[i][j] = cut_slater[i][j] * cut_slater[i][j];
cut_slatersq[j][i] = cut_slatersq[i][j];
a0[j][i] = a0[i][j];
gamma[j][i] = gamma[i][j];
sigma[j][i] = sigma[i][j];
scale[j][i] = scale[i][j];
cut_dpd[j][i] = cut_dpd[i][j];
cut_slater[j][i] = cut_slater[i][j];
cut[j][i] = cut[i][j];
cut_dpdsq[j][i] = cut_dpdsq[i][j];
cut_slatersq[j][i] = cut_slatersq[i][j];
//return cut[i][j];
return MAX(cut_dpd[i][j], cut_slater[i][j]);
return MAX(cut_dpd[i][j], sqrt(cut_slatersq[i][j]));
}
/* ----------------------------------------------------------------------
@ -385,7 +378,7 @@ void PairDPDCoulSlaterLong::write_restart(FILE *fp)
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
@ -393,10 +386,8 @@ void PairDPDCoulSlaterLong::write_restart(FILE *fp)
fwrite(&gamma[i][j],sizeof(double),1,fp);
fwrite(&cut[i][j],sizeof(double),1,fp);
fwrite(&cut_dpd[i][j],sizeof(double),1,fp);
fwrite(&cut_dpdsq[i][j],sizeof(double),1,fp);
fwrite(&cut_slater[i][j],sizeof(double),1,fp);
fwrite(&cut_slatersq[i][j],sizeof(double),1,fp);
fwrite(&scale[i][j],sizeof(double),1,fp);
}
}
}
}
@ -413,7 +404,7 @@ void PairDPDCoulSlaterLong::read_restart(FILE *fp)
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
@ -422,16 +413,15 @@ void PairDPDCoulSlaterLong::read_restart(FILE *fp)
utils::sfread(FLERR,&a0[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&gamma[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &scale[i][j],sizeof(double),1,fp, nullptr, error);
utils::sfread(FLERR,&cut_dpd[i][j],sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_slatersq[i][j],sizeof(double),1,fp,nullptr,error);
}
MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&gamma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_dpd[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_dpdsq[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_slater[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_slatersq[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world);
}
}
}
}
@ -445,14 +435,8 @@ void PairDPDCoulSlaterLong::write_restart_settings(FILE *fp)
fwrite(&temperature,sizeof(double),1,fp);
fwrite(&cut_global,sizeof(double),1,fp);
fwrite(&seed,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&cut_coul,sizeof(double),1,fp);
fwrite(&cut_dpd,sizeof(double),1,fp);
fwrite(&cut_dpdsq,sizeof(double),1,fp);
fwrite(&cut_slater,sizeof(double),1,fp);
fwrite(&cut_slatersq,sizeof(double),1,fp);
fwrite(&lamda,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&cut_coul,sizeof(double),1,fp);
}
/* ----------------------------------------------------------------------
@ -465,22 +449,14 @@ void PairDPDCoulSlaterLong::read_restart_settings(FILE *fp)
utils::sfread(FLERR,&temperature,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&seed,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_coul,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_dpd,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_dpdsq,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_slater,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &cut_slatersq,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &lamda,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR, &offset_flag,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&lamda,sizeof(double),1,fp,nullptr,error);
utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error);
}
MPI_Bcast(&temperature,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&seed,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
// initialize Marsaglia RNG with processor-unique seed
// same seed that pair_style command initially specified
@ -507,7 +483,8 @@ void PairDPDCoulSlaterLong::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g %g\n",i,j,a0[i][j],gamma[i][j],cut[i][j]);
fprintf(fp,"%d %d %g %g %s %g\n",i,j,a0[i][j],gamma[i][j],
(cut_slatersq[i][j] == 0.0) ? "yes" : "no", cut_dpd[i][j]);
}
/* ---------------------------------------------------------------------- */
@ -557,17 +534,11 @@ double PairDPDCoulSlaterLong::single(int i, int j, int itype, int jtype, double
void *PairDPDCoulSlaterLong::extract(const char *str, int &dim)
{
if (strcmp(str,"cut_coul") == 0) {
dim = 0;
return (void *) &cut_coul;
}
if (strcmp(str,"lamda") == 0) {
dim = 0;
return (void *) &lamda;
}
if (strcmp(str,"scale") == 0) {
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
if (strcmp(str,"lamda") == 0) return (void *) &lamda;
dim = 2;
return (void *) scale;
}
if (strcmp(str,"a0") == 0) return (void *) a0;
if (strcmp(str,"gamma") == 0) return (void *) gamma;
return nullptr;
}

View File

@ -47,15 +47,13 @@ class PairDPDCoulSlaterLong : public Pair {
double special_sqrt[4];
int seed;
double **cut;
double **cut_dpd, **cut_dpdsq;
double **cut_slater, **cut_slatersq;
double **cut_dpd, **cut_dpdsq, **cut_slatersq;
double **a0, **gamma;
double **sigma;
class RanMars *random;
double cut_coul, qdist;
double lamda;
double g_ewald;
double **scale;
virtual void allocate();
};

View File

@ -39,26 +39,32 @@ using namespace EwaldConst;
// External functions from cuda library for atom decomposition
int dpd_coul_slater_long_gpu_init(const int ntypes, double **cutsq, double **host_a0, double **host_gamma,
double **host_sigma, double **host_cut_dpd, double **host_cut_dpdsq, double **host_cut_slatersq,
double **host_scale, double *special_lj, const int inum,
const int nall, const int max_nbors, const int maxspecial, const double cell_size,
int &gpu_mode, FILE *screen, double *host_special_coul,
const double qqrd2e, const double g_ewald, const double lamda);
int dpd_coul_slater_long_gpu_init(const int ntypes, double **cutsq, double **host_a0,
double **host_gamma, double **host_sigma, double **host_cut_dpd,
double **host_cut_dpdsq, double **host_cut_slatersq,
double *special_lj, const int inum, const int nall,
const int max_nbors, const int maxspecial,
const double cell_size, int &gpu_mode, FILE *screen,
double *host_special_coul, const double qqrd2e,
const double g_ewald, const double lamda);
void dpd_coul_slater_long_gpu_clear();
int **dpd_coul_slater_long_gpu_compute_n(const int ago, const int inum_full, const int nall, double **host_x,
int *host_type, double *sublo, double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag, const bool eatom,
const bool vatom, int &host_start, int **ilist, int **jnum,
const double cpu_time, bool &success, double **host_v,
const double dtinvsqrt, const int seed, const int timestep, double *boxlo,
int **dpd_coul_slater_long_gpu_compute_n(const int ago, const int inum_full, const int nall,
double **host_x, int *host_type, double *sublo,
double *subhi, tagint *tag, int **nspecial,
tagint **special, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
int **ilist, int **jnum, const double cpu_time,
bool &success, double **host_v, const double dtinvsqrt,
const int seed, const int timestep, double *boxlo,
double *prd);
void dpd_coul_slater_long_gpu_compute(const int ago, const int inum_full, const int nall,
double **host_x, int *host_type, int *ilist, int *numj,
int **firstneigh, const bool eflag, const bool vflag,
const bool eatom, const bool vatom, int &host_start,
const double cpu_time, bool &success, tagint *tag,
double **host_v, const double dtinvsqrt, const int seed,
const int timestep, const int nlocal, double *boxlo,
double *prd);
void dpd_coul_slater_long_gpu_compute(const int ago, const int inum_full, const int nall, double **host_x,
int *host_type, int *ilist, int *numj, int **firstneigh, const bool eflag,
const bool vflag, const bool eatom, const bool vatom, int &host_start,
const double cpu_time, bool &success, tagint *tag, double **host_v,
const double dtinvsqrt, const int seed, const int timestep, const int nlocal,
double *boxlo, double *prd);
void dpd_coul_slater_long_gpu_get_extra_data(double *host_q);
@ -316,8 +322,7 @@ void PairDPDCoulSlaterLongGPU::init_style()
int mnf = 5e-2 * neighbor->oneatom;
int success =
dpd_coul_slater_long_gpu_init(atom->ntypes + 1, cutsq, a0, gamma, sigma,
cut_dpd, cut_dpdsq, cut_slatersq, scale,
force->special_lj, atom->nlocal,
cut_dpd, cut_dpdsq, cut_slatersq, force->special_lj, atom->nlocal,
atom->nlocal + atom->nghost, mnf, maxspecial, cell_size, gpu_mode, screen,
force->special_coul, force->qqrd2e, g_ewald, lamda);
GPU_EXTRA::check_flag(success, error, world);
@ -444,7 +449,7 @@ void PairDPDCoulSlaterLongGPU::cpu_compute(int start, int inum, int eflag, int /
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
prefactor = qqrd2e * qtmp*q[j]/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor*(1-slater_term);
forcecoul *= r2inv;