pair_atm.cpp 11.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author: Sergey Lishchuk
------------------------------------------------------------------------- */

Axel Kohlmeyer's avatar
Axel Kohlmeyer committed
18
19
#include <cmath>
#include "pair_atm.h"
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"

using namespace LAMMPS_NS;

static const char cite_atm_package[] =
  "ATM package:\n\n"
  "@Article{Lishchuk:2012:164501,\n"
  " author = {S. V. Lishchuk},\n"
  " title = {Role of three-body interactions in formation of bulk viscosity in liquid argon},\n"
  " journal = {J.~Chem.~Phys.},\n"
  " year =    2012,\n"
  " volume =  136,\n"
  " pages =   {164501}\n"
  "}\n\n";

/* ---------------------------------------------------------------------- */

PairATM::PairATM(LAMMPS *lmp) : Pair(lmp)
{
  if (lmp->citeme) lmp->citeme->add(cite_atm_package);

  single_enable = 0;
50
51
  restartinfo = 1;
  one_coeff = 0;
52
53
54
55
56
57
58
59
60
61
  manybody_flag = 1;
}

/* ----------------------------------------------------------------------
   check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */

PairATM::~PairATM()
{
  if (copymode) return;
62

63
64
65
  if (allocated) {
    memory->destroy(setflag);
    memory->destroy(cutsq);
66
67

    memory->destroy(nu);
68
69
70
71
72
73
74
75
76
77
78
  }
}

/* ----------------------------------------------------------------------
   workhorse routine that computes pairwise interactions
------------------------------------------------------------------------- */

void PairATM::compute(int eflag, int vflag)
{
  int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
  double xi,yi,zi,evdwl;
79
  double rij2,rik2,rjk2;
80
  double rij[3],rik[3],rjk[3],fj[3],fk[3];
81
  double nu_local;
82
83
84
85
86
87
88
89
90
91
  int *ilist,*jlist,*numneigh,**firstneigh;

  evdwl = 0.0;
  if (eflag || vflag) ev_setup(eflag,vflag);
  else evflag = vflag_fdotr = 0;

  double **x = atom->x;
  double **f = atom->f;
  int *type = atom->type;

92
  double cutoff_squared = cut_global*cut_global;
93
94
95
  double triple = cut_triple*cut_triple*cut_triple;
  double cutoff_triple_sixth = triple*triple;

96
97
98
99
100
  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;

101
102
103
104
  int count1 = 0;
  int count2 = 0;
  int count3 = 0;

105
106
107
108
109
110
111
112
113
114
  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    xi = x[i][0];
    yi = x[i][1];
    zi = x[i][2];

    jlist = firstneigh[i];
    jnum = numneigh[i];
    jnumm1 = jnum - 1;

115
116
117
118
    //    for (jj = 0; jj < jnumm1; jj++) {
    // replace with this line:
    for (jj = 0; jj < jnum; jj++) {

119
120
121
      j = jlist[jj];
      j &= NEIGHMASK;
      rij[0] = x[j][0] - xi;
122
      if (rij[0] < 0.0) continue;
123
      rij[1] = x[j][1] - yi;
124
      if (rij[0] == 0.0 and rij[1] < 0.0) continue;
125
      rij[2] = x[j][2] - zi;
126
      if (rij[0] == 0.0 and rij[1] == 0.0 and rij[2] < 0.0) continue;
127
      rij2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
128
      count1++;
129
      if (rij2 > cutoff_squared) continue;
130
131
132
133
134
135
      count2++;

      //for (kk = jj+1; kk < jnum; kk++) {
      // replace with these two lines:
      for (kk = 0; kk < jnum; kk++) {
        if (kk == jj) continue;
136
137
138
139
140

        k = jlist[kk];
        k &= NEIGHMASK;

        rjk[0] = x[k][0] - x[j][0];
141
        if (rjk[0] < 0.0) continue;
142
        rjk[1] = x[k][1] - x[j][1];
143
        if (rjk[0] == 0.0 and rjk[1] < 0.0) continue;
144
        rjk[2] = x[k][2] - x[j][2];
145
        if (rjk[0] == 0.0 and rjk[1] == 0.0 and rjk[2] < 0.0) continue;
146
        rjk2 = rjk[0]*rjk[0] + rjk[1]*rjk[1] + rjk[2]*rjk[2];
147
        if (rjk2 > cutoff_squared) continue;
148

149
150
151
152
        rik[0] = x[k][0] - xi;
        rik[1] = x[k][1] - yi;
        rik[2] = x[k][2] - zi;
        rik2 = rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2];
153
        if (rik2 > cutoff_squared) continue;
154

155
156
        double r6 = rij2*rjk2*rik2;
        if (r6 > cutoff_triple_sixth) continue;
157

158
159
        nu_local = nu[type[i]][type[j]][type[k]];
        if (nu_local == 0.0) continue;
160

161
        count3++;
162
        interaction_ddd(nu_local,
163
                        r6,rij2,rik2,rjk2,rij,rik,rjk,fj,fk,eflag,evdwl);
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

        f[i][0] -= fj[0] + fk[0];
        f[i][1] -= fj[1] + fk[1];
        f[i][2] -= fj[2] + fk[2];
        f[j][0] += fj[0];
        f[j][1] += fj[1];
        f[j][2] += fj[2];
        f[k][0] += fk[0];
        f[k][1] += fk[1];
        f[k][2] += fk[2];

        if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,rij,rik);
      }
    }
  }

180
181
182
183
184
185
186
187
188
  int count = count1;
  MPI_Allreduce(&count,&count1,1,MPI_INT,MPI_SUM,world);
  count = count2;
  MPI_Allreduce(&count,&count2,1,MPI_INT,MPI_SUM,world);
  count = count3;
  MPI_Allreduce(&count,&count3,1,MPI_INT,MPI_SUM,world);
  printf("FORCE %g %d %d %d\n",cutoff_squared,count1,count2,count3);


189
190
191
  if (vflag_fdotr) virial_fdotr_compute();
}

192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
/* ---------------------------------------------------------------------- */

void PairATM::allocate()
{
  allocated = 1;
  int n = atom->ntypes;

  memory->create(setflag,n+1,n+1,"pair:setflag");
  for (int i = 1; i <= n; i++)
    for (int j = i; j <= n; j++)
      setflag[i][j] = 0;

  memory->create(cutsq,n+1,n+1,"pair:cutsq");

  memory->create(nu,n+1,n+1,n+1,"pair:nu");

  // initialize all nu values to 0.0

  for (int i = 1; i <= n; i++)
    for (int j = 1; j <= n; j++)
      for (int k = 1; k <= n; k++)
        nu[i][j][k] = 0.0;
}

216
/* ----------------------------------------------------------------------
217
   global settings
218
219
220
221
------------------------------------------------------------------------- */

void PairATM::settings(int narg, char **arg)
{
222
  if (narg != 2) error->all(FLERR,"Illegal pair_style command");
223

224
  cut_global = force->numeric(FLERR,arg[0]);
225
  cut_triple = force->numeric(FLERR,arg[1]);
226
227
228
}

/* ----------------------------------------------------------------------
229
   set coefficients for one I,J,K type triplet
230
231
232
233
------------------------------------------------------------------------- */

void PairATM::coeff(int narg, char **arg)
{
234
  if (narg != 4) error->all(FLERR,"Incorrect args for pair coefficients");
235
236
237
238
239
240
241
242
243
244
245
246
247
  if (!allocated) allocate();

  int ilo,ihi,jlo,jhi,klo,khi;
  force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
  force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
  force->bounds(FLERR,arg[2],atom->ntypes,klo,khi);

  double nu_one = force->numeric(FLERR,arg[3]);

  int count = 0;
  for (int i = ilo; i <= ihi; i++) {
    for (int j = MAX(jlo,i); j<=jhi; j++) {
      for (int k = MAX(klo,j); k<=khi; k++) {
248
        nu[i][j][k] = nu_one;
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
        count++;
      }
      setflag[i][j] = 1;
    }
  }

  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}

/* ----------------------------------------------------------------------
   init specific to this pair style
------------------------------------------------------------------------- */

void PairATM::init_style()
{
264
265
266
  if (force->newton_pair == 0)
    error->all(FLERR,"Pair style ATM requires newton pair on");

267
  // need a full neighbor list
268

269
270
271
272
273
274
  int irequest = neighbor->request(this,instance_me);
  neighbor->requests[irequest]->half = 0;
  neighbor->requests[irequest]->full = 1;
}

/* ----------------------------------------------------------------------
275
276
   init for one i,j type pair and corresponding j,i
   also for all k type permutations
277
278
279
280
281
------------------------------------------------------------------------- */

double PairATM::init_one(int i, int j)
{
  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
282
283
284
285
286
287
288
289

  // set all 6 symmetric permutations of I,J,K types to same nu value

  int ntypes = atom->ntypes;
  for (int k = j; k <= ntypes; k++)
    nu[i][k][j] = nu[j][i][k] = nu[j][k][i] = nu[k][i][j] = nu[k][j][i] = 
      nu[i][j][k];

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
  return cut_global;
}

/* ----------------------------------------------------------------------
   proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairATM::write_restart(FILE *fp)
{
  write_restart_settings(fp);

  int i,j,k;
  for (i = 1; i <= atom->ntypes; i++) {
    for (j = i; j <= atom->ntypes; j++) {
      fwrite(&setflag[i][j],sizeof(int),1,fp);
305
      if (setflag[i][j]) 
306
        for (k = j; k <= atom->ntypes; k++) 
307
          fwrite(&nu[i][j][k],sizeof(double),1,fp);
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
    }
  }
}

/* ----------------------------------------------------------------------
   proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairATM::read_restart(FILE *fp)
{
  read_restart_settings(fp);
  allocate();

  int i,j,k;
  int me = comm->me;
  for (i = 1; i <= atom->ntypes; i++) {
    for (j = i; j <= atom->ntypes; j++) {
      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
327
      if (setflag[i][j]) for (k = j; k <= atom->ntypes; k++) {
328
329
330
331
332
333
334
335
336
337
338
339
340
341
        if (me == 0) fread(&nu[i][j][k],sizeof(double),1,fp);
        MPI_Bcast(&nu[i][j][k],1,MPI_DOUBLE,0,world);
      }
    }
  }
}

/* ----------------------------------------------------------------------
   proc 0 writes to restart file
------------------------------------------------------------------------- */

void PairATM::write_restart_settings(FILE *fp)
{
  fwrite(&cut_global,sizeof(double),1,fp);
342
  fwrite(&cut_triple,sizeof(double),1,fp);
343
344
345
346
347
348
349
350
351
}

/* ----------------------------------------------------------------------
   proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */

void PairATM::read_restart_settings(FILE *fp)
{
  int me = comm->me;
352
353
354
355
  if (me == 0) {
    fread(&cut_global,sizeof(double),1,fp);
    fread(&cut_triple,sizeof(double),1,fp);
  }
356
  MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
357
  MPI_Bcast(&cut_triple,1,MPI_DOUBLE,0,world);
358
359
360
361
362
363
}

/* ----------------------------------------------------------------------
   Axilrod-Teller-Muto (dipole-dipole-dipole) potential
------------------------------------------------------------------------- */

364
void PairATM::interaction_ddd(double nu, double r6,
365
                              double rij2, double rik2, double rjk2,
366
367
                              double *rij, double *rik, double *rjk,
                              double *fj, double *fk, int eflag, double &eng)
368
{
369
  double r5inv,rri,rrj,rrk,rrr;
370
371
372
373
374
  r5inv = nu / (r6*r6*sqrt(r6));
  rri = rik[0]*rij[0] + rik[1]*rij[1] + rik[2]*rij[2];
  rrj = rij[0]*rjk[0] + rij[1]*rjk[1] + rij[2]*rjk[2];
  rrk = rjk[0]*rik[0] + rjk[1]*rik[1] + rjk[2]*rik[2];
  rrr = 5.0*rri*rrj*rrk;
375
376
377
378
  for (int i = 0; i < 3; i++) {
    fj[i] = rrj*(rrk - rri)*rik[i] - 
      (rrk*rri - rjk2*rik2 + rrr/rij2) * rij[i] + 
      (rrk*rri - rik2*rij2 + rrr/rjk2) * rjk[i];
379
    fj[i] *= 3.0*r5inv;
380
381
382
    fk[i] = rrk*(rri + rrj)*rij[i] + 
      (rri*rrj + rik2*rij2 - rrr/rjk2) * rjk[i] + 
      (rri*rrj + rij2*rjk2 - rrr/rik2) * rik[i];
383
384
385
386
    fk[i] *= 3.0*r5inv;
  }
  if (eflag) eng = (r6 - 0.6*rrr)*r5inv;
}