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 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 #include "pair_lj_class2_coul_cut_kokkos.h"
16 #include <cmath>
17 #include <cstring>
18 #include "kokkos.h"
19 #include "atom_kokkos.h"
20 #include "comm.h"
21 #include "force.h"
22 #include "neighbor.h"
23 #include "neigh_list.h"
24 #include "neigh_request.h"
25 #include "update.h"
26 #include "respa.h"
27 #include "memory_kokkos.h"
28 #include "error.h"
29 #include "atom_masks.h"
30
31 using namespace LAMMPS_NS;
32
33 #define KOKKOS_CUDA_MAX_THREADS 256
34 #define KOKKOS_CUDA_MIN_BLOCKS 8
35
36 /* ---------------------------------------------------------------------- */
37
38 template<class DeviceType>
PairLJClass2CoulCutKokkos(LAMMPS * lmp)39 PairLJClass2CoulCutKokkos<DeviceType>::PairLJClass2CoulCutKokkos(LAMMPS *lmp):PairLJClass2CoulCut(lmp)
40 {
41 respa_enable = 0;
42
43 kokkosable = 1;
44 atomKK = (AtomKokkos *) atom;
45 execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
46 datamask_read = X_MASK | F_MASK | TYPE_MASK | Q_MASK | ENERGY_MASK | VIRIAL_MASK;
47 datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK;
48 }
49
50 /* ---------------------------------------------------------------------- */
51
52 template<class DeviceType>
~PairLJClass2CoulCutKokkos()53 PairLJClass2CoulCutKokkos<DeviceType>::~PairLJClass2CoulCutKokkos()
54 {
55 if (copymode) return;
56
57 if (allocated) {
58 memoryKK->destroy_kokkos(k_eatom,eatom);
59 memoryKK->destroy_kokkos(k_vatom,vatom);
60 memoryKK->destroy_kokkos(k_cutsq,cutsq);
61 memoryKK->destroy_kokkos(k_cut_ljsq,cut_ljsq);
62 memoryKK->destroy_kokkos(k_cut_coulsq,cut_coulsq);
63 }
64 }
65
66 /* ---------------------------------------------------------------------- */
67
68 template<class DeviceType>
compute(int eflag_in,int vflag_in)69 void PairLJClass2CoulCutKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
70 {
71 eflag = eflag_in;
72 vflag = vflag_in;
73
74 if (neighflag == FULL) no_virial_fdotr_compute = 1;
75
76 ev_init(eflag,vflag,0);
77
78 // reallocate per-atom arrays if necessary
79
80 if (eflag_atom) {
81 memoryKK->destroy_kokkos(k_eatom,eatom);
82 memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom");
83 d_eatom = k_eatom.view<DeviceType>();
84 }
85 if (vflag_atom) {
86 memoryKK->destroy_kokkos(k_vatom,vatom);
87 memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
88 d_vatom = k_vatom.view<DeviceType>();
89 }
90
91 atomKK->sync(execution_space,datamask_read);
92 k_cutsq.template sync<DeviceType>();
93 k_cut_ljsq.template sync<DeviceType>();
94 k_cut_coulsq.template sync<DeviceType>();
95 k_params.template sync<DeviceType>();
96 if (eflag || vflag) atomKK->modified(execution_space,datamask_modify);
97 else atomKK->modified(execution_space,F_MASK);
98
99 x = atomKK->k_x.view<DeviceType>();
100 c_x = atomKK->k_x.view<DeviceType>();
101 f = atomKK->k_f.view<DeviceType>();
102 q = atomKK->k_q.view<DeviceType>();
103 type = atomKK->k_type.view<DeviceType>();
104 nlocal = atom->nlocal;
105 nall = atom->nlocal + atom->nghost;
106 special_lj[0] = force->special_lj[0];
107 special_lj[1] = force->special_lj[1];
108 special_lj[2] = force->special_lj[2];
109 special_lj[3] = force->special_lj[3];
110 special_coul[0] = force->special_coul[0];
111 special_coul[1] = force->special_coul[1];
112 special_coul[2] = force->special_coul[2];
113 special_coul[3] = force->special_coul[3];
114 qqrd2e = force->qqrd2e;
115 newton_pair = force->newton_pair;
116
117 // loop over neighbors of my atoms
118
119 copymode = 1;
120
121 EV_FLOAT ev = pair_compute<PairLJClass2CoulCutKokkos<DeviceType>,void >
122 (this,(NeighListKokkos<DeviceType>*)list);
123
124 if (eflag) {
125 eng_vdwl += ev.evdwl;
126 eng_coul += ev.ecoul;
127 }
128 if (vflag_global) {
129 virial[0] += ev.v[0];
130 virial[1] += ev.v[1];
131 virial[2] += ev.v[2];
132 virial[3] += ev.v[3];
133 virial[4] += ev.v[4];
134 virial[5] += ev.v[5];
135 }
136
137 if (eflag_atom) {
138 k_eatom.template modify<DeviceType>();
139 k_eatom.template sync<LMPHostType>();
140 }
141
142 if (vflag_atom) {
143 k_vatom.template modify<DeviceType>();
144 k_vatom.template sync<LMPHostType>();
145 }
146
147 if (vflag_fdotr) pair_virial_fdotr_compute(this);
148
149 copymode = 0;
150 }
151
152 /* ----------------------------------------------------------------------
153 compute LJ 12-6 pair force between atoms i and j
154 ---------------------------------------------------------------------- */
155 template<class DeviceType>
156 template<bool STACKPARAMS, class Specialisation>
157 KOKKOS_INLINE_FUNCTION
158 F_FLOAT PairLJClass2CoulCutKokkos<DeviceType>::
compute_fpair(const F_FLOAT & rsq,const int & i,const int & j,const int & itype,const int & jtype) const159 compute_fpair(const F_FLOAT& rsq, const int& i, const int&j,
160 const int& itype, const int& jtype) const {
161 (void) i;
162 (void) j;
163 const F_FLOAT r2inv = 1.0/rsq;
164 const F_FLOAT rinv = sqrt(r2inv);
165 const F_FLOAT r3inv = r2inv*rinv;
166 const F_FLOAT r6inv = r3inv*r3inv;
167
168 const F_FLOAT forcelj = r6inv *
169 ((STACKPARAMS?m_params[itype][jtype].lj1:params(itype,jtype).lj1)*r3inv -
170 (STACKPARAMS?m_params[itype][jtype].lj2:params(itype,jtype).lj2));
171
172 return forcelj*r2inv;
173 }
174
175 /* ----------------------------------------------------------------------
176 compute coulomb pair force between atoms i and j
177 ---------------------------------------------------------------------- */
178 template<class DeviceType>
179 template<bool STACKPARAMS, class Specialisation>
180 KOKKOS_INLINE_FUNCTION
181 F_FLOAT PairLJClass2CoulCutKokkos<DeviceType>::
compute_fcoul(const F_FLOAT & rsq,const int &,const int & j,const int &,const int &,const F_FLOAT & factor_coul,const F_FLOAT & qtmp) const182 compute_fcoul(const F_FLOAT& rsq, const int& /*i*/, const int&j,
183 const int& /*itype*/, const int& /*jtype*/,
184 const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const {
185 const F_FLOAT r2inv = 1.0/rsq;
186 const F_FLOAT rinv = sqrt(r2inv);
187 F_FLOAT forcecoul;
188
189 forcecoul = qqrd2e*qtmp*q(j) *rinv;
190
191 return factor_coul*forcecoul*r2inv;
192 }
193
194 /* ----------------------------------------------------------------------
195 compute LJ 12-6 pair potential energy between atoms i and j
196 ---------------------------------------------------------------------- */
197 template<class DeviceType>
198 template<bool STACKPARAMS, class Specialisation>
199 KOKKOS_INLINE_FUNCTION
200 F_FLOAT PairLJClass2CoulCutKokkos<DeviceType>::
compute_evdwl(const F_FLOAT & rsq,const int & i,const int & j,const int & itype,const int & jtype) const201 compute_evdwl(const F_FLOAT& rsq, const int& i, const int&j,
202 const int& itype, const int& jtype) const {
203 (void) i;
204 (void) j;
205 const F_FLOAT r2inv = 1.0/rsq;
206 const F_FLOAT rinv = sqrt(r2inv);
207 const F_FLOAT r3inv = r2inv*rinv;
208 const F_FLOAT r6inv = r3inv*r3inv;
209
210 return r6inv*((STACKPARAMS?m_params[itype][jtype].lj3:params(itype,jtype).lj3)*r3inv -
211 (STACKPARAMS?m_params[itype][jtype].lj4:params(itype,jtype).lj4)) -
212 (STACKPARAMS?m_params[itype][jtype].offset:params(itype,jtype).offset);
213 }
214
215 /* ----------------------------------------------------------------------
216 compute coulomb pair potential energy between atoms i and j
217 ---------------------------------------------------------------------- */
218 template<class DeviceType>
219 template<bool STACKPARAMS, class Specialisation>
220 KOKKOS_INLINE_FUNCTION
221 F_FLOAT PairLJClass2CoulCutKokkos<DeviceType>::
compute_ecoul(const F_FLOAT & rsq,const int &,const int & j,const int &,const int &,const F_FLOAT & factor_coul,const F_FLOAT & qtmp) const222 compute_ecoul(const F_FLOAT& rsq, const int& /*i*/, const int&j,
223 const int& /*itype*/, const int& /*jtype*/,
224 const F_FLOAT& factor_coul, const F_FLOAT& qtmp) const {
225 const F_FLOAT r2inv = 1.0/rsq;
226 const F_FLOAT rinv = sqrt(r2inv);
227
228 return factor_coul*qqrd2e*qtmp*q(j)*rinv;
229
230 }
231
232 /* ----------------------------------------------------------------------
233 allocate all arrays
234 ------------------------------------------------------------------------- */
235
236 template<class DeviceType>
allocate()237 void PairLJClass2CoulCutKokkos<DeviceType>::allocate()
238 {
239 PairLJClass2CoulCut::allocate();
240
241 int n = atom->ntypes;
242 memory->destroy(cutsq);
243 memoryKK->create_kokkos(k_cutsq,cutsq,n+1,n+1,"pair:cutsq");
244 d_cutsq = k_cutsq.template view<DeviceType>();
245 memory->destroy(cut_ljsq);
246 memoryKK->create_kokkos(k_cut_ljsq,cut_ljsq,n+1,n+1,"pair:cut_ljsq");
247 d_cut_ljsq = k_cut_ljsq.template view<DeviceType>();
248 memory->destroy(cut_coulsq);
249 memoryKK->create_kokkos(k_cut_coulsq,cut_coulsq,n+1,n+1,"pair:cut_coulsq");
250 d_cut_coulsq = k_cut_coulsq.template view<DeviceType>();
251 k_params = Kokkos::DualView<params_lj_coul**,Kokkos::LayoutRight,DeviceType>("PairLJClass2CoulCut::params",n+1,n+1);
252 params = k_params.template view<DeviceType>();
253 }
254
255 /* ----------------------------------------------------------------------
256 global settings
257 ------------------------------------------------------------------------- */
258
259 template<class DeviceType>
settings(int narg,char ** arg)260 void PairLJClass2CoulCutKokkos<DeviceType>::settings(int narg, char **arg)
261 {
262 if (narg > 2) error->all(FLERR,"Illegal pair_style command");
263
264 PairLJClass2CoulCut::settings(1,arg);
265 }
266
267 /* ----------------------------------------------------------------------
268 init specific to this pair style
269 ------------------------------------------------------------------------- */
270
271 template<class DeviceType>
init_style()272 void PairLJClass2CoulCutKokkos<DeviceType>::init_style()
273 {
274 PairLJClass2CoulCut::init_style();
275
276 // error if rRESPA with inner levels
277
278 if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
279 int respa = 0;
280 if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
281 if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
282 if (respa)
283 error->all(FLERR,"Cannot use Kokkos pair style with rRESPA inner/middle");
284 }
285
286 // irequest = neigh request made by parent class
287
288 neighflag = lmp->kokkos->neighflag;
289 int irequest = neighbor->nrequest - 1;
290
291 neighbor->requests[irequest]->
292 kokkos_host = std::is_same<DeviceType,LMPHostType>::value &&
293 !std::is_same<DeviceType,LMPDeviceType>::value;
294 neighbor->requests[irequest]->
295 kokkos_device = std::is_same<DeviceType,LMPDeviceType>::value;
296
297 if (neighflag == FULL) {
298 neighbor->requests[irequest]->full = 1;
299 neighbor->requests[irequest]->half = 0;
300 } else if (neighflag == HALF || neighflag == HALFTHREAD) {
301 neighbor->requests[irequest]->full = 0;
302 neighbor->requests[irequest]->half = 1;
303 } else {
304 error->all(FLERR,"Cannot use chosen neighbor list style with lj/class2/coul/cut/kk");
305 }
306 }
307
308 /* ----------------------------------------------------------------------
309 init for one type pair i,j and corresponding j,i
310 ------------------------------------------------------------------------- */
311
312 template<class DeviceType>
init_one(int i,int j)313 double PairLJClass2CoulCutKokkos<DeviceType>::init_one(int i, int j)
314 {
315 double cutone = PairLJClass2CoulCut::init_one(i,j);
316 double cut_ljsqm = cut_ljsq[i][j];
317 double cut_coulsqm = cut_coulsq[i][j];
318
319 k_params.h_view(i,j).lj1 = lj1[i][j];
320 k_params.h_view(i,j).lj2 = lj2[i][j];
321 k_params.h_view(i,j).lj3 = lj3[i][j];
322 k_params.h_view(i,j).lj4 = lj4[i][j];
323 k_params.h_view(i,j).offset = offset[i][j];
324 k_params.h_view(i,j).cut_ljsq = cut_ljsqm;
325 k_params.h_view(i,j).cut_coulsq = cut_coulsqm;
326
327 k_params.h_view(j,i) = k_params.h_view(i,j);
328 if (i<MAX_TYPES_STACKPARAMS+1 && j<MAX_TYPES_STACKPARAMS+1) {
329 m_params[i][j] = m_params[j][i] = k_params.h_view(i,j);
330 m_cutsq[j][i] = m_cutsq[i][j] = cutone*cutone;
331 m_cut_ljsq[j][i] = m_cut_ljsq[i][j] = cut_ljsqm;
332 m_cut_coulsq[j][i] = m_cut_coulsq[i][j] = cut_coulsqm;
333 }
334
335 k_cutsq.h_view(i,j) = k_cutsq.h_view(j,i) = cutone*cutone;
336 k_cutsq.template modify<LMPHostType>();
337 k_cut_ljsq.h_view(i,j) = k_cut_ljsq.h_view(j,i) = cut_ljsqm;
338 k_cut_ljsq.template modify<LMPHostType>();
339 k_cut_coulsq.h_view(i,j) = k_cut_coulsq.h_view(j,i) = cut_coulsqm;
340 k_cut_coulsq.template modify<LMPHostType>();
341 k_params.template modify<LMPHostType>();
342
343 return cutone;
344 }
345
346
347 namespace LAMMPS_NS {
348 template class PairLJClass2CoulCutKokkos<LMPDeviceType>;
349 #ifdef LMP_KOKKOS_GPU
350 template class PairLJClass2CoulCutKokkos<LMPHostType>;
351 #endif
352 }
353
354