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_zbl_omp.h"
17
18 #include "atom.h"
19 #include "comm.h"
20 #include "force.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
29 /* ---------------------------------------------------------------------- */
30
PairZBLOMP(LAMMPS * lmp)31 PairZBLOMP::PairZBLOMP(LAMMPS *lmp) :
32 PairZBL(lmp), ThrOMP(lmp, THR_PAIR)
33 {
34 suffix_flag |= Suffix::OMP;
35 respa_enable = 0;
36 }
37
38 /* ---------------------------------------------------------------------- */
39
compute(int eflag,int vflag)40 void PairZBLOMP::compute(int eflag, int vflag)
41 {
42 ev_init(eflag,vflag);
43
44 const int nall = atom->nlocal + atom->nghost;
45 const int nthreads = comm->nthreads;
46 const int inum = list->inum;
47
48 #if defined(_OPENMP)
49 #pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag,vflag)
50 #endif
51 {
52 int ifrom, ito, tid;
53
54 loop_setup_thr(ifrom, ito, tid, inum, nthreads);
55 ThrData *thr = fix->get_thr(tid);
56 thr->timer(Timer::START);
57 ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
58
59 if (evflag) {
60 if (eflag) {
61 if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
62 else eval<1,1,0>(ifrom, ito, thr);
63 } else {
64 if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
65 else eval<1,0,0>(ifrom, ito, thr);
66 }
67 } else {
68 if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
69 else eval<0,0,0>(ifrom, ito, thr);
70 }
71 thr->timer(Timer::PAIR);
72 reduce_thr(this, eflag, vflag, thr);
73 } // end of omp parallel region
74 }
75
76 template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
eval(int iifrom,int iito,ThrData * const thr)77 void PairZBLOMP::eval(int iifrom, int iito, ThrData * const thr)
78 {
79 const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
80 dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
81 const int * _noalias const type = atom->type;
82 const int * _noalias const ilist = list->ilist;
83 const int * _noalias const numneigh = list->numneigh;
84 const int * const * const firstneigh = list->firstneigh;
85
86 double xtmp,ytmp,ztmp,delx,dely,delz,fxtmp,fytmp,fztmp;
87 double rsq,t,fswitch,eswitch,evdwl,fpair;
88
89 const int nlocal = atom->nlocal;
90 int j,jj,jnum,jtype;
91
92 evdwl = 0.0;
93
94 // loop over neighbors of my atoms
95
96 for (int ii = iifrom; ii < iito; ++ii) {
97 const int i = ilist[ii];
98 const int itype = type[i];
99 const int * _noalias const jlist = firstneigh[i];
100 const double * _noalias const sw1i = sw1[itype];
101 const double * _noalias const sw2i = sw2[itype];
102 const double * _noalias const sw3i = sw3[itype];
103 const double * _noalias const sw4i = sw4[itype];
104 const double * _noalias const sw5i = sw5[itype];
105
106 xtmp = x[i].x;
107 ytmp = x[i].y;
108 ztmp = x[i].z;
109 jnum = numneigh[i];
110 fxtmp=fytmp=fztmp=0.0;
111
112 for (jj = 0; jj < jnum; jj++) {
113 j = jlist[jj];
114 j &= NEIGHMASK;
115
116 delx = xtmp - x[j].x;
117 dely = ytmp - x[j].y;
118 delz = ztmp - x[j].z;
119 rsq = delx*delx + dely*dely + delz*delz;
120 jtype = type[j];
121
122 if (rsq < cut_globalsq) {
123 const double r = sqrt(rsq);
124 fpair = dzbldr(r, itype, jtype);
125
126 if (r > cut_inner) {
127 t = r - cut_inner;
128 fswitch = t*t *
129 (sw1i[jtype] + sw2i[jtype]*t);
130 fpair += fswitch;
131 }
132
133 fpair *= -1.0/r;
134 fxtmp += delx*fpair;
135 fytmp += dely*fpair;
136 fztmp += delz*fpair;
137
138 if (NEWTON_PAIR || j < nlocal) {
139 f[j].x -= delx*fpair;
140 f[j].y -= dely*fpair;
141 f[j].z -= delz*fpair;
142 }
143
144 if (EFLAG) {
145 evdwl = e_zbl(r, itype, jtype);
146 evdwl += sw5i[jtype];
147 if (r > cut_inner) {
148 eswitch = t*t*t *
149 (sw3i[jtype] + sw4i[jtype]*t);
150 evdwl += eswitch;
151 }
152 }
153
154 if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
155 evdwl,0.0,fpair,delx,dely,delz,thr);
156 }
157 }
158 f[i].x += fxtmp;
159 f[i].y += fytmp;
160 f[i].z += fztmp;
161 }
162 }
163
164 /* ---------------------------------------------------------------------- */
165
memory_usage()166 double PairZBLOMP::memory_usage()
167 {
168 double bytes = memory_usage_thr();
169 bytes += PairZBL::memory_usage();
170
171 return bytes;
172 }
173