From f82096c46ccf7e284c6f3d1c3b8af2eaa2bd4e1d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 11 Jul 2021 00:55:22 -0400 Subject: [PATCH] correct nofdotr handling for tersoff styles --- src/MANYBODY/pair_tersoff.cpp | 12 +++++------ src/OPENMP/pair_tersoff_mod_c_omp.cpp | 12 +++++------ src/OPENMP/pair_tersoff_mod_omp.cpp | 12 +++++------ src/OPENMP/pair_tersoff_omp.cpp | 12 +++++------ src/OPENMP/pair_tersoff_table_omp.cpp | 8 +++---- src/OPENMP/thr_omp.cpp | 31 ++++++++++++++++++--------- src/OPENMP/thr_omp.h | 4 ++-- src/pair.cpp | 23 ++++++++++++++------ 8 files changed, 68 insertions(+), 46 deletions(-) diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 33682ce3fa..f21c3bffca 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -88,10 +88,10 @@ void PairTersoff::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(); + if (vflag_either) eval<1,1,1,1>(); else eval<1,1,1,0>(); } else { - if (vflag_atom) eval<1,1,0,1>(); + if (vflag_either) eval<1,1,0,1>(); else eval<1,1,0,0>(); } } else eval<1,0,0,0>(); @@ -100,17 +100,17 @@ void PairTersoff::compute(int eflag, int vflag) if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(); + if (vflag_either) eval<0,1,1,1>(); else eval<0,1,1,0>(); } else { - if (vflag_atom) eval<0,1,0,1>(); + if (vflag_either) eval<0,1,0,1>(); else eval<0,1,0,0>(); } } else eval<0,0,0,0>(); } } -template +template void PairTersoff::eval() { int i,j,k,ii,jj,kk,inum,jnum; @@ -315,7 +315,7 @@ void PairTersoff::eval() f[k][1] += fk[1]; f[k][2] += fk[2]; - if (VFLAG_ATOM) v_tally3(i,j,k,fj,fk,delr1,delr2); + if (VFLAG_EITHER) v_tally3(i,j,k,fj,fk,delr1,delr2); } f[j][0] += fjxtmp; f[j][1] += fjytmp; diff --git a/src/OPENMP/pair_tersoff_mod_c_omp.cpp b/src/OPENMP/pair_tersoff_mod_c_omp.cpp index 0d66e83b00..cd59e06673 100644 --- a/src/OPENMP/pair_tersoff_mod_c_omp.cpp +++ b/src/OPENMP/pair_tersoff_mod_c_omp.cpp @@ -60,20 +60,20 @@ void PairTersoffMODCOMP::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr); else eval<1,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr); else eval<1,1,0,0>(ifrom, ito, thr); } } else eval<1,0,0,0>(ifrom, ito, thr); } else { if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr); else eval<0,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr); else eval<0,1,0,0>(ifrom, ito, thr); } } else eval<0,0,0,0>(ifrom, ito, thr); @@ -84,7 +84,7 @@ void PairTersoffMODCOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum; @@ -282,7 +282,7 @@ void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr) f[k].y += fk[1]; f[k].z += fk[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j].x += fjxtmp; f[j].y += fjytmp; diff --git a/src/OPENMP/pair_tersoff_mod_omp.cpp b/src/OPENMP/pair_tersoff_mod_omp.cpp index b3e3b89283..de19aa3872 100644 --- a/src/OPENMP/pair_tersoff_mod_omp.cpp +++ b/src/OPENMP/pair_tersoff_mod_omp.cpp @@ -60,20 +60,20 @@ void PairTersoffMODOMP::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr); else eval<1,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr); else eval<1,1,0,0>(ifrom, ito, thr); } } else eval<1,0,0,0>(ifrom, ito, thr); } else { if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr); else eval<0,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr); else eval<0,1,0,0>(ifrom, ito, thr); } } else eval<0,0,0,0>(ifrom, ito, thr); @@ -84,7 +84,7 @@ void PairTersoffMODOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum; @@ -282,7 +282,7 @@ void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr) f[k].y += fk[1]; f[k].z += fk[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j].x += fjxtmp; f[j].y += fjytmp; diff --git a/src/OPENMP/pair_tersoff_omp.cpp b/src/OPENMP/pair_tersoff_omp.cpp index 6bf48792eb..2a3cfb5d54 100644 --- a/src/OPENMP/pair_tersoff_omp.cpp +++ b/src/OPENMP/pair_tersoff_omp.cpp @@ -61,10 +61,10 @@ void PairTersoffOMP::compute(int eflag, int vflag) if (shift_flag) { if (evflag) { if (eflag) { - if (vflag_atom) eval<1,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr); else eval<1,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<1,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr); else eval<1,1,0,0>(ifrom, ito, thr); } } else eval<1,0,0,0>(ifrom, ito, thr); @@ -73,10 +73,10 @@ void PairTersoffOMP::compute(int eflag, int vflag) if (evflag) { if (eflag) { - if (vflag_atom) eval<0,1,1,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr); else eval<0,1,1,0>(ifrom, ito, thr); } else { - if (vflag_atom) eval<0,1,0,1>(ifrom, ito, thr); + if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr); else eval<0,1,0,0>(ifrom, ito, thr); } } else eval<0,0,0,0>(ifrom, ito, thr); @@ -87,7 +87,7 @@ void PairTersoffOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jj,kk,jnum,maxshort_thr; @@ -293,7 +293,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr) f[k].y += fk[1]; f[k].z += fk[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr); } f[j].x += fjxtmp; f[j].y += fjytmp; diff --git a/src/OPENMP/pair_tersoff_table_omp.cpp b/src/OPENMP/pair_tersoff_table_omp.cpp index a91f361f79..c5fc60ea90 100644 --- a/src/OPENMP/pair_tersoff_table_omp.cpp +++ b/src/OPENMP/pair_tersoff_table_omp.cpp @@ -80,7 +80,7 @@ void PairTersoffTableOMP::compute(int eflag, int vflag) ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr); if (evflag) - if (vflag_atom) eval<1,1>(ifrom, ito, thr); + if (vflag_either) eval<1,1>(ifrom, ito, thr); else eval<1,0>(ifrom, ito, thr); else eval<0,0>(ifrom, ito, thr); @@ -89,7 +89,7 @@ void PairTersoffTableOMP::compute(int eflag, int vflag) } // end of omp parallel region } -template +template void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) { int i,j,k,ii,jnum; @@ -427,7 +427,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) fytmp += f_ij[1] + f_ik[1]; fztmp += f_ij[2] + f_ik[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); } // second loop over neighbors of atom i except j, forces and virial only - part 2/2 @@ -493,7 +493,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) fytmp += f_ij[1] + f_ik[1]; fztmp += f_ij[2] + f_ik[2]; - if (VFLAG_ATOM) v_tally3_thr(i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); + if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,f_ij,f_ik,dr_ij,dr_ik,thr); } } // loop on J diff --git a/src/OPENMP/thr_omp.cpp b/src/OPENMP/thr_omp.cpp index 4685e116be..867623361c 100644 --- a/src/OPENMP/thr_omp.cpp +++ b/src/OPENMP/thr_omp.cpp @@ -1486,23 +1486,34 @@ void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair, called by AIREBO and Tersoff potential, newton_pair is always on ------------------------------------------------------------------------- */ -void ThrOMP::v_tally3_thr(const int i, const int j, const int k, +void ThrOMP::v_tally3_thr(Pair *pair, const int i, const int j, const int k, const double * const fi, const double * const fj, const double * const drik, const double * const drjk, ThrData * const thr) { double v[6]; - v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); - v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); - v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); - v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); - v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); - v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); + v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]); + v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]); + v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]); + v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]); + v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]); + v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]); - v_tally(thr->vatom_pair[i],v); - v_tally(thr->vatom_pair[j],v); - v_tally(thr->vatom_pair[k],v); + if (pair->vflag_global) v_tally(thr->virial_pair,v); + + if (pair->vflag_atom) { + v[0] *= THIRD; + v[1] *= THIRD; + v[2] *= THIRD; + v[3] *= THIRD; + v[4] *= THIRD; + v[5] *= THIRD; + + v_tally(thr->vatom_pair[i],v); + v_tally(thr->vatom_pair[j],v); + v_tally(thr->vatom_pair[k],v); + } } /* ---------------------------------------------------------------------- diff --git a/src/OPENMP/thr_omp.h b/src/OPENMP/thr_omp.h index 77835122e0..22fdb36f38 100644 --- a/src/OPENMP/thr_omp.h +++ b/src/OPENMP/thr_omp.h @@ -138,6 +138,8 @@ class ThrOMP { void ev_tally3_thr(Pair *const, const int, const int, const int, const double, const double, const double *const, const double *const, const double *const, const double *const, ThrData *const); + void v_tally3_thr(Pair *const, const int, const int, const int, const double *const, + const double *const, const double *const, const double *const, ThrData *const); void ev_tally4_thr(Pair *const, const int, const int, const int, const int, const double, const double *const, const double *const, const double *const, const double *const, const double *const, const double *const, ThrData *const); @@ -170,8 +172,6 @@ class ThrOMP { // style independent versions void v_tally2_thr(const int, const int, const double, const double *const, ThrData *const); - void v_tally3_thr(const int, const int, const int, const double *const, const double *const, - const double *const, const double *const, ThrData *const); void v_tally4_thr(const int, const int, const int, const int, const double *const, const double *const, const double *const, const double *const, const double *const, const double *const, ThrData *const); diff --git a/src/pair.cpp b/src/pair.cpp index 5448a42e1c..f27d2ac370 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1529,12 +1529,23 @@ void Pair::v_tally3(int i, int j, int k, v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); - vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; - vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; - vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; - vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; - vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; - vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + if (vflag_global) { + virial[0] += 3.0*v[0]; + virial[1] += 3.0*v[1]; + virial[2] += 3.0*v[2]; + virial[3] += 3.0*v[3]; + virial[4] += 3.0*v[4]; + virial[5] += 3.0*v[5]; + } + + if (vflag_atom) { + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; + vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + } } /* ----------------------------------------------------------------------