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