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_vashishta_table_omp.h"
17
18 #include "atom.h"
19 #include "comm.h"
20 #include "memory.h"
21 #include "neigh_list.h"
22 #include "suffix.h"
23
24 #include "omp_compat.h"
25 using namespace LAMMPS_NS;
26
27 /* ---------------------------------------------------------------------- */
28
PairVashishtaTableOMP(LAMMPS * lmp)29 PairVashishtaTableOMP::PairVashishtaTableOMP(LAMMPS *lmp) :
30 PairVashishtaTable(lmp), ThrOMP(lmp, THR_PAIR)
31 {
32 suffix_flag |= Suffix::OMP;
33 respa_enable = 0;
34 }
35
36 /* ---------------------------------------------------------------------- */
37
compute(int eflag,int vflag)38 void PairVashishtaTableOMP::compute(int eflag, int vflag)
39 {
40 ev_init(eflag,vflag);
41
42 const int nall = atom->nlocal + atom->nghost;
43 const int nthreads = comm->nthreads;
44 const int inum = list->inum;
45
46 #if defined(_OPENMP)
47 #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
48 #endif
49 {
50 int ifrom, ito, tid;
51
52 loop_setup_thr(ifrom, ito, tid, inum, nthreads);
53 ThrData *thr = fix->get_thr(tid);
54 thr->timer(Timer::START);
55 ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
56
57 if (evflag) {
58 if (eflag) {
59 eval<1,1>(ifrom, ito, thr);
60 } else {
61 eval<1,0>(ifrom, ito, thr);
62 }
63 } else eval<0,0>(ifrom, ito, thr);
64
65 thr->timer(Timer::PAIR);
66 reduce_thr(this, eflag, vflag, thr);
67 } // end of omp parallel region
68 }
69
70 template <int EVFLAG, int EFLAG>
eval(int iifrom,int iito,ThrData * const thr)71 void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
72 {
73 int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
74 tagint itag,jtag;
75 int itype,jtype,ktype,ijparam,ikparam,ijkparam;
76 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
77 double rsq,rsq1,rsq2;
78 double delr1[3],delr2[3],fj[3],fk[3];
79 int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
80
81 evdwl = 0.0;
82
83 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
84 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
85 const tagint * _noalias const tag = atom->tag;
86 const int * _noalias const type = atom->type;
87 const int nlocal = atom->nlocal;
88 const double cutshortsq = r0max*r0max;
89
90 ilist = list->ilist;
91 numneigh = list->numneigh;
92 firstneigh = list->firstneigh;
93 maxshort_thr = maxshort;
94 memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
95
96 double fxtmp,fytmp,fztmp;
97
98 // loop over full neighbor list of my atoms
99
100 for (ii = iifrom; ii < iito; ++ii) {
101
102 i = ilist[ii];
103 itag = tag[i];
104 itype = map[type[i]];
105 xtmp = x[i].x;
106 ytmp = x[i].y;
107 ztmp = x[i].z;
108 fxtmp = fytmp = fztmp = 0.0;
109
110 // two-body interactions, skip half of them
111
112 jlist = firstneigh[i];
113 jnum = numneigh[i];
114 int numshort = 0;
115
116 for (jj = 0; jj < jnum; jj++) {
117 j = jlist[jj];
118 j &= NEIGHMASK;
119
120 delx = xtmp - x[j].x;
121 dely = ytmp - x[j].y;
122 delz = ztmp - x[j].z;
123 rsq = delx*delx + dely*dely + delz*delz;
124
125 if (rsq < cutshortsq) {
126 neighshort_thr[numshort++] = j;
127 if (numshort >= maxshort_thr) {
128 maxshort_thr += maxshort_thr/2;
129 memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
130 }
131 }
132
133 jtag = tag[j];
134 if (itag > jtag) {
135 if ((itag+jtag) % 2 == 0) continue;
136 } else if (itag < jtag) {
137 if ((itag+jtag) % 2 == 1) continue;
138 } else {
139 if (x[j].z < ztmp) continue;
140 if (x[j].z == ztmp && x[j].y < ytmp) continue;
141 if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
142 }
143
144 jtype = map[type[j]];
145 ijparam = elem3param[itype][jtype][jtype];
146 if (rsq >= params[ijparam].cutsq) continue;
147
148 twobody_table(params[ijparam],rsq,fpair,EFLAG,evdwl);
149
150 fxtmp += delx*fpair;
151 fytmp += dely*fpair;
152 fztmp += delz*fpair;
153 f[j].x -= delx*fpair;
154 f[j].y -= dely*fpair;
155 f[j].z -= delz*fpair;
156
157 if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
158 evdwl,0.0,fpair,delx,dely,delz,thr);
159 }
160
161 jnumm1 = numshort - 1;
162
163 for (jj = 0; jj < jnumm1; jj++) {
164 j = neighshort_thr[jj];
165 jtype = map[type[j]];
166 ijparam = elem3param[itype][jtype][jtype];
167 delr1[0] = x[j].x - xtmp;
168 delr1[1] = x[j].y - ytmp;
169 delr1[2] = x[j].z - ztmp;
170 rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
171 if (rsq1 >= params[ijparam].cutsq2) continue;
172
173 double fjxtmp,fjytmp,fjztmp;
174 fjxtmp = fjytmp = fjztmp = 0.0;
175
176 for (kk = jj+1; kk < numshort; kk++) {
177 k = neighshort_thr[kk];
178 ktype = map[type[k]];
179 ikparam = elem3param[itype][ktype][ktype];
180 ijkparam = elem3param[itype][jtype][ktype];
181
182 delr2[0] = x[k].x - xtmp;
183 delr2[1] = x[k].y - ytmp;
184 delr2[2] = x[k].z - ztmp;
185 rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
186 if (rsq2 >= params[ikparam].cutsq2) continue;
187
188 threebody(¶ms[ijparam],¶ms[ikparam],¶ms[ijkparam],
189 rsq1,rsq2,delr1,delr2,fj,fk,EFLAG,evdwl);
190
191 fxtmp -= fj[0] + fk[0];
192 fytmp -= fj[1] + fk[1];
193 fztmp -= fj[2] + fk[2];
194 fjxtmp += fj[0];
195 fjytmp += fj[1];
196 fjztmp += fj[2];
197 f[k].x += fk[0];
198 f[k].y += fk[1];
199 f[k].z += fk[2];
200
201 if (EVFLAG) ev_tally3_thr(this,i,j,k,evdwl,0.0,fj,fk,delr1,delr2,thr);
202 }
203 f[j].x += fjxtmp;
204 f[j].y += fjytmp;
205 f[j].z += fjztmp;
206 }
207 f[i].x += fxtmp;
208 f[i].y += fytmp;
209 f[i].z += fztmp;
210 }
211 memory->destroy(neighshort_thr);
212 }
213
214 /* ---------------------------------------------------------------------- */
215
memory_usage()216 double PairVashishtaTableOMP::memory_usage()
217 {
218 double bytes = memory_usage_thr();
219 bytes += PairVashishtaTable::memory_usage();
220
221 return bytes;
222 }
223