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