diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index da608724e3..e6fa767f0d 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -50,193 +50,287 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int // Compute forces atom i elti = fmap[type[i]]; + if (elti < 0) return; - if (elti >= 0) { - xitmp = x[i][0]; - yitmp = x[i][1]; - zitmp = x[i][2]; + xitmp = x[i][0]; + yitmp = x[i][1]; + zitmp = x[i][2]; - // Treat each pair - for (jn = 0; jn < numneigh; jn++) { - j = firstneigh[jn]; - eltj = fmap[type[j]]; + // Treat each pair + for (jn = 0; jn < numneigh; jn++) { + j = firstneigh[jn]; + eltj = fmap[type[j]]; - if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) { + if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) { - sij = scrfcn[fnoffset + jn] * fcpair[fnoffset + jn]; - 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; + sij = scrfcn[fnoffset + jn] * fcpair[fnoffset + jn]; + 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; - // Compute phi and phip - ind = this->eltind[elti][eltj]; - pp = rij * this->rdrar; - kk = (int)pp; - kk = std::min(kk, this->nrar - 2); - pp = pp - kk; - pp = std::min(pp, 1.0); - phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + - this->phirar[ind][kk]; - phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk]; - recip = 1.0 / r; + // Compute phi and phip + ind = this->eltind[elti][eltj]; + pp = rij * this->rdrar; + kk = (int)pp; + kk = std::min(kk, this->nrar - 2); + pp = pp - kk; + pp = std::min(pp, 1.0); + phi = ((this->phirar3[ind][kk] * pp + this->phirar2[ind][kk]) * pp + this->phirar1[ind][kk]) * pp + + this->phirar[ind][kk]; + phip = (this->phirar6[ind][kk] * pp + this->phirar5[ind][kk]) * pp + this->phirar4[ind][kk]; + recip = 1.0 / r; - if (eflag_either != 0) { - if (eflag_global != 0) - *eng_vdwl = *eng_vdwl + phi * sij; - if (eflag_atom != 0) { - eatom[i] = eatom[i] + 0.5 * phi * sij; - eatom[j] = eatom[j] + 0.5 * phi * sij; + if (eflag_either != 0) { + if (eflag_global != 0) + *eng_vdwl = *eng_vdwl + phi * sij; + if (eflag_atom != 0) { + eatom[i] = eatom[i] + 0.5 * phi * sij; + eatom[j] = eatom[j] + 0.5 * phi * sij; + } + } + + // write(1,*) "force_meamf: phi: ",phi + // write(1,*) "force_meamf: phip: ",phip + + // Compute pair densities and derivatives + invrei = 1.0 / this->re_meam[elti][elti]; + ai = rij * invrei - 1.0; + ro0i = this->rho0_meam[elti]; + rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai); + drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i; + rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai); + drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i; + rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai); + drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i; + rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai); + drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i; + + if (elti != eltj) { + invrej = 1.0 / this->re_meam[eltj][eltj]; + aj = rij * invrej - 1.0; + ro0j = this->rho0_meam[eltj]; + rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj); + drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j; + rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj); + drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j; + rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj); + drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j; + rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj); + drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j; + } else { + rhoa0j = rhoa0i; + drhoa0j = drhoa0i; + rhoa1j = rhoa1i; + drhoa1j = drhoa1i; + rhoa2j = rhoa2i; + drhoa2j = drhoa2i; + rhoa3j = rhoa3i; + drhoa3j = drhoa3i; + } + + const double t1mi = this->t1_meam[elti]; + const double t2mi = this->t2_meam[elti]; + const double t3mi = this->t3_meam[elti]; + const double t1mj = this->t1_meam[eltj]; + const double t2mj = this->t2_meam[eltj]; + const double t3mj = this->t3_meam[eltj]; + + if (this->ialloy == 1) { + rhoa1j *= t1mj; + rhoa2j *= t2mj; + rhoa3j *= t3mj; + rhoa1i *= t1mi; + rhoa2i *= t2mi; + rhoa3i *= t3mi; + drhoa1j *= t1mj; + drhoa2j *= t2mj; + drhoa3j *= t3mj; + drhoa1i *= t1mi; + drhoa2i *= t2mi; + drhoa3i *= t3mi; + } + + nv2 = 0; + nv3 = 0; + arg1i1 = 0.0; + arg1j1 = 0.0; + arg1i2 = 0.0; + arg1j2 = 0.0; + arg1i3 = 0.0; + arg1j3 = 0.0; + arg3i3 = 0.0; + arg3j3 = 0.0; + 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 + arho3[i][nv3] * arg; + arg1j3 = arg1j3 - arho3[j][nv3] * arg; + nv3 = nv3 + 1; } + arg = delij[n] * delij[p] * this->v2D[nv2]; + arg1i2 = arg1i2 + arho2[i][nv2] * arg; + arg1j2 = arg1j2 + arho2[j][nv2] * arg; + nv2 = nv2 + 1; } + 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]; + } - // write(1,*) "force_meamf: phi: ",phi - // write(1,*) "force_meamf: phip: ",phip + // rho0 terms + drho0dr1 = drhoa0j * sij; + drho0dr2 = drhoa0i * sij; - // Compute pair densities and derivatives - invrei = 1.0 / this->re_meam[elti][elti]; - ai = rij * invrei - 1.0; - ro0i = this->rho0_meam[elti]; - rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai); - drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i; - rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai); - drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i; - rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai); - drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i; - rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai); - drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i; + // rho1 terms + a1 = 2 * sij / rij; + drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1; + drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1; + a1 = 2.0 * sij / rij; + for (m = 0; m < 3; m++) { + drho1drm1[m] = a1 * rhoa1j * arho1[i][m]; + drho1drm2[m] = -a1 * rhoa1i * arho1[j][m]; + } - if (elti != eltj) { - invrej = 1.0 / this->re_meam[eltj][eltj]; - aj = rij * invrej - 1.0; - ro0j = this->rho0_meam[eltj]; - rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj); - drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j; - rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj); - drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j; - rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj); - drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j; - rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj); - drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j; - } else { - rhoa0j = rhoa0i; - drhoa0j = drhoa0i; - rhoa1j = rhoa1i; - drhoa1j = drhoa1i; - rhoa2j = rhoa2i; - drhoa2j = drhoa2i; - rhoa3j = rhoa3i; - drhoa3j = drhoa3i; - } - - if (this->ialloy == 1) { - rhoa1j = rhoa1j * this->t1_meam[eltj]; - rhoa2j = rhoa2j * this->t2_meam[eltj]; - rhoa3j = rhoa3j * this->t3_meam[eltj]; - rhoa1i = rhoa1i * this->t1_meam[elti]; - rhoa2i = rhoa2i * this->t2_meam[elti]; - rhoa3i = rhoa3i * this->t3_meam[elti]; - drhoa1j = drhoa1j * this->t1_meam[eltj]; - drhoa2j = drhoa2j * this->t2_meam[eltj]; - drhoa3j = drhoa3j * this->t3_meam[eltj]; - drhoa1i = drhoa1i * this->t1_meam[elti]; - drhoa2i = drhoa2i * this->t2_meam[elti]; - drhoa3i = drhoa3i * this->t3_meam[elti]; + // rho2 terms + a2 = 2 * sij / rij2; + drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - 2.0 / 3.0 * arho2b[i] * drhoa2j * sij; + drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - 2.0 / 3.0 * arho2b[j] * drhoa2i * sij; + a2 = 4 * sij / rij2; + for (m = 0; m < 3; m++) { + drho2drm1[m] = 0.0; + drho2drm2[m] = 0.0; + 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]; + } + // rho3 terms + rij3 = rij * rij2; + a3 = 2 * sij / rij3; + a3a = 6.0 / 5.0 * sij / rij; + drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - a3a * (drhoa3j - rhoa3j / rij) * arg3i3; + drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - a3a * (drhoa3i - rhoa3i / rij) * arg3j3; + a3 = 6 * sij / rij3; + a3a = 6 * sij / (5 * rij); + for (m = 0; m < 3; m++) { + drho3drm1[m] = 0.0; + drho3drm2[m] = 0.0; nv2 = 0; - nv3 = 0; - arg1i1 = 0.0; - arg1j1 = 0.0; - arg1i2 = 0.0; - arg1j2 = 0.0; - arg1i3 = 0.0; - arg1j3 = 0.0; - arg3i3 = 0.0; - arg3j3 = 0.0; 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 + arho3[i][nv3] * arg; - arg1j3 = arg1j3 - arho3[j][nv3] * arg; - nv3 = nv3 + 1; - } arg = delij[n] * delij[p] * this->v2D[nv2]; - arg1i2 = arg1i2 + arho2[i][nv2] * arg; - arg1j2 = arg1j2 + arho2[j][nv2] * 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; } - 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]; } + drho3drm1[m] = (a3 * drho3drm1[m] - a3a * arho3b[i][m]) * rhoa3j; + drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i; + } - // rho0 terms - drho0dr1 = drhoa0j * sij; - drho0dr2 = drhoa0i * sij; + // Compute derivatives of weighting functions t wrt rij + t1i = t_ave[i][0]; + t2i = t_ave[i][1]; + t3i = t_ave[i][2]; + t1j = t_ave[j][0]; + t2j = t_ave[j][1]; + t3j = t_ave[j][2]; - // rho1 terms - a1 = 2 * sij / rij; - drho1dr1 = a1 * (drhoa1j - rhoa1j / rij) * arg1i1; - drho1dr2 = a1 * (drhoa1i - rhoa1i / rij) * arg1j1; - a1 = 2.0 * sij / rij; - for (m = 0; m < 3; m++) { - drho1drm1[m] = a1 * rhoa1j * arho1[i][m]; - drho1drm2[m] = -a1 * rhoa1i * arho1[j][m]; - } + if (this->ialloy == 1) { - // rho2 terms - a2 = 2 * sij / rij2; - drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - 2.0 / 3.0 * arho2b[i] * drhoa2j * sij; - drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - 2.0 / 3.0 * arho2b[j] * drhoa2i * sij; - a2 = 4 * sij / rij2; - for (m = 0; m < 3; m++) { - drho2drm1[m] = 0.0; - drho2drm2[m] = 0.0; - 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]; - } + a1i = 0.0; + a1j = 0.0; + a2i = 0.0; + a2j = 0.0; + a3i = 0.0; + a3j = 0.0; + if (!iszero(tsq_ave[i][0])) + a1i = drhoa0j * sij / tsq_ave[i][0]; + if (!iszero(tsq_ave[j][0])) + a1j = drhoa0i * sij / tsq_ave[j][0]; + if (!iszero(tsq_ave[i][1])) + a2i = drhoa0j * sij / tsq_ave[i][1]; + if (!iszero(tsq_ave[j][1])) + a2j = drhoa0i * sij / tsq_ave[j][1]; + if (!iszero(tsq_ave[i][2])) + a3i = drhoa0j * sij / tsq_ave[i][2]; + if (!iszero(tsq_ave[j][2])) + a3j = drhoa0i * sij / tsq_ave[j][2]; - // rho3 terms - rij3 = rij * rij2; - a3 = 2 * sij / rij3; - a3a = 6.0 / 5.0 * sij / rij; - drho3dr1 = a3 * (drhoa3j - 3 * rhoa3j / rij) * arg1i3 - a3a * (drhoa3j - rhoa3j / rij) * arg3i3; - drho3dr2 = a3 * (drhoa3i - 3 * rhoa3i / rij) * arg1j3 - a3a * (drhoa3i - rhoa3i / rij) * arg3j3; - a3 = 6 * sij / rij3; - a3a = 6 * sij / (5 * rij); - for (m = 0; m < 3; m++) { - drho3drm1[m] = 0.0; - drho3drm2[m] = 0.0; - 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] + 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 * arho3b[i][m]) * rhoa3j; - drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * arho3b[j][m]) * rhoa3i; - } + dt1dr1 = a1i * (t1mj - t1i * t1mj*t1mj); + dt1dr2 = a1j * (t1mi - t1j * t1mi*t1mi); + dt2dr1 = a2i * (t2mj - t2i * t2mj*t2mj); + dt2dr2 = a2j * (t2mi - t2j * t2mi*t2mi); + dt3dr1 = a3i * (t3mj - t3i * t3mj*t3mj); + dt3dr2 = a3j * (t3mi - t3j * t3mi*t3mi); - // Compute derivatives of weighting functions t wrt rij - t1i = t_ave[i][0]; - t2i = t_ave[i][1]; - t3i = t_ave[i][2]; - t1j = t_ave[j][0]; - t2j = t_ave[j][1]; - t3j = t_ave[j][2]; + } else if (this->ialloy == 2) { + + dt1dr1 = 0.0; + dt1dr2 = 0.0; + dt2dr1 = 0.0; + dt2dr2 = 0.0; + dt3dr1 = 0.0; + dt3dr2 = 0.0; + + } else { + + ai = 0.0; + if (!iszero(rho0[i])) + ai = drhoa0j * sij / rho0[i]; + aj = 0.0; + if (!iszero(rho0[j])) + aj = drhoa0i * sij / rho0[j]; + + dt1dr1 = ai * (t1mj - t1i); + dt1dr2 = aj * (t1mi - t1j); + dt2dr1 = ai * (t2mj - t2i); + dt2dr2 = aj * (t2mi - t2j); + dt3dr1 = ai * (t3mj - t3i); + dt3dr2 = aj * (t3mi - t3j); + } + + // Compute derivatives of total density wrt rij, sij and rij(3) + get_shpfcn(this->lattce_meam[elti][elti], shpi); + get_shpfcn(this->lattce_meam[eltj][eltj], shpj); + drhodr1 = dgamma1[i] * drho0dr1 + + dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 + + dt3dr1 * rho3[i] + t3i * drho3dr1) - + dgamma3[i] * (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[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2); + for (m = 0; m < 3; m++) { + drhodrm1[m] = 0.0; + drhodrm2[m] = 0.0; + drhodrm1[m] = dgamma2[i] * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]); + drhodrm2[m] = dgamma2[j] * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]); + } + + // Compute derivatives wrt sij, but only if necessary + if (!iszero(dscrfcn[fnoffset + jn])) { + drho0ds1 = rhoa0j; + drho0ds2 = rhoa0i; + a1 = 2.0 / rij; + drho1ds1 = a1 * rhoa1j * arg1i1; + drho1ds2 = a1 * rhoa1i * arg1j1; + a2 = 2.0 / rij2; + drho2ds1 = a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j; + drho2ds2 = a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arho2b[j] * rhoa2i; + a3 = 2.0 / rij3; + a3a = 6.0 / (5.0 * rij); + drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3; + drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3; if (this->ialloy == 1) { @@ -247,281 +341,193 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int a3i = 0.0; a3j = 0.0; if (!iszero(tsq_ave[i][0])) - a1i = drhoa0j * sij / tsq_ave[i][0]; + a1i = rhoa0j / tsq_ave[i][0]; if (!iszero(tsq_ave[j][0])) - a1j = drhoa0i * sij / tsq_ave[j][0]; + a1j = rhoa0i / tsq_ave[j][0]; if (!iszero(tsq_ave[i][1])) - a2i = drhoa0j * sij / tsq_ave[i][1]; + a2i = rhoa0j / tsq_ave[i][1]; if (!iszero(tsq_ave[j][1])) - a2j = drhoa0i * sij / tsq_ave[j][1]; + a2j = rhoa0i / tsq_ave[j][1]; if (!iszero(tsq_ave[i][2])) - a3i = drhoa0j * sij / tsq_ave[i][2]; + a3i = rhoa0j / tsq_ave[i][2]; if (!iszero(tsq_ave[j][2])) - a3j = drhoa0i * sij / tsq_ave[j][2]; + a3j = rhoa0i / tsq_ave[j][2]; - dt1dr1 = a1i * (this->t1_meam[eltj] - t1i * pow(this->t1_meam[eltj], 2)); - dt1dr2 = a1j * (this->t1_meam[elti] - t1j * pow(this->t1_meam[elti], 2)); - dt2dr1 = a2i * (this->t2_meam[eltj] - t2i * pow(this->t2_meam[eltj], 2)); - dt2dr2 = a2j * (this->t2_meam[elti] - t2j * pow(this->t2_meam[elti], 2)); - dt3dr1 = a3i * (this->t3_meam[eltj] - t3i * pow(this->t3_meam[eltj], 2)); - dt3dr2 = a3j * (this->t3_meam[elti] - t3j * pow(this->t3_meam[elti], 2)); + dt1ds1 = a1i * (t1mj - t1i * pow(t1mj, 2)); + dt1ds2 = a1j * (t1mi - t1j * pow(t1mi, 2)); + dt2ds1 = a2i * (t2mj - t2i * pow(t2mj, 2)); + dt2ds2 = a2j * (t2mi - t2j * pow(t2mi, 2)); + dt3ds1 = a3i * (t3mj - t3i * pow(t3mj, 2)); + dt3ds2 = a3j * (t3mi - t3j * pow(t3mi, 2)); } else if (this->ialloy == 2) { - dt1dr1 = 0.0; - dt1dr2 = 0.0; - dt2dr1 = 0.0; - dt2dr2 = 0.0; - dt3dr1 = 0.0; - dt3dr2 = 0.0; + dt1ds1 = 0.0; + dt1ds2 = 0.0; + dt2ds1 = 0.0; + dt2ds2 = 0.0; + dt3ds1 = 0.0; + dt3ds2 = 0.0; } else { ai = 0.0; if (!iszero(rho0[i])) - ai = drhoa0j * sij / rho0[i]; + ai = rhoa0j / rho0[i]; aj = 0.0; if (!iszero(rho0[j])) - aj = drhoa0i * sij / rho0[j]; + aj = rhoa0i / rho0[j]; - dt1dr1 = ai * (this->t1_meam[eltj] - t1i); - dt1dr2 = aj * (this->t1_meam[elti] - t1j); - dt2dr1 = ai * (this->t2_meam[eltj] - t2i); - dt2dr2 = aj * (this->t2_meam[elti] - t2j); - dt3dr1 = ai * (this->t3_meam[eltj] - t3i); - dt3dr2 = aj * (this->t3_meam[elti] - t3j); + dt1ds1 = ai * (t1mj - t1i); + dt1ds2 = aj * (t1mi - t1j); + dt2ds1 = ai * (t2mj - t2i); + dt2ds2 = aj * (t2mi - t2j); + dt3ds1 = ai * (t3mj - t3i); + dt3ds2 = aj * (t3mi - t3j); } - // Compute derivatives of total density wrt rij, sij and rij(3) - get_shpfcn(this->lattce_meam[elti][elti], shpi); - get_shpfcn(this->lattce_meam[eltj][eltj], shpj); - drhodr1 = dgamma1[i] * drho0dr1 + - dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 + - dt3dr1 * rho3[i] + t3i * drho3dr1) - - dgamma3[i] * (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[0] * dt1dr2 + shpj[1] * dt2dr2 + shpj[2] * dt3dr2); - for (m = 0; m < 3; m++) { - drhodrm1[m] = 0.0; - drhodrm2[m] = 0.0; - drhodrm1[m] = dgamma2[i] * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]); - drhodrm2[m] = dgamma2[j] * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]); + drhods1 = dgamma1[i] * drho0ds1 + + dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 + + dt3ds1 * rho3[i] + t3i * drho3ds1) - + dgamma3[i] * (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[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2); + } + + // Compute derivatives of energy wrt rij, sij and rij[3] + dUdrij = phip * sij + frhop[i] * drhodr1 + frhop[j] * drhodr2; + dUdsij = 0.0; + if (!iszero(dscrfcn[fnoffset + jn])) { + dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2; + } + for (m = 0; m < 3; m++) { + dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m]; + } + + // Add the part of the force due to dUdrij and dUdsij + + force = dUdrij * recip + dUdsij * dscrfcn[fnoffset + jn]; + for (m = 0; m < 3; m++) { + forcem = delij[m] * force + dUdrijm[m]; + 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]; + 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.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]; } + } - // Compute derivatives wrt sij, but only if necessary - if (!iszero(dscrfcn[fnoffset + jn])) { - drho0ds1 = rhoa0j; - drho0ds2 = rhoa0i; - a1 = 2.0 / rij; - drho1ds1 = a1 * rhoa1j * arg1i1; - drho1ds2 = a1 * rhoa1i * arg1j1; - a2 = 2.0 / rij2; - drho2ds1 = a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j; - drho2ds2 = a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arho2b[j] * rhoa2i; - a3 = 2.0 / rij3; - a3a = 6.0 / (5.0 * rij); - drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3; - drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3; + // Now compute forces on other atoms k due to change in sij - if (this->ialloy == 1) { + if (iszero(sij) || iszero(sij - 1.0)) continue; //: cont jn loop - a1i = 0.0; - a1j = 0.0; - a2i = 0.0; - a2j = 0.0; - a3i = 0.0; - a3j = 0.0; - if (!iszero(tsq_ave[i][0])) - a1i = rhoa0j / tsq_ave[i][0]; - if (!iszero(tsq_ave[j][0])) - a1j = rhoa0i / tsq_ave[j][0]; - if (!iszero(tsq_ave[i][1])) - a2i = rhoa0j / tsq_ave[i][1]; - if (!iszero(tsq_ave[j][1])) - a2j = rhoa0i / tsq_ave[j][1]; - if (!iszero(tsq_ave[i][2])) - a3i = rhoa0j / tsq_ave[i][2]; - if (!iszero(tsq_ave[j][2])) - a3j = rhoa0i / tsq_ave[j][2]; + double dxik(0), dyik(0), dzik(0); + double dxjk(0), dyjk(0), dzjk(0); - dt1ds1 = a1i * (this->t1_meam[eltj] - t1i * pow(this->t1_meam[eltj], 2)); - dt1ds2 = a1j * (this->t1_meam[elti] - t1j * pow(this->t1_meam[elti], 2)); - dt2ds1 = a2i * (this->t2_meam[eltj] - t2i * pow(this->t2_meam[eltj], 2)); - dt2ds2 = a2j * (this->t2_meam[elti] - t2j * pow(this->t2_meam[elti], 2)); - dt3ds1 = a3i * (this->t3_meam[eltj] - t3i * pow(this->t3_meam[eltj], 2)); - dt3ds2 = a3j * (this->t3_meam[elti] - t3j * pow(this->t3_meam[elti], 2)); + for (kn = 0; kn < numneigh_full; kn++) { + k = firstneigh_full[kn]; + eltk = fmap[type[k]]; + if (k != j && eltk >= 0) { + double xik, xjk, cikj, sikj, dfc, a; + double dCikj1, dCikj2; + double delc, rik2, rjk2; - } else if (this->ialloy == 2) { + 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]; - dt1ds1 = 0.0; - dt1ds2 = 0.0; - dt2ds1 = 0.0; - dt2ds2 = 0.0; - dt3ds1 = 0.0; - dt3ds2 = 0.0; - - } else { - - ai = 0.0; - if (!iszero(rho0[i])) - ai = rhoa0j / rho0[i]; - aj = 0.0; - if (!iszero(rho0[j])) - aj = rhoa0i / rho0[j]; - - dt1ds1 = ai * (this->t1_meam[eltj] - t1i); - dt1ds2 = aj * (this->t1_meam[elti] - t1j); - dt2ds1 = ai * (this->t2_meam[eltj] - t2i); - dt2ds2 = aj * (this->t2_meam[elti] - t2j); - dt3ds1 = ai * (this->t3_meam[eltj] - t3i); - dt3ds2 = aj * (this->t3_meam[elti] - t3j); - } - - drhods1 = dgamma1[i] * drho0ds1 + - dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 + - dt3ds1 * rho3[i] + t3i * drho3ds1) - - dgamma3[i] * (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[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2); - } - - // Compute derivatives of energy wrt rij, sij and rij[3] - dUdrij = phip * sij + frhop[i] * drhodr1 + frhop[j] * drhodr2; - dUdsij = 0.0; - if (!iszero(dscrfcn[fnoffset + jn])) { - dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2; - } - for (m = 0; m < 3; m++) { - dUdrijm[m] = frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m]; - } - - // Add the part of the force due to dUdrij and dUdsij - - force = dUdrij * recip + dUdsij * dscrfcn[fnoffset + jn]; - for (m = 0; m < 3; m++) { - forcem = delij[m] * force + dUdrijm[m]; - 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]; - 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.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 - - if (iszero(sij) || iszero(sij - 1.0)) - continue; //: cont jn loop - for (kn = 0; kn < numneigh_full; kn++) { - k = firstneigh_full[kn]; - eltk = fmap[type[k]]; - if (k != j && eltk >= 0) { - 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; - } + 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; - for (m = 0; m < 3; m++) { - delik[m] = x[k][m] - x[i][m]; - deljk[m] = x[k][m] - x[j][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]; - } + if (!iszero(dsij1) || !iszero(dsij2)) { + force1 = dUdsij * dsij1; + force2 = dUdsij * dsij2; - // Tabulate per-atom virial as symmetrized stress tensor + f[i][0] += force1 * dxik; + f[i][1] += force1 * dyik; + f[i][2] += force1 * dzik; + f[j][0] += force2 * dxjk; + f[j][1] += force2 * dyjk; + f[j][2] += force2 * dzjk; + f[k][0] -= force1 * dxik + force2 * dxjk; + f[k][1] -= force1 * dyik + force2 * dyjk; + f[k][2] -= force1 * dzik + force2 * dzjk; - if (vflag_atom != 0) { - fi[0] = force1 * delik[0]; - fi[1] = force1 * delik[1]; - fi[2] = force1 * delik[2]; - fj[0] = force2 * deljk[0]; - fj[1] = force2 * deljk[1]; - fj[2] = force2 * deljk[2]; - 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] = -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]); + // Tabulate per-atom virial as symmetrized stress tensor - 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]; - } + if (vflag_atom != 0) { + fi[0] = force1 * dxik; + fi[1] = force1 * dyik; + fi[2] = force1 * dzik; + fj[0] = force2 * dxjk; + fj[1] = force2 * dyjk; + fj[2] = force2 * dzjk; + v[0] = -third * (dxik * fi[0] + dxjk * fj[0]); + v[1] = -third * (dyik * fi[1] + dyjk * fj[1]); + v[2] = -third * (dzik * fi[2] + dzjk * fj[2]); + v[3] = -sixth * (dxik * fi[1] + dxjk * fj[1] + dyik * fi[0] + dyjk * fj[0]); + v[4] = -sixth * (dxik * fi[2] + dxjk * fj[2] + dzik * fi[0] + dzjk * fj[0]); + v[5] = -sixth * (dyik * fi[2] + dyjk * fj[2] + dzik * fi[1] + dzjk * fj[1]); + + 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]; } } } - // end of k loop } + // end of k loop } } - // end of j loop } - - // else if elti<0, this is not a meam atom + // end of j loop } }