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_yukawa.h"
16 
17 #include <cmath>
18 #include "atom.h"
19 #include "force.h"
20 #include "comm.h"
21 #include "neigh_list.h"
22 #include "memory.h"
23 #include "error.h"
24 
25 
26 using namespace LAMMPS_NS;
27 
28 /* ---------------------------------------------------------------------- */
29 
PairYukawa(LAMMPS * lmp)30 PairYukawa::PairYukawa(LAMMPS *lmp) : Pair(lmp)
31 {
32   writedata = 1;
33 }
34 
35 /* ---------------------------------------------------------------------- */
36 
~PairYukawa()37 PairYukawa::~PairYukawa()
38 {
39   if (copymode) return;
40 
41   if (allocated) {
42     memory->destroy(setflag);
43     memory->destroy(cutsq);
44 
45     memory->destroy(rad);
46     memory->destroy(cut);
47     memory->destroy(a);
48     memory->destroy(offset);
49   }
50 }
51 
52 /* ---------------------------------------------------------------------- */
53 
compute(int eflag,int vflag)54 void PairYukawa::compute(int eflag, int vflag)
55 {
56   int i,j,ii,jj,inum,jnum,itype,jtype;
57   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
58   double rsq,r2inv,r,rinv,screening,forceyukawa,factor;
59   int *ilist,*jlist,*numneigh,**firstneigh;
60 
61   evdwl = 0.0;
62   ev_init(eflag,vflag);
63 
64   double **x = atom->x;
65   double **f = atom->f;
66   int *type = atom->type;
67   int nlocal = atom->nlocal;
68   double *special_lj = force->special_lj;
69   int newton_pair = force->newton_pair;
70 
71   inum = list->inum;
72   ilist = list->ilist;
73   numneigh = list->numneigh;
74   firstneigh = list->firstneigh;
75 
76   // loop over neighbors of my atoms
77 
78   for (ii = 0; ii < inum; ii++) {
79     i = ilist[ii];
80     xtmp = x[i][0];
81     ytmp = x[i][1];
82     ztmp = x[i][2];
83     itype = type[i];
84     jlist = firstneigh[i];
85     jnum = numneigh[i];
86 
87     for (jj = 0; jj < jnum; jj++) {
88       j = jlist[jj];
89       factor = special_lj[sbmask(j)];
90       j &= NEIGHMASK;
91 
92       delx = xtmp - x[j][0];
93       dely = ytmp - x[j][1];
94       delz = ztmp - x[j][2];
95       rsq = delx*delx + dely*dely + delz*delz;
96       jtype = type[j];
97 
98       if (rsq < cutsq[itype][jtype]) {
99         r2inv = 1.0/rsq;
100         r = sqrt(rsq);
101         rinv = 1.0/r;
102         screening = exp(-kappa*r);
103         forceyukawa = a[itype][jtype] * screening * (kappa + rinv);
104 
105         fpair = factor*forceyukawa * r2inv;
106 
107         f[i][0] += delx*fpair;
108         f[i][1] += dely*fpair;
109         f[i][2] += delz*fpair;
110         if (newton_pair || j < nlocal) {
111           f[j][0] -= delx*fpair;
112           f[j][1] -= dely*fpair;
113           f[j][2] -= delz*fpair;
114         }
115 
116         if (eflag) {
117           evdwl = a[itype][jtype] * screening * rinv - offset[itype][jtype];
118           evdwl *= factor;
119         }
120 
121         if (evflag) ev_tally(i,j,nlocal,newton_pair,
122                              evdwl,0.0,fpair,delx,dely,delz);
123       }
124     }
125   }
126 
127   if (vflag_fdotr) virial_fdotr_compute();
128 }
129 
130 /* ----------------------------------------------------------------------
131    allocate all arrays
132 ------------------------------------------------------------------------- */
133 
allocate()134 void PairYukawa::allocate()
135 {
136   allocated = 1;
137   int n = atom->ntypes;
138 
139   memory->create(setflag,n+1,n+1,"pair:setflag");
140   for (int i = 1; i <= n; i++)
141     for (int j = i; j <= n; j++)
142       setflag[i][j] = 0;
143 
144   memory->create(cutsq,n+1,n+1,"pair:cutsq");
145   memory->create(rad,n+1,"pair:rad");
146   memory->create(cut,n+1,n+1,"pair:cut");
147   memory->create(a,n+1,n+1,"pair:a");
148   memory->create(offset,n+1,n+1,"pair:offset");
149 }
150 
151 /* ----------------------------------------------------------------------
152    global settings
153 ------------------------------------------------------------------------- */
154 
settings(int narg,char ** arg)155 void PairYukawa::settings(int narg, char **arg)
156 {
157   if (narg != 2) error->all(FLERR,"Illegal pair_style command");
158 
159   kappa = utils::numeric(FLERR,arg[0],false,lmp);
160   cut_global = utils::numeric(FLERR,arg[1],false,lmp);
161 
162   // reset cutoffs that have been explicitly set
163 
164   if (allocated) {
165     int i,j;
166     for (i = 1; i <= atom->ntypes; i++)
167       for (j = i; j <= atom->ntypes; j++)
168         if (setflag[i][j]) cut[i][j] = cut_global;
169   }
170 }
171 
172 /* ----------------------------------------------------------------------
173    set coeffs for one or more type pairs
174 ------------------------------------------------------------------------- */
175 
coeff(int narg,char ** arg)176 void PairYukawa::coeff(int narg, char **arg)
177 {
178   if (narg < 3 || narg > 4)
179     error->all(FLERR,"Incorrect args for pair coefficients");
180   if (!allocated) allocate();
181 
182   int ilo,ihi,jlo,jhi;
183   utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
184   utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
185 
186   double a_one = utils::numeric(FLERR,arg[2],false,lmp);
187 
188   double cut_one = cut_global;
189   if (narg == 4) cut_one = utils::numeric(FLERR,arg[3],false,lmp);
190 
191   int count = 0;
192   for (int i = ilo; i <= ihi; i++) {
193     for (int j = MAX(jlo,i); j <= jhi; j++) {
194       a[i][j] = a_one;
195       cut[i][j] = cut_one;
196       setflag[i][j] = 1;
197       count++;
198     }
199   }
200 
201   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
202 }
203 
204 /* ----------------------------------------------------------------------
205    init for one type pair i,j and corresponding j,i
206 ------------------------------------------------------------------------- */
207 
init_one(int i,int j)208 double PairYukawa::init_one(int i, int j)
209 {
210   if (setflag[i][j] == 0) {
211     a[i][j] = mix_energy(a[i][i],a[j][j],1.0,1.0);
212     cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
213   }
214 
215   if (offset_flag && (cut[i][j] > 0.0)) {
216     double screening = exp(-kappa * cut[i][j]);
217     offset[i][j] = a[i][j] * screening / cut[i][j];
218   } else offset[i][j] = 0.0;
219 
220   a[j][i] = a[i][j];
221   offset[j][i] = offset[i][j];
222 
223   return cut[i][j];
224 }
225 
226 /* ----------------------------------------------------------------------
227    proc 0 writes to restart file
228 ------------------------------------------------------------------------- */
229 
write_restart(FILE * fp)230 void PairYukawa::write_restart(FILE *fp)
231 {
232   write_restart_settings(fp);
233 
234   int i,j;
235   for (i = 1; i <= atom->ntypes; i++)
236     for (j = i; j <= atom->ntypes; j++) {
237       fwrite(&setflag[i][j],sizeof(int),1,fp);
238       if (setflag[i][j]) {
239         fwrite(&a[i][j],sizeof(double),1,fp);
240         fwrite(&cut[i][j],sizeof(double),1,fp);
241       }
242     }
243 }
244 
245 /* ----------------------------------------------------------------------
246    proc 0 reads from restart file, bcasts
247 ------------------------------------------------------------------------- */
248 
read_restart(FILE * fp)249 void PairYukawa::read_restart(FILE *fp)
250 {
251   read_restart_settings(fp);
252 
253   allocate();
254 
255   int i,j;
256   int me = comm->me;
257   for (i = 1; i <= atom->ntypes; i++)
258     for (j = i; j <= atom->ntypes; j++) {
259       if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
260       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
261       if (setflag[i][j]) {
262         if (me == 0) {
263           utils::sfread(FLERR,&a[i][j],sizeof(double),1,fp,nullptr,error);
264           utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
265         }
266         MPI_Bcast(&a[i][j],1,MPI_DOUBLE,0,world);
267         MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
268       }
269     }
270 }
271 
272 /* ----------------------------------------------------------------------
273    proc 0 writes to restart file
274 ------------------------------------------------------------------------- */
275 
write_restart_settings(FILE * fp)276 void PairYukawa::write_restart_settings(FILE *fp)
277 {
278   fwrite(&kappa,sizeof(double),1,fp);
279   fwrite(&cut_global,sizeof(double),1,fp);
280   fwrite(&offset_flag,sizeof(int),1,fp);
281   fwrite(&mix_flag,sizeof(int),1,fp);
282 }
283 
284 /* ----------------------------------------------------------------------
285    proc 0 reads from restart file, bcasts
286 ------------------------------------------------------------------------- */
287 
read_restart_settings(FILE * fp)288 void PairYukawa::read_restart_settings(FILE *fp)
289 {
290   if (comm->me == 0) {
291     utils::sfread(FLERR,&kappa,sizeof(double),1,fp,nullptr,error);
292     utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
293     utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
294     utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
295   }
296   MPI_Bcast(&kappa,1,MPI_DOUBLE,0,world);
297   MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
298   MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
299   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
300 }
301 
302 /* ----------------------------------------------------------------------
303    proc 0 writes to data file
304 ------------------------------------------------------------------------- */
305 
write_data(FILE * fp)306 void PairYukawa::write_data(FILE *fp)
307 {
308   for (int i = 1; i <= atom->ntypes; i++)
309     fprintf(fp,"%d %g\n",i,a[i][i]);
310 }
311 
312 /* ----------------------------------------------------------------------
313    proc 0 writes all pairs to data file
314 ------------------------------------------------------------------------- */
315 
write_data_all(FILE * fp)316 void PairYukawa::write_data_all(FILE *fp)
317 {
318   for (int i = 1; i <= atom->ntypes; i++)
319     for (int j = i; j <= atom->ntypes; j++)
320       fprintf(fp,"%d %d %g %g\n",i,j,a[i][j],cut[i][j]);
321 }
322 
323 /* ---------------------------------------------------------------------- */
324 
single(int,int,int itype,int jtype,double rsq,double,double factor_lj,double & fforce)325 double PairYukawa::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
326                           double /*factor_coul*/, double factor_lj,
327                           double &fforce)
328 {
329   double r2inv,r,rinv,screening,forceyukawa,phi;
330 
331   r2inv = 1.0/rsq;
332   r = sqrt(rsq);
333   rinv = 1.0/r;
334   screening = exp(-kappa*r);
335   forceyukawa = a[itype][jtype] * screening * (kappa + rinv);
336   fforce = factor_lj*forceyukawa * r2inv;
337 
338   phi = a[itype][jtype] * screening * rinv - offset[itype][jtype];
339   return factor_lj*phi;
340 }
341