Change indexing of remaining variables and locals (8fca667e) · Commits · 郑智淋 / lammps

src/USER-MEAMC/meam.h

+2 −15

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

		int nr, nrar;
		double dr, rdrar;
		@@ -181,18 +181,5 @@ public:
		arr[__i][__j][__k] = v; \
		}

		/*
		Fortran Array Semantics in C.
		- Stack-Allocated and global arrays are 1-based, declared as foo[N+1] and simply ignoring the first element
		- Multi-Dimensional MUST be declared in reverse, so that the order when accessing is the same as in Fortran
		- arrays that are passed externally via the meam_* functions must use the arr*v() functions below
		(or be used with 0-based indexing)
		*/

		// access data with same index as used in fortran (1-based)
		#define arr1v(ptr, i) ptr[i - 1]
		#define arr2v(ptr, i, j) ptr[j - 1][i - 1]


		};
		#endif

src/USER-MEAMC/meam_dens_final.cpp

+15 −15

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

		// Complete the calculation of density
		@@ -38,19 +38,19 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
		rho1[i] = 0.0;
		rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i];
		rho3[i] = 0.0;
		for (m = 1; m <= 3; m++) {
		rho1[i] = rho1[i] + arr2v(arho1, m, i+1) * arr2v(arho1, m, i+1);
		rho3[i] = rho3[i] - 3.0 / 5.0 * arr2v(arho3b, m, i+1) * arr2v(arho3b, m, i+1);
		for (m = 0; m < 3; m++) {
		rho1[i] = rho1[i] + arho1[i][m] * arho1[i][m];
		rho3[i] = rho3[i] - 3.0 / 5.0 * arho3b[i][m] * arho3b[i][m];
		}
		for (m = 1; m <= 6; m++) {
		for (m = 0; m < 6; m++) {
		rho2[i] =
		rho2[i] +
		this->v2D[m] * arr2v(arho2, m, i+1) * arr2v(arho2, m, i+1);
		this->v2D[m] * arho2[i][m] * arho2[i][m];
		}
		for (m = 1; m <= 10; m++) {
		for (m = 0; m < 10; m++) {
		rho3[i] =
		rho3[i] +
		this->v3D[m] * arr2v(arho3, m, i+1) * arr2v(arho3, m, i+1);
		this->v3D[m] * arho3[i][m] * arho3[i][m];
		}

		if (rho0[i] > 0.0) {
		@@ -88,13 +88,13 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
		dGbar = 0.0;
		} else {
		if (this->mix_ref_t == 1) {
		gam = (t_ave[i][0] * shp[1] + t_ave[i][1] * shp[2] +
		t_ave[i][2] * shp[3]) /
		gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] +
		t_ave[i][2] * shp[2]) /
		(Z * Z);
		} else {
		gam = (this->t1_meam[elti] * shp[1] +
		this->t2_meam[elti] * shp[2] +
		this->t3_meam[elti] * shp[3]) /
		gam = (this->t1_meam[elti] * shp[0] +
		this->t2_meam[elti] * shp[1] +
		this->t3_meam[elti] * shp[2]) /
		(Z * Z);
		}
		G_gam(gam, this->ibar_meam[elti], &Gbar, errorflag);
		@@ -106,8 +106,8 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global,
		Gbar = 1.0;
		dGbar = 0.0;
		} else {
		gam = (t_ave[i][0] * shp[1] + t_ave[i][1] * shp[2] +
		t_ave[i][2] * shp[3]) /
		gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] +
		t_ave[i][2] * shp[2]) /
		(Z * Z);
		dG_gam(gam, this->ibar_meam[elti], &Gbar, &dGbar);
		}

src/USER-MEAMC/meam_dens_init.cpp

+18 −22

Original line number	Diff line number	Diff line
		@@ -246,7 +246,7 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
		{
		int jn, j, m, n, p, elti, eltj;
		int nv2, nv3;
		double xtmp, ytmp, ztmp, delij[3 + 1], rij2, rij, sij;
		double xtmp, ytmp, ztmp, delij[3], rij2, rij, sij;
		double ai, aj, rhoa0j, rhoa1j, rhoa2j, rhoa3j, A1j, A2j, A3j;
		// double G,Gbar,gam,shp[3+1];
		double ro0i, ro0j;
		@@ -260,10 +260,10 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
		if (!iszero(scrfcn[jn])) {
		j = firstneigh[jn];
		sij = scrfcn[jn] * fcpair[jn];
		delij[1] = x[j][0] - xtmp;
		delij[2] = x[j][1] - ytmp;
		delij[3] = x[j][2] - ztmp;
		rij2 = delij[1] * delij[1] + delij[2] * delij[2] + delij[3] * delij[3];
		delij[0] = x[j][0] - xtmp;
		delij[1] = x[j][1] - ytmp;
		delij[2] = x[j][2] - ztmp;
		rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
		if (rij2 < this->cutforcesq) {
		eltj = fmap[type[j]];
		rij = sqrt(rij2);
		@@ -315,24 +315,20 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x,
		A1i = rhoa1i / rij;
		A2i = rhoa2i / rij2;
		A3i = rhoa3i / (rij2 * rij);
		nv2 = 1;
		nv3 = 1;
		for (m = 1; m <= 3; m++) {
		arr2v(arho1, m, i+1) = arr2v(arho1, m, i+1) + A1j * delij[m];
		arr2v(arho1, m, j+1) = arr2v(arho1, m, j+1) - A1i * delij[m];
		arr2v(arho3b, m, i+1) = arr2v(arho3b, m, i+1) + rhoa3j * delij[m] / rij;
		arr2v(arho3b, m, j+1) = arr2v(arho3b, m, j+1) - rhoa3i * delij[m] / rij;
		for (n = m; n <= 3; n++) {
		arr2v(arho2, nv2, i+1) =
		arr2v(arho2, nv2, i+1) + A2j * delij[m] * delij[n];
		arr2v(arho2, nv2, j+1) =
		arr2v(arho2, nv2, j+1) + A2i * delij[m] * delij[n];
		nv2 = 0;
		nv3 = 0;
		for (m = 0; m < 3; m++) {
		arho1[i][m] = arho1[i][m] + A1j * delij[m];
		arho1[j][m] = arho1[j][m] - A1i * delij[m];
		arho3b[i][m] = arho3b[i][m] + rhoa3j * delij[m] / rij;
		arho3b[j][m] = arho3b[j][m] - rhoa3i * delij[m] / rij;
		for (n = m; n < 3; n++) {
		arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n];
		arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n];
		nv2 = nv2 + 1;
		for (p = n; p <= 3; p++) {
		arr2v(arho3, nv3, i+1) =
		arr2v(arho3, nv3, i+1) + A3j * delij[m] * delij[n] * delij[p];
		arr2v(arho3, nv3, j+1) =
		arr2v(arho3, nv3, j+1) - A3i * delij[m] * delij[n] * delij[p];
		for (p = n; p < 3; p++) {
		arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p];
		arho3[j][nv3] = arho3[j][nv3] - A3i * delij[m] * delij[n] * delij[p];
		nv3 = nv3 + 1;
		}
		}

src/USER-MEAMC/meam_force.cpp

+80 −105

Original line number	Diff line number	Diff line
		@@ -35,15 +35,15 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		{
		int j, jn, k, kn, kk, m, n, p, q;
		int nv2, nv3, elti, eltj, eltk, ind;
		double xitmp, yitmp, zitmp, delij[3 + 1], rij2, rij, rij3;
		double delik[3 + 1], deljk[3 + 1], v[6 + 1], fi[3 + 1], fj[3 + 1];
		double xitmp, yitmp, zitmp, delij[3], rij2, rij, rij3;
		double delik[3], deljk[3], v[6], fi[3], fj[3];
		double third, sixth;
		double pp, dUdrij, dUdsij, dUdrijm[3 + 1], force, forcem;
		double pp, dUdrij, dUdsij, dUdrijm[3], force, forcem;
		double r, recip, phi, phip;
		double sij;
		double a1, a1i, a1j, a2, a2i, a2j;
		double a3i, a3j;
		double shpi[3 + 1], shpj[3 + 1];
		double shpi[3], shpj[3];
		double ai, aj, ro0i, ro0j, invrei, invrej;
		double rhoa0j, drhoa0j, rhoa0i, drhoa0i;
		double rhoa1j, drhoa1j, rhoa1i, drhoa1i;
		@@ -51,15 +51,15 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		double a3, a3a, rhoa3j, drhoa3j, rhoa3i, drhoa3i;
		double drho0dr1, drho0dr2, drho0ds1, drho0ds2;
		double drho1dr1, drho1dr2, drho1ds1, drho1ds2;
		double drho1drm1[3 + 1], drho1drm2[3 + 1];
		double drho1drm1[3], drho1drm2[3];
		double drho2dr1, drho2dr2, drho2ds1, drho2ds2;
		double drho2drm1[3 + 1], drho2drm2[3 + 1];
		double drho2drm1[3], drho2drm2[3];
		double drho3dr1, drho3dr2, drho3ds1, drho3ds2;
		double drho3drm1[3 + 1], drho3drm2[3 + 1];
		double drho3drm1[3], drho3drm2[3];
		double dt1dr1, dt1dr2, dt1ds1, dt1ds2;
		double dt2dr1, dt2dr2, dt2ds1, dt2ds2;
		double dt3dr1, dt3dr2, dt3ds1, dt3ds2;
		double drhodr1, drhodr2, drhods1, drhods2, drhodrm1[3 + 1], drhodrm2[3 + 1];
		double drhodr1, drhodr2, drhods1, drhods2, drhodrm1[3], drhodrm2[3];
		double arg;
		double arg1i1, arg1j1, arg1i2, arg1j2, arg1i3, arg1j3, arg3i3, arg3j3;
		double dsij1, dsij2, force1, force2;
		@@ -86,10 +86,10 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) {

		sij = scrfcn[fnoffset + jn] * fcpair[fnoffset + jn];
		delij[1] = x[j][0] - xitmp;
		delij[2] = x[j][1] - yitmp;
		delij[3] = x[j][2] - zitmp;
		rij2 = delij[1] * delij[1] + delij[2] * delij[2] + delij[3] * delij[3];
		delij[0] = x[j][0] - xitmp;
		delij[1] = x[j][1] - yitmp;
		delij[2] = x[j][2] - zitmp;
		rij2 = delij[0] * delij[0] + delij[1] * delij[1] + delij[2] * delij[2];
		if (rij2 < this->cutforcesq) {
		rij = sqrt(rij2);
		r = rij;
		@@ -176,8 +176,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		drhoa3i = drhoa3i * this->t3_meam[elti];
		}

		nv2 = 1;
		nv3 = 1;
		nv2 = 0;
		nv3 = 0;
		arg1i1 = 0.0;
		arg1j1 = 0.0;
		arg1i2 = 0.0;
		@@ -186,23 +186,23 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		arg1j3 = 0.0;
		arg3i3 = 0.0;
		arg3j3 = 0.0;
		for (n = 1; n <= 3; n++) {
		for (p = n; p <= 3; p++) {
		for (q = p; q <= 3; q++) {
		for (n = 0; n < 3; n++) {
		for (p = n; p < 3; p++) {
		for (q = p; q < 3; q++) {
		arg = delij[n] * delij[p] * delij[q] * this->v3D[nv3];
		arg1i3 = arg1i3 + arr2v(arho3, nv3, i+1) * arg;
		arg1j3 = arg1j3 - arr2v(arho3, nv3, j+1) * arg;
		arg1i3 = arg1i3 + arho3[i][nv3] * arg;
		arg1j3 = arg1j3 - arho3[j][nv3] * arg;
		nv3 = nv3 + 1;
		}
		arg = delij[n] * delij[p] * this->v2D[nv2];
		arg1i2 = arg1i2 + arr2v(arho2, nv2, i+1) * arg;
		arg1j2 = arg1j2 + arr2v(arho2, nv2, j+1) * arg;
		arg1i2 = arg1i2 + arho2[i][nv2] * arg;
		arg1j2 = arg1j2 + arho2[j][nv2] * arg;
		nv2 = nv2 + 1;
		}
		arg1i1 = arg1i1 + arr2v(arho1, n, i+1) * delij[n];
		arg1j1 = arg1j1 - arr2v(arho1, n, j+1) * delij[n];
		arg3i3 = arg3i3 + arr2v(arho3b, n, i+1) * delij[n];
		arg3j3 = arg3j3 - arr2v(arho3b, n, j+1) * delij[n];
		arg1i1 = arg1i1 + arho1[i][n] * delij[n];
		arg1j1 = arg1j1 - arho1[j][n] * delij[n];
		arg3i3 = arg3i3 + arho3b[i][n] * delij[n];
		arg3j3 = arg3j3 - arho3b[j][n] * delij[n];
		}

		// rho0 terms
		@@ -214,9 +214,9 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1;
		drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1;
		a1 = 2.0 * sij / rij;
		for (m = 1; m <= 3; m++) {
		drho1drm1[m] = a1 * rhoa1j * arr2v(arho1, m, i+1);
		drho1drm2[m] = -a1 * rhoa1i * arr2v(arho1, m, j+1);
		for (m = 0; m < 3; m++) {
		drho1drm1[m] = a1 * rhoa1j * arho1[i][m];
		drho1drm2[m] = -a1 * rhoa1i * arho1[j][m];
		}

		// rho2 terms
		@@ -226,14 +226,12 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 -
		2.0 / 3.0 * arho2b[j] * drhoa2i * sij;
		a2 = 4 * sij / rij2;
		for (m = 1; m <= 3; m++) {
		for (m = 0; m < 3; m++) {
		drho2drm1[m] = 0.0;
		drho2drm2[m] = 0.0;
		for (n = 1; n <= 3; n++) {
		drho2drm1[m] = drho2drm1[m] +
		arr2v(arho2, this->vind2D[m][n], i+1) * delij[n];
		drho2drm2[m] = drho2drm2[m] -
		arr2v(arho2, this->vind2D[m][n], j+1) * delij[n];
		for (n = 0; n < 3; n++) {
		drho2drm1[m] = drho2drm1[m] + arho2[i][this->vind2D[m][n]] * delij[n];
		drho2drm2[m] = drho2drm2[m] - arho2[j][this->vind2D[m][n]] * delij[n];
		}
		drho2drm1[m] = a2 * rhoa2j * drho2drm1[m];
		drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m];
		@@ -249,24 +247,22 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		a3a * (drhoa3i - rhoa3i / rij) * arg3j3;
		a3 = 6 * sij / rij3;
		a3a = 6 * sij / (5 * rij);
		for (m = 1; m <= 3; m++) {
		for (m = 0; m < 3; m++) {
		drho3drm1[m] = 0.0;
		drho3drm2[m] = 0.0;
		nv2 = 1;
		for (n = 1; n <= 3; n++) {
		for (p = n; p <= 3; p++) {
		nv2 = 0;
		for (n = 0; n < 3; n++) {
		for (p = n; p < 3; p++) {
		arg = delij[n] * delij[p] * this->v2D[nv2];
		drho3drm1[m] = drho3drm1[m] +
		arr2v(arho3, this->vind3D[m][n][p], i+1) * arg;
		drho3drm2[m] = drho3drm2[m] +
		arr2v(arho3, this->vind3D[m][n][p], j+1) * arg;
		drho3drm1[m] = drho3drm1[m] + arho3[i][this->vind3D[m][n][p]] * arg;
		drho3drm2[m] = drho3drm2[m] + arho3[j][this->vind3D[m][n][p]] * arg;
		nv2 = nv2 + 1;
		}
		}
		drho3drm1[m] =
		(a3 * drho3drm1[m] - a3a * arr2v(arho3b, m, i+1)) * rhoa3j;
		(a3 * drho3drm1[m] - a3a * arho3b[i][m]) * rhoa3j;
		drho3drm2[m] =
		(-a3 * drho3drm2[m] + a3a * arr2v(arho3b, m, j+1)) * rhoa3i;
		(-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i;
		}

		// Compute derivatives of weighting functions t wrt rij
		@@ -346,15 +342,15 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		dt2dr1 * rho2[i] + t2i * drho2dr1 +
		dt3dr1 * rho3[i] + t3i * drho3dr1) -
		dgamma3[i] *
		(shpi[1] * dt1dr1 + shpi[2] * dt2dr1 + shpi[3] * dt3dr1);
		(shpi[0] * dt1dr1 + shpi[1] * dt2dr1 + shpi[2] * dt3dr1);
		drhodr2 =
		dgamma1[j] * drho0dr2 +
		dgamma2[j] * (dt1dr2 * rho1[j] + t1j * drho1dr2 +
		dt2dr2 * rho2[j] + t2j * drho2dr2 +
		dt3dr2 * rho3[j] + t3j * drho3dr2) -
		dgamma3[j] *
		(shpj[1] * dt1dr2 + shpj[2] * dt2dr2 + shpj[3] * dt3dr2);
		for (m = 1; m <= 3; m++) {
		(shpj[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2);
		for (m = 0; m < 3; m++) {
		drhodrm1[m] = 0.0;
		drhodrm2[m] = 0.0;
		drhodrm1[m] =
		@@ -442,14 +438,14 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		dt2ds1 * rho2[i] + t2i * drho2ds1 +
		dt3ds1 * rho3[i] + t3i * drho3ds1) -
		dgamma3[i] *
		(shpi[1] * dt1ds1 + shpi[2] * dt2ds1 + shpi[3] * dt3ds1);
		(shpi[0] * dt1ds1 + shpi[1] * dt2ds1 + shpi[2] * dt3ds1);
		drhods2 =
		dgamma1[j] * drho0ds2 +
		dgamma2[j] * (dt1ds2 * rho1[j] + t1j * drho1ds2 +
		dt2ds2 * rho2[j] + t2j * drho2ds2 +
		dt3ds2 * rho3[j] + t3j * drho3ds2) -
		dgamma3[j] *
		(shpj[1] * dt1ds2 + shpj[2] * dt2ds2 + shpj[3] * dt3ds2);
		(shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
		}

		// Compute derivatives of energy wrt rij, sij and rij[3]
		@@ -458,7 +454,7 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		if (!iszero(dscrfcn[fnoffset + jn])) {
		dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2;
		}
		for (m = 1; m <= 3; m++) {
		for (m = 0; m < 3; m++) {
		dUdrijm[m] =
		frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m];
		}
		@@ -466,37 +462,29 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		// Add the part of the force due to dUdrij and dUdsij

		force = dUdrij * recip + dUdsij * dscrfcn[fnoffset + jn];
		for (m = 1; m <= 3; m++) {
		for (m = 0; m < 3; m++) {
		forcem = delij[m] * force + dUdrijm[m];
		f[i][m-1] = f[i][m-1] + forcem;
		f[j][m-1] = f[j][m-1] - forcem;
		f[i][m] = f[i][m] + forcem;
		f[j][m] = f[j][m] - forcem;
		}

		// Tabulate per-atom virial as symmetrized stress tensor

		if (vflag_atom != 0) {
		fi[0] = delij[0] * force + dUdrijm[0];
		fi[1] = delij[1] * force + dUdrijm[1];
		fi[2] = delij[2] * force + dUdrijm[2];
		fi[3] = delij[3] * force + dUdrijm[3];
		v[0] = -0.5 * (delij[0] * fi[0]);
		v[1] = -0.5 * (delij[1] * fi[1]);
		v[2] = -0.5 * (delij[2] * fi[2]);
		v[3] = -0.5 * (delij[3] * fi[3]);
		v[4] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]);
		v[5] = -0.25 * (delij[1] * fi[3] + delij[3] * fi[1]);
		v[6] = -0.25 * (delij[2] * fi[3] + delij[3] * fi[2]);

		vatom[i][0] = vatom[i][0] + v[1];
		vatom[i][1] = vatom[i][1] + v[2];
		vatom[i][2] = vatom[i][2] + v[3];
		vatom[i][3] = vatom[i][3] + v[4];
		vatom[i][4] = vatom[i][4] + v[5];
		vatom[i][5] = vatom[i][5] + v[6];
		vatom[j][0] = vatom[j][0] + v[1];
		vatom[j][1] = vatom[j][1] + v[2];
		vatom[j][2] = vatom[j][2] + v[3];
		vatom[j][3] = vatom[j][3] + v[4];
		vatom[j][4] = vatom[j][4] + v[5];
		vatom[j][5] = vatom[j][5] + v[6];
		v[3] = -0.25 * (delij[0] * fi[1] + delij[1] * fi[0]);
		v[4] = -0.25 * (delij[0] * fi[2] + delij[2] * fi[0]);
		v[5] = -0.25 * (delij[1] * fi[2] + delij[2] * fi[1]);

		for (m = 0; m<6; m++) {
		vatom[i][m] = vatom[i][m] + v[m];
		vatom[j][m] = vatom[j][m] + v[m];
		}
		}

		// Now compute forces on other atoms k due to change in sij
		@@ -512,53 +500,40 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global,
		if (!iszero(dsij1) \|\| !iszero(dsij2)) {
		force1 = dUdsij * dsij1;
		force2 = dUdsij * dsij2;
		for (m = 1; m <= 3; m++) {
		delik[m] = arr2v(x, m, k+1) - arr2v(x, m, i+1);
		deljk[m] = arr2v(x, m, k+1) - arr2v(x, m, j+1);
		for (m = 0; m < 3; m++) {
		delik[m] = x[k][m] - x[i][m];
		deljk[m] = x[k][m] - x[j][m];
		}
		for (m = 1; m <= 3; m++) {
		arr2v(f, m, i+1) = arr2v(f, m, i+1) + force1 * delik[m];
		arr2v(f, m, j+1) = arr2v(f, m, j+1) + force2 * deljk[m];
		arr2v(f, m, k+1) = arr2v(f, m, k+1) - force1 * delik[m] - force2 * deljk[m];
		for (m = 0; m < 3; m++) {
		f[i][m] = f[i][m] + force1 * delik[m];
		f[j][m] = f[j][m] + force2 * deljk[m];
		f[k][m] = f[k][m] - force1 * delik[m] - force2 * deljk[m];
		}

		// Tabulate per-atom virial as symmetrized stress tensor

		if (vflag_atom != 0) {
		fi[0] = force1 * delik[0];
		fi[1] = force1 * delik[1];
		fi[2] = force1 * delik[2];
		fi[3] = force1 * delik[3];
		fj[0] = force2 * deljk[0];
		fj[1] = force2 * deljk[1];
		fj[2] = force2 * deljk[2];
		fj[3] = force2 * deljk[3];
		v[0] = -third * (delik[0] * fi[0] + deljk[0] * fj[0]);
		v[1] = -third * (delik[1] * fi[1] + deljk[1] * fj[1]);
		v[2] = -third * (delik[2] * fi[2] + deljk[2] * fj[2]);
		v[3] = -third * (delik[3] * fi[3] + deljk[3] * fj[3]);
		v[4] = -sixth * (delik[1] * fi[2] + deljk[1] * fj[2] +
		v[3] = -sixth * (delik[0] * fi[1] + deljk[0] * fj[1] +
		delik[1] * fi[0] + deljk[1] * fj[0]);
		v[4] = -sixth * (delik[0] * fi[2] + deljk[0] * fj[2] +
		delik[2] * fi[0] + deljk[2] * fj[0]);
		v[5] = -sixth * (delik[1] * fi[2] + deljk[1] * fj[2] +
		delik[2] * fi[1] + deljk[2] * fj[1]);
		v[5] = -sixth * (delik[1] * fi[3] + deljk[1] * fj[3] +
		delik[3] * fi[1] + deljk[3] * fj[1]);
		v[6] = -sixth * (delik[2] * fi[3] + deljk[2] * fj[3] +
		delik[3] * fi[2] + deljk[3] * fj[2]);

		vatom[i][0] = vatom[i][0] + v[1];
		vatom[i][1] = vatom[i][1] + v[2];
		vatom[i][2] = vatom[i][2] + v[3];
		vatom[i][3] = vatom[i][3] + v[4];
		vatom[i][4] = vatom[i][4] + v[5];
		vatom[i][5] = vatom[i][5] + v[6];
		vatom[j][0] = vatom[j][0] + v[1];
		vatom[j][1] = vatom[j][1] + v[2];
		vatom[j][2] = vatom[j][2] + v[3];
		vatom[j][3] = vatom[j][3] + v[4];
		vatom[j][4] = vatom[j][4] + v[5];
		vatom[j][5] = vatom[j][5] + v[6];
		vatom[k][0] = vatom[k][0] + v[1];
		vatom[k][1] = vatom[k][1] + v[2];
		vatom[k][2] = vatom[k][2] + v[3];
		vatom[k][3] = vatom[k][3] + v[4];
		vatom[k][4] = vatom[k][4] + v[5];
		vatom[k][5] = vatom[k][5] + v[6];

		for (m = 0; m<6; m++) {
		vatom[i][m] = vatom[i][m] + v[m];
		vatom[j][m] = vatom[j][m] + v[m];
		vatom[k][m] = vatom[k][m] + v[m];
		}
		}
		}
		}

src/USER-MEAMC/meam_setup_done.cpp

+39 −39

Original line number	Diff line number	Diff line
		@@ -34,14 +34,14 @@ MEAM::meam_setup_done(double* cutmax)
		alloyparams();

		// indices and factors for Voight notation
		nv2 = 1;
		nv3 = 1;
		for (m = 1; m <= 3; m++) {
		for (n = m; n <= 3; n++) {
		nv2 = 0;
		nv3 = 0;
		for (m = 0; m < 3; m++) {
		for (n = m; n < 3; n++) {
		this->vind2D[m][n] = nv2;
		this->vind2D[n][m] = nv2;
		nv2 = nv2 + 1;
		for (p = n; p <= 3; p++) {
		for (p = n; p < 3; p++) {
		this->vind3D[m][n][p] = nv3;
		this->vind3D[m][p][n] = nv3;
		this->vind3D[n][m][p] = nv3;
		@@ -53,23 +53,23 @@ MEAM::meam_setup_done(double* cutmax)
		}
		}

		this->v2D[1] = 1;
		this->v2D[0] = 1;
		this->v2D[1] = 2;
		this->v2D[2] = 2;
		this->v2D[3] = 2;
		this->v2D[4] = 1;
		this->v2D[5] = 2;
		this->v2D[6] = 1;
		this->v2D[3] = 1;
		this->v2D[4] = 2;
		this->v2D[5] = 1;

		this->v3D[1] = 1;
		this->v3D[0] = 1;
		this->v3D[1] = 3;
		this->v3D[2] = 3;
		this->v3D[3] = 3;
		this->v3D[4] = 3;
		this->v3D[5] = 6;
		this->v3D[6] = 3;
		this->v3D[7] = 1;
		this->v3D[4] = 6;
		this->v3D[5] = 3;
		this->v3D[6] = 1;
		this->v3D[7] = 3;
		this->v3D[8] = 3;
		this->v3D[9] = 3;
		this->v3D[10] = 1;
		this->v3D[9] = 1;

		nv2 = 0;
		for (m = 0; m < this->neltypes; m++) {
		@@ -379,7 +379,7 @@ MEAM::phi_meam(double r, int a, int b)
		{
		/unused:double a1,a2,a12;/
		double t11av, t21av, t31av, t12av, t22av, t32av;
		double G1, G2, s1[3 + 1], s2[3 + 1] /,s12[3+1]/, rho0_1, rho0_2;
		double G1, G2, s1[3], s2[3], rho0_1, rho0_2;
		double Gam1, Gam2, Z1, Z2;
		double rhobar1, rhobar2, F1, F2;
		double rho01, rho11, rho21, rho31;
		@@ -470,14 +470,14 @@ MEAM::phi_meam(double r, int a, int b)
		G1 = 1.0;
		else {
		get_shpfcn(s1, this->lattce_meam[a][a]);
		Gam1 = (s1[1] * t11av + s1[2] * t21av + s1[3] * t31av) / (Z1 * Z1);
		Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1);
		G_gam(Gam1, this->ibar_meam[a], &G1, &errorflag);
		}
		if (this->ibar_meam[b] <= 0)
		G2 = 1.0;
		else {
		get_shpfcn(s2, this->lattce_meam[b][b]);
		Gam2 = (s2[1] * t12av + s2[2] * t22av + s2[3] * t32av) / (Z2 * Z2);
		Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2);
		G_gam(Gam2, this->ibar_meam[b], &G2, &errorflag);
		}
		rho0_1 = this->rho0_meam[a] * Z1 * G1;
		@@ -585,7 +585,7 @@ void
		MEAM::compute_reference_density(void)
		{
		int a, Z, Z2, errorflag;
		double gam, Gbar, shp[3 + 1];
		double gam, Gbar, shp[3];
		double rho0, rho0_2nn, arat, scrn;

		// loop over element types
		@@ -595,8 +595,8 @@ MEAM::compute_reference_density(void)
		Gbar = 1.0;
		else {
		get_shpfcn(shp, this->lattce_meam[a][a]);
		gam = (this->t1_meam[a] * shp[1] + this->t2_meam[a] * shp[2] +
		this->t3_meam[a] * shp[3]) /
		gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] +
		this->t3_meam[a] * shp[2]) /
		(Z * Z);
		G_gam(gam, this->ibar_meam[a], &Gbar, &errorflag);
		}
		@@ -627,24 +627,24 @@ void
		MEAM::get_shpfcn(double* s /* s(3) */, lattice_t latt)
		{
		if (latt == FCC \|\| latt == BCC \|\| latt == B1 \|\| latt == B2) {
		s[0] = 0.0;
		s[1] = 0.0;
		s[2] = 0.0;
		s[3] = 0.0;
		} else if (latt == HCP) {
		s[0] = 0.0;
		s[1] = 0.0;
		s[2] = 0.0;
		s[3] = 1.0 / 3.0;
		s[2] = 1.0 / 3.0;
		} else if (latt == DIA) {
		s[0] = 0.0;
		s[1] = 0.0;
		s[2] = 0.0;
		s[3] = 32.0 / 9.0;
		s[2] = 32.0 / 9.0;
		} else if (latt == DIM) {
		s[1] = 1.0;
		s[2] = 2.0 / 3.0;
		s[0] = 1.0;
		s[1] = 2.0 / 3.0;
		// s(3) = 1.d0
		s[3] = 0.40;
		s[2] = 0.40;
		} else {
		s[1] = 0.0;
		s[0] = 0.0;
		// call error('Lattice not defined in get_shpfcn.')
		}
		}
		@@ -804,7 +804,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
		double* rho32)
		{
		double a1, a2;
		double s[3 + 1];
		double s[3];
		lattice_t lat;
		int Zij1nn, Zij2nn;
		double rhoa01nn, rhoa02nn;
		@@ -859,12 +859,12 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double*
		get_shpfcn(s, DIM);
		*rho01 = rhoa02;
		*rho02 = rhoa01;
		rho11 = s[1] rhoa12 * rhoa12;
		rho12 = s[1] rhoa11 * rhoa11;
		rho21 = s[2] rhoa22 * rhoa22;
		rho22 = s[2] rhoa21 * rhoa21;
		rho31 = s[3] rhoa32 * rhoa32;
		rho32 = s[3] rhoa31 * rhoa31;
		rho11 = s[0] rhoa12 * rhoa12;
		rho12 = s[0] rhoa11 * rhoa11;
		rho21 = s[1] rhoa22 * rhoa22;
		rho22 = s[1] rhoa21 * rhoa21;
		rho31 = s[2] rhoa32 * rhoa32;
		rho32 = s[2] rhoa31 * rhoa31;
		} else if (lat == C11) {
		*rho01 = rhoa01;
		*rho02 = rhoa02;

Admin message