From dadd1c8b4d6e2e277666bb3f36715352f375f435 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20H=C3=BCtter?= Date: Thu, 22 Jun 2017 19:02:14 +0200 Subject: [PATCH] Remove neigh_f2c/c2f, related cleanup - neighbour lists now use C indexing - removed many arr*v() macros - removed some unneccessary pointers - minor reformatting --- src/USER-MEAMC/meam.h | 28 +-- src/USER-MEAMC/meam_dens_final.cpp | 108 +++++------ src/USER-MEAMC/meam_dens_init.cpp | 192 +++++++++--------- src/USER-MEAMC/meam_force.cpp | 300 ++++++++++++++--------------- src/USER-MEAMC/pair_meamc.cpp | 62 +----- src/USER-MEAMC/pair_meamc.h | 2 - 6 files changed, 310 insertions(+), 382 deletions(-) diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 2ac487483e..7f210eac18 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -141,20 +141,24 @@ public: void interpolate_meam(int); double compute_phi(double, int, int); public: - void meam_setup_global(int, lattice_t*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); - void meam_setup_param(int, double, int, int*, int*); - void meam_setup_done(double*); - void meam_dens_setup(int, int, int); - void meam_dens_init(int* i, int* ntype, int* type, int* fmap, double** x, - int* numneigh, int* firstneigh, int* numneigh_full, + void meam_setup_global(int nelt, lattice_t* lat, double* z, int* ielement, double* atwt, + double* alpha, double* b0, double* b1, double* b2, + double* b3, double* alat, double* esub, double* asub, + double* t0, double* t1, double* t2, double* t3, + double* rozero, int* ibar); + void meam_setup_param(int which, double value, int nindex, int* index /*index(3)*/, int* errorflag); + void meam_setup_done(double* cutmax); + void meam_dens_setup(int atom_nmax, int nall, int n_neigh); + void meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, + int numneigh, int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset, int* errorflag); - void meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, - int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype, + void meam_dens_final(int nlocal, int eflag_either, int eflag_global, + int eflag_atom, double* eng_vdwl, double* eatom, int ntype, int* type, int* fmap, int* errorflag); - void meam_force(int* iptr, int* eflag_either, int* eflag_global, - int* eflag_atom, int* vflag_atom, double* eng_vdwl, double* eatom, - int* ntype, int* type, int* fmap, double** x, int* numneigh, - int* firstneigh, int* numneigh_full, int* firstneigh_full, + void meam_force(int i, int eflag_either, int eflag_global, + int eflag_atom, int vflag_atom, double* eng_vdwl, double* eatom, + int ntype, int* type, int* fmap, double** x, int numneigh, + int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom, int* errorflag); }; diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index 716cd8cada..b753542db2 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -21,8 +21,8 @@ using namespace LAMMPS_NS; // void -MEAM::meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, - int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype, +MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, + int eflag_atom, double* eng_vdwl, double* eatom, int ntype, int* type, int* fmap, int* errorflag) { int i, elti; @@ -32,56 +32,54 @@ MEAM::meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, // Complete the calculation of density - for (i = 1; i <= *nlocal; i++) { - elti = fmap[arr1v(type, i)]; + for (i = 0; i < nlocal; i++) { + elti = fmap[type[i]]; if (elti >= 0) { - arr1v(rho1, i) = 0.0; - arr1v(rho2, i) = -1.0 / 3.0 * arr1v(arho2b, i) * arr1v(arho2b, i); - arr1v(rho3, i) = 0.0; + rho1[i] = 0.0; + rho2[i] = -1.0 / 3.0 * arho2b[i] * arho2b[i]; + rho3[i] = 0.0; for (m = 1; m <= 3; m++) { - arr1v(rho1, i) = - arr1v(rho1, i) + arr2v(arho1, m, i) * arr2v(arho1, m, i); - arr1v(rho3, i) = arr1v(rho3, i) - - 3.0 / 5.0 * arr2v(arho3b, m, i) * arr2v(arho3b, m, i); + 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 = 1; m <= 6; m++) { - arr1v(rho2, i) = - arr1v(rho2, i) + - this->v2D[m] * arr2v(arho2, m, i) * arr2v(arho2, m, i); + rho2[i] = + rho2[i] + + this->v2D[m] * arr2v(arho2, m, i+1) * arr2v(arho2, m, i+1); } for (m = 1; m <= 10; m++) { - arr1v(rho3, i) = - arr1v(rho3, i) + - this->v3D[m] * arr2v(arho3, m, i) * arr2v(arho3, m, i); + rho3[i] = + rho3[i] + + this->v3D[m] * arr2v(arho3, m, i+1) * arr2v(arho3, m, i+1); } - if (arr1v(rho0, i) > 0.0) { + if (rho0[i] > 0.0) { if (this->ialloy == 1) { - arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) / arr2v(tsq_ave, 1, i); - arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) / arr2v(tsq_ave, 2, i); - arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr2v(tsq_ave, 3, i); + t_ave[i][0] = t_ave[i][0] / tsq_ave[i][0]; + t_ave[i][1] = t_ave[i][1] / tsq_ave[i][1]; + t_ave[i][2] = t_ave[i][2] / tsq_ave[i][2]; } else if (this->ialloy == 2) { - arr2v(t_ave, 1, i) = this->t1_meam[elti]; - arr2v(t_ave, 2, i) = this->t2_meam[elti]; - arr2v(t_ave, 3, i) = this->t3_meam[elti]; + t_ave[i][0] = this->t1_meam[elti]; + t_ave[i][1] = this->t2_meam[elti]; + t_ave[i][2] = this->t3_meam[elti]; } else { - arr2v(t_ave, 1, i) = arr2v(t_ave, 1, i) / arr1v(rho0, i); - arr2v(t_ave, 2, i) = arr2v(t_ave, 2, i) / arr1v(rho0, i); - arr2v(t_ave, 3, i) = arr2v(t_ave, 3, i) / arr1v(rho0, i); + t_ave[i][0] = t_ave[i][0] / rho0[i]; + t_ave[i][1] = t_ave[i][1] / rho0[i]; + t_ave[i][2] = t_ave[i][2] / rho0[i]; } } - arr1v(gamma, i) = arr2v(t_ave, 1, i) * arr1v(rho1, i) + - arr2v(t_ave, 2, i) * arr1v(rho2, i) + - arr2v(t_ave, 3, i) * arr1v(rho3, i); + gamma[i] = t_ave[i][0] * rho1[i] + + t_ave[i][1] * rho2[i] + + t_ave[i][2] * rho3[i]; - if (arr1v(rho0, i) > 0.0) { - arr1v(gamma, i) = arr1v(gamma, i) / (arr1v(rho0, i) * arr1v(rho0, i)); + if (rho0[i] > 0.0) { + gamma[i] = gamma[i] / (rho0[i] * rho0[i]); } Z = this->Z_meam[elti]; - G_gam(arr1v(gamma, i), this->ibar_meam[elti], &G, errorflag); + G_gam(gamma[i], this->ibar_meam[elti], &G, errorflag); if (*errorflag != 0) return; get_shpfcn(shp, this->lattce_meam[elti][elti]); @@ -90,8 +88,8 @@ MEAM::meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, dGbar = 0.0; } else { if (this->mix_ref_t == 1) { - gam = (arr2v(t_ave, 1, i) * shp[1] + arr2v(t_ave, 2, i) * shp[2] + - arr2v(t_ave, 3, i) * shp[3]) / + gam = (t_ave[i][0] * shp[1] + t_ave[i][1] * shp[2] + + t_ave[i][2] * shp[3]) / (Z * Z); } else { gam = (this->t1_meam[elti] * shp[1] + @@ -101,15 +99,15 @@ MEAM::meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, } G_gam(gam, this->ibar_meam[elti], &Gbar, errorflag); } - arr1v(rho, i) = arr1v(rho0, i) * G; + rho[i] = rho0[i] * G; if (this->mix_ref_t == 1) { if (this->ibar_meam[elti] <= 0) { Gbar = 1.0; dGbar = 0.0; } else { - gam = (arr2v(t_ave, 1, i) * shp[1] + arr2v(t_ave, 2, i) * shp[2] + - arr2v(t_ave, 3, i) * shp[3]) / + gam = (t_ave[i][0] * shp[1] + t_ave[i][1] * shp[2] + + t_ave[i][2] * shp[3]) / (Z * Z); dG_gam(gam, this->ibar_meam[elti], &Gbar, &dGbar); } @@ -121,57 +119,57 @@ MEAM::meam_dens_final(int* nlocal, int* eflag_either, int* eflag_global, rho_bkgd = this->rho_ref_meam[elti]; } } - rhob = arr1v(rho, i) / rho_bkgd; + rhob = rho[i] / rho_bkgd; denom = 1.0 / rho_bkgd; - dG_gam(arr1v(gamma, i), this->ibar_meam[elti], &G, &dG); + dG_gam(gamma[i], this->ibar_meam[elti], &G, &dG); - arr1v(dgamma1, i) = (G - 2 * dG * arr1v(gamma, i)) * denom; + dgamma1[i] = (G - 2 * dG * gamma[i]) * denom; - if (!iszero(arr1v(rho0, i))) { - arr1v(dgamma2, i) = (dG / arr1v(rho0, i)) * denom; + if (!iszero(rho0[i])) { + dgamma2[i] = (dG / rho0[i]) * denom; } else { - arr1v(dgamma2, i) = 0.0; + dgamma2[i] = 0.0; } // dgamma3 is nonzero only if we are using the "mixed" rule for // computing t in the reference system (which is not correct, but // included for backward compatibility if (this->mix_ref_t == 1) { - arr1v(dgamma3, i) = arr1v(rho0, i) * G * dGbar / (Gbar * Z * Z) * denom; + dgamma3[i] = rho0[i] * G * dGbar / (Gbar * Z * Z) * denom; } else { - arr1v(dgamma3, i) = 0.0; + dgamma3[i] = 0.0; } B = this->A_meam[elti] * this->Ec_meam[elti][elti]; if (!iszero(rhob)) { if (this->emb_lin_neg == 1 && rhob <= 0) { - arr1v(frhop, i) = -B; + frhop[i] = -B; } else { - arr1v(frhop, i) = B * (log(rhob) + 1.0); + frhop[i] = B * (log(rhob) + 1.0); } - if (*eflag_either != 0) { - if (*eflag_global != 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); } } - if (*eflag_atom != 0) { + if (eflag_atom != 0) { if (this->emb_lin_neg == 1 && rhob <= 0) { - arr1v(eatom, i) = arr1v(eatom, i) - B * rhob; + eatom[i] = eatom[i] - B * rhob; } else { - arr1v(eatom, i) = arr1v(eatom, i) + B * rhob * log(rhob); + eatom[i] = eatom[i] + B * rhob * log(rhob); } } } } else { if (this->emb_lin_neg == 1) { - arr1v(frhop, i) = -B; + frhop[i] = -B; } else { - arr1v(frhop, i) = B; + frhop[i] = B; } } } diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index e530a5feb1..932346a19a 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -101,18 +101,18 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh) // void -MEAM::meam_dens_init(int* i, int* ntype, int* type, int* fmap, double** x, - int* numneigh, int* firstneigh, int* numneigh_full, +MEAM::meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, + int numneigh, int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset, int* errorflag) { *errorflag = 0; // Compute screening function and derivatives - getscreen(*i, &scrfcn[fnoffset], &dscrfcn[fnoffset], &fcpair[fnoffset], x, *numneigh, firstneigh, - *numneigh_full, firstneigh_full, *ntype, type, fmap); + getscreen(i, &scrfcn[fnoffset], &dscrfcn[fnoffset], &fcpair[fnoffset], x, numneigh, firstneigh, + numneigh_full, firstneigh_full, ntype, type, fmap); // Calculate intermediate density terms to be communicated - calc_rho1(*i, *ntype, type, fmap, x, *numneigh, firstneigh, &scrfcn[fnoffset], &fcpair[fnoffset]); + calc_rho1(i, ntype, type, fmap, x, numneigh, firstneigh, &scrfcn[fnoffset], &fcpair[fnoffset]); } // ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -135,24 +135,24 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double rnorm, fc, dfc, drinv; drinv = 1.0 / this->delr_meam; - elti = fmap[arr1v(type, i)]; + elti = fmap[type[i]]; if (elti >= 0) { - xitmp = arr2v(x, 1, i); - yitmp = arr2v(x, 2, i); - zitmp = arr2v(x, 3, i); + xitmp = x[i][0]; + yitmp = x[i][1]; + zitmp = x[i][2]; - for (jn = 1; jn <= numneigh; jn++) { - j = arr1v(firstneigh, jn); + for (jn = 0; jn < numneigh; jn++) { + j = firstneigh[jn]; - eltj = fmap[arr1v(type, j)]; + eltj = fmap[type[j]]; if (eltj >= 0) { // First compute screening function itself, sij - xjtmp = arr2v(x, 1, j); - yjtmp = arr2v(x, 2, j); - zjtmp = arr2v(x, 3, j); + xjtmp = x[j][0]; + yjtmp = x[j][1]; + zjtmp = x[j][2]; delxij = xjtmp - xitmp; delyij = yjtmp - yitmp; delzij = zjtmp - zitmp; @@ -171,21 +171,21 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, } // Now compute derivatives - arr1v(dscrfcn, jn) = 0.0; + dscrfcn[jn] = 0.0; sfcij = sij * fcij; if (iszero(sfcij) || iszero(sfcij - 1.0)) goto LABEL_100; rbound = this->ebound_meam[elti][eltj] * rij2; - for (kn = 1; kn <= numneigh_full; kn++) { - k = arr1v(firstneigh_full, kn); + for (kn = 0; kn < numneigh_full; kn++) { + k = firstneigh_full[kn]; if (k == j) continue; - eltk = fmap[arr1v(type, k)]; + eltk = fmap[type[k]]; if (eltk < 0) continue; - xktmp = arr2v(x, 1, k); - yktmp = arr2v(x, 2, k); - zktmp = arr2v(x, 3, k); + xktmp = x[k][0]; + yktmp = x[k][1]; + zktmp = x[k][2]; delxjk = xktmp - xjtmp; delyjk = yktmp - yjtmp; delzjk = zktmp - zjtmp; @@ -223,16 +223,16 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, dfcut(cikj, &sikj, &dfikj); coef1 = dfikj / (delc * sikj); dCfunc(rij2, rik2, rjk2, &dCikj); - arr1v(dscrfcn, jn) = arr1v(dscrfcn, jn) + coef1 * dCikj; + dscrfcn[jn] = dscrfcn[jn] + coef1 * dCikj; } } coef1 = sfcij; coef2 = sij * dfcij / rij; - arr1v(dscrfcn, jn) = arr1v(dscrfcn, jn) * coef1 - coef2; + dscrfcn[jn] = dscrfcn[jn] * coef1 - coef2; LABEL_100: - arr1v(scrfcn, jn) = sij; - arr1v(fcpair, jn) = fcij; + scrfcn[jn] = sij; + fcpair[jn] = fcij; } } } @@ -252,20 +252,20 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, double ro0i, ro0j; double rhoa0i, rhoa1i, rhoa2i, rhoa3i, A1i, A2i, A3i; - elti = fmap[arr1v(type, i)]; - xtmp = arr2v(x, 1, i); - ytmp = arr2v(x, 2, i); - ztmp = arr2v(x, 3, i); - for (jn = 1; jn <= numneigh; jn++) { - if (!iszero(arr1v(scrfcn, jn))) { - j = arr1v(firstneigh, jn); - sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); - delij[1] = arr2v(x, 1, j) - xtmp; - delij[2] = arr2v(x, 2, j) - ytmp; - delij[3] = arr2v(x, 3, j) - ztmp; + elti = fmap[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + for (jn = 0; jn < numneigh; jn++) { + 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]; if (rij2 < this->cutforcesq) { - eltj = fmap[arr1v(type, j)]; + eltj = fmap[type[j]]; rij = sqrt(rij2); ai = rij / this->re_meam[elti][elti] - 1.0; aj = rij / this->re_meam[eltj][eltj] - 1.0; @@ -287,45 +287,27 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, rhoa2i = rhoa2i * this->t2_meam[elti]; rhoa3i = rhoa3i * this->t3_meam[elti]; } - arr1v(rho0, i) = arr1v(rho0, i) + rhoa0j; - arr1v(rho0, j) = arr1v(rho0, j) + rhoa0i; + rho0[i] = rho0[i] + rhoa0j; + rho0[j] = rho0[j] + rhoa0i; // For ialloy = 2, use single-element value (not average) if (this->ialloy != 2) { - arr2v(t_ave, 1, i) = - arr2v(t_ave, 1, i) + this->t1_meam[eltj] * rhoa0j; - arr2v(t_ave, 2, i) = - arr2v(t_ave, 2, i) + this->t2_meam[eltj] * rhoa0j; - arr2v(t_ave, 3, i) = - arr2v(t_ave, 3, i) + this->t3_meam[eltj] * rhoa0j; - arr2v(t_ave, 1, j) = - arr2v(t_ave, 1, j) + this->t1_meam[elti] * rhoa0i; - arr2v(t_ave, 2, j) = - arr2v(t_ave, 2, j) + this->t2_meam[elti] * rhoa0i; - arr2v(t_ave, 3, j) = - arr2v(t_ave, 3, j) + this->t3_meam[elti] * rhoa0i; + t_ave[i][0] = t_ave[i][0] + this->t1_meam[eltj] * rhoa0j; + t_ave[i][1] = t_ave[i][1] + this->t2_meam[eltj] * rhoa0j; + t_ave[i][2] = t_ave[i][2] + this->t3_meam[eltj] * rhoa0j; + t_ave[j][0] = t_ave[j][0] + this->t1_meam[elti] * rhoa0i; + t_ave[j][1] = t_ave[j][1] + this->t2_meam[elti] * rhoa0i; + t_ave[j][2] = t_ave[j][2] + this->t3_meam[elti] * rhoa0i; } if (this->ialloy == 1) { - arr2v(tsq_ave, 1, i) = - arr2v(tsq_ave, 1, i) + - this->t1_meam[eltj] * this->t1_meam[eltj] * rhoa0j; - arr2v(tsq_ave, 2, i) = - arr2v(tsq_ave, 2, i) + - this->t2_meam[eltj] * this->t2_meam[eltj] * rhoa0j; - arr2v(tsq_ave, 3, i) = - arr2v(tsq_ave, 3, i) + - this->t3_meam[eltj] * this->t3_meam[eltj] * rhoa0j; - arr2v(tsq_ave, 1, j) = - arr2v(tsq_ave, 1, j) + - this->t1_meam[elti] * this->t1_meam[elti] * rhoa0i; - arr2v(tsq_ave, 2, j) = - arr2v(tsq_ave, 2, j) + - this->t2_meam[elti] * this->t2_meam[elti] * rhoa0i; - arr2v(tsq_ave, 3, j) = - arr2v(tsq_ave, 3, j) + - this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i; + tsq_ave[i][0] = tsq_ave[i][0] + this->t1_meam[eltj] * this->t1_meam[eltj] * rhoa0j; + tsq_ave[i][1] = tsq_ave[i][1] + this->t2_meam[eltj] * this->t2_meam[eltj] * rhoa0j; + tsq_ave[i][2] = tsq_ave[i][2] + this->t3_meam[eltj] * this->t3_meam[eltj] * rhoa0j; + tsq_ave[j][0] = tsq_ave[j][0] + this->t1_meam[elti] * this->t1_meam[elti] * rhoa0i; + tsq_ave[j][1] = tsq_ave[j][1] + this->t2_meam[elti] * this->t2_meam[elti] * rhoa0i; + tsq_ave[j][2] = tsq_ave[j][2] + this->t3_meam[elti] * this->t3_meam[elti] * rhoa0i; } - arr1v(arho2b, i) = arr1v(arho2b, i) + rhoa2j; - arr1v(arho2b, j) = arr1v(arho2b, j) + rhoa2i; + arho2b[i] = arho2b[i] + rhoa2j; + arho2b[j] = arho2b[j] + rhoa2i; A1j = rhoa1j / rij; A2j = rhoa2j / rij2; @@ -336,21 +318,21 @@ MEAM::calc_rho1(int i, int ntype, int* type, int* fmap, double** x, nv2 = 1; nv3 = 1; for (m = 1; m <= 3; m++) { - arr2v(arho1, m, i) = arr2v(arho1, m, i) + A1j * delij[m]; - arr2v(arho1, m, j) = arr2v(arho1, m, j) - A1i * delij[m]; - arr2v(arho3b, m, i) = arr2v(arho3b, m, i) + rhoa3j * delij[m] / rij; - arr2v(arho3b, m, j) = arr2v(arho3b, m, j) - rhoa3i * delij[m] / rij; + 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) = - arr2v(arho2, nv2, i) + A2j * delij[m] * delij[n]; - arr2v(arho2, nv2, j) = - arr2v(arho2, nv2, j) + A2i * delij[m] * delij[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 = nv2 + 1; for (p = n; p <= 3; p++) { - arr2v(arho3, nv3, i) = - arr2v(arho3, nv3, i) + A3j * delij[m] * delij[n] * delij[p]; - arr2v(arho3, nv3, j) = - arr2v(arho3, nv3, j) - A3i * delij[m] * delij[n] * delij[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]; nv3 = nv3 + 1; } } @@ -379,26 +361,26 @@ MEAM::screen(int i, int j, double** x, double rijsq, double* sij, double Cmax, Cmin, rbound; *sij = 1.0; - elti = fmap[arr1v(type, i)]; - eltj = fmap[arr1v(type, j)]; + elti = fmap[type[i]]; + eltj = fmap[type[j]]; // if rjksq > ebound*rijsq, atom k is definitely outside the ellipse rbound = this->ebound_meam[elti][eltj] * rijsq; - for (nk = 1; nk <= numneigh_full; nk++) { - k = arr1v(firstneigh_full, nk); - eltk = fmap[arr1v(type, k)]; + for (nk = 0; nk < numneigh_full; nk++) { + k = firstneigh_full[nk]; + eltk = fmap[type[k]]; if (k == j) continue; - delxjk = arr2v(x, 1, k) - arr2v(x, 1, j); - delyjk = arr2v(x, 2, k) - arr2v(x, 2, j); - delzjk = arr2v(x, 3, k) - arr2v(x, 3, j); + delxjk = x[k][0] - x[j][0]; + delyjk = x[k][1] - x[j][1]; + delzjk = x[k][2] - x[j][2]; rjksq = delxjk * delxjk + delyjk * delyjk + delzjk * delzjk; if (rjksq > rbound) continue; - delxik = arr2v(x, 1, k) - arr2v(x, 1, i); - delyik = arr2v(x, 2, k) - arr2v(x, 2, i); - delzik = arr2v(x, 3, k) - arr2v(x, 3, i); + delxik = x[k][0] - x[i][0]; + delyik = x[k][1] - x[i][1]; + delzik = x[k][2] - x[i][2]; riksq = delxik * delxik + delyik * delyik + delzik * delzik; if (riksq > rbound) continue; @@ -450,10 +432,10 @@ MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2, double rbound, delc, sij, xik, xjk, cikj, sikj, dfc, a; double Cmax, Cmin, dCikj1, dCikj2; - sij = arr1v(scrfcn, jn) * arr1v(fcpair, jn); - elti = fmap[arr1v(type, i)]; - eltj = fmap[arr1v(type, j)]; - eltk = fmap[arr1v(type, k)]; + 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]; @@ -462,14 +444,14 @@ MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2, if (!iszero(sij) && !iszero(sij - 1.0)) { rbound = rij2 * this->ebound_meam[elti][eltj]; delc = Cmax - Cmin; - dxjk = arr2v(x, 1, k) - arr2v(x, 1, j); - dyjk = arr2v(x, 2, k) - arr2v(x, 2, j); - dzjk = arr2v(x, 3, k) - arr2v(x, 3, j); + 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 = arr2v(x, 1, k) - arr2v(x, 1, i); - dyik = arr2v(x, 2, k) - arr2v(x, 2, i); - dzik = arr2v(x, 3, k) - arr2v(x, 3, i); + 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; diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index 674f766c44..f29021dd4f 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -26,14 +26,14 @@ using namespace LAMMPS_NS; // void -MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, - int* eflag_atom, int* vflag_atom, double* eng_vdwl, double* eatom, - int* ntype, int* type, int* fmap, double** x, int* numneigh, - int* firstneigh, int* numneigh_full, int* firstneigh_full, +MEAM::meam_force(int i, int eflag_either, int eflag_global, + int eflag_atom, int vflag_atom, double* eng_vdwl, double* eatom, + int ntype, int* type, int* fmap, double** x, int numneigh, + int* firstneigh, int numneigh_full, int* firstneigh_full, int fnoffset, double** f, double** vatom, int* errorflag) { - int i, j, jn, k, kn, kk, m, n, p, q; + 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]; @@ -69,29 +69,26 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, third = 1.0 / 3.0; sixth = 1.0 / 6.0; - //: aliased - i = *iptr; - // Compute forces atom i - elti = fmap[arr1v(type, i)]; + elti = fmap[type[i]]; if (elti >= 0) { - xitmp = arr2v(x, 1, i); - yitmp = arr2v(x, 2, i); - zitmp = arr2v(x, 3, i); + xitmp = x[i][0]; + yitmp = x[i][1]; + zitmp = x[i][2]; // Treat each pair - for (jn = 1; jn <= *numneigh; jn++) { - j = arr1v(firstneigh, jn); - eltj = fmap[arr1v(type, j)]; + for (jn = 0; jn < numneigh; jn++) { + j = firstneigh[jn]; + eltj = fmap[type[j]]; - if (!iszero(arr1v(scrfcn, fnoffset + jn)) && eltj >= 0) { + if (!iszero(scrfcn[fnoffset + jn]) && eltj >= 0) { - sij = arr1v(scrfcn, fnoffset + jn) * arr1v(fcpair, fnoffset + jn); - delij[1] = arr2v(x, 1, j) - xitmp; - delij[2] = arr2v(x, 2, j) - yitmp; - delij[3] = arr2v(x, 3, j) - zitmp; + 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]; if (rij2 < this->cutforcesq) { rij = sqrt(rij2); @@ -116,12 +113,12 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, this->phirar4[ind][kk]; recip = 1.0 / r; - if (*eflag_either != 0) { - if (*eflag_global != 0) + if (eflag_either != 0) { + if (eflag_global != 0) *eng_vdwl = *eng_vdwl + phi * sij; - if (*eflag_atom != 0) { - arr1v(eatom, i) = arr1v(eatom, i) + 0.5 * phi * sij; - arr1v(eatom, j) = arr1v(eatom, j) + 0.5 * phi * sij; + if (eflag_atom != 0) { + eatom[i] = eatom[i] + 0.5 * phi * sij; + eatom[j] = eatom[j] + 0.5 * phi * sij; } } @@ -193,19 +190,19 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, 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) * arg; - arg1j3 = arg1j3 - arr2v(arho3, nv3, j) * arg; + arg1i3 = arg1i3 + arr2v(arho3, nv3, i+1) * arg; + arg1j3 = arg1j3 - arr2v(arho3, nv3, j+1) * arg; nv3 = nv3 + 1; } arg = delij[n] * delij[p] * this->v2D[nv2]; - arg1i2 = arg1i2 + arr2v(arho2, nv2, i) * arg; - arg1j2 = arg1j2 + arr2v(arho2, nv2, j) * arg; + arg1i2 = arg1i2 + arr2v(arho2, nv2, i+1) * arg; + arg1j2 = arg1j2 + arr2v(arho2, nv2, j+1) * arg; nv2 = nv2 + 1; } - arg1i1 = arg1i1 + arr2v(arho1, n, i) * delij[n]; - arg1j1 = arg1j1 - arr2v(arho1, n, j) * delij[n]; - arg3i3 = arg3i3 + arr2v(arho3b, n, i) * delij[n]; - arg3j3 = arg3j3 - arr2v(arho3b, n, j) * delij[n]; + 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]; } // rho0 terms @@ -218,25 +215,25 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, 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); - drho1drm2[m] = -a1 * rhoa1i * arr2v(arho1, m, j); + drho1drm1[m] = a1 * rhoa1j * arr2v(arho1, m, i+1); + drho1drm2[m] = -a1 * rhoa1i * arr2v(arho1, m, j+1); } // rho2 terms a2 = 2 * sij / rij2; drho2dr1 = a2 * (drhoa2j - 2 * rhoa2j / rij) * arg1i2 - - 2.0 / 3.0 * arr1v(arho2b, i) * drhoa2j * sij; + 2.0 / 3.0 * arho2b[i] * drhoa2j * sij; drho2dr2 = a2 * (drhoa2i - 2 * rhoa2i / rij) * arg1j2 - - 2.0 / 3.0 * arr1v(arho2b, j) * drhoa2i * sij; + 2.0 / 3.0 * arho2b[j] * drhoa2i * sij; a2 = 4 * sij / rij2; for (m = 1; 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) * delij[n]; + arr2v(arho2, this->vind2D[m][n], i+1) * delij[n]; drho2drm2[m] = drho2drm2[m] - - arr2v(arho2, this->vind2D[m][n], j) * delij[n]; + arr2v(arho2, this->vind2D[m][n], j+1) * delij[n]; } drho2drm1[m] = a2 * rhoa2j * drho2drm1[m]; drho2drm2[m] = -a2 * rhoa2i * drho2drm2[m]; @@ -260,25 +257,25 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, 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) * arg; + arr2v(arho3, this->vind3D[m][n][p], i+1) * arg; drho3drm2[m] = drho3drm2[m] + - arr2v(arho3, this->vind3D[m][n][p], j) * arg; + arr2v(arho3, this->vind3D[m][n][p], j+1) * arg; nv2 = nv2 + 1; } } drho3drm1[m] = - (a3 * drho3drm1[m] - a3a * arr2v(arho3b, m, i)) * rhoa3j; + (a3 * drho3drm1[m] - a3a * arr2v(arho3b, m, i+1)) * rhoa3j; drho3drm2[m] = - (-a3 * drho3drm2[m] + a3a * arr2v(arho3b, m, j)) * rhoa3i; + (-a3 * drho3drm2[m] + a3a * arr2v(arho3b, m, j+1)) * rhoa3i; } // Compute derivatives of weighting functions t wrt rij - t1i = arr2v(t_ave, 1, i); - t2i = arr2v(t_ave, 2, i); - t3i = arr2v(t_ave, 3, i); - t1j = arr2v(t_ave, 1, j); - t2j = arr2v(t_ave, 2, j); - t3j = arr2v(t_ave, 3, j); + 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]; if (this->ialloy == 1) { @@ -288,18 +285,18 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, a2j = 0.0; a3i = 0.0; a3j = 0.0; - if (!iszero(arr2v(tsq_ave, 1, i))) - a1i = drhoa0j * sij / arr2v(tsq_ave, 1, i); - if (!iszero(arr2v(tsq_ave, 1, j))) - a1j = drhoa0i * sij / arr2v(tsq_ave, 1, j); - if (!iszero(arr2v(tsq_ave, 2, i))) - a2i = drhoa0j * sij / arr2v(tsq_ave, 2, i); - if (!iszero(arr2v(tsq_ave, 2, j))) - a2j = drhoa0i * sij / arr2v(tsq_ave, 2, j); - if (!iszero(arr2v(tsq_ave, 3, i))) - a3i = drhoa0j * sij / arr2v(tsq_ave, 3, i); - if (!iszero(arr2v(tsq_ave, 3, j))) - a3j = drhoa0i * sij / arr2v(tsq_ave, 3, j); + 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]; dt1dr1 = a1i * (this->t1_meam[eltj] - t1i * pow(this->t1_meam[eltj], 2)); @@ -326,11 +323,11 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, } else { ai = 0.0; - if (!iszero(arr1v(rho0, i))) - ai = drhoa0j * sij / arr1v(rho0, i); + if (!iszero(rho0[i])) + ai = drhoa0j * sij / rho0[i]; aj = 0.0; - if (!iszero(arr1v(rho0, j))) - aj = drhoa0i * sij / arr1v(rho0, j); + if (!iszero(rho0[j])) + aj = drhoa0i * sij / rho0[j]; dt1dr1 = ai * (this->t1_meam[eltj] - t1i); dt1dr2 = aj * (this->t1_meam[elti] - t1j); @@ -344,32 +341,32 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, get_shpfcn(shpi, this->lattce_meam[elti][elti]); get_shpfcn(shpj, this->lattce_meam[eltj][eltj]); drhodr1 = - arr1v(dgamma1, i) * drho0dr1 + - arr1v(dgamma2, i) * (dt1dr1 * arr1v(rho1, i) + t1i * drho1dr1 + - dt2dr1 * arr1v(rho2, i) + t2i * drho2dr1 + - dt3dr1 * arr1v(rho3, i) + t3i * drho3dr1) - - arr1v(dgamma3, i) * + dgamma1[i] * drho0dr1 + + dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + + dt2dr1 * rho2[i] + t2i * drho2dr1 + + dt3dr1 * rho3[i] + t3i * drho3dr1) - + dgamma3[i] * (shpi[1] * dt1dr1 + shpi[2] * dt2dr1 + shpi[3] * dt3dr1); drhodr2 = - arr1v(dgamma1, j) * drho0dr2 + - arr1v(dgamma2, j) * (dt1dr2 * arr1v(rho1, j) + t1j * drho1dr2 + - dt2dr2 * arr1v(rho2, j) + t2j * drho2dr2 + - dt3dr2 * arr1v(rho3, j) + t3j * drho3dr2) - - arr1v(dgamma3, j) * + 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++) { drhodrm1[m] = 0.0; drhodrm2[m] = 0.0; drhodrm1[m] = - arr1v(dgamma2, i) * + dgamma2[i] * (t1i * drho1drm1[m] + t2i * drho2drm1[m] + t3i * drho3drm1[m]); drhodrm2[m] = - arr1v(dgamma2, j) * + dgamma2[j] * (t1j * drho1drm2[m] + t2j * drho2drm2[m] + t3j * drho3drm2[m]); } // Compute derivatives wrt sij, but only if necessary - if (!iszero(arr1v(dscrfcn, fnoffset + jn))) { + if (!iszero(dscrfcn[fnoffset + jn])) { drho0ds1 = rhoa0j; drho0ds2 = rhoa0i; a1 = 2.0 / rij; @@ -377,9 +374,9 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, drho1ds2 = a1 * rhoa1i * arg1j1; a2 = 2.0 / rij2; drho2ds1 = - a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arr1v(arho2b, i) * rhoa2j; + a2 * rhoa2j * arg1i2 - 2.0 / 3.0 * arho2b[i] * rhoa2j; drho2ds2 = - a2 * rhoa2i * arg1j2 - 2.0 / 3.0 * arr1v(arho2b, j) * rhoa2i; + 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; @@ -393,18 +390,12 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, a2j = 0.0; a3i = 0.0; a3j = 0.0; - if (!iszero(arr2v(tsq_ave, 1, i))) - a1i = rhoa0j / arr2v(tsq_ave, 1, i); - if (!iszero(arr2v(tsq_ave, 1, j))) - a1j = rhoa0i / arr2v(tsq_ave, 1, j); - if (!iszero(arr2v(tsq_ave, 2, i))) - a2i = rhoa0j / arr2v(tsq_ave, 2, i); - if (!iszero(arr2v(tsq_ave, 2, j))) - a2j = rhoa0i / arr2v(tsq_ave, 2, j); - if (!iszero(arr2v(tsq_ave, 3, i))) - a3i = rhoa0j / arr2v(tsq_ave, 3, i); - if (!iszero(arr2v(tsq_ave, 3, j))) - a3j = rhoa0i / arr2v(tsq_ave, 3, j); + 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]; dt1ds1 = a1i * (this->t1_meam[eltj] - t1i * pow(this->t1_meam[eltj], 2)); @@ -431,11 +422,11 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, } else { ai = 0.0; - if (!iszero(arr1v(rho0, i))) - ai = rhoa0j / arr1v(rho0, i); + if (!iszero(rho0[i])) + ai = rhoa0j / rho0[i]; aj = 0.0; - if (!iszero(arr1v(rho0, j))) - aj = rhoa0i / arr1v(rho0, j); + if (!iszero(rho0[j])) + aj = rhoa0i / rho0[j]; dt1ds1 = ai * (this->t1_meam[eltj] - t1i); dt1ds2 = aj * (this->t1_meam[elti] - t1j); @@ -446,44 +437,44 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, } drhods1 = - arr1v(dgamma1, i) * drho0ds1 + - arr1v(dgamma2, i) * (dt1ds1 * arr1v(rho1, i) + t1i * drho1ds1 + - dt2ds1 * arr1v(rho2, i) + t2i * drho2ds1 + - dt3ds1 * arr1v(rho3, i) + t3i * drho3ds1) - - arr1v(dgamma3, i) * + dgamma1[i] * drho0ds1 + + dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + + dt2ds1 * rho2[i] + t2i * drho2ds1 + + dt3ds1 * rho3[i] + t3i * drho3ds1) - + dgamma3[i] * (shpi[1] * dt1ds1 + shpi[2] * dt2ds1 + shpi[3] * dt3ds1); drhods2 = - arr1v(dgamma1, j) * drho0ds2 + - arr1v(dgamma2, j) * (dt1ds2 * arr1v(rho1, j) + t1j * drho1ds2 + - dt2ds2 * arr1v(rho2, j) + t2j * drho2ds2 + - dt3ds2 * arr1v(rho3, j) + t3j * drho3ds2) - - arr1v(dgamma3, j) * + 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); } // Compute derivatives of energy wrt rij, sij and rij[3] - dUdrij = phip * sij + arr1v(frhop, i) * drhodr1 + arr1v(frhop, j) * drhodr2; + dUdrij = phip * sij + frhop[i] * drhodr1 + frhop[j] * drhodr2; dUdsij = 0.0; - if (!iszero(arr1v(dscrfcn, fnoffset + jn))) { - dUdsij = phi + arr1v(frhop, i) * drhods1 + arr1v(frhop, j) * drhods2; + if (!iszero(dscrfcn[fnoffset + jn])) { + dUdsij = phi + frhop[i] * drhods1 + frhop[j] * drhods2; } for (m = 1; m <= 3; m++) { dUdrijm[m] = - arr1v(frhop, i) * drhodrm1[m] + arr1v(frhop, j) * drhodrm2[m]; + frhop[i] * drhodrm1[m] + frhop[j] * drhodrm2[m]; } // Add the part of the force due to dUdrij and dUdsij - force = dUdrij * recip + dUdsij * arr1v(dscrfcn, fnoffset + jn); + force = dUdrij * recip + dUdsij * dscrfcn[fnoffset + jn]; for (m = 1; m <= 3; m++) { forcem = delij[m] * force + dUdrijm[m]; - arr2v(f, m, i) = arr2v(f, m, i) + forcem; - arr2v(f, m, j) = arr2v(f, m, j) - forcem; + f[i][m-1] = f[i][m-1] + forcem; + f[j][m-1] = f[j][m-1] - forcem; } // Tabulate per-atom virial as symmetrized stress tensor - if (*vflag_atom != 0) { + if (vflag_atom != 0) { fi[1] = delij[1] * force + dUdrijm[1]; fi[2] = delij[2] * force + dUdrijm[2]; fi[3] = delij[3] * force + dUdrijm[3]; @@ -494,47 +485,46 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, v[5] = -0.25 * (delij[1] * fi[3] + delij[3] * fi[1]); v[6] = -0.25 * (delij[2] * fi[3] + delij[3] * fi[2]); - arr2v(vatom, 1, i) = arr2v(vatom, 1, i) + v[1]; - arr2v(vatom, 2, i) = arr2v(vatom, 2, i) + v[2]; - arr2v(vatom, 3, i) = arr2v(vatom, 3, i) + v[3]; - arr2v(vatom, 4, i) = arr2v(vatom, 4, i) + v[4]; - arr2v(vatom, 5, i) = arr2v(vatom, 5, i) + v[5]; - arr2v(vatom, 6, i) = arr2v(vatom, 6, i) + v[6]; - arr2v(vatom, 1, j) = arr2v(vatom, 1, j) + v[1]; - arr2v(vatom, 2, j) = arr2v(vatom, 2, j) + v[2]; - arr2v(vatom, 3, j) = arr2v(vatom, 3, j) + v[3]; - arr2v(vatom, 4, j) = arr2v(vatom, 4, j) + v[4]; - arr2v(vatom, 5, j) = arr2v(vatom, 5, j) + v[5]; - arr2v(vatom, 6, j) = arr2v(vatom, 6, j) + v[6]; + 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]; } // 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 = 1; kn <= *numneigh_full; kn++) { - k = arr1v(firstneigh_full, kn); - eltk = fmap[arr1v(type, k)]; + for (kn = 0; kn < numneigh_full; kn++) { + k = firstneigh_full[kn]; + eltk = fmap[type[k]]; if (k != j && eltk >= 0) { - dsij(i, j, k, jn, *numneigh, rij2, &dsij1, &dsij2, *ntype, + dsij(i, j, k, jn, numneigh, rij2, &dsij1, &dsij2, ntype, type, fmap, x, &scrfcn[fnoffset], &fcpair[fnoffset]); if (!iszero(dsij1) || !iszero(dsij2)) { force1 = dUdsij * dsij1; force2 = dUdsij * dsij2; for (m = 1; m <= 3; m++) { - delik[m] = arr2v(x, m, k) - arr2v(x, m, i); - deljk[m] = arr2v(x, m, k) - arr2v(x, m, j); + 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 = 1; m <= 3; m++) { - arr2v(f, m, i) = arr2v(f, m, i) + force1 * delik[m]; - arr2v(f, m, j) = arr2v(f, m, j) + force2 * deljk[m]; - arr2v(f, m, k) = - arr2v(f, m, k) - force1 * delik[m] - force2 * deljk[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]; } // Tabulate per-atom virial as symmetrized stress tensor - if (*vflag_atom != 0) { + if (vflag_atom != 0) { fi[1] = force1 * delik[1]; fi[2] = force1 * delik[2]; fi[3] = force1 * delik[3]; @@ -551,24 +541,24 @@ MEAM::meam_force(int* iptr, int* eflag_either, int* eflag_global, v[6] = -sixth * (delik[2] * fi[3] + deljk[2] * fj[3] + delik[3] * fi[2] + deljk[3] * fj[2]); - arr2v(vatom, 1, i) = arr2v(vatom, 1, i) + v[1]; - arr2v(vatom, 2, i) = arr2v(vatom, 2, i) + v[2]; - arr2v(vatom, 3, i) = arr2v(vatom, 3, i) + v[3]; - arr2v(vatom, 4, i) = arr2v(vatom, 4, i) + v[4]; - arr2v(vatom, 5, i) = arr2v(vatom, 5, i) + v[5]; - arr2v(vatom, 6, i) = arr2v(vatom, 6, i) + v[6]; - arr2v(vatom, 1, j) = arr2v(vatom, 1, j) + v[1]; - arr2v(vatom, 2, j) = arr2v(vatom, 2, j) + v[2]; - arr2v(vatom, 3, j) = arr2v(vatom, 3, j) + v[3]; - arr2v(vatom, 4, j) = arr2v(vatom, 4, j) + v[4]; - arr2v(vatom, 5, j) = arr2v(vatom, 5, j) + v[5]; - arr2v(vatom, 6, j) = arr2v(vatom, 6, j) + v[6]; - arr2v(vatom, 1, k) = arr2v(vatom, 1, k) + v[1]; - arr2v(vatom, 2, k) = arr2v(vatom, 2, k) + v[2]; - arr2v(vatom, 3, k) = arr2v(vatom, 3, k) + v[3]; - arr2v(vatom, 4, k) = arr2v(vatom, 4, k) + v[4]; - arr2v(vatom, 5, k) = arr2v(vatom, 5, k) + v[5]; - arr2v(vatom, 6, k) = arr2v(vatom, 6, k) + v[6]; + 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]; } } } diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index bf096c1a1e..e54f346589 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -129,24 +129,17 @@ void PairMEAMC::compute(int eflag, int vflag) int *type = atom->type; int ntype = atom->ntypes; - // change neighbor list indices to Fortran indexing - - neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half); - neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full); - // 3 stages of MEAM calculation // loop over my atoms followed by communication - int ifort; int offset = 0; errorflag = 0; for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; - ifort = i+1; - meam_inst->meam_dens_init(&ifort,&ntype,type,map,x, - &numneigh_half[i],firstneigh_half[i], - &numneigh_full[i],firstneigh_full[i], + meam_inst->meam_dens_init(i,ntype,type,map,x, + numneigh_half[i],firstneigh_half[i], + numneigh_full[i],firstneigh_full[i], offset, &errorflag); if (errorflag) { @@ -159,8 +152,8 @@ void PairMEAMC::compute(int eflag, int vflag) comm->reverse_comm_pair(this); - meam_inst->meam_dens_final(&nlocal,&eflag_either,&eflag_global,&eflag_atom, - &eng_vdwl,eatom,&ntype,type,map, + meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom, + &eng_vdwl,eatom,ntype,type,map, &errorflag); if (errorflag) { char str[128]; @@ -181,11 +174,10 @@ void PairMEAMC::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; - ifort = i+1; - meam_inst->meam_force(&ifort,&eflag_either,&eflag_global,&eflag_atom, - &vflag_atom,&eng_vdwl,eatom,&ntype,type,map,x, - &numneigh_half[i],firstneigh_half[i], - &numneigh_full[i],firstneigh_full[i], + meam_inst->meam_force(i,eflag_either,eflag_global,eflag_atom, + vflag_atom,&eng_vdwl,eatom,ntype,type,map,x, + numneigh_half[i],firstneigh_half[i], + numneigh_full[i],firstneigh_full[i], offset, f,vptr,&errorflag); if (errorflag) { @@ -196,11 +188,6 @@ void PairMEAMC::compute(int eflag, int vflag) offset += numneigh_half[i]; } - // change neighbor list indices back to C indexing - - neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); - neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); - if (vflag_fdotr) virial_fdotr_compute(); } @@ -806,34 +793,3 @@ void PairMEAMC::neigh_strip(int inum, int *ilist, for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK; } } - -/* ---------------------------------------------------------------------- - toggle neighbor list indices between zero- and one-based values - needed for access by MEAM Fortran library -------------------------------------------------------------------------- */ - -void PairMEAMC::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh) -{ - int i,j,ii,jnum; - int *jlist; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - for (j = 0; j < jnum; j++) jlist[j]--; - } -} - -void PairMEAMC::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh) -{ - int i,j,ii,jnum; - int *jlist; - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - for (j = 0; j < jnum; j++) jlist[j]++; - } -} diff --git a/src/USER-MEAMC/pair_meamc.h b/src/USER-MEAMC/pair_meamc.h index fc31cd55c6..476a70dd04 100644 --- a/src/USER-MEAMC/pair_meamc.h +++ b/src/USER-MEAMC/pair_meamc.h @@ -54,8 +54,6 @@ class PairMEAMC : public Pair { void allocate(); void read_files(char *, char *); void neigh_strip(int, int *, int *, int **); - void neigh_f2c(int, int *, int *, int **); - void neigh_c2f(int, int *, int *, int **); }; }