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: Chuanfu Luo (luochuanfu@gmail.com)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_lj96_cut.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "error.h"
24 #include "force.h"
25 #include "math_const.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "neigh_request.h"
29 #include "neighbor.h"
30 #include "respa.h"
31 #include "update.h"
32
33 #include <cmath>
34
35 using namespace LAMMPS_NS;
36 using namespace MathConst;
37
38 /* ---------------------------------------------------------------------- */
39
PairLJ96Cut(LAMMPS * lmp)40 PairLJ96Cut::PairLJ96Cut(LAMMPS *lmp) : Pair(lmp)
41 {
42 respa_enable = 1;
43 writedata = 1;
44 cut_respa = nullptr;
45 }
46
47 /* ---------------------------------------------------------------------- */
48
~PairLJ96Cut()49 PairLJ96Cut::~PairLJ96Cut()
50 {
51 if (allocated) {
52 memory->destroy(setflag);
53 memory->destroy(cutsq);
54
55 memory->destroy(cut);
56 memory->destroy(epsilon);
57 memory->destroy(sigma);
58 memory->destroy(lj1);
59 memory->destroy(lj2);
60 memory->destroy(lj3);
61 memory->destroy(lj4);
62 memory->destroy(offset);
63 }
64 }
65
66 /* ---------------------------------------------------------------------- */
67
compute(int eflag,int vflag)68 void PairLJ96Cut::compute(int eflag, int vflag)
69 {
70 int i,j,ii,jj,inum,jnum,itype,jtype;
71 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
72 double rsq,r2inv,r3inv,r6inv,forcelj,factor_lj;
73 int *ilist,*jlist,*numneigh,**firstneigh;
74
75 evdwl = 0.0;
76 ev_init(eflag,vflag);
77
78 double **x = atom->x;
79 double **f = atom->f;
80 int *type = atom->type;
81 int nlocal = atom->nlocal;
82 double *special_lj = force->special_lj;
83 int newton_pair = force->newton_pair;
84
85 inum = list->inum;
86 ilist = list->ilist;
87 numneigh = list->numneigh;
88 firstneigh = list->firstneigh;
89
90 // loop over neighbors of my atoms
91
92 for (ii = 0; ii < inum; ii++) {
93 i = ilist[ii];
94 xtmp = x[i][0];
95 ytmp = x[i][1];
96 ztmp = x[i][2];
97 itype = type[i];
98 jlist = firstneigh[i];
99 jnum = numneigh[i];
100
101 for (jj = 0; jj < jnum; jj++) {
102 j = jlist[jj];
103 factor_lj = special_lj[sbmask(j)];
104 j &= NEIGHMASK;
105
106 delx = xtmp - x[j][0];
107 dely = ytmp - x[j][1];
108 delz = ztmp - x[j][2];
109 rsq = delx*delx + dely*dely + delz*delz;
110 jtype = type[j];
111
112 if (rsq < cutsq[itype][jtype]) {
113 r2inv = 1.0/rsq;
114 r6inv = r2inv*r2inv*r2inv;
115 r3inv = sqrt(r6inv);
116 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
117 fpair = factor_lj*forcelj*r2inv;
118
119 f[i][0] += delx*fpair;
120 f[i][1] += dely*fpair;
121 f[i][2] += delz*fpair;
122 if (newton_pair || j < nlocal) {
123 f[j][0] -= delx*fpair;
124 f[j][1] -= dely*fpair;
125 f[j][2] -= delz*fpair;
126 }
127
128 if (eflag) {
129 evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
130 offset[itype][jtype];
131 evdwl *= factor_lj;
132 }
133
134 if (evflag) ev_tally(i,j,nlocal,newton_pair,
135 evdwl,0.0,fpair,delx,dely,delz);
136 }
137 }
138 }
139
140 if (vflag_fdotr) virial_fdotr_compute();
141 }
142
143 /* ---------------------------------------------------------------------- */
144
compute_inner()145 void PairLJ96Cut::compute_inner()
146 {
147 int i,j,ii,jj,inum,jnum,itype,jtype;
148 double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
149 double rsq,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw;
150 int *ilist,*jlist,*numneigh,**firstneigh;
151
152 double **x = atom->x;
153 double **f = atom->f;
154 int *type = atom->type;
155 int nlocal = atom->nlocal;
156 double *special_lj = force->special_lj;
157 int newton_pair = force->newton_pair;
158
159 inum = list->inum_inner;
160 ilist = list->ilist_inner;
161 numneigh = list->numneigh_inner;
162 firstneigh = list->firstneigh_inner;
163
164 double cut_out_on = cut_respa[0];
165 double cut_out_off = cut_respa[1];
166
167 double cut_out_diff = cut_out_off - cut_out_on;
168 double cut_out_on_sq = cut_out_on*cut_out_on;
169 double cut_out_off_sq = cut_out_off*cut_out_off;
170
171 // loop over neighbors of my atoms
172
173 for (ii = 0; ii < inum; ii++) {
174 i = ilist[ii];
175 xtmp = x[i][0];
176 ytmp = x[i][1];
177 ztmp = x[i][2];
178 itype = type[i];
179 jlist = firstneigh[i];
180 jnum = numneigh[i];
181
182 for (jj = 0; jj < jnum; jj++) {
183 j = jlist[jj];
184 factor_lj = special_lj[sbmask(j)];
185 j &= NEIGHMASK;
186
187 delx = xtmp - x[j][0];
188 dely = ytmp - x[j][1];
189 delz = ztmp - x[j][2];
190 rsq = delx*delx + dely*dely + delz*delz;
191
192 if (rsq < cut_out_off_sq) {
193 r2inv = 1.0/rsq;
194 r6inv = r2inv*r2inv*r2inv;
195 r3inv = sqrt(r6inv);
196 jtype = type[j];
197 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
198 fpair = factor_lj*forcelj*r2inv;
199 if (rsq > cut_out_on_sq) {
200 rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
201 fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw);
202 }
203
204 f[i][0] += delx*fpair;
205 f[i][1] += dely*fpair;
206 f[i][2] += delz*fpair;
207 if (newton_pair || j < nlocal) {
208 f[j][0] -= delx*fpair;
209 f[j][1] -= dely*fpair;
210 f[j][2] -= delz*fpair;
211 }
212 }
213 }
214 }
215 }
216
217 /* ---------------------------------------------------------------------- */
218
compute_middle()219 void PairLJ96Cut::compute_middle()
220 {
221 int i,j,ii,jj,inum,jnum,itype,jtype;
222 double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
223 double rsq,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw;
224 int *ilist,*jlist,*numneigh,**firstneigh;
225
226 double **x = atom->x;
227 double **f = atom->f;
228 int *type = atom->type;
229 int nlocal = atom->nlocal;
230 double *special_lj = force->special_lj;
231 int newton_pair = force->newton_pair;
232
233 inum = list->inum_middle;
234 ilist = list->ilist_middle;
235 numneigh = list->numneigh_middle;
236 firstneigh = list->firstneigh_middle;
237
238 double cut_in_off = cut_respa[0];
239 double cut_in_on = cut_respa[1];
240 double cut_out_on = cut_respa[2];
241 double cut_out_off = cut_respa[3];
242
243 double cut_in_diff = cut_in_on - cut_in_off;
244 double cut_out_diff = cut_out_off - cut_out_on;
245 double cut_in_off_sq = cut_in_off*cut_in_off;
246 double cut_in_on_sq = cut_in_on*cut_in_on;
247 double cut_out_on_sq = cut_out_on*cut_out_on;
248 double cut_out_off_sq = cut_out_off*cut_out_off;
249
250 // loop over neighbors of my atoms
251
252 for (ii = 0; ii < inum; ii++) {
253 i = ilist[ii];
254 xtmp = x[i][0];
255 ytmp = x[i][1];
256 ztmp = x[i][2];
257 itype = type[i];
258 jlist = firstneigh[i];
259 jnum = numneigh[i];
260
261 for (jj = 0; jj < jnum; jj++) {
262 j = jlist[jj];
263 factor_lj = special_lj[sbmask(j)];
264 j &= NEIGHMASK;
265
266 delx = xtmp - x[j][0];
267 dely = ytmp - x[j][1];
268 delz = ztmp - x[j][2];
269 rsq = delx*delx + dely*dely + delz*delz;
270
271 if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
272 r2inv = 1.0/rsq;
273 r6inv = r2inv*r2inv*r2inv;
274 r3inv = sqrt(r6inv);
275 jtype = type[j];
276 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
277 fpair = factor_lj*forcelj*r2inv;
278 if (rsq < cut_in_on_sq) {
279 rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
280 fpair *= rsw*rsw*(3.0 - 2.0*rsw);
281 }
282 if (rsq > cut_out_on_sq) {
283 rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
284 fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
285 }
286
287 f[i][0] += delx*fpair;
288 f[i][1] += dely*fpair;
289 f[i][2] += delz*fpair;
290 if (newton_pair || j < nlocal) {
291 f[j][0] -= delx*fpair;
292 f[j][1] -= dely*fpair;
293 f[j][2] -= delz*fpair;
294 }
295 }
296 }
297 }
298 }
299
300 /* ---------------------------------------------------------------------- */
301
compute_outer(int eflag,int vflag)302 void PairLJ96Cut::compute_outer(int eflag, int vflag)
303 {
304 int i,j,ii,jj,inum,jnum,itype,jtype;
305 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
306 double rsq,r2inv,r3inv,r6inv,forcelj,factor_lj,rsw;
307 int *ilist,*jlist,*numneigh,**firstneigh;
308
309 evdwl = 0.0;
310 ev_init(eflag,vflag);
311
312 double **x = atom->x;
313 double **f = atom->f;
314 int *type = atom->type;
315 int nlocal = atom->nlocal;
316 double *special_lj = force->special_lj;
317 int newton_pair = force->newton_pair;
318
319 inum = list->inum;
320 ilist = list->ilist;
321 numneigh = list->numneigh;
322 firstneigh = list->firstneigh;
323
324 double cut_in_off = cut_respa[2];
325 double cut_in_on = cut_respa[3];
326
327 double cut_in_diff = cut_in_on - cut_in_off;
328 double cut_in_off_sq = cut_in_off*cut_in_off;
329 double cut_in_on_sq = cut_in_on*cut_in_on;
330
331 // loop over neighbors of my atoms
332
333 for (ii = 0; ii < inum; ii++) {
334 i = ilist[ii];
335 xtmp = x[i][0];
336 ytmp = x[i][1];
337 ztmp = x[i][2];
338 itype = type[i];
339 jlist = firstneigh[i];
340 jnum = numneigh[i];
341
342 for (jj = 0; jj < jnum; jj++) {
343 j = jlist[jj];
344 factor_lj = special_lj[sbmask(j)];
345 j &= NEIGHMASK;
346
347 delx = xtmp - x[j][0];
348 dely = ytmp - x[j][1];
349 delz = ztmp - x[j][2];
350 rsq = delx*delx + dely*dely + delz*delz;
351 jtype = type[j];
352
353 if (rsq < cutsq[itype][jtype]) {
354 if (rsq > cut_in_off_sq) {
355 r2inv = 1.0/rsq;
356 r6inv = r2inv*r2inv*r2inv;
357 r3inv = sqrt(r6inv);
358 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
359 fpair = factor_lj*forcelj*r2inv;
360 if (rsq < cut_in_on_sq) {
361 rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
362 fpair *= rsw*rsw*(3.0 - 2.0*rsw);
363 }
364
365 f[i][0] += delx*fpair;
366 f[i][1] += dely*fpair;
367 f[i][2] += delz*fpair;
368 if (newton_pair || j < nlocal) {
369 f[j][0] -= delx*fpair;
370 f[j][1] -= dely*fpair;
371 f[j][2] -= delz*fpair;
372 }
373 }
374
375 if (eflag) {
376 r2inv = 1.0/rsq;
377 r6inv = r2inv*r2inv*r2inv;
378 r3inv = sqrt(r6inv);
379 evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) - offset[itype][jtype];
380 evdwl *= factor_lj;
381 }
382
383 if (vflag) {
384 if (rsq <= cut_in_off_sq) {
385 r2inv = 1.0/rsq;
386 r6inv = r2inv*r2inv*r2inv;
387 r3inv = sqrt(r6inv);
388 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
389 fpair = factor_lj*forcelj*r2inv;
390 } else if (rsq < cut_in_on_sq)
391 fpair = factor_lj*forcelj*r2inv;
392 }
393
394 if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
395 }
396 }
397 }
398 }
399
400 /* ----------------------------------------------------------------------
401 allocate all arrays
402 ------------------------------------------------------------------------- */
403
allocate()404 void PairLJ96Cut::allocate()
405 {
406 allocated = 1;
407 int n = atom->ntypes;
408
409 memory->create(setflag,n+1,n+1,"pair:setflag");
410 for (int i = 1; i <= n; i++)
411 for (int j = i; j <= n; j++)
412 setflag[i][j] = 0;
413
414 memory->create(cutsq,n+1,n+1,"pair:cutsq");
415
416 memory->create(cut,n+1,n+1,"pair:cut");
417 memory->create(epsilon,n+1,n+1,"pair:epsilon");
418 memory->create(sigma,n+1,n+1,"pair:sigma");
419 memory->create(lj1,n+1,n+1,"pair:lj1");
420 memory->create(lj2,n+1,n+1,"pair:lj2");
421 memory->create(lj3,n+1,n+1,"pair:lj3");
422 memory->create(lj4,n+1,n+1,"pair:lj4");
423 memory->create(offset,n+1,n+1,"pair:offset");
424 }
425
426 /* ----------------------------------------------------------------------
427 global settings
428 ------------------------------------------------------------------------- */
429
settings(int narg,char ** arg)430 void PairLJ96Cut::settings(int narg, char **arg)
431 {
432 if (narg != 1) error->all(FLERR,"Illegal pair_style command");
433
434 cut_global = utils::numeric(FLERR,arg[0],false,lmp);
435
436 // reset cutoffs that have been explicitly set
437
438 if (allocated) {
439 int i,j;
440 for (i = 1; i <= atom->ntypes; i++)
441 for (j = i; j <= atom->ntypes; j++)
442 if (setflag[i][j]) cut[i][j] = cut_global;
443 }
444 }
445
446 /* ----------------------------------------------------------------------
447 set coeffs for one or more type pairs
448 ------------------------------------------------------------------------- */
449
coeff(int narg,char ** arg)450 void PairLJ96Cut::coeff(int narg, char **arg)
451 {
452 if (narg < 4 || narg > 5)
453 error->all(FLERR,"Incorrect args for pair coefficients");
454 if (!allocated) allocate();
455
456 int ilo,ihi,jlo,jhi;
457 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
458 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
459
460 double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
461 double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
462
463 double cut_one = cut_global;
464 if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp);
465
466 int count = 0;
467 for (int i = ilo; i <= ihi; i++) {
468 for (int j = MAX(jlo,i); j <= jhi; j++) {
469 epsilon[i][j] = epsilon_one;
470 sigma[i][j] = sigma_one;
471 cut[i][j] = cut_one;
472 setflag[i][j] = 1;
473 count++;
474 }
475 }
476
477 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
478 }
479
480 /* ----------------------------------------------------------------------
481 init specific to this pair style
482 ------------------------------------------------------------------------- */
483
init_style()484 void PairLJ96Cut::init_style()
485 {
486 // request regular or rRESPA neighbor list
487
488 int irequest;
489 int respa = 0;
490
491 if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
492 if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
493 if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
494 }
495
496 irequest = neighbor->request(this,instance_me);
497
498 if (respa >= 1) {
499 neighbor->requests[irequest]->respaouter = 1;
500 neighbor->requests[irequest]->respainner = 1;
501 }
502 if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
503
504 // set rRESPA cutoffs
505
506 if (utils::strmatch(update->integrate_style,"^respa") &&
507 ((Respa *) update->integrate)->level_inner >= 0)
508 cut_respa = ((Respa *) update->integrate)->cutoff;
509 else cut_respa = nullptr;
510 }
511
512 /* ----------------------------------------------------------------------
513 init for one type pair i,j and corresponding j,i
514 ------------------------------------------------------------------------- */
515
init_one(int i,int j)516 double PairLJ96Cut::init_one(int i, int j)
517 {
518 if (setflag[i][j] == 0) {
519 epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
520 sigma[i][i],sigma[j][j]);
521 sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
522 cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
523 }
524
525 lj1[i][j] = 36.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
526 lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
527 lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],9.0);
528 lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
529
530 if (offset_flag && (cut[i][j] > 0.0)) {
531 double ratio = sigma[i][j] / cut[i][j];
532 offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,9.0) - pow(ratio,6.0));
533 } else offset[i][j] = 0.0;
534
535 lj1[j][i] = lj1[i][j];
536 lj2[j][i] = lj2[i][j];
537 lj3[j][i] = lj3[i][j];
538 lj4[j][i] = lj4[i][j];
539 offset[j][i] = offset[i][j];
540
541 // check interior rRESPA cutoff
542
543 if (cut_respa && cut[i][j] < cut_respa[3])
544 error->all(FLERR,"Pair cutoff < Respa interior cutoff");
545
546 // compute I,J contribution to long-range tail correction
547 // count total # of atoms of type I and J via Allreduce
548
549 if (tail_flag) {
550 int *type = atom->type;
551 int nlocal = atom->nlocal;
552
553 double count[2],all[2];
554 count[0] = count[1] = 0.0;
555 for (int k = 0; k < nlocal; k++) {
556 if (type[k] == i) count[0] += 1.0;
557 if (type[k] == j) count[1] += 1.0;
558 }
559 MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
560
561 double sig3 = sigma[i][j]*sigma[i][j]*sigma[i][j];
562 double sig6 = sig3*sig3;
563 double rc3 = cut[i][j]*cut[i][j]*cut[i][j];
564 double rc6 = rc3*rc3;
565
566 etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
567 sig6 * (sig3 - 2.0*rc3) / (6.0*rc6);
568 ptail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
569 sig6 * (3.0*sig3 - 4.0*rc3) / (6.0*rc6);
570 }
571
572 return cut[i][j];
573 }
574
575 /* ----------------------------------------------------------------------
576 proc 0 writes to restart file
577 ------------------------------------------------------------------------- */
578
write_restart(FILE * fp)579 void PairLJ96Cut::write_restart(FILE *fp)
580 {
581 write_restart_settings(fp);
582
583 int i,j;
584 for (i = 1; i <= atom->ntypes; i++)
585 for (j = i; j <= atom->ntypes; j++) {
586 fwrite(&setflag[i][j],sizeof(int),1,fp);
587 if (setflag[i][j]) {
588 fwrite(&epsilon[i][j],sizeof(double),1,fp);
589 fwrite(&sigma[i][j],sizeof(double),1,fp);
590 fwrite(&cut[i][j],sizeof(double),1,fp);
591 }
592 }
593 }
594
595 /* ----------------------------------------------------------------------
596 proc 0 reads from restart file, bcasts
597 ------------------------------------------------------------------------- */
598
read_restart(FILE * fp)599 void PairLJ96Cut::read_restart(FILE *fp)
600 {
601 read_restart_settings(fp);
602 allocate();
603
604 int i,j;
605 int me = comm->me;
606 for (i = 1; i <= atom->ntypes; i++)
607 for (j = i; j <= atom->ntypes; j++) {
608 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
609 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
610 if (setflag[i][j]) {
611 if (me == 0) {
612 utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error);
613 utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error);
614 utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
615 }
616 MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
617 MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
618 MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
619 }
620 }
621 }
622
623 /* ----------------------------------------------------------------------
624 proc 0 writes to restart file
625 ------------------------------------------------------------------------- */
626
write_restart_settings(FILE * fp)627 void PairLJ96Cut::write_restart_settings(FILE *fp)
628 {
629 fwrite(&cut_global,sizeof(double),1,fp);
630 fwrite(&offset_flag,sizeof(int),1,fp);
631 fwrite(&mix_flag,sizeof(int),1,fp);
632 fwrite(&tail_flag,sizeof(int),1,fp);
633 }
634
635 /* ----------------------------------------------------------------------
636 proc 0 reads from restart file, bcasts
637 ------------------------------------------------------------------------- */
638
read_restart_settings(FILE * fp)639 void PairLJ96Cut::read_restart_settings(FILE *fp)
640 {
641 int me = comm->me;
642 if (me == 0) {
643 utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
644 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
645 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
646 utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error);
647 }
648 MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
649 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
650 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
651 MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
652 }
653
654 /* ----------------------------------------------------------------------
655 proc 0 writes to data file
656 ------------------------------------------------------------------------- */
657
write_data(FILE * fp)658 void PairLJ96Cut::write_data(FILE *fp)
659 {
660 for (int i = 1; i <= atom->ntypes; i++)
661 fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]);
662 }
663
664 /* ----------------------------------------------------------------------
665 proc 0 writes all pairs to data file
666 ------------------------------------------------------------------------- */
667
write_data_all(FILE * fp)668 void PairLJ96Cut::write_data_all(FILE *fp)
669 {
670 for (int i = 1; i <= atom->ntypes; i++)
671 for (int j = i; j <= atom->ntypes; j++)
672 fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]);
673 }
674
675 /* ---------------------------------------------------------------------- */
676
single(int,int,int itype,int jtype,double rsq,double,double factor_lj,double & fforce)677 double PairLJ96Cut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
678 double /*factor_coul*/, double factor_lj,
679 double &fforce)
680 {
681 double r2inv,r3inv,r6inv,forcelj,philj;
682
683 r2inv = 1.0/rsq;
684 r6inv = r2inv*r2inv*r2inv;
685 r3inv = sqrt(r6inv);
686 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
687 fforce = factor_lj*forcelj*r2inv;
688
689 philj = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
690 offset[itype][jtype];
691 return factor_lj*philj;
692 }
693