revised implementation of inclusion of embedding energy in PairEAM::single() (e0646b73) · Commits · 郑智淋 / lammps

src/MANYBODY/pair_eam.cpp

+8 −7

Original line number	Diff line number	Diff line
		@@ -43,7 +43,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
		nmax = 0;
		rho = NULL;
		fp = NULL;
		need_embed = NULL;
		count_embed = NULL;
		map = NULL;
		type2frho = NULL;

		@@ -78,7 +78,7 @@ PairEAM::~PairEAM()

		memory->destroy(rho);
		memory->destroy(fp);
		memory->destroy(need_embed);
		memory->destroy(count_embed);

		if (allocated) {
		memory->destroy(setflag);
		@@ -153,11 +153,11 @@ void PairEAM::compute(int eflag, int vflag)
		if (atom->nmax > nmax) {
		memory->destroy(rho);
		memory->destroy(fp);
		memory->destroy(need_embed);
		memory->destroy(count_embed);
		nmax = atom->nmax;
		memory->create(rho,nmax,"pair:rho");
		memory->create(fp,nmax,"pair:fp");
		memory->create(need_embed,nmax,"pair:need_embed");
		memory->create(count_embed,nmax,"pair:count_embed");
		}

		double **x = atom->x;
		@@ -234,7 +234,7 @@ void PairEAM::compute(int eflag, int vflag)
		p = MIN(p,1.0);
		coeff = frho_spline[type2frho[type[i]]][m];
		fp[i] = (coeff[0]p + coeff[1])p + coeff[2];
		need_embed[i] = 1;
		count_embed[i] = 0;
		if (eflag) {
		phi = ((coeff[3]p + coeff[4])p + coeff[5])*p + coeff[6];
		if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
		@@ -271,6 +271,7 @@ void PairEAM::compute(int eflag, int vflag)
		rsq = delxdelx + delydely + delz*delz;

		if (rsq < cutforcesq) {
		++count_embed[i];
		jtype = type[j];
		r = sqrt(rsq);
		p = r*rdr + 1.0;
		@@ -807,7 +808,7 @@ double PairEAM::single(int i, int j, int itype, int jtype,
		double r,p,rhoip,rhojp,z2,z2p,recip,phi,phip,psip;
		double *coeff;

		if (need_embed[i]) {
		if (count_embed[i] > 0) {
		p = rho[i]*rdrho + 1.0;
		m = static_cast<int> (p);
		m = MAX(1,MIN(m,nrho-1));
		@@ -816,7 +817,7 @@ double PairEAM::single(int i, int j, int itype, int jtype,
		coeff = frho_spline[type2frho[itype]][m];
		phi = ((coeff[3]p + coeff[4])p + coeff[5])*p + coeff[6];
		if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
		need_embed[i] = 0;
		phi *= 1.0/static_cast<double>(count_embed[i]);
		} else phi = 0.0;

		r = sqrt(rsq);

src/MANYBODY/pair_eam.h

+1 −1

Original line number	Diff line number	Diff line
		@@ -70,7 +70,7 @@ class PairEAM : public Pair {
		// per-atom arrays

		double rho,fp;
		int *need_embed;
		int *count_embed;

		// potentials as file data

src/OPT/pair_eam_opt.cpp

+8 −4

Original line number	Diff line number	Diff line
		@@ -80,11 +80,13 @@ void PairEAMOpt::eval()
		// grow energy array if necessary

		if (atom->nmax > nmax) {
		memory->sfree(rho);
		memory->sfree(fp);
		memory->destroy(rho);
		memory->destroy(fp);
		memory->destroy(count_embed);
		nmax = atom->nmax;
		rho = (double ) memory->smalloc(nmaxsizeof(double),"pair:rho");
		fp = (double ) memory->smalloc(nmaxsizeof(double),"pair:fp");
		memory->create(rho,nmax,"pair:rho");
		memory->create(fp,nmax,"pair:fp");
		memory->create(count_embed,nmax,"pair:count_embed");
		}

		double** _noalias x = atom->x;
		@@ -238,6 +240,7 @@ void PairEAMOpt::eval()
		++m;
		coeff = frho_spline[type2frho[type[i]]][m];
		fp[i] = (coeff[0]p + coeff[1])p + coeff[2];
		count_embed[i] = 0;
		if (EFLAG) {
		double phi = ((coeff[3]p + coeff[4])p + coeff[5])*p + coeff[6];
		if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
		@@ -280,6 +283,7 @@ void PairEAMOpt::eval()
		double rsq = delxdelx + delydely + delz*delz;

		if (rsq < tmp_cutforcesq) {
		++count_embed[i];
		jtype = type[j] - 1;
		double r = sqrt(rsq);
		double rhoip,rhojp,z2,z2p;

src/USER-OMP/pair_eam_omp.cpp

+4 −0

Original line number	Diff line number	Diff line
		@@ -51,9 +51,11 @@ void PairEAMOMP::compute(int eflag, int vflag)
		if (atom->nmax > nmax) {
		memory->destroy(rho);
		memory->destroy(fp);
		memory->destroy(count_embed);
		nmax = atom->nmax;
		memory->create(rho,nthreads*nmax,"pair:rho");
		memory->create(fp,nmax,"pair:fp");
		memory->create(count_embed,nmax,"pair:count_embed");
		}

		#if defined(_OPENMP)
		@@ -198,6 +200,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
		p = MIN(p,1.0);
		coeff = frho_spline[type2frho[type[i]]][m];
		fp[i] = (coeff[0]p + coeff[1])p + coeff[2];
		count_embed[i] = 0;
		if (EFLAG) {
		phi = ((coeff[3]p + coeff[4])p + coeff[5])*p + coeff[6];
		if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
		@@ -243,6 +246,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
		rsq = delxdelx + delydely + delz*delz;

		if (rsq < cutforcesq) {
		++count_embed[i];
		jtype = type[j];
		r = sqrt(rsq);
		p = r*rdr + 1.0;

Admin message