1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 /* -------------------------------------------------------------------------
43 Thanks to Chris Stoltz (P&G) for providing
44 a Fortran version of the MC integrator
45 ------------------------------------------------------------------------- */
46
47 #include <cmath>
48 #include <stdio.h>
49 #include <stdlib.h>
50 #include <string.h>
51 #include "fix_template_multiplespheres.h"
52 #include "fix_property_atom.h"
53 #include "math_extra.h"
54 #include "math_extra_liggghts.h"
55 #include "vector_liggghts.h"
56 #include "atom.h"
57 #include "atom_vec.h"
58 #include "domain.h"
59 #include "modify.h"
60 #include "comm.h"
61 #include "force.h"
62 #include "update.h"
63 #include "output.h"
64 #include "memory.h"
65 #include "error.h"
66 #include "random_mars.h"
67 #include "random_park.h"
68 #include "fix_rigid.h"
69 #include "particleToInsert.h"
70 #include "input_multisphere.h"
71
72 using namespace LAMMPS_NS;
73 using namespace LMP_PROBABILITY_NS;
74
75 #define LARGE 1e8
76 #define EPSILON 1.0e-7
77 #define N_SHUFFLE_BOUND 200
78
79 /* ---------------------------------------------------------------------- */
80
FixTemplateMultiplespheres(LAMMPS * lmp,int narg,char ** arg)81 FixTemplateMultiplespheres::FixTemplateMultiplespheres(LAMMPS *lmp, int narg, char **arg) :
82 FixTemplateSphere(lmp, narg, arg),
83 scale_fact(1.0)
84 {
85 if(pdf_density->rand_style() != RANDOM_CONSTANT) error->all(FLERR,"Fix particletemplate/multiplespheres currently only supports constant density");
86 if(pdf_radius) error->fix_error(FLERR,this,"currently does not support keyword 'radius'");
87 if(domain->dimension != 3) error->fix_error(FLERR,this,"only supports 3D simulations");
88
89 // parse number of spheres
90 if (strcmp(arg[iarg++],"nspheres") != 0) error->fix_error(FLERR,this,"expecting argument 'nspheres'");
91 nspheres = atoi(arg[iarg++]);
92 if(nspheres < 2) error->fix_error(FLERR,this,"illegal number of spheres");
93
94 // allocate arrays
95 memory->create(x_sphere,nspheres,3,"FixTemplateMultiplespheres:x_sphere");
96 r_sphere = new double[nspheres];
97 atom_type_sphere = 0;
98
99 // re-create pti with correct nspheres
100 delete pti;
101 pti = new ParticleToInsert(lmp,nspheres);
102
103 bonded = false;
104 fix_bond_random_id = 0;
105
106 for (int i = 0; i < 3; i++) {
107 x_min[i] = LARGE;
108 x_max[i] = -LARGE;
109 }
110
111 overlap_slightly = no_overlap = false;
112
113 // parse args
114
115 if (narg < iarg+4) error->fix_error(FLERR,this,"not enough arguments");
116 if (strcmp(arg[iarg++],"ntry") != 0) error->fix_error(FLERR,this,"need 'ntry' to be defined");
117 ntry = static_cast<int>(atoi(arg[iarg++]));
118
119 bool spheres_read = false;
120
121 bool hasargs = true;
122 while (iarg < narg && hasargs)
123 {
124 hasargs = false;
125
126 if ((strcmp(arg[iarg],"spheres") == 0) || (strcmp(arg[iarg],"spheres_different_types") == 0))
127 {
128 bool different_type = false;
129 if(strcmp(arg[iarg],"spheres_different_types") == 0)
130 different_type= true;
131
132 hasargs = true;
133 spheres_read = true;
134 iarg++;
135
136 if (strcmp(arg[iarg],"file") == 0)
137 {
138 iarg++;
139 if (narg < iarg+3) error->fix_error(FLERR,this,"not enough arguments for 'file'");
140
141 if(different_type)
142 atom_type_sphere = new int[nspheres];
143
144 char *clmp_filename = arg[iarg++];
145
146 if (strcmp(arg[iarg++],"scale") != 0) error->fix_error(FLERR,this,"you have to specify a scale factor");
147 scale_fact = atof(arg[iarg++]);
148 if(scale_fact<=0.) error->fix_error(FLERR,this,"scale factor must be >0");
149
150 // allocate input class, try to open file, read data from file
151 InputMultisphere *myclmp_input = new InputMultisphere(lmp,0,NULL);
152 myclmp_input->clmpfile(clmp_filename,x_sphere,r_sphere,atom_type_sphere,nspheres);
153 delete myclmp_input;
154
155 for(int i = 0; i < nspheres; i++)
156 {
157 if(r_sphere[i] <= 0.) error->fix_error(FLERR,this,"radius must be > 0");
158 if(different_type)
159 {
160 r_sphere[i] *= (scale_fact*force->cg(atom_type_sphere[i]));
161 vectorScalarMult3D(x_sphere[i],scale_fact*force->cg(atom_type_sphere[i]));
162 }
163 else
164 {
165 r_sphere[i] *= (scale_fact*force->cg(atom_type));
166 vectorScalarMult3D(x_sphere[i],scale_fact*force->cg(atom_type));
167 }
168 }
169
170 // calculate bounding box
171 for(int i = 0; i < nspheres; i++)
172 {
173 for(int j=0;j<3;j++)
174 {
175 if (x_sphere[i][j]-r_sphere[i]<x_min[j]) x_min[j] = x_sphere[i][j]-r_sphere[i];
176 if (x_sphere[i][j]+r_sphere[i]>x_max[j]) x_max[j] = x_sphere[i][j]+r_sphere[i];
177 }
178 }
179 }
180 else
181 {
182 if (narg < iarg + 4*nspheres) error->fix_error(FLERR,this,"not enough arguments");
183
184 if(different_type)
185 error->fix_error(FLERR,this,"have to use keyword 'file' with option 'spheres_different_type'");
186
187 //read sphere r and coos, determine min and max
188 for(int i = 0; i < nspheres; i++)
189 {
190 if(different_type)
191 r_sphere[i] = atof(arg[iarg+3])*force->cg(atom_type_sphere[i]);
192 else
193 r_sphere[i] = atof(arg[iarg+3])*force->cg(atom_type);
194 if(r_sphere[i] <= 0.) error->fix_error(FLERR,this,"radius must be >0");
195 for(int j = 0; j < 3; j++)
196 {
197 if(different_type)
198 x_sphere[i][j] = atof(arg[iarg+j])*force->cg(atom_type_sphere[i]);
199 else
200 x_sphere[i][j] = atof(arg[iarg+j])*force->cg(atom_type);
201 if (x_sphere[i][j]-r_sphere[i]<x_min[j]) x_min[j]=x_sphere[i][j]-r_sphere[i];
202 if (x_sphere[i][j]+r_sphere[i]>x_max[j]) x_max[j]=x_sphere[i][j]+r_sphere[i];
203 }
204 iarg+=4;
205 }
206 }
207 }
208 else if(strcmp(arg[iarg],"bonded") == 0)
209 {
210 if (narg < iarg+2)
211 error->fix_error(FLERR,this,"not enough arguments for 'bonded'");
212 if(0 == strcmp(arg[iarg+1],"yes"))
213 bonded = true;
214 else if(0 == strcmp(arg[iarg+1],"no"))
215 bonded = false;
216 else
217 error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'bonded'");
218 iarg+=2;
219 }
220 else if(strcmp(style,"particletemplate/multiplespheres") == 0)
221 error->fix_error(FLERR,this,"unknown keyword");
222 }
223
224 if(!spheres_read) error->fix_error(FLERR,this,"need to define spheres for the template");
225
226 if(comm->me == 0 && screen) fprintf(screen,"Calculating the properties of the given template.\n Depending on ntry, this may take a while...\n");
227
228 if(ntry < 1e3) error->fix_error(FLERR,this,"ntry is too low");
229 if(comm->me == 0 && ntry < 1e5) error->warning(FLERR,"fix particletemplate/multisphere: ntry is very low");
230
231 }
232
233 /* ---------------------------------------------------------------------- */
234
~FixTemplateMultiplespheres()235 FixTemplateMultiplespheres::~FixTemplateMultiplespheres()
236 {
237 memory->destroy(x_sphere);
238 delete []r_sphere;
239 if(atom_type_sphere) delete []atom_type_sphere;
240 }
241
242 /* ---------------------------------------------------------------------- */
243
post_create()244 void FixTemplateMultiplespheres::post_create()
245 {
246 // calculate bounding sphere and center of mass
247 // also transforms sphere coordinates so that com = 0/0/0
248
249 calc_bounding_sphere();
250 calc_center_of_mass();
251
252 // check amount of overlap
253 // needed for some functionalities
254 check_overlap();
255
256 if(0 == strcmp(style,"particletemplate/multiplespheres"))
257 print_info();
258
259 if(bonded && !fix_bond_random_id)
260 {
261
262 fix_bond_random_id = static_cast<FixPropertyAtom*>(modify->find_fix_property("bond_random_id","property/atom","scalar",0,0,this->style,false));
263
264 if(!fix_bond_random_id)
265 {
266 const char *fixarg[] = {
267 "bond_random_id", // fix id
268 "all", // fix group
269 "property/atom", // fix style: property/atom
270 "bond_random_id", // property name
271 "scalar", // 1 vector per particle
272 "yes", // restart
273 "yes", // communicate ghost
274 "no", // communicate rev
275 "-1."
276 };
277 fix_bond_random_id = modify->add_fix_property_atom(9,const_cast<char**>(fixarg),style);
278
279 }
280 }
281 }
282
283 /* ---------------------------------------------------------------------- */
284
print_info()285 void FixTemplateMultiplespheres::print_info()
286 {
287 if (logfile)
288 {
289 fprintf(logfile,"Finished calculating properties of template\n");
290 fprintf(logfile," mass = %e, radius of bounding sphere = %e, radius of equivalent sphere = %e\n",mass_expect,r_bound,r_equiv);
291 fprintf(logfile," center of mass = %e, %e, %e\n",0.,0.,0.);
292 fprintf(logfile," center of bounding sphere in global coords = %e, %e, %e\n",x_bound[0],x_bound[1],x_bound[2]);
293 }
294 }
295
296 /* ----------------------------------------------------------------------*/
297
maxtype()298 int FixTemplateMultiplespheres::maxtype()
299 {
300 if(!atom_type_sphere)
301 return atom_type;
302 return vectorMaxN(atom_type_sphere,nspheres);
303 }
304
305 /* ----------------------------------------------------------------------*/
306
mintype()307 int FixTemplateMultiplespheres::mintype()
308 {
309 if(!atom_type_sphere)
310 return atom_type;
311 return vectorMinN(atom_type_sphere,nspheres);
312 }
313
314 /* ----------------------------------------------------------------------
315 calc bounding sphere with iterative procedure
316
317 do this multiple times, randomizing point order every time, see
318 http://hacksoflife.blogspot.com/2009/01/randomized-bounding-spheres.html
319 choose optimal result at the end - this gives linear run-time
320 ------------------------------------------------------------------------- */
321
calc_bounding_sphere()322 void FixTemplateMultiplespheres::calc_bounding_sphere()
323 {
324 r_bound = LARGE;
325 int *visited = new int[nspheres];
326 double d[3],dist;
327
328 for(int shuffle = 0; shuffle < N_SHUFFLE_BOUND; shuffle ++)
329 {
330 for(int i = 0; i < nspheres; i++) visited[i] = 0;
331
332 int isphere = -1;
333 int nvisited = 0;
334 double x_bound_temp[3],rbound_temp;
335
336 while(isphere < 0 || visited[isphere] || isphere >= nspheres )
337 isphere = static_cast<int>(random_mc->uniform()*nspheres);
338
339 nvisited++;
340 visited[isphere] = 1;
341
342 vectorCopy3D(x_sphere[isphere],x_bound_temp);
343 rbound_temp = r_sphere[isphere];
344
345 while(nvisited < nspheres)
346 {
347 while(isphere < 0 || visited[isphere] || isphere >= nspheres )
348 isphere = static_cast<int>(random_mc->uniform()*nspheres);
349
350 nvisited++;
351 visited[isphere] = 1;
352
353 vectorSubtract3D(x_sphere[isphere],x_bound_temp,d);
354 dist = vectorMag3D(d);
355
356 // do nothing if sphere is completely contained in bounding sphere
357 // if not contained, shift and extend bounding sphere
358 if(dist + r_sphere[isphere] > rbound_temp)
359 {
360 double fact = (dist + r_sphere[isphere] - rbound_temp) / (2. * dist);
361 vectorScalarMult3D(d,fact);
362 vectorAdd3D(x_bound_temp,d,x_bound_temp);
363 rbound_temp += vectorMag3D(d);
364 }
365
366 }
367 if(rbound_temp < r_bound)
368 {
369 r_bound = rbound_temp;
370 vectorCopy3D(x_bound_temp,x_bound);
371 }
372
373 }
374 delete []visited;
375
376 // do a coarse check on the validity of the bounding sphere calc
377 for(int i = 0; i < nspheres; i++)
378 {
379 double temp[3];
380 vectorSubtract3D(x_bound,x_sphere[i],temp);
381 if(vectorMag3D(temp) > r_bound) error->fix_error(FLERR,this,"Bounding sphere calculation for template failed");
382 }
383
384 }
385
386 /* ----------------------------------------------------------------------
387 check overlap of spheres against each other
388 ------------------------------------------------------------------------- */
389
check_overlap()390 void FixTemplateMultiplespheres::check_overlap()
391 {
392
393 double overlap_slightly_min = 1.0000001;
394 double overlap_slightly_max = 1.01;
395
396 overlap_slightly = true;
397 no_overlap = true;
398
399 double dist;
400 bool *overlap_slightly_vec = new bool[nspheres];
401 vectorInitializeN(overlap_slightly_vec,nspheres,false);
402
403 for(int i = 0; i < nspheres; i++)
404 {
405 for(int j = i+1; j < nspheres; j++)
406 {
407 dist = pointDistance(x_sphere[i],x_sphere[j]);
408
409 if(dist < (r_sphere[i] + r_sphere[j]))
410 no_overlap = false;
411
412 if( (dist > overlap_slightly_min*(r_sphere[i] + r_sphere[j])) &&
413 (dist < overlap_slightly_max*(r_sphere[i] + r_sphere[j]))
414 )
415 {
416 overlap_slightly_vec[i] = true;
417 overlap_slightly_vec[j] = true;
418 }
419 }
420 }
421
422 // if any sphere has no slight overlap, result is false
423 for(int i = 0; i < nspheres; i++)
424 if(!overlap_slightly_vec[i])
425 overlap_slightly = false;
426
427 delete []overlap_slightly_vec;
428 }
429
430 /* ----------------------------------------------------------------------
431 sqr distance from x_sphere[j] to xtest
432 ------------------------------------------------------------------------- */
433
dist_sqr(int j,double * xtest)434 double FixTemplateMultiplespheres::dist_sqr(int j, double *xtest)
435 {
436 double dSqr = 0.;
437 dSqr += (xtest[0]-x_sphere[j][0])*(xtest[0]-x_sphere[j][0]);
438 dSqr += (xtest[1]-x_sphere[j][1])*(xtest[1]-x_sphere[j][1]);
439 dSqr += (xtest[2]-x_sphere[j][2])*(xtest[2]-x_sphere[j][2]);
440 return dSqr;
441 }
442
443 /* ----------------------------------------------------------------------
444 generate random point in bbox
445 ------------------------------------------------------------------------- */
446
generate_xtry(double * x_try)447 void FixTemplateMultiplespheres::generate_xtry(double *x_try)
448 {
449 x_try[0] = x_min[0]+(x_max[0]-x_min[0])*random_mc->uniform();
450 x_try[1] = x_min[1]+(x_max[1]-x_min[1])*random_mc->uniform();
451 x_try[2] = x_min[2]+(x_max[2]-x_min[2])*random_mc->uniform();
452 }
453
454 /* ----------------------------------------------------------------------
455 calc center of mass
456 ------------------------------------------------------------------------- */
457
calc_center_of_mass()458 void FixTemplateMultiplespheres::calc_center_of_mass()
459 {
460 // mc integration, calc volume and com, mass weight
461 int nsuccess = 0;
462
463 double x_try[3],xcm[3],dist_j_sqr;
464
465 vectorZeroize3D(xcm);
466
467 bool alreadyChecked = false;
468 for(int i = 0; i < ntry; i++)
469 {
470 generate_xtry(x_try);
471
472 alreadyChecked = false;
473 for(int j = 0; j < nspheres; j++)
474 {
475 dist_j_sqr = dist_sqr(j,x_try);
476
477 // only count once if contained in multiple spheres
478 if (alreadyChecked) break;
479 if(dist_j_sqr < r_sphere[j]*r_sphere[j])
480 {
481 xcm[0] = (xcm[0]*static_cast<double>(nsuccess)+x_try[0])/static_cast<double>(nsuccess+1);
482 xcm[1] = (xcm[1]*static_cast<double>(nsuccess)+x_try[1])/static_cast<double>(nsuccess+1);
483 xcm[2] = (xcm[2]*static_cast<double>(nsuccess)+x_try[2])/static_cast<double>(nsuccess+1);
484 nsuccess++;
485 alreadyChecked = true;
486 }
487 }
488 }
489
490 // expectancy values
491 volume_expect = static_cast<double>(nsuccess)/static_cast<double>(ntry)*(x_max[0]-x_min[0])*(x_max[1]-x_min[1])*(x_max[2]-x_min[2]);
492 mass_expect = volume_expect*expectancy(pdf_density);
493 r_equiv = pow(6.*mass_expect/(8.*expectancy(pdf_density)*M_PI),1./3.);
494
495 // transform into a system with center of mass=0/0/0
496
497 for(int i = 0; i < nspheres; i++)
498 vectorSubtract3D(x_sphere[i],xcm,x_sphere[i]);
499
500 vectorSubtract3D(x_min,xcm,x_min);
501 vectorSubtract3D(x_max,xcm,x_max);
502 vectorSubtract3D(x_bound,xcm,x_bound);
503
504 }
505
506 /* ----------------------------------------------------------------------*/
507
max_r_bound()508 double FixTemplateMultiplespheres::max_r_bound()
509 {
510 return r_bound;
511 }
512
513 /* ----------------------------------------------------------------------*/
514
min_rad()515 double FixTemplateMultiplespheres::min_rad()
516 {
517 double rmin = 100000000.;
518
519 for(int j = 0; j < nspheres; j++)
520 if(rmin > r_sphere[j]) rmin = r_sphere[j];
521
522 return rmin;
523 }
524
525 /* ----------------------------------------------------------------------*/
526
max_rad()527 double FixTemplateMultiplespheres::max_rad()
528 {
529 double rmax = 0.;
530
531 for(int j = 0; j < nspheres; j++)
532 if(rmax < r_sphere[j]) rmax = r_sphere[j];
533
534 return rmax;
535 }
536
537 /* ----------------------------------------------------------------------*/
538
number_spheres()539 int FixTemplateMultiplespheres::number_spheres()
540 {
541 return nspheres;
542 }
543
544 /* ----------------------------------------------------------------------*/
545
randomize_single()546 void FixTemplateMultiplespheres::randomize_single()
547 {
548
549 pti->nparticles = nspheres;
550 pti->density_ins = expectancy(pdf_density);
551 pti->volume_ins = volume_expect;
552 pti->mass_ins = mass_expect;
553 pti->r_bound_ins = r_bound;
554 vectorCopy3D(x_bound,pti->x_bound_ins);
555 pti->atom_type = atom_type;
556 if(atom_type_sphere)
557 {
558 vectorCopyN(atom_type_sphere,pti->atom_type_vector,nspheres);
559 pti->atom_type_vector_flag = true;
560 }
561
562 for(int j = 0; j < nspheres; j++)
563 {
564 pti->radius_ins[j] = r_sphere[j];
565 vectorCopy3D(x_sphere[j],pti->x_ins[j]);
566 }
567
568 vectorZeroize3D(pti->v_ins);
569 vectorZeroize3D(pti->omega_ins);
570
571 pti->groupbit = groupbit;
572 }
573
574 /* ----------------------------------------------------------------------*/
575
init_ptilist(int n_random_max,const bool enforce_single,FixPropertyAtom * const fix_release)576 void FixTemplateMultiplespheres::init_ptilist(int n_random_max, const bool enforce_single, FixPropertyAtom * const fix_release)
577 {
578 if(pti_list) error->all(FLERR,"invalid FixTemplateSphere::init_list()");
579 n_pti_max = n_random_max;
580 pti_list = (ParticleToInsert**) memory->smalloc(n_pti_max*sizeof(ParticleToInsert*),"pti_list");
581 for(int i = 0; i < n_pti_max; i++)
582 pti_list[i] = new ParticleToInsert(lmp, enforce_single ? 1 : nspheres, fix_release);
583 }
584
585 /* ----------------------------------------------------------------------*/
586
randomize_ptilist(int n_random,int distribution_groupbit,int distorder)587 void FixTemplateMultiplespheres::randomize_ptilist(int n_random,int distribution_groupbit,int distorder)
588 {
589
590 for(int i = 0; i < n_random; i++)
591 {
592 ParticleToInsert *pti = pti_list[i];
593
594 pti->density_ins = expectancy(pdf_density);
595 pti->volume_ins = volume_expect;
596 pti->mass_ins = mass_expect;
597 pti->r_bound_ins = r_bound;
598 vectorCopy3D(x_bound,pti->x_bound_ins);
599 pti->atom_type = atom_type;
600 if(atom_type_sphere)
601 {
602 vectorCopyN(atom_type_sphere,pti->atom_type_vector,nspheres);
603 pti->atom_type_vector_flag = true;
604 }
605
606 for(int j = 0; j < nspheres; j++)
607 {
608 pti->radius_ins[j] = r_sphere[j];
609 vectorCopy3D(x_sphere[j],pti->x_ins[j]);
610 }
611
612 vectorZeroize3D(pti->v_ins);
613 vectorZeroize3D(pti->omega_ins);
614
615 pti->groupbit = groupbit | distribution_groupbit;
616
617 pti_list[i]->distorder = distorder;
618
619 if(bonded)
620 {
621 if (!pti_list[i]->fix_property)
622 {
623 pti_list[i]->fix_property = new FixPropertyAtom*[1];
624 if (pti_list[i]->fix_property_value)
625 error->one(FLERR, "Internal error (fix property pti list)");
626 pti_list[i]->fix_property_value = new double*[1];
627 pti_list[i]->fix_property_value[0] = new double[1];
628 if (pti_list[i]->fix_property_nentry)
629 error->one(FLERR, "Internal error (fix property pti list)");
630 pti_list[i]->fix_property_nentry = new int[1];
631 }
632 pti_list[i]->fix_property[0] = fix_bond_random_id;
633
634 pti_list[i]->fix_property_value[0][0] = static_cast<double>(update->ntimestep)+random_insertion->uniform();
635 pti_list[i]->n_fix_property = 1;
636 pti_list[i]->fix_property_nentry[0] = 1;
637 }
638 }
639 }
640
641 /* ----------------------------------------------------------------------*/
642
direct_set_ptlist(const int i,const void * const data,const int distribution_groupbit,const int distorder)643 void FixTemplateMultiplespheres::direct_set_ptlist(const int i, const void * const data, const int distribution_groupbit, const int distorder)
644 {
645 const PARTICLE_PACKING::MultipleSphere * const ms = static_cast<const PARTICLE_PACKING::MultipleSphere * const>(data);
646 pti_list[i]->atom_type = ms->get_type();
647 const double radius = ms->get_radius();
648
649 pti_list[i]->radius_ins[0] = radius;
650 pti_list[i]->density_ins = ms->get_density();
651 pti_list[i]->volume_ins = radius*radius*radius*4.1887902047863909;
652 pti_list[i]->mass_ins = pti_list[i]->density_ins*pti_list[i]->volume_ins;
653 pti_list[i]->id_ins = ms->get_id();
654
655 // set fix_property_atom
656 if (pti_list[i]->fix_property && ms->n_fix_properties() != (unsigned int)pti_list[i]->n_fix_property)
657 error->one(FLERR, "Inconsistent fix_property count");
658 if (pti_list[i]->fix_property_value)
659 {
660 if (!pti_list[i]->fix_property_nentry)
661 error->one(FLERR, "Nentry not available");
662 for (int j = 0; j < pti_list[i]->n_fix_property; j++)
663 {
664 if (ms->fix_property_nentries(j) != (unsigned int)pti_list[i]->fix_property_nentry[j])
665 error->one(FLERR, "Inconsistent fix property entries");
666 }
667 }
668
669 if (ms->n_fix_properties() > 0)
670 {
671 const int n = ms->n_fix_properties();
672 pti_list[i]->n_fix_property = n;
673 if (!pti_list[i]->fix_property)
674 pti_list[i]->fix_property = new FixPropertyAtom*[n];
675 const bool create_fpv = !pti_list[i]->fix_property_value;
676 if (create_fpv)
677 pti_list[i]->fix_property_value = new double*[n];
678 if (!pti_list[i]->fix_property_nentry)
679 pti_list[i]->fix_property_nentry = new int[n];
680 bool found_bonded = false;
681 for (int j = 0; j < n; j++)
682 {
683 pti_list[i]->fix_property[j] = ms->get_fix_property(j);
684 const int m = ms->fix_property_nentries(j);
685 if (create_fpv)
686 pti_list[i]->fix_property_value[j] = new double[m];
687 pti_list[i]->fix_property_nentry[j] = m;
688 for (int k = 0; k < m; k++)
689 pti_list[i]->fix_property_value[j][k] = ms->fix_property_value(j, k);
690 if (pti_list[i]->fix_property[j] == fix_bond_random_id)
691 {
692 found_bonded = true;
693
694 pti_list[i]->fix_property_value[j][0] += static_cast<double>(update->ntimestep);
695 }
696 }
697 if (bonded && !found_bonded)
698 error->one(FLERR, "Bond random id could not be found");
699 }
700
701 // init insertion position
702 vectorZeroize3D(pti_list[i]->x_ins[0]);
703 vectorZeroize3D(pti_list[i]->v_ins);
704 vectorZeroize3D(pti_list[i]->omega_ins);
705
706 pti_list[i]->groupbit = groupbit | distribution_groupbit;
707
708 pti_list[i]->distorder = distorder;
709 pti_list[i]->distorder = distorder;
710 }
711
712 /* ----------------------------------------------------------------------*/
713
generate_hash()714 unsigned int FixTemplateMultiplespheres::generate_hash()
715 {
716 unsigned int hash = 0;
717 unsigned int start = seed_orig*123457; // it's magic
718 if (atom_type_sphere)
719 {
720 for (int i = 0; i < nspheres; i++)
721 add_hash_value(atom_type_sphere[i], start, hash);
722 }
723 else
724 add_hash_value(atom_type, start, hash);
725 add_hash_value(nspheres, start, hash);
726 for (int i = 0; i < nspheres; i++)
727 add_hash_value(r_sphere[i], start, hash);
728 add_hash_value(pdf_density->rand_style(), start, hash);
729 add_hash_value(expectancy(pdf_density), start, hash);
730 add_hash_value(cubic_expectancy(pdf_density), start, hash);
731 add_hash_value(bonded ? 1 : 0, start, hash);
732 return hash;
733 }
734