Factor out embedding function, make sure documented logic for emb_lin_neg is obeyed (2f72b9f1) · Commits · 郑智淋 / lammps

src/USER-MEAMC/meam.h

+1 −0

Original line number	Diff line number	Diff line
		@@ -187,6 +187,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]);

+8 −29

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;
		}
		}
		}

+22 −0

Original line number	Diff line number	Diff line
		@@ -143,6 +143,28 @@ 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) {
		dF = AEc * (1.0 + log(rhobar));
		return AEc * rhobar * log(rhobar);
		} 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
		//

+5 −17

Original line number	Diff line number	Diff line
		@@ -319,7 +319,7 @@ MEAM::phi_meam(double r, int a, int b)
		double t11av, t21av, t31av, t12av, t22av, t32av;
		double G1, G2, s1[3], s2[3], rho0_1, rho0_2;
		double Gam1, Gam2, Z1, Z2;
		double rhobar1, rhobar2, F1, F2;
		double rhobar1, rhobar2, F1, F2, dF;
		double rho01, rho11, rho21, rho31;
		double rho02, rho12, rho22, rho32;
		double scalfac, phiaa, phibb;
		@@ -447,22 +447,10 @@ MEAM::phi_meam(double r, int a, int b)
		}

		// compute embedding functions, eqn I.5
		if (iszero(rhobar1))
		F1 = 0.0;
		else {
		if (this->emb_lin_neg == 1 && rhobar1 <= 0)
		F1 = -this->A_meam[a] * this->Ec_meam[a][a] * rhobar1;
		else
		F1 = this->A_meam[a] * this->Ec_meam[a][a] * rhobar1 * log(rhobar1);
		}
		if (iszero(rhobar2))
		F2 = 0.0;
		else {
		if (this->emb_lin_neg == 1 && rhobar2 <= 0)
		F2 = -this->A_meam[b] * this->Ec_meam[b][b] * rhobar2;
		else
		F2 = this->A_meam[b] * this->Ec_meam[b][b] * rhobar2 * log(rhobar2);
		}

		F1 = embedding(this->A_meam[a], this->Ec_meam[a][a], rhobar1, dF);
		F2 = embedding(this->A_meam[b], this->Ec_meam[b][b], rhobar2, dF);


		// compute Rose function, I.16
		Eu = erose(r, this->re_meam[a][b], this->alpha_meam[a][b], this->Ec_meam[a][b], this->repuls_meam[a][b],