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 This style is a simplified re-implementation of the CG/CMM pair style
15 ------------------------------------------------------------------------- */
16
17 #include "pair_lj_sdk_omp.h"
18 #include "lj_sdk_common.h"
19
20 #include "atom.h"
21 #include "comm.h"
22 #include "force.h"
23 #include "neigh_list.h"
24 #include "suffix.h"
25
26 #include <cmath>
27
28 #include "omp_compat.h"
29 using namespace LAMMPS_NS;
30 using namespace LJSDKParms;
31
32 /* ---------------------------------------------------------------------- */
33
PairLJSDKOMP(LAMMPS * lmp)34 PairLJSDKOMP::PairLJSDKOMP(LAMMPS *lmp) :
35 PairLJSDK(lmp), ThrOMP(lmp, THR_PAIR)
36 {
37 suffix_flag |= Suffix::OMP;
38 respa_enable = 0;
39 }
40
41 /* ---------------------------------------------------------------------- */
42
compute(int eflag,int vflag)43 void PairLJSDKOMP::compute(int eflag, int vflag)
44 {
45 ev_init(eflag,vflag);
46
47 const int nall = atom->nlocal + atom->nghost;
48 const int nthreads = comm->nthreads;
49 const int inum = list->inum;
50
51 #if defined(_OPENMP)
52 #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
53 #endif
54 {
55 int ifrom, ito, tid;
56
57 loop_setup_thr(ifrom, ito, tid, inum, nthreads);
58 ThrData *thr = fix->get_thr(tid);
59 thr->timer(Timer::START);
60 ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
61
62 if (evflag) {
63 if (eflag) {
64 if (force->newton_pair) eval_thr<1,1,1>(ifrom, ito, thr);
65 else eval_thr<1,1,0>(ifrom, ito, thr);
66 } else {
67 if (force->newton_pair) eval_thr<1,0,1>(ifrom, ito, thr);
68 else eval_thr<1,0,0>(ifrom, ito, thr);
69 }
70 } else {
71 if (force->newton_pair) eval_thr<0,0,1>(ifrom, ito, thr);
72 else eval_thr<0,0,0>(ifrom, ito, thr);
73 }
74
75 thr->timer(Timer::PAIR);
76 reduce_thr(this, eflag, vflag, thr);
77 } // end of omp parallel region
78 }
79
80 /* ---------------------------------------------------------------------- */
81
82 template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
eval_thr(int iifrom,int iito,ThrData * const thr)83 void PairLJSDKOMP::eval_thr(int iifrom, int iito, ThrData * const thr)
84 {
85 int i,j,ii,jj,jtype;
86 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
87 double rsq,r2inv,forcelj,factor_lj;
88
89 evdwl = 0.0;
90
91 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
92 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
93 const int * _noalias const type = atom->type;
94 const int nlocal = atom->nlocal;
95 const double * _noalias const special_lj = force->special_lj;
96 double fxtmp,fytmp,fztmp;
97
98 const int * const ilist = list->ilist;
99 const int * const numneigh = list->numneigh;
100 const int * const * const firstneigh = list->firstneigh;
101
102 // loop over neighbors of my atoms
103
104 for (ii = iifrom; ii < iito; ++ii) {
105
106 i = ilist[ii];
107 xtmp = x[i].x;
108 ytmp = x[i].y;
109 ztmp = x[i].z;
110 fxtmp=fytmp=fztmp=0.0;
111
112 const int itype = type[i];
113 const int * const jlist = firstneigh[i];
114 const int jnum = numneigh[i];
115
116 for (jj = 0; jj < jnum; jj++) {
117 j = jlist[jj];
118 factor_lj = special_lj[sbmask(j)];
119 j &= NEIGHMASK;
120
121 delx = xtmp - x[j].x;
122 dely = ytmp - x[j].y;
123 delz = ztmp - x[j].z;
124 rsq = delx*delx + dely*dely + delz*delz;
125 jtype = type[j];
126
127 if (rsq < cutsq[itype][jtype]) {
128 r2inv = 1.0/rsq;
129 const int ljt = lj_type[itype][jtype];
130
131 if (ljt == LJ12_4) {
132 const double r4inv=r2inv*r2inv;
133 forcelj = r4inv*(lj1[itype][jtype]*r4inv*r4inv
134 - lj2[itype][jtype]);
135
136 if (EFLAG)
137 evdwl = r4inv*(lj3[itype][jtype]*r4inv*r4inv
138 - lj4[itype][jtype]) - offset[itype][jtype];
139
140 } else if (ljt == LJ9_6) {
141 const double r3inv = r2inv*sqrt(r2inv);
142 const double r6inv = r3inv*r3inv;
143 forcelj = r6inv*(lj1[itype][jtype]*r3inv
144 - lj2[itype][jtype]);
145 if (EFLAG)
146 evdwl = r6inv*(lj3[itype][jtype]*r3inv
147 - lj4[itype][jtype]) - offset[itype][jtype];
148
149 } else if (ljt == LJ12_6) {
150 const double r6inv = r2inv*r2inv*r2inv;
151 forcelj = r6inv*(lj1[itype][jtype]*r6inv
152 - lj2[itype][jtype]);
153 if (EFLAG)
154 evdwl = r6inv*(lj3[itype][jtype]*r6inv
155 - lj4[itype][jtype]) - offset[itype][jtype];
156 } else continue;
157
158 fpair = factor_lj*forcelj*r2inv;
159
160 fxtmp += delx*fpair;
161 fytmp += dely*fpair;
162 fztmp += delz*fpair;
163 if (NEWTON_PAIR || j < nlocal) {
164 f[j].x -= delx*fpair;
165 f[j].y -= dely*fpair;
166 f[j].z -= delz*fpair;
167 }
168
169 if (EFLAG) evdwl *= factor_lj;
170 if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
171 evdwl,0.0,fpair,delx,dely,delz,thr);
172 }
173 }
174 f[i].x += fxtmp;
175 f[i].y += fytmp;
176 f[i].z += fztmp;
177 }
178 }
179
180 /* ---------------------------------------------------------------------- */
181
memory_usage()182 double PairLJSDKOMP::memory_usage()
183 {
184 double bytes = memory_usage_thr();
185 bytes += PairLJSDK::memory_usage();
186
187 return bytes;
188 }
189