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 authors: Amit Kumar and Michael Bybee (UIUC)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_brownian.h"
20
21 #include <cmath>
22 #include <cstring>
23 #include "atom.h"
24 #include "comm.h"
25 #include "force.h"
26 #include "neighbor.h"
27 #include "neigh_list.h"
28 #include "domain.h"
29 #include "update.h"
30 #include "modify.h"
31 #include "fix.h"
32 #include "fix_wall.h"
33 #include "input.h"
34 #include "variable.h"
35 #include "random_mars.h"
36 #include "math_const.h"
37 #include "math_special.h"
38 #include "memory.h"
39 #include "error.h"
40
41
42 using namespace LAMMPS_NS;
43 using namespace MathConst;
44 using namespace MathSpecial;
45
46 // same as fix_wall.cpp
47
48 enum{EDGE,CONSTANT,VARIABLE};
49
50 /* ---------------------------------------------------------------------- */
51
PairBrownian(LAMMPS * lmp)52 PairBrownian::PairBrownian(LAMMPS *lmp) : Pair(lmp)
53 {
54 single_enable = 0;
55 random = nullptr;
56 }
57
58 /* ---------------------------------------------------------------------- */
59
~PairBrownian()60 PairBrownian::~PairBrownian()
61 {
62 if (allocated) {
63 memory->destroy(setflag);
64 memory->destroy(cutsq);
65
66 memory->destroy(cut);
67 memory->destroy(cut_inner);
68 }
69 delete random;
70 }
71
72 /* ---------------------------------------------------------------------- */
73
compute(int eflag,int vflag)74 void PairBrownian::compute(int eflag, int vflag)
75 {
76 int i,j,ii,jj,inum,jnum,itype,jtype;
77 double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,tx,ty,tz;
78 double rsq,r,h_sep,radi;
79 int *ilist,*jlist,*numneigh,**firstneigh;
80
81 ev_init(eflag,vflag);
82
83 double **x = atom->x;
84 double **f = atom->f;
85 double **torque = atom->torque;
86 double *radius = atom->radius;
87 int *type = atom->type;
88 int nlocal = atom->nlocal;
89 int newton_pair = force->newton_pair;
90
91 double vxmu2f = force->vxmu2f;
92 double randr;
93 double prethermostat;
94 double xl[3],a_sq,a_sh,a_pu,Fbmag;
95 double p1[3],p2[3],p3[3];
96 int overlaps = 0;
97
98 // This section of code adjusts R0/RT0/RS0 if necessary due to changes
99 // in the volume fraction as a result of fix deform or moving walls
100
101 double dims[3], wallcoord;
102 if (flagVF) // Flag for volume fraction corrections
103 if (flagdeform || flagwall == 2) { // Possible changes in volume fraction
104 if (flagdeform && !flagwall)
105 for (j = 0; j < 3; j++)
106 dims[j] = domain->prd[j];
107 else if (flagwall == 2 || (flagdeform && flagwall == 1)) {
108 double wallhi[3], walllo[3];
109 for (int j = 0; j < 3; j++) {
110 wallhi[j] = domain->prd[j];
111 walllo[j] = 0;
112 }
113 for (int m = 0; m < wallfix->nwall; m++) {
114 int dim = wallfix->wallwhich[m] / 2;
115 int side = wallfix->wallwhich[m] % 2;
116 if (wallfix->xstyle[m] == VARIABLE) {
117 wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
118 }
119 else wallcoord = wallfix->coord0[m];
120 if (side == 0) walllo[dim] = wallcoord;
121 else wallhi[dim] = wallcoord;
122 }
123 for (int j = 0; j < 3; j++)
124 dims[j] = wallhi[j] - walllo[j];
125 }
126 double vol_T = dims[0]*dims[1]*dims[2];
127 double vol_f = vol_P/vol_T;
128 if (flaglog == 0) {
129 R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
130 RT0 = 8*MY_PI*mu*cube(rad);
131 //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.33*vol_f + 2.80*vol_f*vol_f);
132 } else {
133 R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
134 RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
135 //RS0 = 20.0/3.0*MY_PI*mu*pow(rad,3)*(1.0 + 3.64*vol_f - 6.95*vol_f*vol_f);
136 }
137 }
138
139 // scale factor for Brownian moments
140
141 prethermostat = sqrt(24.0*force->boltz*t_target/update->dt);
142 prethermostat *= sqrt(force->vxmu2f/force->ftm2v/force->mvv2e);
143
144 inum = list->inum;
145 ilist = list->ilist;
146 numneigh = list->numneigh;
147 firstneigh = list->firstneigh;
148
149 for (ii = 0; ii < inum; ii++) {
150 i = ilist[ii];
151 xtmp = x[i][0];
152 ytmp = x[i][1];
153 ztmp = x[i][2];
154 itype = type[i];
155 radi = radius[i];
156 jlist = firstneigh[i];
157 jnum = numneigh[i];
158
159 // FLD contribution to force and torque due to isotropic terms
160
161 if (flagfld) {
162 f[i][0] += prethermostat*sqrt(R0)*(random->uniform()-0.5);
163 f[i][1] += prethermostat*sqrt(R0)*(random->uniform()-0.5);
164 f[i][2] += prethermostat*sqrt(R0)*(random->uniform()-0.5);
165 if (flaglog) {
166 torque[i][0] += prethermostat*sqrt(RT0)*(random->uniform()-0.5);
167 torque[i][1] += prethermostat*sqrt(RT0)*(random->uniform()-0.5);
168 torque[i][2] += prethermostat*sqrt(RT0)*(random->uniform()-0.5);
169 }
170 }
171
172 if (!flagHI) continue;
173
174 for (jj = 0; jj < jnum; jj++) {
175 j = jlist[jj];
176 j &= NEIGHMASK;
177
178 delx = xtmp - x[j][0];
179 dely = ytmp - x[j][1];
180 delz = ztmp - x[j][2];
181 rsq = delx*delx + dely*dely + delz*delz;
182 jtype = type[j];
183
184 if (rsq < cutsq[itype][jtype]) {
185 r = sqrt(rsq);
186
187 // scalar resistances a_sq and a_sh
188
189 h_sep = r - 2.0*radi;
190
191 // check for overlaps
192
193 if (h_sep < 0.0) overlaps++;
194
195 // if less than minimum gap, use minimum gap instead
196
197 if (r < cut_inner[itype][jtype])
198 h_sep = cut_inner[itype][jtype] - 2.0*radi;
199
200 // scale h_sep by radi
201
202 h_sep = h_sep/radi;
203
204 // scalar resistances
205
206 if (flaglog) {
207 a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep + 9.0/40.0*log(1.0/h_sep));
208 a_sh = 6.0*MY_PI*mu*radi*(1.0/6.0*log(1.0/h_sep));
209 a_pu = 8.0*MY_PI*mu*cube(radi)*(3.0/160.0*log(1.0/h_sep));
210 } else
211 a_sq = 6.0*MY_PI*mu*radi*(1.0/4.0/h_sep);
212
213 // generate the Pairwise Brownian Force: a_sq
214
215 Fbmag = prethermostat*sqrt(a_sq);
216
217 // generate a random number
218
219 randr = random->uniform()-0.5;
220
221 // contribution due to Brownian motion
222
223 fx = Fbmag*randr*delx/r;
224 fy = Fbmag*randr*dely/r;
225 fz = Fbmag*randr*delz/r;
226
227 // add terms due to a_sh
228
229 if (flaglog) {
230
231 // generate two orthogonal vectors to the line of centers
232
233 p1[0] = delx/r; p1[1] = dely/r; p1[2] = delz/r;
234 set_3_orthogonal_vectors(p1,p2,p3);
235
236 // magnitude
237
238 Fbmag = prethermostat*sqrt(a_sh);
239
240 // force in each of the two directions
241
242 randr = random->uniform()-0.5;
243 fx += Fbmag*randr*p2[0];
244 fy += Fbmag*randr*p2[1];
245 fz += Fbmag*randr*p2[2];
246
247 randr = random->uniform()-0.5;
248 fx += Fbmag*randr*p3[0];
249 fy += Fbmag*randr*p3[1];
250 fz += Fbmag*randr*p3[2];
251 }
252
253 // scale forces to appropriate units
254
255 fx = vxmu2f*fx;
256 fy = vxmu2f*fy;
257 fz = vxmu2f*fz;
258
259 // sum to total force
260
261 f[i][0] -= fx;
262 f[i][1] -= fy;
263 f[i][2] -= fz;
264
265 if (newton_pair || j < nlocal) {
266 //randr = random->uniform()-0.5;
267 //fx = Fbmag*randr*delx/r;
268 //fy = Fbmag*randr*dely/r;
269 //fz = Fbmag*randr*delz/r;
270
271 f[j][0] += fx;
272 f[j][1] += fy;
273 f[j][2] += fz;
274 }
275
276 // torque due to the Brownian Force
277
278 if (flaglog) {
279
280 // location of the point of closest approach on I from its center
281
282 xl[0] = -delx/r*radi;
283 xl[1] = -dely/r*radi;
284 xl[2] = -delz/r*radi;
285
286 // torque = xl_cross_F
287
288 tx = xl[1]*fz - xl[2]*fy;
289 ty = xl[2]*fx - xl[0]*fz;
290 tz = xl[0]*fy - xl[1]*fx;
291
292 // torque is same on both particles
293
294 torque[i][0] -= tx;
295 torque[i][1] -= ty;
296 torque[i][2] -= tz;
297
298 if (newton_pair || j < nlocal) {
299 torque[j][0] -= tx;
300 torque[j][1] -= ty;
301 torque[j][2] -= tz;
302 }
303
304 // torque due to a_pu
305
306 Fbmag = prethermostat*sqrt(a_pu);
307
308 // force in each direction
309
310 randr = random->uniform()-0.5;
311 tx = Fbmag*randr*p2[0];
312 ty = Fbmag*randr*p2[1];
313 tz = Fbmag*randr*p2[2];
314
315 randr = random->uniform()-0.5;
316 tx += Fbmag*randr*p3[0];
317 ty += Fbmag*randr*p3[1];
318 tz += Fbmag*randr*p3[2];
319
320 // torque has opposite sign on two particles
321
322 torque[i][0] -= tx;
323 torque[i][1] -= ty;
324 torque[i][2] -= tz;
325
326 if (newton_pair || j < nlocal) {
327 torque[j][0] += tx;
328 torque[j][1] += ty;
329 torque[j][2] += tz;
330 }
331 }
332
333 if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
334 0.0,0.0,-fx,-fy,-fz,delx,dely,delz);
335 }
336 }
337 }
338
339 int print_overlaps = 0;
340 if (print_overlaps && overlaps)
341 printf("Number of overlaps=%d\n",overlaps);
342
343 if (vflag_fdotr) virial_fdotr_compute();
344 }
345
346 /* ----------------------------------------------------------------------
347 allocate all arrays
348 ------------------------------------------------------------------------- */
349
allocate()350 void PairBrownian::allocate()
351 {
352 allocated = 1;
353 int n = atom->ntypes;
354
355 memory->create(setflag,n+1,n+1,"pair:setflag");
356 for (int i = 1; i <= n; i++)
357 for (int j = i; j <= n; j++)
358 setflag[i][j] = 0;
359
360 memory->create(cutsq,n+1,n+1,"pair:cutsq");
361
362 memory->create(cut,n+1,n+1,"pair:cut");
363 memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
364 }
365
366 /* ----------------------------------------------------------------------
367 global settings
368 ------------------------------------------------------------------------- */
369
settings(int narg,char ** arg)370 void PairBrownian::settings(int narg, char **arg)
371 {
372 if (narg != 7 && narg != 9) error->all(FLERR,"Illegal pair_style command");
373
374 mu = utils::numeric(FLERR,arg[0],false,lmp);
375 flaglog = utils::inumeric(FLERR,arg[1],false,lmp);
376 flagfld = utils::inumeric(FLERR,arg[2],false,lmp);
377 cut_inner_global = utils::numeric(FLERR,arg[3],false,lmp);
378 cut_global = utils::numeric(FLERR,arg[4],false,lmp);
379 t_target = utils::numeric(FLERR,arg[5],false,lmp);
380 seed = utils::inumeric(FLERR,arg[6],false,lmp);
381
382 flagHI = flagVF = 1;
383 if (narg == 9) {
384 flagHI = utils::inumeric(FLERR,arg[7],false,lmp);
385 flagVF = utils::inumeric(FLERR,arg[8],false,lmp);
386 }
387
388 if (flaglog == 1 && flagHI == 0) {
389 error->warning(FLERR,"Cannot include log terms without 1/r terms; "
390 "setting flagHI to 1");
391 flagHI = 1;
392 }
393
394 // initialize Marsaglia RNG with processor-unique seed
395
396 delete random;
397 random = new RanMars(lmp,seed + comm->me);
398
399 // reset cutoffs that have been explicitly set
400
401 if (allocated) {
402 for (int i = 1; i <= atom->ntypes; i++)
403 for (int j = i; j <= atom->ntypes; j++)
404 if (setflag[i][j]) {
405 cut_inner[i][j] = cut_inner_global;
406 cut[i][j] = cut_global;
407 }
408 }
409 }
410
411 /* ----------------------------------------------------------------------
412 set coeffs for one or more type pairs
413 ------------------------------------------------------------------------- */
414
coeff(int narg,char ** arg)415 void PairBrownian::coeff(int narg, char **arg)
416 {
417 if (narg != 2 && narg != 4)
418 error->all(FLERR,"Incorrect args for pair coefficients");
419
420 if (!allocated) allocate();
421
422 int ilo,ihi,jlo,jhi;
423 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
424 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
425
426 double cut_inner_one = cut_inner_global;
427 double cut_one = cut_global;
428
429 if (narg == 4) {
430 cut_inner_one = utils::numeric(FLERR,arg[2],false,lmp);
431 cut_one = utils::numeric(FLERR,arg[3],false,lmp);
432 }
433
434 int count = 0;
435 for (int i = ilo; i <= ihi; i++)
436 for (int j = MAX(jlo,i); j <= jhi; j++) {
437 cut_inner[i][j] = cut_inner_one;
438 cut[i][j] = cut_one;
439 setflag[i][j] = 1;
440 count++;
441 }
442
443 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
444 }
445
446 /* ----------------------------------------------------------------------
447 init specific to this pair style
448 ------------------------------------------------------------------------- */
449
init_style()450 void PairBrownian::init_style()
451 {
452 if (!atom->sphere_flag)
453 error->all(FLERR,"Pair brownian requires atom style sphere");
454
455 // if newton off, forces between atoms ij will be double computed
456 // using different random numbers
457
458 if (force->newton_pair == 0 && comm->me == 0)
459 error->warning(FLERR,
460 "Pair brownian needs newton pair on for "
461 "momentum conservation");
462
463 neighbor->request(this,instance_me);
464
465 // insure all particles are finite-size
466 // for pair hybrid, should limit test to types using the pair style
467
468 double *radius = atom->radius;
469 int nlocal = atom->nlocal;
470
471 for (int i = 0; i < nlocal; i++)
472 if (radius[i] == 0.0)
473 error->one(FLERR,"Pair brownian requires extended particles");
474
475 // require monodisperse system with same radii for all types
476
477 double radtype;
478 for (int i = 1; i <= atom->ntypes; i++) {
479 if (!atom->radius_consistency(i,radtype))
480 error->all(FLERR,"Pair brownian requires monodisperse particles");
481 if (i > 1 && radtype != rad)
482 error->all(FLERR,"Pair brownian requires monodisperse particles");
483 rad = radtype;
484 }
485
486 // set the isotropic constants that depend on the volume fraction
487 // vol_T = total volume
488 // check for fix deform, if exists it must use "remap v"
489 // If box will change volume, set appropriate flag so that volume
490 // and v.f. corrections are re-calculated at every step.
491 //
492 // If available volume is different from box volume
493 // due to walls, set volume appropriately; if walls will
494 // move, set appropriate flag so that volume and v.f. corrections
495 // are re-calculated at every step.
496
497 flagdeform = flagwall = 0;
498 for (int i = 0; i < modify->nfix; i++) {
499 if (strcmp(modify->fix[i]->style,"deform") == 0)
500 flagdeform = 1;
501 else if (strstr(modify->fix[i]->style,"wall") != nullptr) {
502 if (flagwall)
503 error->all(FLERR,
504 "Cannot use multiple fix wall commands with pair brownian");
505 flagwall = 1; // Walls exist
506 wallfix = (FixWall *) modify->fix[i];
507 if (wallfix->xflag) flagwall = 2; // Moving walls exist
508 }
509 }
510
511 // set the isotropic constants depending on the volume fraction
512 // vol_T = total volumeshearing = flagdeform = flagwall = 0;
513
514 double vol_T, wallcoord;
515 if (!flagwall) vol_T = domain->xprd*domain->yprd*domain->zprd;
516 else {
517 double wallhi[3], walllo[3];
518 for (int j = 0; j < 3; j++) {
519 wallhi[j] = domain->prd[j];
520 walllo[j] = 0;
521 }
522 for (int m = 0; m < wallfix->nwall; m++) {
523 int dim = wallfix->wallwhich[m] / 2;
524 int side = wallfix->wallwhich[m] % 2;
525 if (wallfix->xstyle[m] == VARIABLE) {
526 wallfix->xindex[m] = input->variable->find(wallfix->xstr[m]);
527 // Since fix->wall->init happens after pair->init_style
528 wallcoord = input->variable->compute_equal(wallfix->xindex[m]);
529 }
530
531 else wallcoord = wallfix->coord0[m];
532
533 if (side == 0) walllo[dim] = wallcoord;
534 else wallhi[dim] = wallcoord;
535 }
536 vol_T = (wallhi[0] - walllo[0]) * (wallhi[1] - walllo[1]) *
537 (wallhi[2] - walllo[2]);
538 }
539
540 // vol_P = volume of particles, assuming mono-dispersity
541 // vol_f = volume fraction
542
543 vol_P = atom->natoms*(4.0/3.0)*MY_PI*cube(rad);
544
545 double vol_f = vol_P/vol_T;
546
547 // set isotropic constants
548 if (!flagVF) vol_f = 0;
549
550 if (flaglog == 0) {
551 R0 = 6*MY_PI*mu*rad*(1.0 + 2.16*vol_f);
552 RT0 = 8*MY_PI*mu*cube(rad); // not actually needed
553 } else {
554 R0 = 6*MY_PI*mu*rad*(1.0 + 2.725*vol_f - 6.583*vol_f*vol_f);
555 RT0 = 8*MY_PI*mu*cube(rad)*(1.0 + 0.749*vol_f - 2.469*vol_f*vol_f);
556 }
557 }
558
559 /* ----------------------------------------------------------------------
560 init for one type pair i,j and corresponding j,i
561 ------------------------------------------------------------------------- */
562
init_one(int i,int j)563 double PairBrownian::init_one(int i, int j)
564 {
565 if (setflag[i][j] == 0) {
566 cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
567 cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
568 }
569
570 cut_inner[j][i] = cut_inner[i][j];
571
572 return cut[i][j];
573 }
574
575 /* ----------------------------------------------------------------------
576 proc 0 writes to restart file
577 ------------------------------------------------------------------------- */
578
write_restart(FILE * fp)579 void PairBrownian::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(&cut_inner[i][j],sizeof(double),1,fp);
589 fwrite(&cut[i][j],sizeof(double),1,fp);
590 }
591 }
592 }
593
594 /* ----------------------------------------------------------------------
595 proc 0 reads from restart file, bcasts
596 ------------------------------------------------------------------------- */
597
read_restart(FILE * fp)598 void PairBrownian::read_restart(FILE *fp)
599 {
600 read_restart_settings(fp);
601 allocate();
602
603 int i,j;
604 int me = comm->me;
605 for (i = 1; i <= atom->ntypes; i++)
606 for (j = i; j <= atom->ntypes; j++) {
607 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
608 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
609 if (setflag[i][j]) {
610 if (me == 0) {
611 utils::sfread(FLERR,&cut_inner[i][j],sizeof(double),1,fp,nullptr,error);
612 utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
613 }
614 MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
615 MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
616 }
617 }
618 }
619
620 /* ----------------------------------------------------------------------
621 proc 0 writes to restart file
622 ------------------------------------------------------------------------- */
623
write_restart_settings(FILE * fp)624 void PairBrownian::write_restart_settings(FILE *fp)
625 {
626 fwrite(&mu,sizeof(double),1,fp);
627 fwrite(&flaglog,sizeof(int),1,fp);
628 fwrite(&flagfld,sizeof(int),1,fp);
629 fwrite(&cut_inner_global,sizeof(double),1,fp);
630 fwrite(&cut_global,sizeof(double),1,fp);
631 fwrite(&t_target,sizeof(double),1,fp);
632 fwrite(&seed,sizeof(int),1,fp);
633 fwrite(&offset_flag,sizeof(int),1,fp);
634 fwrite(&mix_flag,sizeof(int),1,fp);
635 fwrite(&flagHI,sizeof(int),1,fp);
636 fwrite(&flagVF,sizeof(int),1,fp);
637 }
638
639 /* ----------------------------------------------------------------------
640 proc 0 reads from restart file, bcasts
641 ------------------------------------------------------------------------- */
642
read_restart_settings(FILE * fp)643 void PairBrownian::read_restart_settings(FILE *fp)
644 {
645 int me = comm->me;
646 if (me == 0) {
647 utils::sfread(FLERR,&mu,sizeof(double),1,fp,nullptr,error);
648 utils::sfread(FLERR,&flaglog,sizeof(int),1,fp,nullptr,error);
649 utils::sfread(FLERR,&flagfld,sizeof(int),1,fp,nullptr,error);
650 utils::sfread(FLERR,&cut_inner_global,sizeof(double),1,fp,nullptr,error);
651 utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
652 utils::sfread(FLERR,&t_target, sizeof(double),1,fp,nullptr,error);
653 utils::sfread(FLERR,&seed, sizeof(int),1,fp,nullptr,error);
654 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
655 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
656 utils::sfread(FLERR,&flagHI,sizeof(int),1,fp,nullptr,error);
657 utils::sfread(FLERR,&flagVF,sizeof(int),1,fp,nullptr,error);
658 }
659 MPI_Bcast(&mu,1,MPI_DOUBLE,0,world);
660 MPI_Bcast(&flaglog,1,MPI_INT,0,world);
661 MPI_Bcast(&flagfld,1,MPI_INT,0,world);
662 MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
663 MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
664 MPI_Bcast(&t_target,1,MPI_DOUBLE,0,world);
665 MPI_Bcast(&seed,1,MPI_INT,0,world);
666 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
667 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
668 MPI_Bcast(&flagHI,1,MPI_INT,0,world);
669 MPI_Bcast(&flagVF,1,MPI_INT,0,world);
670
671 // additional setup based on restart parameters
672
673 delete random;
674 random = new RanMars(lmp,seed + comm->me);
675 }
676
677 /* ----------------------------------------------------------------------*/
678
set_3_orthogonal_vectors(double p1[3],double p2[3],double p3[3])679 void PairBrownian::set_3_orthogonal_vectors(double p1[3],
680 double p2[3], double p3[3])
681 {
682 double norm;
683 int ix,iy,iz;
684
685 // find the index of maximum magnitude and store it in iz
686
687 if (fabs(p1[0]) > fabs(p1[1])) {
688 iz=0;
689 ix=1;
690 iy=2;
691 } else {
692 iz=1;
693 ix=2;
694 iy=0;
695 }
696
697 if (iz==0) {
698 if (fabs(p1[0]) < fabs(p1[2])) {
699 iz = 2;
700 ix = 0;
701 iy = 1;
702 }
703 } else {
704 if (fabs(p1[1]) < fabs(p1[2])) {
705 iz = 2;
706 ix = 0;
707 iy = 1;
708 }
709 }
710
711 // set p2 arbitrarily such that it's orthogonal to p1
712
713 p2[ix]=1.0;
714 p2[iy]=1.0;
715 p2[iz] = -(p1[ix]*p2[ix] + p1[iy]*p2[iy])/p1[iz];
716
717 // normalize p2
718
719 norm = sqrt(p2[0]*p2[0] + p2[1]*p2[1] + p2[2]*p2[2]);
720
721 p2[0] = p2[0]/norm;
722 p2[1] = p2[1]/norm;
723 p2[2] = p2[2]/norm;
724
725 // Set p3 by taking the cross product p3=p2xp1
726
727 p3[0] = p1[1]*p2[2] - p1[2]*p2[1];
728 p3[1] = p1[2]*p2[0] - p1[0]*p2[2];
729 p3[2] = p1[0]*p2[1] - p1[1]*p2[0];
730 }
731