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: Ray Shan (Sandia, tnshan@sandia.gov)
17 Oleg Sergeev (VNIIA, sergeev@vniia.ru)
18 ------------------------------------------------------------------------- */
19
20 #include "fix_reaxff_species.h"
21
22 #include "atom.h"
23 #include "comm.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "fix_ave_atom.h"
27 #include "force.h"
28 #include "memory.h"
29 #include "modify.h"
30 #include "neigh_list.h"
31 #include "neighbor.h"
32 #include "update.h"
33
34 #include "pair_reaxff.h"
35 #include "reaxff_defs.h"
36
37 #include <cstring>
38
39 using namespace LAMMPS_NS;
40 using namespace FixConst;
41
42 /* ---------------------------------------------------------------------- */
43
FixReaxFFSpecies(LAMMPS * lmp,int narg,char ** arg)44 FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
45 Fix(lmp, narg, arg)
46 {
47 if (narg < 7) error->all(FLERR,"Illegal fix reaxff/species command");
48
49 force_reneighbor = 0;
50
51 vector_flag = 1;
52 size_vector = 2;
53 extvector = 0;
54
55 peratom_flag = 1;
56 size_peratom_cols = 0;
57 peratom_freq = 1;
58
59 nvalid = -1;
60
61 MPI_Comm_rank(world,&me);
62 MPI_Comm_size(world,&nprocs);
63 ntypes = atom->ntypes;
64
65 nevery = atoi(arg[3]);
66 nrepeat = atoi(arg[4]);
67 global_freq = nfreq = atoi(arg[5]);
68
69 comm_forward = 4;
70
71 if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
72 error->all(FLERR,"Illegal fix reaxff/species command");
73 if (nfreq % nevery || nrepeat*nevery > nfreq)
74 error->all(FLERR,"Illegal fix reaxff/species command");
75
76 // Neighbor lists must stay unchanged during averaging of bonds,
77 // but may be updated when no averaging is performed.
78
79 int rene_flag = 0;
80 if (nevery * nrepeat != 1 && (nfreq % neighbor->every != 0 || neighbor->every < nevery * nrepeat)) {
81 int newneighborevery = nevery * nrepeat;
82 while (nfreq % newneighborevery != 0 && newneighborevery <= nfreq / 2)
83 newneighborevery++;
84
85 if (nfreq % newneighborevery != 0)
86 newneighborevery = nfreq;
87
88 neighbor->every = newneighborevery;
89 rene_flag = 1;
90 }
91
92 if (nevery * nrepeat != 1 && (neighbor->delay != 0 || neighbor->dist_check != 0)) {
93 neighbor->delay = 0;
94 neighbor->dist_check = 0;
95 rene_flag = 1;
96 }
97
98 if (me == 0 && rene_flag) {
99 error->warning(FLERR,"Resetting reneighboring criteria for fix reaxff/species");
100 }
101
102 tmparg = nullptr;
103 memory->create(tmparg,4,4,"reaxff/species:tmparg");
104 strcpy(tmparg[0],arg[3]);
105 strcpy(tmparg[1],arg[4]);
106 strcpy(tmparg[2],arg[5]);
107
108 if (me == 0) {
109 char *suffix = strrchr(arg[6],'.');
110 if (suffix && strcmp(suffix,".gz") == 0) {
111 #ifdef LAMMPS_GZIP
112 auto gzip = fmt::format("gzip -6 > {}",arg[6]);
113 #ifdef _WIN32
114 fp = _popen(gzip.c_str(),"wb");
115 #else
116 fp = popen(gzip.c_str(),"w");
117 #endif
118 #else
119 error->one(FLERR,"Cannot open gzipped file");
120 #endif
121 } else fp = fopen(arg[6],"w");
122
123 if (!fp)
124 error->one(FLERR,fmt::format("Cannot open fix reaxff/species file {}: "
125 "{}",arg[6],utils::getsyserror()));
126 }
127
128 x0 = nullptr;
129 clusterID = nullptr;
130
131 int ntmp = 1;
132 memory->create(x0,ntmp,"reaxff/species:x0");
133 memory->create(clusterID,ntmp,"reaxff/species:clusterID");
134 vector_atom = clusterID;
135
136 BOCut = nullptr;
137 Name = nullptr;
138 MolName = nullptr;
139 MolType = nullptr;
140 NMol = nullptr;
141 nd = nullptr;
142 molmap = nullptr;
143
144 nmax = 0;
145 setupflag = 0;
146
147 // set default bond order cutoff
148 int n, i, j, itype, jtype;
149 double bo_cut;
150 bg_cut = 0.30;
151 n = ntypes+1;
152 memory->create(BOCut,n,n,"reaxff/species:BOCut");
153 for (i = 1; i < n; i ++)
154 for (j = 1; j < n; j ++)
155 BOCut[i][j] = bg_cut;
156
157 // optional args
158 eletype = nullptr;
159 ele = filepos = nullptr;
160 eleflag = posflag = padflag = 0;
161
162 singlepos_opened = multipos_opened = 0;
163 multipos = 0;
164 posfreq = 0;
165
166 int iarg = 7;
167 while (iarg < narg) {
168
169 // set BO cutoff
170 if (strcmp(arg[iarg],"cutoff") == 0) {
171 if (iarg+4 > narg) error->all(FLERR,"Illegal fix reaxff/species command");
172 itype = atoi(arg[iarg+1]);
173 jtype = atoi(arg[iarg+2]);
174 bo_cut = atof(arg[iarg+3]);
175 if (itype > ntypes || jtype > ntypes)
176 error->all(FLERR,"Illegal fix reaxff/species command");
177 if (itype <= 0 || jtype <= 0)
178 error->all(FLERR,"Illegal fix reaxff/species command");
179 if (bo_cut > 1.0 || bo_cut < 0.0)
180 error->all(FLERR,"Illegal fix reaxff/species command");
181
182 BOCut[itype][jtype] = bo_cut;
183 BOCut[jtype][itype] = bo_cut;
184 iarg += 4;
185
186 // modify element type names
187 } else if (strcmp(arg[iarg],"element") == 0) {
188 if (iarg+ntypes+1 > narg) error->all(FLERR,"Illegal fix reaxff/species command");
189
190 eletype = (char**) malloc(ntypes*sizeof(char*));
191 int len;
192 for (int i = 0; i < ntypes; i ++) {
193 len = strlen(arg[iarg+1+i])+1;
194 eletype[i] = (char*) malloc(len*sizeof(char));
195 strcpy(eletype[i],arg[iarg+1+i]);
196 }
197 eleflag = 1;
198 iarg += ntypes + 1;
199
200 // position of molecules
201 } else if (strcmp(arg[iarg],"position") == 0) {
202 if (iarg+3 > narg) error->all(FLERR,"Illegal fix reaxff/species command");
203 posflag = 1;
204 posfreq = atoi(arg[iarg+1]);
205 if (posfreq < nfreq || (posfreq%nfreq != 0))
206 error->all(FLERR,"Illegal fix reaxff/species command");
207
208 filepos = new char[255];
209 strcpy(filepos,arg[iarg+2]);
210 if (strchr(filepos,'*')) {
211 multipos = 1;
212 } else {
213 if (me == 0) {
214 pos = fopen(filepos, "w");
215 if (pos == nullptr) error->one(FLERR,"Cannot open fix reaxff/species position file");
216 }
217 singlepos_opened = 1;
218 multipos = 0;
219 }
220 iarg += 3;
221 } else error->all(FLERR,"Illegal fix reaxff/species command");
222 }
223
224 if (!eleflag) {
225 memory->create(ele,ntypes+1,"reaxff/species:ele");
226 ele[0]='C';
227 if (ntypes > 1)
228 ele[1]='H';
229 if (ntypes > 2)
230 ele[2]='O';
231 if (ntypes > 3)
232 ele[3]='N';
233 }
234
235 vector_nmole = 0;
236 vector_nspec = 0;
237
238 }
239
240 /* ---------------------------------------------------------------------- */
241
~FixReaxFFSpecies()242 FixReaxFFSpecies::~FixReaxFFSpecies()
243 {
244 memory->destroy(ele);
245 memory->destroy(BOCut);
246 memory->destroy(clusterID);
247 memory->destroy(x0);
248
249 memory->destroy(nd);
250 memory->destroy(Name);
251 memory->destroy(NMol);
252 memory->destroy(MolType);
253 memory->destroy(MolName);
254 memory->destroy(tmparg);
255
256 if (filepos)
257 delete [] filepos;
258
259 if (me == 0) fclose(fp);
260 if (me == 0 && posflag && multipos_opened) fclose(pos);
261
262 modify->delete_compute("SPECATOM");
263 modify->delete_fix("SPECBOND");
264 }
265
266 /* ---------------------------------------------------------------------- */
267
setmask()268 int FixReaxFFSpecies::setmask()
269 {
270 int mask = 0;
271 mask |= POST_INTEGRATE;
272 return mask;
273 }
274
275 /* ---------------------------------------------------------------------- */
276
setup(int)277 void FixReaxFFSpecies::setup(int /*vflag*/)
278 {
279 ntotal = static_cast<int> (atom->natoms);
280 if (Name == nullptr)
281 memory->create(Name,ntypes,"reaxff/species:Name");
282
283 post_integrate();
284 }
285
286 /* ---------------------------------------------------------------------- */
287
init()288 void FixReaxFFSpecies::init()
289 {
290 if (atom->tag_enable == 0)
291 error->all(FLERR,"Cannot use fix reaxff/species unless atoms have IDs");
292
293 reaxff = (PairReaxFF *) force->pair_match("^reax..",0);
294 if (reaxff == nullptr) error->all(FLERR,"Cannot use fix reaxff/species without "
295 "pair_style reaxff, reaxff/kk, or reaxff/omp");
296
297 reaxff->fixspecies_flag = 1;
298
299 // reset next output timestep if not yet set or timestep has been reset
300 if (nvalid != update->ntimestep)
301 nvalid = update->ntimestep+nfreq;
302
303 // check if this fix has been called twice
304 int count = 0;
305 for (int i = 0; i < modify->nfix; i++)
306 if (strcmp(modify->fix[i]->style,"reaxff/species") == 0) count++;
307 if (count > 1 && comm->me == 0)
308 error->warning(FLERR,"More than one fix reaxff/species");
309
310 if (!setupflag) {
311 // create a compute to store properties
312 modify->add_compute("SPECATOM all SPEC/ATOM q x y z vx vy vz abo01 abo02 abo03 abo04 "
313 "abo05 abo06 abo07 abo08 abo09 abo10 abo11 abo12 abo13 abo14 "
314 "abo15 abo16 abo17 abo18 abo19 abo20 abo21 abo22 abo23 abo24");
315
316 // create a fix to point to fix_ave_atom for averaging stored properties
317 auto fixcmd = fmt::format("SPECBOND all ave/atom {} {} {}",tmparg[0],tmparg[1],tmparg[2]);
318 for (int i = 1; i < 32; ++i) fixcmd += " c_SPECATOM[" + std::to_string(i) + "]";
319 f_SPECBOND = (FixAveAtom *) modify->add_fix(fixcmd);
320 setupflag = 1;
321 }
322 }
323
324 /* ---------------------------------------------------------------------- */
325
init_list(int,NeighList * ptr)326 void FixReaxFFSpecies::init_list(int /*id*/, NeighList *ptr)
327 {
328 list = ptr;
329 }
330
331 /* ---------------------------------------------------------------------- */
332
post_integrate()333 void FixReaxFFSpecies::post_integrate()
334 {
335 Output_ReaxFF_Bonds(update->ntimestep,fp);
336 if (me == 0) fflush(fp);
337 }
338
339 /* ---------------------------------------------------------------------- */
340
Output_ReaxFF_Bonds(bigint ntimestep,FILE *)341 void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/)
342
343 {
344 int Nmole, Nspec;
345
346 // point to fix_ave_atom
347 f_SPECBOND->end_of_step();
348
349 if (ntimestep != nvalid) return;
350
351 nlocal = atom->nlocal;
352
353 if (atom->nmax > nmax) {
354 nmax = atom->nmax;
355 memory->destroy(x0);
356 memory->destroy(clusterID);
357 memory->create(x0,nmax,"reaxff/species:x0");
358 memory->create(clusterID,nmax,"reaxff/species:clusterID");
359 vector_atom = clusterID;
360 }
361
362 for (int i = 0; i < nmax; i++) {
363 x0[i].x = x0[i].y = x0[i].z = 0.0;
364 }
365
366 Nmole = Nspec = 0;
367
368 FindMolecule();
369
370 SortMolecule (Nmole);
371
372 FindSpecies(Nmole, Nspec);
373
374 vector_nmole = Nmole;
375 vector_nspec = Nspec;
376
377 if (me == 0 && ntimestep >= 0)
378 WriteFormulas (Nmole, Nspec);
379
380 if (posflag && ((ntimestep)%posfreq==0)) {
381 WritePos(Nmole, Nspec);
382 if (me == 0) fflush(pos);
383 }
384
385 nvalid += nfreq;
386 }
387
388 /* ---------------------------------------------------------------------- */
389
chAnchor(AtomCoord in1,AtomCoord in2)390 AtomCoord FixReaxFFSpecies::chAnchor(AtomCoord in1, AtomCoord in2)
391 {
392 if (in1.x < in2.x)
393 return in1;
394 else if (in1.x == in2.x) {
395 if (in1.y < in2.y)
396 return in1;
397 else if (in1.y == in2.y) {
398 if (in1.z < in2.z)
399 return in1;
400 }
401 }
402 return in2;
403 }
404
405 /* ---------------------------------------------------------------------- */
406
FindMolecule()407 void FixReaxFFSpecies::FindMolecule ()
408 {
409 int i,j,ii,jj,inum,itype,jtype,loop,looptot;
410 int change,done,anychange;
411 int *mask = atom->mask;
412 int *ilist;
413 double bo_tmp,bo_cut;
414 double **spec_atom = f_SPECBOND->array_atom;
415
416 inum = reaxff->list->inum;
417 ilist = reaxff->list->ilist;
418
419 for (ii = 0; ii < inum; ii++) {
420 i = ilist[ii];
421 if (mask[i] & groupbit) {
422 clusterID[i] = atom->tag[i];
423 x0[i].x = spec_atom[i][1];
424 x0[i].y = spec_atom[i][2];
425 x0[i].z = spec_atom[i][3];
426 }
427 else clusterID[i] = 0.0;
428 }
429
430 loop = 0;
431 while (1) {
432 comm->forward_comm_fix(this);
433 loop ++;
434
435 change = 0;
436 while (1) {
437 done = 1;
438
439 for (ii = 0; ii < inum; ii++) {
440 i = ilist[ii];
441 if (!(mask[i] & groupbit)) continue;
442
443 itype = atom->type[i];
444
445 for (jj = 0; jj < MAXSPECBOND; jj++) {
446 j = reaxff->tmpid[i][jj];
447
448 if ((j == 0) || (j < i)) continue;
449 if (!(mask[j] & groupbit)) continue;
450
451 if (clusterID[i] == clusterID[j]
452 && x0[i].x == x0[j].x
453 && x0[i].y == x0[j].y
454 && x0[i].z == x0[j].z) continue;
455
456 jtype = atom->type[j];
457 bo_cut = BOCut[itype][jtype];
458 bo_tmp = spec_atom[i][jj+7];
459
460 if (bo_tmp > bo_cut) {
461 clusterID[i] = clusterID[j] = MIN(clusterID[i], clusterID[j]);
462 x0[i] = x0[j] = chAnchor(x0[i], x0[j]);
463 done = 0;
464 }
465 }
466 }
467 if (!done) change = 1;
468 if (done) break;
469 }
470 MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world);
471 if (!anychange) break;
472
473 MPI_Allreduce(&loop,&looptot,1,MPI_INT,MPI_SUM,world);
474 if (looptot >= 400*nprocs) break;
475
476 }
477 }
478
479 /* ---------------------------------------------------------------------- */
480
SortMolecule(int & Nmole)481 void FixReaxFFSpecies::SortMolecule(int &Nmole)
482 {
483 memory->destroy(molmap);
484 molmap = nullptr;
485
486 int n, idlo, idhi;
487 int *mask =atom->mask;
488 int lo = ntotal;
489 int hi = -ntotal;
490 int flag = 0;
491 for (n = 0; n < nlocal; n++) {
492 if (!(mask[n] & groupbit)) continue;
493 if (clusterID[n] == 0.0) flag = 1;
494 lo = MIN(lo,nint(clusterID[n]));
495 hi = MAX(hi,nint(clusterID[n]));
496 }
497 int flagall;
498 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
499 if (flagall && me == 0)
500 error->warning(FLERR,"Atom with cluster ID = 0 included in "
501 "fix reaxff/species group");
502 MPI_Allreduce(&lo,&idlo,1,MPI_INT,MPI_MIN,world);
503 MPI_Allreduce(&hi,&idhi,1,MPI_INT,MPI_MAX,world);
504 if (idlo == ntotal)
505 if (me == 0)
506 error->warning(FLERR,"Atom with cluster ID = maxmol "
507 "included in fix reaxff/species group");
508
509 int nlen = idhi - idlo + 1;
510 memory->create(molmap,nlen,"reaxff/species:molmap");
511 for (n = 0; n < nlen; n++) molmap[n] = 0;
512
513 for (n = 0; n < nlocal; n++) {
514 if (!(mask[n] & groupbit)) continue;
515 molmap[nint(clusterID[n])-idlo] = 1;
516 }
517
518 int *molmapall;
519 memory->create(molmapall,nlen,"reaxff/species:molmapall");
520 MPI_Allreduce(molmap,molmapall,nlen,MPI_INT,MPI_MAX,world);
521
522 Nmole = 0;
523 for (n = 0; n < nlen; n++) {
524 if (molmapall[n]) molmap[n] = Nmole++;
525 else molmap[n] = -1;
526 }
527 memory->destroy(molmapall);
528
529 flag = 0;
530 for (n = 0; n < nlocal; n++) {
531 if (mask[n] & groupbit) continue;
532 if (nint(clusterID[n]) < idlo || nint(clusterID[n]) > idhi) continue;
533 if (molmap[nint(clusterID[n])-idlo] >= 0) flag = 1;
534 }
535
536 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
537 if (flagall && comm->me == 0)
538 error->warning(FLERR,"One or more cluster has atoms not in group");
539
540 for (n = 0; n < nlocal; n++) {
541 if (!(mask[n] & groupbit)) continue;
542 clusterID[n] = molmap[nint(clusterID[n])-idlo] + 1;
543 }
544
545 memory->destroy(molmap);
546 molmap = nullptr;
547
548 }
549
550 /* ---------------------------------------------------------------------- */
551
FindSpecies(int Nmole,int & Nspec)552 void FixReaxFFSpecies::FindSpecies(int Nmole, int &Nspec)
553 {
554 int k, l, m, n, itype, cid;
555 int flag_identity, flag_mol, flag_spec;
556 int flag_tmp;
557 int *mask =atom->mask;
558 int *Nameall, *NMolall;
559
560 memory->destroy(MolName);
561 MolName = nullptr;
562 memory->create(MolName,Nmole*(ntypes+1),"reaxff/species:MolName");
563
564 memory->destroy(NMol);
565 NMol = nullptr;
566 memory->create(NMol,Nmole,"reaxff/species:NMol");
567 for (m = 0; m < Nmole; m ++)
568 NMol[m] = 1;
569
570 memory->create(Nameall,ntypes,"reaxff/species:Nameall");
571 memory->create(NMolall,Nmole,"reaxff/species:NMolall");
572
573 for (m = 1, Nspec = 0; m <= Nmole; m ++) {
574 for (n = 0; n < ntypes; n ++) Name[n] = 0;
575 for (n = 0, flag_mol = 0; n < nlocal; n ++) {
576 if (!(mask[n] & groupbit)) continue;
577 cid = nint(clusterID[n]);
578 if (cid == m) {
579 itype = atom->type[n]-1;
580 Name[itype] ++;
581 flag_mol = 1;
582 }
583 }
584 MPI_Allreduce(&flag_mol,&flag_tmp,1,MPI_INT,MPI_MAX,world);
585 flag_mol = flag_tmp;
586
587 MPI_Allreduce(Name,Nameall,ntypes,MPI_INT,MPI_SUM,world);
588 for (n = 0; n < ntypes; n++) Name[n] = Nameall[n];
589
590 if (flag_mol == 1) {
591 flag_identity = 1;
592 for (k = 0; k < Nspec; k ++) {
593 flag_spec=0;
594 for (l = 0; l < ntypes; l ++)
595 if (MolName[ntypes*k+l] != Name[l]) flag_spec = 1;
596 if (flag_spec == 0) NMol[k] ++;
597 flag_identity *= flag_spec;
598 }
599 if (Nspec == 0 || flag_identity == 1) {
600 for (l = 0; l < ntypes; l ++)
601 MolName[ntypes*Nspec+l] = Name[l];
602 Nspec ++;
603 }
604 }
605 }
606 memory->destroy(NMolall);
607 memory->destroy(Nameall);
608
609 memory->destroy(nd);
610 nd = nullptr;
611 memory->create(nd,Nspec,"reaxff/species:nd");
612
613 memory->destroy(MolType);
614 MolType = nullptr;
615 memory->create(MolType,Nspec*(ntypes+2),"reaxff/species:MolType");
616 }
617
618 /* ---------------------------------------------------------------------- */
619
CheckExistence(int id,int ntypes)620 int FixReaxFFSpecies::CheckExistence(int id, int ntypes)
621 {
622 int i, j, molid, flag;
623
624 for (i = 0; i < Nmoltype; i ++) {
625 flag = 0;
626 for (j = 0; j < ntypes; j ++) {
627 molid = MolType[ntypes * i + j];
628 if (molid != MolName[ntypes * id + j]) flag = 1;
629 }
630 if (flag == 0) return i;
631 }
632 for (i = 0; i < ntypes; i ++)
633 MolType[ntypes * Nmoltype + i] = MolName[ntypes *id + i];
634
635 Nmoltype ++;
636 return Nmoltype - 1;
637 }
638
639 /* ---------------------------------------------------------------------- */
640
WriteFormulas(int Nmole,int Nspec)641 void FixReaxFFSpecies::WriteFormulas(int Nmole, int Nspec)
642 {
643 int i, j, itemp;
644 bigint ntimestep = update->ntimestep;
645
646 fprintf(fp,"# Timestep No_Moles No_Specs ");
647
648 Nmoltype = 0;
649
650 for (i = 0; i < Nspec; i ++)
651 nd[i] = CheckExistence(i, ntypes);
652
653 for (i = 0; i < Nmoltype; i ++) {
654 for (j = 0;j < ntypes; j ++) {
655 itemp = MolType[ntypes * i + j];
656 if (itemp != 0) {
657 if (eletype) fprintf(fp,"%s",eletype[j]);
658 else fprintf(fp,"%c",ele[j]);
659 if (itemp != 1) fprintf(fp,"%d",itemp);
660 }
661 }
662 fprintf(fp,"\t");
663 }
664 fprintf(fp,"\n");
665
666 fprintf(fp,BIGINT_FORMAT,ntimestep);
667 fprintf(fp,"%11d%11d\t",Nmole,Nspec);
668
669 for (i = 0; i < Nmoltype; i ++)
670 fprintf(fp," %d\t",NMol[i]);
671 fprintf(fp,"\n");
672
673 }
674
675 /* ---------------------------------------------------------------------- */
676
OpenPos()677 void FixReaxFFSpecies::OpenPos()
678 {
679 char *filecurrent;
680 bigint ntimestep = update->ntimestep;
681
682 filecurrent = (char*) malloc((strlen(filepos)+16)*sizeof(char));
683 char *ptr = strchr(filepos,'*');
684 *ptr = '\0';
685 if (padflag == 0)
686 sprintf(filecurrent,"%s" BIGINT_FORMAT "%s",
687 filepos,ntimestep,ptr+1);
688 else {
689 char bif[8],pad[16];
690 strcpy(bif,BIGINT_FORMAT);
691 sprintf(pad,"%%s%%0%d%s%%s",padflag,&bif[1]);
692 sprintf(filecurrent,pad,filepos,ntimestep,ptr+1);
693 }
694 *ptr = '*';
695
696 if (me == 0) {
697 pos = fopen(filecurrent, "w");
698 if (pos == nullptr) error->one(FLERR,"Cannot open fix reaxff/species position file");
699 } else pos = nullptr;
700 multipos_opened = 1;
701
702 free(filecurrent);
703 }
704
705 /* ---------------------------------------------------------------------- */
706
WritePos(int Nmole,int Nspec)707 void FixReaxFFSpecies::WritePos(int Nmole, int Nspec)
708 {
709 int i, itype, cid;
710 int count, count_tmp, m, n, k;
711 int *Nameall;
712 int *mask =atom->mask;
713 double avq, avq_tmp, avx[3], avx_tmp, box[3], halfbox[3];
714 double **spec_atom = f_SPECBOND->array_atom;
715
716 if (multipos) OpenPos();
717
718 box[0] = domain->boxhi[0] - domain->boxlo[0];
719 box[1] = domain->boxhi[1] - domain->boxlo[1];
720 box[2] = domain->boxhi[2] - domain->boxlo[2];
721
722 for (int j = 0; j < 3; j++)
723 halfbox[j] = box[j] / 2;
724
725 if (me == 0) {
726 fprintf(pos,"Timestep " BIGINT_FORMAT " NMole %d NSpec %d xlo %f "
727 "xhi %f ylo %f yhi %f zlo %f zhi %f\n",
728 update->ntimestep,Nmole, Nspec,
729 domain->boxlo[0],domain->boxhi[0],
730 domain->boxlo[1],domain->boxhi[1],
731 domain->boxlo[2],domain->boxhi[2]);
732
733 fprintf(pos,"ID\tAtom_Count\tType\tAve_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n");
734 }
735
736 Nameall = nullptr;
737 memory->create(Nameall,ntypes,"reaxff/species:Nameall");
738
739 for (m = 1; m <= Nmole; m ++) {
740
741 count = 0;
742 avq = 0.0;
743 for (n = 0; n < 3; n++)
744 avx[n] = 0.0;
745 for (n = 0; n < ntypes; n ++)
746 Name[n] = 0;
747
748 for (i = 0; i < nlocal; i ++) {
749 if (!(mask[i] & groupbit)) continue;
750 cid = nint(clusterID[i]);
751 if (cid == m) {
752 itype = atom->type[i]-1;
753 Name[itype] ++;
754 count ++;
755 avq += spec_atom[i][0];
756 if ((x0[i].x - spec_atom[i][1]) > halfbox[0])
757 spec_atom[i][1] += box[0];
758 if ((spec_atom[i][1] - x0[i].x) > halfbox[0])
759 spec_atom[i][1] -= box[0];
760 if ((x0[i].y - spec_atom[i][2]) > halfbox[1])
761 spec_atom[i][2] += box[1];
762 if ((spec_atom[i][2] - x0[i].y) > halfbox[1])
763 spec_atom[i][2] -= box[1];
764 if ((x0[i].z - spec_atom[i][3]) > halfbox[2])
765 spec_atom[i][3] += box[2];
766 if ((spec_atom[i][3] - x0[i].z) > halfbox[2])
767 spec_atom[i][3] -= box[2];
768 for (n = 0; n < 3; n++)
769 avx[n] += spec_atom[i][n+1];
770 }
771 }
772
773 avq_tmp = 0.0;
774 MPI_Allreduce(&avq,&avq_tmp,1,MPI_DOUBLE,MPI_SUM,world);
775 avq = avq_tmp;
776
777 for (n = 0; n < 3; n++) {
778 avx_tmp = 0.0;
779 MPI_Reduce(&avx[n],&avx_tmp,1,MPI_DOUBLE,MPI_SUM,0,world);
780 avx[n] = avx_tmp;
781 }
782
783 MPI_Reduce(&count,&count_tmp,1,MPI_INT,MPI_SUM,0,world);
784 count = count_tmp;
785
786 MPI_Reduce(Name,Nameall,ntypes,MPI_INT,MPI_SUM,0,world);
787 for (n = 0; n < ntypes; n++) Name[n] = Nameall[n];
788
789 if (me == 0) {
790 fprintf(pos,"%d\t%d\t",m,count);
791 for (n = 0; n < ntypes; n++) {
792 if (Name[n] != 0) {
793 if (eletype) fprintf(pos,"%s",eletype[n]);
794 else fprintf(pos,"%c",ele[n]);
795 if (Name[n] != 1) fprintf(pos,"%d",Name[n]);
796 }
797 }
798 if (count > 0) {
799 avq /= count;
800 for (k = 0; k < 3; k++) {
801 avx[k] /= count;
802 if (avx[k] >= domain->boxhi[k])
803 avx[k] -= box[k];
804 if (avx[k] < domain->boxlo[k])
805 avx[k] += box[k];
806
807 avx[k] -= domain->boxlo[k];
808 avx[k] /= box[k];
809 }
810 fprintf(pos,"\t%.8f \t%.8f \t%.8f \t%.8f",
811 avq,avx[0],avx[1],avx[2]);
812 }
813 fprintf(pos,"\n");
814 }
815 }
816 if (me == 0 && !multipos) fprintf(pos,"#\n");
817 memory->destroy(Nameall);
818 }
819
820 /* ---------------------------------------------------------------------- */
821
compute_vector(int n)822 double FixReaxFFSpecies::compute_vector(int n)
823 {
824 if (n == 0)
825 return vector_nmole;
826 if (n == 1)
827 return vector_nspec;
828 return 0.0;
829
830 }
831
832 /* ---------------------------------------------------------------------- */
833
nint(const double & r)834 int FixReaxFFSpecies::nint(const double &r)
835 {
836 int i = 0;
837 if (r>0.0) i = static_cast<int>(r+0.5);
838 else if (r<0.0) i = static_cast<int>(r-0.5);
839 return i;
840 }
841
842 /* ---------------------------------------------------------------------- */
843
pack_forward_comm(int n,int * list,double * buf,int,int *)844 int FixReaxFFSpecies::pack_forward_comm(int n, int *list, double *buf,
845 int /*pbc_flag*/, int * /*pbc*/)
846 {
847 int i,j,m;
848
849 m = 0;
850 for (i = 0; i < n; i++) {
851 j = list[i];
852 buf[m] = clusterID[j];
853 buf[m+1] = x0[j].x;
854 buf[m+2] = x0[j].y;
855 buf[m+3] = x0[j].z;
856 m += 4;
857 }
858 return m;
859 }
860
861 /* ---------------------------------------------------------------------- */
862
unpack_forward_comm(int n,int first,double * buf)863 void FixReaxFFSpecies::unpack_forward_comm(int n, int first, double *buf)
864 {
865 int i,m,last;
866
867 m = 0;
868 last = first + n;
869 for (i = first; i < last; i++) {
870 clusterID[i] = buf[m];
871 x0[i].x = buf[m+1];
872 x0[i].y = buf[m+2];
873 x0[i].z = buf[m+3];
874 m += 4;
875 }
876 }
877
878 /* ---------------------------------------------------------------------- */
879
memory_usage()880 double FixReaxFFSpecies::memory_usage()
881 {
882 double bytes;
883
884 bytes = 4*nmax*sizeof(double); // clusterID + x0
885
886 return bytes;
887 }
888
889 /* ---------------------------------------------------------------------- */
890