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 /* ----------------------------------------------------------------------
16 Contributing author: Sergey Lishchuk
17 ------------------------------------------------------------------------- */
18
19 #include "pair_atm.h"
20
21 #include <cmath>
22 #include "atom.h"
23 #include "citeme.h"
24 #include "comm.h"
25 #include "error.h"
26 #include "force.h"
27 #include "memory.h"
28 #include "neigh_list.h"
29 #include "neigh_request.h"
30 #include "neighbor.h"
31
32
33 using namespace LAMMPS_NS;
34
35 static const char cite_atm_package[] =
36 "ATM package:\n\n"
37 "@Article{Lishchuk:2012:164501,\n"
38 " author = {S. V. Lishchuk},\n"
39 " title = {Role of three-body interactions in formation of bulk viscosity in liquid argon},\n"
40 " journal = {J.~Chem.~Phys.},\n"
41 " year = 2012,\n"
42 " volume = 136,\n"
43 " pages = {164501}\n"
44 "}\n\n";
45
46 /* ---------------------------------------------------------------------- */
47
PairATM(LAMMPS * lmp)48 PairATM::PairATM(LAMMPS *lmp) : Pair(lmp)
49 {
50 if (lmp->citeme) lmp->citeme->add(cite_atm_package);
51
52 single_enable = 0;
53 restartinfo = 1;
54 one_coeff = 0;
55 manybody_flag = 1;
56 centroidstressflag = CENTROID_NOTAVAIL;
57 }
58
59 /* ----------------------------------------------------------------------
60 check if allocated, since class can be destructed when incomplete
61 ------------------------------------------------------------------------- */
62
~PairATM()63 PairATM::~PairATM()
64 {
65 if (copymode) return;
66
67 if (allocated) {
68 memory->destroy(setflag);
69 memory->destroy(cutsq);
70
71 memory->destroy(nu);
72 }
73 }
74
75 /* ----------------------------------------------------------------------
76 workhorse routine that computes pairwise interactions
77 ------------------------------------------------------------------------- */
78
compute(int eflag,int vflag)79 void PairATM::compute(int eflag, int vflag)
80 {
81 int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
82 double xi,yi,zi,evdwl;
83 double rij2,rik2,rjk2;
84 double rij[3],rik[3],rjk[3],fj[3],fk[3];
85 double nu_local;
86 int *ilist,*jlist,*numneigh,**firstneigh;
87
88 evdwl = 0.0;
89 ev_init(eflag,vflag);
90
91 double **x = atom->x;
92 double **f = atom->f;
93 int *type = atom->type;
94
95 double cutoff_squared = cut_global*cut_global;
96 double triple = cut_triple*cut_triple*cut_triple;
97 double cutoff_triple_sixth = triple*triple;
98
99 inum = list->inum;
100 ilist = list->ilist;
101 numneigh = list->numneigh;
102 firstneigh = list->firstneigh;
103
104 // triple loop over local atoms and neighbors twice
105 // must compute each IJK triplet interaction exactly once
106 // by proc that owns the triplet atom with smallest x coord
107 // special logic to break ties if multiple atoms have same x or y coords
108 // inner two loops for jj=1,Jnum and kk=jj+1,Jnum insure
109 // the pair of other 2 non-minimum-x atoms is only considered once
110 // triplet geometry criteria for calculation:
111 // each pair distance <= cutoff
112 // produce of 3 pair distances <= cutoff_triple^3
113
114 for (ii = 0; ii < inum; ii++) {
115 i = ilist[ii];
116 xi = x[i][0];
117 yi = x[i][1];
118 zi = x[i][2];
119
120 jlist = firstneigh[i];
121 jnum = numneigh[i];
122 jnumm1 = jnum - 1;
123
124 for (jj = 0; jj < jnumm1; jj++) {
125 j = jlist[jj];
126 j &= NEIGHMASK;
127
128 rij[0] = x[j][0] - xi;
129 if (rij[0] < 0.0) continue;
130 rij[1] = x[j][1] - yi;
131 if (rij[0] == 0.0 and rij[1] < 0.0) continue;
132 rij[2] = x[j][2] - zi;
133 if (rij[0] == 0.0 and rij[1] == 0.0 and rij[2] < 0.0) continue;
134 rij2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
135 if (rij2 > cutoff_squared) continue;
136
137 for (kk = jj+1; kk < jnum; kk++) {
138 k = jlist[kk];
139 k &= NEIGHMASK;
140
141 rik[0] = x[k][0] - xi;
142 if (rik[0] < 0.0) continue;
143 rik[1] = x[k][1] - yi;
144 if (rik[0] == 0.0 and rik[1] < 0.0) continue;
145 rik[2] = x[k][2] - zi;
146 if (rik[0] == 0.0 and rik[1] == 0.0 and rik[2] < 0.0) continue;
147 rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2];
148 if (rik2 > cutoff_squared) continue;
149
150 rjk[0] = x[k][0] - x[j][0];
151 rjk[1] = x[k][1] - x[j][1];
152 rjk[2] = x[k][2] - x[j][2];
153 rjk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2];
154 if (rjk2 > cutoff_squared) continue;
155
156 double r6 = rij2*rjk2*rik2;
157 if (r6 > cutoff_triple_sixth) continue;
158
159 nu_local = nu[type[i]][type[j]][type[k]];
160 if (nu_local == 0.0) continue;
161
162 interaction_ddd(nu_local,
163 r6,rij2,rik2,rjk2,rij,rik,rjk,fj,fk,eflag,evdwl);
164
165 f[i][0] -= fj[0] + fk[0];
166 f[i][1] -= fj[1] + fk[1];
167 f[i][2] -= fj[2] + fk[2];
168 f[j][0] += fj[0];
169 f[j][1] += fj[1];
170 f[j][2] += fj[2];
171 f[k][0] += fk[0];
172 f[k][1] += fk[1];
173 f[k][2] += fk[2];
174
175 if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,rij,rik);
176 }
177 }
178 }
179
180 if (vflag_fdotr) virial_fdotr_compute();
181 }
182
183 /* ---------------------------------------------------------------------- */
184
allocate()185 void PairATM::allocate()
186 {
187 allocated = 1;
188 int n = atom->ntypes;
189
190 memory->create(setflag,n+1,n+1,"pair:setflag");
191 for (int i = 1; i <= n; i++)
192 for (int j = i; j <= n; j++)
193 setflag[i][j] = 0;
194
195 memory->create(cutsq,n+1,n+1,"pair:cutsq");
196
197 memory->create(nu,n+1,n+1,n+1,"pair:nu");
198
199 // initialize all nu values to 0.0
200
201 for (int i = 1; i <= n; i++)
202 for (int j = 1; j <= n; j++)
203 for (int k = 1; k <= n; k++)
204 nu[i][j][k] = 0.0;
205 }
206
207 /* ----------------------------------------------------------------------
208 global settings
209 ------------------------------------------------------------------------- */
210
settings(int narg,char ** arg)211 void PairATM::settings(int narg, char **arg)
212 {
213 if (narg != 2) error->all(FLERR,"Illegal pair_style command");
214
215 cut_global = utils::numeric(FLERR,arg[0],false,lmp);
216 cut_triple = utils::numeric(FLERR,arg[1],false,lmp);
217 }
218
219 /* ----------------------------------------------------------------------
220 set coefficients for one I,J,K type triplet
221 ------------------------------------------------------------------------- */
222
coeff(int narg,char ** arg)223 void PairATM::coeff(int narg, char **arg)
224 {
225 if (narg != 4) error->all(FLERR,"Incorrect args for pair coefficients");
226 if (!allocated) allocate();
227
228 int ilo,ihi,jlo,jhi,klo,khi;
229 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
230 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
231 utils::bounds(FLERR,arg[2],1,atom->ntypes,klo,khi,error);
232
233 double nu_one = utils::numeric(FLERR,arg[3],false,lmp);
234
235 int count = 0;
236 for (int i = ilo; i <= ihi; i++) {
237 for (int j = MAX(jlo,i); j<=jhi; j++) {
238 for (int k = MAX(klo,j); k<=khi; k++) {
239 nu[i][j][k] = nu_one;
240 count++;
241 }
242 setflag[i][j] = 1;
243 }
244 }
245
246 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
247 }
248
249 /* ----------------------------------------------------------------------
250 init specific to this pair style
251 ------------------------------------------------------------------------- */
252
init_style()253 void PairATM::init_style()
254 {
255 if (force->newton_pair == 0)
256 error->all(FLERR,"Pair style ATM requires newton pair on");
257
258 // need a full neighbor list
259
260 int irequest = neighbor->request(this,instance_me);
261 neighbor->requests[irequest]->half = 0;
262 neighbor->requests[irequest]->full = 1;
263 }
264
265 /* ----------------------------------------------------------------------
266 init for one i,j type pair and corresponding j,i
267 also for all k type permutations
268 ------------------------------------------------------------------------- */
269
init_one(int i,int j)270 double PairATM::init_one(int i, int j)
271 {
272 if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
273
274 // set all 6 symmetric permutations of I,J,K types to same nu value
275
276 int ntypes = atom->ntypes;
277 for (int k = j; k <= ntypes; k++)
278 nu[i][k][j] = nu[j][i][k] = nu[j][k][i] = nu[k][i][j] = nu[k][j][i] =
279 nu[i][j][k];
280
281 return cut_global;
282 }
283
284 /* ----------------------------------------------------------------------
285 proc 0 writes to restart file
286 ------------------------------------------------------------------------- */
287
write_restart(FILE * fp)288 void PairATM::write_restart(FILE *fp)
289 {
290 write_restart_settings(fp);
291
292 int i,j,k;
293 for (i = 1; i <= atom->ntypes; i++) {
294 for (j = i; j <= atom->ntypes; j++) {
295 fwrite(&setflag[i][j],sizeof(int),1,fp);
296 if (setflag[i][j])
297 for (k = j; k <= atom->ntypes; k++)
298 fwrite(&nu[i][j][k],sizeof(double),1,fp);
299 }
300 }
301 }
302
303 /* ----------------------------------------------------------------------
304 proc 0 reads from restart file, bcasts
305 ------------------------------------------------------------------------- */
306
read_restart(FILE * fp)307 void PairATM::read_restart(FILE *fp)
308 {
309 read_restart_settings(fp);
310 allocate();
311
312 int i,j,k;
313 int me = comm->me;
314 for (i = 1; i <= atom->ntypes; i++) {
315 for (j = i; j <= atom->ntypes; j++) {
316 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
317 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
318 if (setflag[i][j]) for (k = j; k <= atom->ntypes; k++) {
319 if (me == 0) utils::sfread(FLERR,&nu[i][j][k],sizeof(double),1,fp,nullptr,error);
320 MPI_Bcast(&nu[i][j][k],1,MPI_DOUBLE,0,world);
321 }
322 }
323 }
324 }
325
326 /* ----------------------------------------------------------------------
327 proc 0 writes to restart file
328 ------------------------------------------------------------------------- */
329
write_restart_settings(FILE * fp)330 void PairATM::write_restart_settings(FILE *fp)
331 {
332 fwrite(&cut_global,sizeof(double),1,fp);
333 fwrite(&cut_triple,sizeof(double),1,fp);
334 }
335
336 /* ----------------------------------------------------------------------
337 proc 0 reads from restart file, bcasts
338 ------------------------------------------------------------------------- */
339
read_restart_settings(FILE * fp)340 void PairATM::read_restart_settings(FILE *fp)
341 {
342 int me = comm->me;
343 if (me == 0) {
344 utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
345 utils::sfread(FLERR,&cut_triple,sizeof(double),1,fp,nullptr,error);
346 }
347 MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
348 MPI_Bcast(&cut_triple,1,MPI_DOUBLE,0,world);
349 }
350
351 /* ----------------------------------------------------------------------
352 Axilrod-Teller-Muto (dipole-dipole-dipole) potential
353 ------------------------------------------------------------------------- */
354
interaction_ddd(double nu,double r6,double rij2,double rik2,double rjk2,double * rij,double * rik,double * rjk,double * fj,double * fk,int eflag,double & eng)355 void PairATM::interaction_ddd(double nu, double r6,
356 double rij2, double rik2, double rjk2,
357 double *rij, double *rik, double *rjk,
358 double *fj, double *fk, int eflag, double &eng)
359 {
360 double r5inv,rri,rrj,rrk,rrr;
361 r5inv = nu / (r6*r6*sqrt(r6));
362 rri = rik[0]*rij[0] + rik[1]*rij[1] + rik[2]*rij[2];
363 rrj = rij[0]*rjk[0] + rij[1]*rjk[1] + rij[2]*rjk[2];
364 rrk = rjk[0]*rik[0] + rjk[1]*rik[1] + rjk[2]*rik[2];
365 rrr = 5.0*rri*rrj*rrk;
366 for (int i = 0; i < 3; i++) {
367 fj[i] = rrj*(rrk - rri)*rik[i] -
368 (rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] +
369 (rrk*rri - rik2*rij2 + rrr/rjk2) * rjk[i];
370 fj[i] *= 3.0*r5inv;
371 fk[i] = rrk*(rri + rrj)*rij[i] +
372 (rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] +
373 (rri*rrj + rij2*rjk2 - rrr/rik2) * rik[i];
374 fk[i] *= 3.0*r5inv;
375 }
376 if (eflag) eng = (r6 - 0.6*rrr)*r5inv;
377 }
378