Merge pull request #1551 from martok/meamc-embedding (691fc357) · Commits · 郑智淋 / lammps

doc/src/pair_meamc.txt

+2 −1

Original line number	Diff line number	Diff line
		@@ -147,7 +147,8 @@ asub = "A" parameter for MEAM (see e.g. "(Baskes)"_#Baskes) :pre

		The alpha, b0, b1, b2, b3, t0, t1, t2, t3 parameters correspond to the
		standard MEAM parameters in the literature "(Baskes)"_#Baskes (the b
		parameters are the standard beta parameters). The rozero parameter is
		parameters are the standard beta parameters). Note that only parameters
		normalized to t0 = 1.0 are supported. The rozero parameter is
		an element-dependent density scaling that weights the reference
		background density (see e.g. equation 4.5 in "(Gullet)"_#Gullet) and
		is typically 1.0 for single-element systems. The ibar parameter

src/USER-MEAMC/meam.h

+5 −2

Original line number	Diff line number	Diff line
		@@ -93,8 +93,9 @@ private:
		int augt1, ialloy, mix_ref_t, erose_form;
		int emb_lin_neg, bkgd_dyn;
		double gsmooth_factor;
		int vind2D[3][3], vind3D[3][3][3];
		int v2D[6], v3D[10];

		int vind2D[3][3], vind3D[3][3][3]; // x-y-z to Voigt-like index
		int v2D[6], v3D[10]; // multiplicity of Voigt index (i.e. [1] -> xy+yx = 2

		int nr, nrar;
		double dr, rdrar;
		@@ -121,6 +122,7 @@ protected:
		else if (xi <= 0.0)
		return 0.0;
		else {
		// ( 1.d0 - (1.d0 - xi)4 )2, but with better codegen
		a = 1.0 - xi;
		a = a; a = a;
		a = 1.0 - a;
		@@ -187,6 +189,7 @@ protected:
		double G_gam(const double gamma, const int ibar, int &errorflag) const;
		double dG_gam(const double gamma, const int ibar, double &dG) const;
		static double zbl(const double r, const int z1, const int z2);
		double embedding(const double A, const double Ec, const double rhobar, double& dF) const;
		static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form);

		static void get_shpfcn(const lattice_t latt, double (&s)[3]);

src/USER-MEAMC/meam_dens_final.cpp

+7 −28

Original line number	Diff line number	Diff line
		@@ -10,7 +10,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
		int i, elti;
		int m;
		double rhob, G, dG, Gbar, dGbar, gam, shp[3], Z;
		double B, denom, rho_bkgd;
		double denom, rho_bkgd, Fl;

		// Complete the calculation of density

		@@ -111,35 +111,14 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
		dgamma3[i] = 0.0;
		}

		B = this->A_meam[elti] * this->Ec_meam[elti][elti];
		Fl = embedding(this->A_meam[elti], this->Ec_meam[elti][elti], rhob, frhop[i]);

		if (!iszero(rhob)) {
		if (this->emb_lin_neg == 1 && rhob <= 0) {
		frhop[i] = -B;
		} else {
		frhop[i] = B * (log(rhob) + 1.0);
		}
		if (eflag_either != 0) {
		if (eflag_global != 0) {
		if (this->emb_lin_neg == 1 && rhob <= 0) {
		eng_vdwl = eng_vdwl - B * rhob;
		} else {
		eng_vdwl = eng_vdwl + B * rhob * log(rhob);
		}
		eng_vdwl = eng_vdwl + Fl;
		}
		if (eflag_atom != 0) {
		if (this->emb_lin_neg == 1 && rhob <= 0) {
		eatom[i] = eatom[i] - B * rhob;
		} else {
		eatom[i] = eatom[i] + B * rhob * log(rhob);
		}
		}
		}
		} else {
		if (this->emb_lin_neg == 1) {
		frhop[i] = -B;
		} else {
		frhop[i] = B;
		eatom[i] = eatom[i] + Fl;
		}
		}
		}

src/USER-MEAMC/meam_dens_init.cpp

+79 −79

Original line number	Diff line number	Diff line
		@@ -127,33 +127,39 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
		delyij = yjtmp - yitmp;
		delzij = zjtmp - zitmp;
		rij2 = delxij * delxij + delyij * delyij + delzij * delzij;
		rij = sqrt(rij2);

		if (rij2 > this->cutforcesq) {
		dscrfcn[jn] = 0.0;
		scrfcn[jn] = 0.0;
		fcpair[jn] = 0.0;
		continue;
		}

		const double rbound = this->ebound_meam[elti][eltj] * rij2;
		if (rij > this->rc_meam) {
		fcij = 0.0;
		dfcij = 0.0;
		sij = 0.0;
		} else {
		rnorm = (this->rc_meam - rij) * drinv;
		rij = sqrt(rij2);
		rnorm = (this->cutforce - rij) * drinv;
		sij = 1.0;

		// if rjk2 > ebound*rijsq, atom k is definitely outside the ellipse
		for (kn = 0; kn < numneigh_full; kn++) {
		k = firstneigh_full[kn];
		if (k == j) continue;
		eltk = fmap[type[k]];
		if (eltk < 0) continue;
		if (k == j) continue;

		delxjk = x[k][0] - xjtmp;
		delyjk = x[k][1] - yjtmp;
		delzjk = x[k][2] - zjtmp;
		xktmp = x[k][0];
		yktmp = x[k][1];
		zktmp = x[k][2];

		delxjk = xktmp - xjtmp;
		delyjk = yktmp - yjtmp;
		delzjk = zktmp - zjtmp;
		rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
		if (rjk2 > rbound) continue;

		delxik = x[k][0] - xitmp;
		delyik = x[k][1] - yitmp;
		delzik = x[k][2] - zitmp;
		delxik = xktmp - xitmp;
		delyik = yktmp - yitmp;
		delzik = zktmp - zitmp;
		rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
		if (rik2 > rbound) continue;

		@@ -185,32 +191,26 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
		fc = dfcut(rnorm, dfc);
		fcij = fc;
		dfcij = dfc * drinv;
		}

		// Now compute derivatives
		dscrfcn[jn] = 0.0;
		sfcij = sij * fcij;
		if (iszero(sfcij) \|\| iszero(sfcij - 1.0))
		goto LABEL_100;

		if (!iszero(sfcij) && !iszero(sfcij - 1.0)) {
		for (kn = 0; kn < numneigh_full; kn++) {
		k = firstneigh_full[kn];
		if (k == j) continue;
		eltk = fmap[type[k]];
		if (eltk < 0) continue;

		xktmp = x[k][0];
		yktmp = x[k][1];
		zktmp = x[k][2];
		delxjk = xktmp - xjtmp;
		delyjk = yktmp - yjtmp;
		delzjk = zktmp - zjtmp;
		delxjk = x[k][0] - xjtmp;
		delyjk = x[k][1] - yjtmp;
		delzjk = x[k][2] - zjtmp;
		rjk2 = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk;
		if (rjk2 > rbound) continue;

		delxik = xktmp - xitmp;
		delyik = yktmp - yitmp;
		delzik = zktmp - zitmp;
		delxik = x[k][0] - xitmp;
		delyik = x[k][1] - yitmp;
		delzik = x[k][2] - zitmp;
		rik2 = delxik * delxik + delyik * delyik + delzik * delzik;
		if (rik2 > rbound) continue;

		@@ -245,8 +245,8 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double**
		coef1 = sfcij;
		coef2 = sij * dfcij / rij;
		dscrfcn[jn] = dscrfcn[jn] * coef1 - coef2;
		}

		LABEL_100:
		scrfcn[jn] = sij;
		fcpair[jn] = fcij;
		}

src/USER-MEAMC/meam_funcs.cpp

+24 −1

Original line number	Diff line number	Diff line
		@@ -52,7 +52,7 @@ MEAM::G_gam(const double gamma, const int ibar, int& errorflag) const
		case 1:
		return MathSpecial::fm_exp(gamma / 2.0);
		case 3:
		return 2.0 / (1.0 + exp(-gamma));
		return 2.0 / (1.0 + MathSpecial::fm_exp(-gamma));
		case -5:
		if ((1.0 + gamma) >= 0) {
		return sqrt(1.0 + gamma);
		@@ -143,6 +143,29 @@ MEAM::zbl(const double r, const int z1, const int z2)
		return result;
		}

		//-----------------------------------------------------------------------------
		// Compute embedding function F(rhobar) and derivative F'(rhobar), eqn I.5
		//
		double
		MEAM::embedding(const double A, const double Ec, const double rhobar, double& dF) const
		{
		const double AEc = A * Ec;

		if (rhobar > 0.0) {
		const double lrb = log(rhobar);
		dF = AEc * (1.0 + lrb);
		return AEc * rhobar * lrb;
		} else {
		if (this->emb_lin_neg == 0) {
		dF = 0.0;
		return 0.0;
		} else {
		dF = - AEc;
		return - AEc * rhobar;
		}
		}
		}

		//-----------------------------------------------------------------------------
		// Compute Rose energy function, I.16
		//

Admin message