another speedup by folding dsij() into meam_force() (060e3297) · Commits · 郑智淋 / lammps

src/USER-MEAMC/meam.h

+0 −2

Original line number	Diff line number	Diff line
		@@ -198,8 +198,6 @@ protected:
		int* firstneigh, int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap);
		void calc_rho1(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,
		double* scrfcn, double* fcpair);
		void dsij(int i, int j, int k, int jn, int numneigh, double rij2, double* dsij1, double* dsij2, int ntype,
		int* type, int* fmap, double** x, double* scrfcn, double* fcpair);

		void alloyparams();
		void compute_pair_meam();

src/USER-MEAMC/meam_dens_init.cpp

+0 −60

Original line number	Diff line number	Diff line
		@@ -354,63 +354,3 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, int numneigh
		}
		}

		// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

		void
		MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2, double* dsij1, double* dsij2, int ntype,
		int* type, int* fmap, double** x, double* scrfcn, double* fcpair)
		{
		// Inputs: i,j,k = id's of 3 atom triplet
		// jn = id of i-j pair
		// rij2 = squared distance between i and j
		// Outputs: dsij1 = deriv. of sij w.r.t. rik
		// dsij2 = deriv. of sij w.r.t. rjk

		int elti, eltj, eltk;
		double rik2, rjk2;

		double dxik, dyik, dzik;
		double dxjk, dyjk, dzjk;
		double rbound, delc, sij, xik, xjk, cikj, sikj, dfc, a;
		double Cmax, Cmin, dCikj1, dCikj2;

		sij = scrfcn[jn] * fcpair[jn];
		elti = fmap[type[i]];
		eltj = fmap[type[j]];
		eltk = fmap[type[k]];
		Cmax = this->Cmax_meam[elti][eltj][eltk];
		Cmin = this->Cmin_meam[elti][eltj][eltk];

		*dsij1 = 0.0;
		*dsij2 = 0.0;
		if (!iszero(sij) && !iszero(sij - 1.0)) {
		rbound = rij2 * this->ebound_meam[elti][eltj];
		delc = Cmax - Cmin;
		dxjk = x[k][0] - x[j][0];
		dyjk = x[k][1] - x[j][1];
		dzjk = x[k][2] - x[j][2];
		rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
		if (rjk2 <= rbound) {
		dxik = x[k][0] - x[i][0];
		dyik = x[k][1] - x[i][1];
		dzik = x[k][2] - x[i][2];
		rik2 = dxik * dxik + dyik * dyik + dzik * dzik;
		if (rik2 <= rbound) {
		xik = rik2 / rij2;
		xjk = rjk2 / rij2;
		a = 1 - (xik - xjk) * (xik - xjk);
		if (!iszero(a)) {
		cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
		if (cikj >= Cmin && cikj <= Cmax) {
		cikj = (cikj - Cmin) / delc;
		sikj = dfcut(cikj, dfc);
		dCfunc2(rij2, rik2, rjk2, dCikj1, dCikj2);
		a = sij / delc * dfc / sikj;
		dsij1 = a dCikj1;
		dsij2 = a dCikj2;
		}
		}
		}
		}
		}
		}

src/USER-MEAMC/meam_force.cpp

+56 −15

Original line number	Diff line number	Diff line
		@@ -435,8 +435,49 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int
		k = firstneigh_full[kn];
		eltk = fmap[type[k]];
		if (k != j && eltk >= 0) {
		dsij(i, j, k, jn, numneigh, rij2, &dsij1, &dsij2, ntype, type, fmap, x, &scrfcn[fnoffset],
		&fcpair[fnoffset]);
		double xik, xjk, cikj, sikj, dfc, a;
		double dCikj1, dCikj2;
		double dxik, dyik, dzik;
		double dxjk, dyjk, dzjk;
		double delc, rik2, rjk2;

		sij = scrfcn[jn+fnoffset] * fcpair[jn+fnoffset];
		const double Cmax = this->Cmax_meam[elti][eltj][eltk];
		const double Cmin = this->Cmin_meam[elti][eltj][eltk];

		dsij1 = 0.0;
		dsij2 = 0.0;
		if (!iszero(sij) && !iszero(sij - 1.0)) {
		const double rbound = rij2 * this->ebound_meam[elti][eltj];
		delc = Cmax - Cmin;
		dxjk = x[k][0] - x[j][0];
		dyjk = x[k][1] - x[j][1];
		dzjk = x[k][2] - x[j][2];
		rjk2 = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk;
		if (rjk2 <= rbound) {
		dxik = x[k][0] - x[i][0];
		dyik = x[k][1] - x[i][1];
		dzik = x[k][2] - x[i][2];
		rik2 = dxik * dxik + dyik * dyik + dzik * dzik;
		if (rik2 <= rbound) {
		xik = rik2 / rij2;
		xjk = rjk2 / rij2;
		a = 1 - (xik - xjk) * (xik - xjk);
		if (!iszero(a)) {
		cikj = (2.0 * (xik + xjk) + a - 2.0) / a;
		if (cikj >= Cmin && cikj <= Cmax) {
		cikj = (cikj - Cmin) / delc;
		sikj = dfcut(cikj, dfc);
		dCfunc2(rij2, rik2, rjk2, dCikj1, dCikj2);
		a = sij / delc * dfc / sikj;
		dsij1 = a * dCikj1;
		dsij2 = a * dCikj2;
		}
		}
		}
		}
		}

		if (!iszero(dsij1) \|\| !iszero(dsij2)) {
		force1 = dUdsij * dsij1;
		force2 = dUdsij * dsij2;

Admin message