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