1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 This software is distributed under the GNU General Public License.
8
9 See the README file in the top-level LAMMPS directory.
10 ------------------------------------------------------------------------- */
11
12 /* ----------------------------------------------------------------------
13 Contributing author: Axel Kohlmeyer (Temple U)
14 ------------------------------------------------------------------------- */
15
16 #include "pair_tersoff_mod_c_omp.h"
17
18 #include "atom.h"
19 #include "comm.h"
20 #include "math_extra.h"
21 #include "neigh_list.h"
22 #include "suffix.h"
23
24 #include <cmath>
25
26 #include "omp_compat.h"
27 using namespace LAMMPS_NS;
28 using namespace MathExtra;
29
30 /* ---------------------------------------------------------------------- */
31
PairTersoffMODCOMP(LAMMPS * lmp)32 PairTersoffMODCOMP::PairTersoffMODCOMP(LAMMPS *lmp) :
33 PairTersoffMODC(lmp), ThrOMP(lmp, THR_PAIR)
34 {
35 suffix_flag |= Suffix::OMP;
36 respa_enable = 0;
37 }
38
39 /* ---------------------------------------------------------------------- */
40
compute(int eflag,int vflag)41 void PairTersoffMODCOMP::compute(int eflag, int vflag)
42 {
43 ev_init(eflag,vflag);
44
45 const int nall = atom->nlocal + atom->nghost;
46 const int nthreads = comm->nthreads;
47 const int inum = list->inum;
48
49 #if defined(_OPENMP)
50 #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
51 #endif
52 {
53 int ifrom, ito, tid;
54
55 loop_setup_thr(ifrom, ito, tid, inum, nthreads);
56 ThrData *thr = fix->get_thr(tid);
57 thr->timer(Timer::START);
58 ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
59
60 if (shift_flag) {
61 if (evflag) {
62 if (eflag) {
63 if (vflag_either) eval<1,1,1,1>(ifrom, ito, thr);
64 else eval<1,1,1,0>(ifrom, ito, thr);
65 } else {
66 if (vflag_either) eval<1,1,0,1>(ifrom, ito, thr);
67 else eval<1,1,0,0>(ifrom, ito, thr);
68 }
69 } else eval<1,0,0,0>(ifrom, ito, thr);
70 } else {
71 if (evflag) {
72 if (eflag) {
73 if (vflag_either) eval<0,1,1,1>(ifrom, ito, thr);
74 else eval<0,1,1,0>(ifrom, ito, thr);
75 } else {
76 if (vflag_either) eval<0,1,0,1>(ifrom, ito, thr);
77 else eval<0,1,0,0>(ifrom, ito, thr);
78 }
79 } else eval<0,0,0,0>(ifrom, ito, thr);
80 }
81
82 thr->timer(Timer::PAIR);
83 reduce_thr(this, eflag, vflag, thr);
84 } // end of omp parallel region
85 }
86
87 template <int SHIFT_FLAG, int EVFLAG, int EFLAG, int VFLAG_EITHER>
eval(int iifrom,int iito,ThrData * const thr)88 void PairTersoffMODCOMP::eval(int iifrom, int iito, ThrData * const thr)
89 {
90 int i,j,k,ii,jj,kk,jnum;
91 tagint itag,jtag;
92 int itype,jtype,ktype,iparam_ij,iparam_ijk;
93 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
94 double fforce;
95 double rsq,rsq1,rsq2;
96 double delr1[3],delr2[3],fi[3],fj[3],fk[3];
97 double r1_hat[3],r2_hat[3];
98 double zeta_ij,prefactor;
99 double forceshiftfac;
100 int *ilist,*jlist,*numneigh,**firstneigh;
101
102 evdwl = 0.0;
103
104 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
105 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
106 const tagint * _noalias const tag = atom->tag;
107 const int * _noalias const type = atom->type;
108 const int nlocal = atom->nlocal;
109
110 ilist = list->ilist;
111 numneigh = list->numneigh;
112 firstneigh = list->firstneigh;
113
114 double fxtmp,fytmp,fztmp;
115
116 // loop over full neighbor list of my atoms
117
118 for (ii = iifrom; ii < iito; ++ii) {
119
120 i = ilist[ii];
121 itag = tag[i];
122 itype = map[type[i]];
123 xtmp = x[i].x;
124 ytmp = x[i].y;
125 ztmp = x[i].z;
126 fxtmp = fytmp = fztmp = 0.0;
127
128 // two-body interactions, skip half of them
129
130 jlist = firstneigh[i];
131 jnum = numneigh[i];
132
133 for (jj = 0; jj < jnum; jj++) {
134 j = jlist[jj];
135 j &= NEIGHMASK;
136 jtag = tag[j];
137
138 if (itag > jtag) {
139 if ((itag+jtag) % 2 == 0) continue;
140 } else if (itag < jtag) {
141 if ((itag+jtag) % 2 == 1) continue;
142 } else {
143 if (x[j].z < ztmp) continue;
144 if (x[j].z == ztmp && x[j].y < ytmp) continue;
145 if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
146 }
147
148 jtype = map[type[j]];
149
150 delx = xtmp - x[j].x;
151 dely = ytmp - x[j].y;
152 delz = ztmp - x[j].z;
153 rsq = delx*delx + dely*dely + delz*delz;
154
155 // shift rsq and store correction for force
156
157 if (shift_flag) {
158 double rsqtmp = rsq + shift*shift + 2*sqrt(rsq)*shift;
159 forceshiftfac = sqrt(rsqtmp/rsq);
160 rsq = rsqtmp;
161 }
162
163 iparam_ij = elem3param[itype][jtype][jtype];
164 if (rsq > params[iparam_ij].cutsq) continue;
165
166 repulsive(¶ms[iparam_ij],rsq,fpair,EFLAG,evdwl);
167
168 // correct force for shift in rsq
169
170 if (shift_flag) fpair *= forceshiftfac;
171
172 fxtmp += delx*fpair;
173 fytmp += dely*fpair;
174 fztmp += delz*fpair;
175 f[j].x -= delx*fpair;
176 f[j].y -= dely*fpair;
177 f[j].z -= delz*fpair;
178
179 if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
180 evdwl,0.0,fpair,delx,dely,delz,thr);
181 }
182
183 // three-body interactions
184 // skip immediately if I-J is not within cutoff
185 double fjxtmp,fjytmp,fjztmp;
186
187 for (jj = 0; jj < jnum; jj++) {
188 j = jlist[jj];
189 j &= NEIGHMASK;
190 jtype = map[type[j]];
191 iparam_ij = elem3param[itype][jtype][jtype];
192
193 delr1[0] = x[j].x - xtmp;
194 delr1[1] = x[j].y - ytmp;
195 delr1[2] = x[j].z - ztmp;
196 rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
197
198 if (shift_flag)
199 rsq1 += shift*shift + 2*sqrt(rsq1)*shift;
200
201 if (rsq1 > params[iparam_ij].cutsq) continue;
202
203 const double r1inv = 1.0/sqrt(dot3(delr1, delr1));
204 scale3(r1inv, delr1, r1_hat);
205
206 // accumulate bondorder zeta for each i-j interaction via loop over k
207
208 fjxtmp = fjytmp = fjztmp = 0.0;
209 zeta_ij = 0.0;
210
211 for (kk = 0; kk < jnum; kk++) {
212 if (jj == kk) continue;
213 k = jlist[kk];
214 k &= NEIGHMASK;
215 ktype = map[type[k]];
216 iparam_ijk = elem3param[itype][jtype][ktype];
217
218 delr2[0] = x[k].x - xtmp;
219 delr2[1] = x[k].y - ytmp;
220 delr2[2] = x[k].z - ztmp;
221 rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
222
223 if (shift_flag)
224 rsq2 += shift*shift + 2*sqrt(rsq2)*shift;
225
226 if (rsq2 > params[iparam_ijk].cutsq) continue;
227
228 const double r2inv = 1.0/sqrt(dot3(delr2, delr2));
229 scale3(r2inv, delr2, r2_hat);
230
231 zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,r1_hat,r2_hat);
232 }
233
234 // pairwise force due to zeta
235
236 force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fforce,prefactor,EFLAG,evdwl);
237
238 fpair = fforce*r1inv;
239
240 fxtmp += delr1[0]*fpair;
241 fytmp += delr1[1]*fpair;
242 fztmp += delr1[2]*fpair;
243 fjxtmp -= delr1[0]*fpair;
244 fjytmp -= delr1[1]*fpair;
245 fjztmp -= delr1[2]*fpair;
246
247 if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,evdwl,0.0,
248 -fpair,-delr1[0],-delr1[1],-delr1[2],thr);
249
250 // attractive term via loop over k
251
252 for (kk = 0; kk < jnum; kk++) {
253 if (jj == kk) continue;
254 k = jlist[kk];
255 k &= NEIGHMASK;
256 ktype = map[type[k]];
257 iparam_ijk = elem3param[itype][jtype][ktype];
258
259 delr2[0] = x[k].x - xtmp;
260 delr2[1] = x[k].y - ytmp;
261 delr2[2] = x[k].z - ztmp;
262 rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
263
264 if (shift_flag)
265 rsq2 += shift*shift + 2*sqrt(rsq2)*shift;
266
267 if (rsq2 > params[iparam_ijk].cutsq) continue;
268
269 const double r2inv = 1.0/sqrt(dot3(delr2, delr2));
270 scale3(r2inv, delr2, r2_hat);
271
272 attractive(¶ms[iparam_ijk],prefactor,
273 rsq1,rsq2,r1_hat,r2_hat,fi,fj,fk);
274
275 fxtmp += fi[0];
276 fytmp += fi[1];
277 fztmp += fi[2];
278 fjxtmp += fj[0];
279 fjytmp += fj[1];
280 fjztmp += fj[2];
281 f[k].x += fk[0];
282 f[k].y += fk[1];
283 f[k].z += fk[2];
284
285 if (VFLAG_EITHER) v_tally3_thr(this,i,j,k,fj,fk,delr1,delr2,thr);
286 }
287 f[j].x += fjxtmp;
288 f[j].y += fjytmp;
289 f[j].z += fjztmp;
290 }
291 f[i].x += fxtmp;
292 f[i].y += fytmp;
293 f[i].z += fztmp;
294 }
295 }
296
297 /* ---------------------------------------------------------------------- */
298
memory_usage()299 double PairTersoffMODCOMP::memory_usage()
300 {
301 double bytes = memory_usage_thr();
302 bytes += PairTersoffMOD::memory_usage();
303
304 return bytes;
305 }
306