jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

OpenMP version added.

This commit is contained in:
DallasTrinkle 2017-05-04 15:08:04 -05:00
parent 754b40cb31
commit f7230006fe
2 changed files with 29 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/USER-MISC/pair_meam_spline.cpp
View File

@ -304,7 +304,7 @@ void PairMEAMSpline::compute(int eflag, int vflag)
double pair_pot_deriv; double pair_pot_deriv;
double pair_pot = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv); double pair_pot = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv);
fpair += pair_pot_deriv; fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient // Divide by r_ij to get forces from gradient

40
src/USER-OMP/pair_meam_spline_omp.cpp
View File

@ -140,19 +140,25 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
nextTwoBodyInfo->tag = j; nextTwoBodyInfo->tag = j;
nextTwoBodyInfo->r = rij; nextTwoBodyInfo->r = rij;
nextTwoBodyInfo->f = f.eval(rij, nextTwoBodyInfo->fprime); // nextTwoBodyInfo->f = f.eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->f = fs[i_to_potl(j)].eval(rij, nextTwoBodyInfo->fprime);
nextTwoBodyInfo->del[0] = jdelx / rij; nextTwoBodyInfo->del[0] = jdelx / rij;
nextTwoBodyInfo->del[1] = jdely / rij; nextTwoBodyInfo->del[1] = jdely / rij;
nextTwoBodyInfo->del[2] = jdelz / rij; nextTwoBodyInfo->del[2] = jdelz / rij;
for(int kk = 0; kk < numBonds; kk++) { for(int kk = 0; kk < numBonds; kk++) {
const MEAM2Body& bondk = myTwoBodyInfo[kk]; const MEAM2Body& bondk = myTwoBodyInfo[kk];
double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] + nextTwoBodyInfo->del[1]*bondk.del[1] + nextTwoBodyInfo->del[2]*bondk.del[2]); double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] +
partial_sum += bondk.f * g.eval(cos_theta); nextTwoBodyInfo->del[1]*bondk.del[1] +
nextTwoBodyInfo->del[2]*bondk.del[2]);
partial_sum += bondk.f * gs[ij_to_potl(j,bondk.tag)].eval(cos_theta);
// double cos_theta = (nextTwoBodyInfo->del[0]*bondk.del[0] + nextTwoBodyInfo->del[1]*bondk.del[1] + nextTwoBodyInfo->del[2]*bondk.del[2]);
// partial_sum += bondk.f * g.eval(cos_theta);
} }
rho_value += nextTwoBodyInfo->f * partial_sum; rho_value += nextTwoBodyInfo->f * partial_sum;
rho_value += rho.eval(rij); // rho_value += rho.eval(rij);
rho_value += rhos[i_to_potl(j)].eval(rij);
numBonds++; numBonds++;
nextTwoBodyInfo++; nextTwoBodyInfo++;
@ -161,7 +167,9 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
// Compute embedding energy and its derivative. // Compute embedding energy and its derivative.
double Uprime_i; double Uprime_i;
double embeddingEnergy = U.eval(rho_value, Uprime_i) - zero_atom_energy; // double embeddingEnergy = U.eval(rho_value, Uprime_i) - zero_atom_energy;
double embeddingEnergy = Us[i_to_potl(i)].eval(rho_value, Uprime_i)
- zero_atom_energies[i_to_potl(i)];
Uprime_thr[i] = Uprime_i; Uprime_thr[i] = Uprime_i;
if (EFLAG) if (EFLAG)
e_tally_thr(this,i,i,nlocal,1/*newton_pair*/,embeddingEnergy,0.0,thr); e_tally_thr(this,i,i,nlocal,1/*newton_pair*/,embeddingEnergy,0.0,thr);
@ -187,7 +195,8 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
+ bondj.del[1]*bondk->del[1] + bondj.del[1]*bondk->del[1]
+ bondj.del[2]*bondk->del[2]); + bondj.del[2]*bondk->del[2]);
double g_prime; double g_prime;
double g_value = g.eval(cos_theta, g_prime); // double g_value = g.eval(cos_theta, g_prime);
double g_value = gs[ij_to_potl(j,bondk->tag)].eval(cos_theta, g_prime);
const double f_rik_prime = bondk->fprime; const double f_rik_prime = bondk->fprime;
const double f_rik = bondk->f; const double f_rik = bondk->f;
@ -292,13 +301,20 @@ void PairMEAMSplineOMP::eval(int iifrom, int iito, ThrData * const thr)
if(rij_sq < cutforcesq) { if(rij_sq < cutforcesq) {
double rij = sqrt(rij_sq); double rij = sqrt(rij_sq);
double rho_prime; // double rho_prime;
rho.eval(rij, rho_prime); // rho.eval(rij, rho_prime);
double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]); // double fpair = rho_prime * (Uprime_values[i] + Uprime_values[j]);
double rho_prime_i,rho_prime_j;
rhos[i_to_potl(i)].eval(rij,rho_prime_i);
rhos[i_to_potl(j)].eval(rij,rho_prime_j);
double fpair = rho_prime_j * Uprime_values[i] + rho_prime_i*Uprime_values[j];
// double pair_pot_deriv;
// double pair_pot = phi.eval(rij, pair_pot_deriv);
double pair_pot_deriv; double pair_pot_deriv;
double pair_pot = phi.eval(rij, pair_pot_deriv); double pair_pot = phis[ij_to_potl(i,j)].eval(rij, pair_pot_deriv);
fpair += pair_pot_deriv;
fpair += pair_pot_deriv;
// Divide by r_ij to get forces from gradient. // Divide by r_ij to get forces from gradient.
fpair /= rij; fpair /= rij;