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