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