1 // Header file for this class
2 #include "LammpsInterface.h"
3 // LAMMPS includes
4 #include "lammps.h"
5 #include "atom.h" // x, v, f
6 #include "atom_vec.h" //for insertion
7 #include "domain.h" // for basing locations on regions
8 #include "region.h" // region bounding box and style
9 #include "force.h" // boltzman constant
10 #include "group.h" // atom masks
11 #include "memory.h" // grow atom information
12 #include "compute.h" // computes
13 #include "compute_pe_atom.h" // computes potential energy per atom
14 #include "compute_stress_atom.h" // computes stress per atom
15 #include "compute_centro_atom.h" // computes centrosymmetry per atom
16 #include "compute_cna_atom.h" // computes common-neighbor-analysis per atom
17 #include "compute_coord_atom.h" // computes coordination number per atom
18 #include "compute_ke_atom.h" // computes kinetic energy per atom
19 #include "modify.h" //
20 #include "neighbor.h" // neighbors
21 #include "neigh_list.h" // neighbor list
22 #include "update.h" // timestepping information
23 #include "pair.h" // pair potentials
24 #include "MANYBODY/pair_eam.h" // pair potentials
25 #include "lattice.h" // lattice parameters
26 #include "bond.h" // bond potentials
27 #include "comm.h" //
28 #include "fix.h"
29 #include "utils.h"
30
31 // ATC includes
32 #include "ATC_Error.h"
33 #include "MatrixLibrary.h"
34 #include "Utility.h"
35 using ATC_Utility::to_string;
36
37 // Other include files
38 #include "mpi.h"
39 #include <cstring>
40 #include <map>
41 #include <typeinfo>
42
43 using std::max;
44 using std::stringstream;
45 using std::copy;
46 using std::map;
47 using std::pair;
48 using std::string;
49 using std::set;
50 using LAMMPS_NS::bigint;
51 using LAMMPS_NS::utils::read_lines_from_file;
52
53 namespace ATC {
54
55 const static double PI = 3.141592653589793238;
56
57 const static int seed_ = 3141592;
58
59 const static int MAX_GROUP_BIT = 2147483647; //4294967295; // pow(2,31)-1;
60
norm(double * v)61 double norm(double * v) {return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); }
62
63 LammpsInterface * LammpsInterface::myInstance_ = nullptr;
64
65 // -----------------------------------------------------------------
66 // instance()
67 // -----------------------------------------------------------------
instance()68 LammpsInterface * LammpsInterface::instance()
69 {
70 if (myInstance_ == nullptr) {
71 myInstance_ = new LammpsInterface();
72 }
73 return myInstance_;
74 }
75
76 // -----------------------------------------------------------------
77 // Destroy()
78 // -----------------------------------------------------------------
Destroy()79 void LammpsInterface::Destroy()
80 {
81 if (myInstance_) delete myInstance_;
82 myInstance_ = nullptr;
83 }
84
85
86 // -----------------------------------------------------------------
87 // constructor
88 // -----------------------------------------------------------------
LammpsInterface()89 LammpsInterface::LammpsInterface()
90 : lammps_(nullptr),
91 fixPointer_(nullptr),
92 commRank_(0),
93 atomPE_(nullptr),
94 refBoxIsSet_(false),
95 random_(nullptr),
96 globalrandom_(nullptr)
97 {
98 }
99
100 // -----------------------------------------------------------------
101 // general interface methods
102 // -----------------------------------------------------------------
103
world() const104 MPI_Comm LammpsInterface::world() const { return lammps_->world; }
set_fix_pointer(LAMMPS_NS::Fix * thisFix)105 void LammpsInterface::set_fix_pointer(LAMMPS_NS::Fix * thisFix) { fixPointer_ = thisFix; }
forward_comm_fix() const106 void LammpsInterface::forward_comm_fix() const { lammps_->comm->forward_comm_fix(fixPointer_); }
comm_borders() const107 void LammpsInterface::comm_borders() const { lammps_->comm->borders(); }
108
109 #ifndef ISOLATE_FE
sparse_allsum(SparseMatrix<double> & toShare) const110 void LammpsInterface::sparse_allsum(SparseMatrix<double> &toShare) const
111 {
112 toShare.compress();
113
114 // initialize MPI information
115 int nProcs;
116 int myRank;
117 MPI_Comm_size(MPI_COMM_WORLD, &nProcs);
118 MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
119
120 int error;
121
122 // get numbers of rows, columns, rowsCRS, and
123 // sizes (number of nonzero elements in matrix)
124 SparseMatInfo *recInfo = new SparseMatInfo[nProcs];
125 SparseMatInfo myInfo;
126 myInfo.rows = toShare.nRows();
127 myInfo.cols = toShare.nCols();
128 myInfo.rowsCRS = toShare.nRowsCRS();
129 myInfo.size = toShare.size();
130
131 error = MPI_Allgather(&myInfo, 4, MPI_INT,
132 recInfo, 4, MPI_INT, lammps_->world);
133 if (error != MPI_SUCCESS) throw ATC_Error("error in sparse_allsum_numrows "+to_string(error));
134
135 // adjust row sendcounts because recRowsCRS is off by one
136 int *rowCounts = new int[nProcs];
137 int *sizeCounts = new int[nProcs];
138 // set up total size of receive buffers for Allgatherv calls
139 int totalRowsCRS = 0;
140 int totalSize = 0;
141 // set up array of displacements for Allgatherv calls
142 int *rowOffsets = new int[nProcs];
143 rowOffsets[0] = 0;
144 int *sizeOffsets = new int[nProcs];
145 sizeOffsets[0] = 0;
146 for (int i = 0; i < nProcs; i++) {
147 // find the total number of entries to share in the mpi calls below
148 rowCounts[i] = recInfo[i].rowsCRS + 1;
149 sizeCounts[i] = recInfo[i].size;
150 totalRowsCRS += rowCounts[i];
151 totalSize += recInfo[i].size;
152 // these already have their 0th slot filled in
153 if (i == 0) continue;
154 rowOffsets[i] = rowOffsets[i-1] + rowCounts[i-1];
155 sizeOffsets[i] = sizeOffsets[i-1] + sizeCounts[i-1];
156 }
157
158 // get actual rows
159 INDEX *rec_ia = new INDEX[totalRowsCRS];
160 if (toShare.size() == 0) {
161 double dummy;
162 error = MPI_Allgatherv(&dummy, 0, MPI_INT,
163 rec_ia, rowCounts, rowOffsets, MPI_INT, lammps_->world);
164 }
165 else
166 error = MPI_Allgatherv(toShare.rows(), rowCounts[myRank], MPI_INT,
167 rec_ia, rowCounts, rowOffsets, MPI_INT, lammps_->world);
168 if (error != MPI_SUCCESS)
169 throw ATC_Error("error in sparse_allsum_rowarray "+to_string(error));
170
171 // get actual cols
172 INDEX *rec_ja = new INDEX[totalSize];
173 error = MPI_Allgatherv(toShare.cols(), sizeCounts[myRank], MPI_INT,
174 rec_ja, sizeCounts, sizeOffsets, MPI_INT, lammps_->world);
175 if (error != MPI_SUCCESS)
176 throw ATC_Error("error in sparse_allsum_colarray "+to_string(error));
177
178 // get the array of values
179 double *rec_vals = new double[totalSize];
180 error = MPI_Allgatherv(toShare.ptr(), sizeCounts[myRank], MPI_DOUBLE,
181 rec_vals, sizeCounts, sizeOffsets, MPI_DOUBLE, lammps_->world);
182 if (error != MPI_SUCCESS)
183 throw ATC_Error("error in sparse_allsum_valarray "+to_string(error));
184
185 INDEX *rec_ia_proc;
186 INDEX *rec_ja_proc;
187 double *rec_vals_proc;
188 for (int i = 0; i < nProcs; i++) {
189 if (myRank != i) {
190 // deallocated when tempMat is deleted since it wraps them
191 rec_ia_proc = new INDEX[rowCounts[i]];
192 rec_ja_proc = new INDEX[sizeCounts[i]];
193 rec_vals_proc = new double[sizeCounts[i]];
194
195 // copy the data passed with MPI into the new spots
196 copy(rec_ia + rowOffsets[i],
197 rec_ia + rowOffsets[i] + rowCounts[i],
198 rec_ia_proc);
199 copy(rec_ja + sizeOffsets[i],
200 rec_ja + sizeOffsets[i] + sizeCounts[i],
201 rec_ja_proc);
202 copy(rec_vals + sizeOffsets[i],
203 rec_vals + sizeOffsets[i] + sizeCounts[i],
204 rec_vals_proc);
205
206 // Does anyone know why we have to declare tempMat here (as well as set it equal to
207 // something) to avoid segfaults? there are still segfaults, but they happen at a much
208 // later stage of the game now (and for less benchmarks overall).
209 SparseMatrix<double> tempMat =
210 SparseMatrix<double>(rec_ia_proc, rec_ja_proc, rec_vals_proc,
211 recInfo[i].size, recInfo[i].rows,
212 recInfo[i].cols, recInfo[i].rowsCRS);
213 toShare += tempMat;
214 }
215 }
216 delete[] rowCounts;
217 delete[] sizeCounts;
218 delete[] rowOffsets;
219 delete[] sizeOffsets;
220
221 delete[] recInfo;
222 delete[] rec_ia;
223 delete[] rec_ja;
224 delete[] rec_vals;
225 }
226 #endif
227
read_file(std::string filename) const228 std::string LammpsInterface::read_file(std::string filename) const
229 {
230 FILE *fp = nullptr;
231 if (! comm_rank()) {
232 fp = fopen(filename.c_str(),"r");
233 if (!fp) throw ATC_Error("can't open file: "+filename);
234 }
235 const int MAXLINE = 256;
236 const int CHUNK = 1024;
237 char * buffer = new char[CHUNK*MAXLINE];
238 std::stringstream s;
239 bool eof = false;
240 while ( ! eof) {
241 eof = read_lines_from_file(fp,1,MAXLINE,buffer,comm_rank(),lammps_->world);
242 s << buffer;
243 }
244 fclose(fp);
245 delete [] buffer;
246 return s.str(); // can't copy the stream itself
247 }
248
249
250 // -----------------------------------------------------------------
251 // atom interface methods
252 // -----------------------------------------------------------------
fix_id() const253 string LammpsInterface::fix_id() const { return string(fixPointer_->id); }
254
nlocal() const255 int LammpsInterface::nlocal() const { return lammps_->atom->nlocal; }
256
nghost() const257 int LammpsInterface::nghost() const { return lammps_->atom->nghost; }
258
atoms_sorted() const259 bool LammpsInterface::atoms_sorted() const
260 {
261 int sortfreq = lammps_->atom->sortfreq;
262 if (sortfreq > 0) { return true; }
263 else { return false;
264 }
265 }
266
267
natoms() const268 bigint LammpsInterface::natoms() const { return lammps_->atom->natoms; }
269
nmax() const270 int LammpsInterface::nmax() const { return lammps_->atom->nmax; }
271
ntypes() const272 int LammpsInterface::ntypes() const { return lammps_->atom->ntypes; }
273
xatom() const274 double ** LammpsInterface::xatom() const { return lammps_->atom->x; }
275
type_to_charge(int atype) const276 int LammpsInterface::type_to_charge(int atype) const {
277 double *q = lammps_->atom->q;
278 if (! q) return 0;
279 int nlocal = lammps_->atom->nlocal;
280 int *type = lammps_->atom->type;
281 double aq = 0.0;
282 for (int i = 0; i < nlocal; i++) {
283 if (type[i] == atype) {
284 aq = q[i];
285 break;
286 }
287 }
288 double pcharge;
289 MPI_Allreduce(&aq,&pcharge,1,MPI_DOUBLE,MPI_MAX,world());
290 double ncharge;
291 MPI_Allreduce(&aq,&ncharge,1,MPI_DOUBLE,MPI_MIN,world());
292 double charge = (pcharge == 0.0) ? ncharge : pcharge;
293 return charge;
294 }
295
296 //const double ** LammpsInterface::xatom() const { return (const double**)(lammps_->atom->x); }
297
vatom() const298 double ** LammpsInterface::vatom() const { return lammps_->atom->v; }
299
fatom() const300 double ** LammpsInterface::fatom() const { return lammps_->atom->f; }
301
atom_mask() const302 const int * LammpsInterface::atom_mask() const { return (const int*)lammps_->atom->mask; }
303
atom_mask()304 int * LammpsInterface::atom_mask() { return lammps_->atom->mask; }
305
atom_type() const306 int * LammpsInterface::atom_type() const { return lammps_->atom->type; }
307
atom_tag() const308 int * LammpsInterface::atom_tag() const { return lammps_->atom->tag; }
309
atom_to_molecule() const310 int * LammpsInterface::atom_to_molecule() const { return lammps_->atom->molecule; }
311
312
num_bond() const313 int * LammpsInterface::num_bond() const { return lammps_->atom->num_bond; }
314
bond_atom() const315 int ** LammpsInterface::bond_atom() const { return lammps_->atom->bond_atom; }
316
image() const317 int * LammpsInterface::image() const { return lammps_->atom->image; }
318
bond_per_atom() const319 int LammpsInterface::bond_per_atom() const { return lammps_->atom->bond_per_atom; }
320
newton_bond() const321 int LammpsInterface::newton_bond() const { return lammps_->force->newton_bond; }
322
local_to_global_map(int global) const323 int LammpsInterface::local_to_global_map(int global) const { return lammps_->atom->map(global); }
324
atom_mass() const325 double * LammpsInterface::atom_mass() const { return lammps_->atom->mass; }
326
atom_mass(int iType) const327 double LammpsInterface::atom_mass(int iType) const { return lammps_->atom->mass[iType]; }
328
atom_rmass() const329 double * LammpsInterface::atom_rmass() const { return lammps_->atom->rmass; }
330
atom_charge() const331 double * LammpsInterface::atom_charge() const { return lammps_->atom->q; }
332
atom_scalar(FundamentalAtomQuantity quantityType) const333 double * LammpsInterface::atom_scalar(FundamentalAtomQuantity quantityType) const
334 {
335 if (quantityType==ATOM_MASS) {
336 if (atom_mass())
337 throw ATC_Error("Atom mass array requested but not defined");
338 return atom_rmass();
339 }
340 else if (quantityType==ATOM_CHARGE) {
341 double * atomCharge = atom_charge();
342 if (!atomCharge)
343 throw ATC_Error("Atom charge array requested but not defined");
344 return atomCharge;
345 }
346 else
347 throw ATC_Error("BAD type requested in atom_scalar");
348 return nullptr;
349 }
350
atom_vector(FundamentalAtomQuantity quantityType) const351 double ** LammpsInterface::atom_vector(FundamentalAtomQuantity quantityType) const
352 {
353 if (quantityType==ATOM_POSITION)
354 return xatom();
355 else if (quantityType==ATOM_VELOCITY)
356 return vatom();
357 else if (quantityType==ATOM_FORCE)
358 return fatom();
359 else
360 throw ATC_Error("BAD type requested in atom_vector");
361 return nullptr;
362 }
363
atom_quantity_ndof(FundamentalAtomQuantity quantityType) const364 int LammpsInterface::atom_quantity_ndof(FundamentalAtomQuantity quantityType) const
365 {
366 if (quantityType==ATOM_MASS || quantityType==ATOM_CHARGE)
367 return 1;
368 else if (quantityType==ATOM_POSITION || quantityType==ATOM_VELOCITY || quantityType==ATOM_FORCE)
369 return dimension();
370 else
371 throw ATC_Error("BAD type requested in atom_quantity_ndof");
372 }
373
atom_quantity_conversion(FundamentalAtomQuantity quantityType) const374 double LammpsInterface::atom_quantity_conversion(FundamentalAtomQuantity quantityType) const
375 {
376 if (quantityType==ATOM_MASS || quantityType==ATOM_CHARGE || quantityType==ATOM_POSITION || quantityType==ATOM_VELOCITY)
377 return 1;
378 else if ( quantityType==ATOM_FORCE)
379 return ftm2v();
380 else
381 throw ATC_Error("BAD type requested in atom_quantity_conversion");
382 }
383
384 // -----------------------------------------------------------------
385 // domain interface methods
386 // -----------------------------------------------------------------
387
dimension() const388 int LammpsInterface::dimension() const { return lammps_->domain->dimension; }
389
nregion() const390 int LammpsInterface::nregion() const { return lammps_->domain->nregion; }
391
box_bounds(double & boxxlo,double & boxxhi,double & boxylo,double & boxyhi,double & boxzlo,double & boxzhi) const392 void LammpsInterface::box_bounds(double & boxxlo, double & boxxhi,
393 double & boxylo, double & boxyhi,
394 double & boxzlo, double &boxzhi) const
395 {
396 if (lammps_->domain->triclinic == 0) {
397 boxxlo = lammps_->domain->boxlo[0];
398 boxxhi = lammps_->domain->boxhi[0];
399 boxylo = lammps_->domain->boxlo[1];
400 boxyhi = lammps_->domain->boxhi[1];
401 boxzlo = lammps_->domain->boxlo[2];
402 boxzhi = lammps_->domain->boxhi[2];
403 }
404 else {
405 boxxlo = lammps_->domain->boxlo_bound[0];
406 boxxhi = lammps_->domain->boxhi_bound[0];
407 boxylo = lammps_->domain->boxlo_bound[1];
408 boxyhi = lammps_->domain->boxhi_bound[1];
409 boxzlo = lammps_->domain->boxlo_bound[2];
410 boxzhi = lammps_->domain->boxhi_bound[2];
411 }
412 }
413
in_box(double * x) const414 bool LammpsInterface::in_box(double * x) const
415 {
416 double xlo,xhi,ylo,yhi,zlo,zhi;
417 box_bounds(xlo,xhi,ylo,yhi,zlo,zhi);
418 if (x[0] >= xlo && x[0] < xhi &&
419 x[1] >= ylo && x[1] < yhi &&
420 x[2] >= zlo && x[2] < zhi)
421 return true;
422 return false;
423 }
424
in_my_processor_box(double * x) const425 bool LammpsInterface::in_my_processor_box(double * x) const
426 {
427 if (x[0] >= lammps_->domain->sublo[0] && x[0] < lammps_->domain->subhi[0] &&
428 x[1] >= lammps_->domain->sublo[1] && x[1] < lammps_->domain->subhi[1] &&
429 x[2] >= lammps_->domain->sublo[2] && x[2] < lammps_->domain->subhi[2])
430 return true;
431 if (! in_box(x))
432 throw ATC_Error("point is in no processors box");
433 return false;
434 }
435
436
sub_bounds(double & subxlo,double & subxhi,double & subylo,double & subyhi,double & subzlo,double & subzhi) const437 void LammpsInterface::sub_bounds(double & subxlo, double & subxhi,
438 double & subylo, double & subyhi,
439 double & subzlo, double & subzhi) const
440 {
441 if (lammps_->domain->triclinic == 0) {
442 subxlo = lammps_->domain->sublo[0];
443 subxhi = lammps_->domain->subhi[0];
444 subylo = lammps_->domain->sublo[1];
445 subyhi = lammps_->domain->subhi[1];
446 subzlo = lammps_->domain->sublo[2];
447 subzhi = lammps_->domain->subhi[2];
448 }
449 else {
450 ATC_Error("Subboxes not accurate when triclinic != 0.");
451 }
452 }
453
xperiodic() const454 int LammpsInterface::xperiodic() const { return lammps_->domain->xperiodic; }
455
yperiodic() const456 int LammpsInterface::yperiodic() const { return lammps_->domain->yperiodic; }
457
zperiodic() const458 int LammpsInterface::zperiodic() const { return lammps_->domain->zperiodic; }
459
nperiodic() const460 int LammpsInterface::nperiodic() const
461 {
462 int nprd = 0;
463 if ( lammps_->domain->xperiodic > 0 ) { nprd++ ; }
464 if ( lammps_->domain->yperiodic > 0 ) { nprd++ ; }
465 if ( lammps_->domain->zperiodic > 0 ) { nprd++ ; }
466 return nprd;
467 }
468
469 // correct posistions for periodic box
periodicity_correction(double * x) const470 void LammpsInterface::periodicity_correction(double * x) const
471 {
472 int* periodicity = lammps_->domain->periodicity;
473 if (!refBoxIsSet_) set_reference_box();
474 for (int m = 0; m < 3; m++) {
475 if ((bool) periodicity[m]) {
476 if (x[m] < lower_[m] || x[m] > upper_[m]) {
477 x[m] -= length_[m]*floor((x[m]-lower_[m])/length_[m]);
478 }
479 if (x[m] < lower_[m] || x[m] > upper_[m]) {
480 throw ATC_Error("periodicity_correction: still out of box bounds");
481 }
482 }
483 }
484 }
485
set_reference_box(void) const486 void LammpsInterface::set_reference_box(void) const
487 {
488 double * hi = lammps_->domain->boxhi;
489 double * lo = lammps_->domain->boxlo;
490 double * len = lammps_->domain->prd;
491 for (int i = 0; i < 3; i++) {
492 upper_[i] = hi[i];
493 lower_[i] = lo[i];
494 length_[i] = len[i];
495 }
496 refBoxIsSet_ = true;
497 }
498
499
domain_xprd() const500 double LammpsInterface::domain_xprd() const { return lammps_->domain->xprd; }
501
domain_yprd() const502 double LammpsInterface::domain_yprd() const { return lammps_->domain->yprd; }
503
domain_zprd() const504 double LammpsInterface::domain_zprd() const { return lammps_->domain->zprd; }
505
domain_volume() const506 double LammpsInterface::domain_volume() const
507 {
508 return (lammps_->domain->xprd)*
509 (lammps_->domain->yprd)*
510 (lammps_->domain->zprd);
511 }
512
domain_xy() const513 double LammpsInterface::domain_xy() const { return lammps_->domain->xy; }
514
domain_xz() const515 double LammpsInterface::domain_xz() const { return lammps_->domain->xz; }
516
domain_yz() const517 double LammpsInterface::domain_yz() const { return lammps_->domain->yz; }
518
domain_triclinic() const519 int LammpsInterface::domain_triclinic() const { return lammps_->domain->triclinic; }
520
box_periodicity(int & xperiodic,int & yperiodic,int & zperiodic) const521 void LammpsInterface::box_periodicity(int & xperiodic,
522 int & yperiodic,
523 int & zperiodic) const
524 {
525 xperiodic = lammps_->domain->xperiodic;
526 yperiodic = lammps_->domain->yperiodic;
527 zperiodic = lammps_->domain->zperiodic;
528 }
529
region_id(const char * regionName) const530 int LammpsInterface::region_id(const char * regionName) const {
531 int nregion = this->nregion();
532 for (int iregion = 0; iregion < nregion; iregion++) {
533 if (strcmp(regionName, region_name(iregion)) == 0) {
534 return iregion;
535 }
536 }
537 throw ATC_Error("Region has not been defined");
538 return -1;
539 }
540
region_bounds(const char * regionName,double & xmin,double & xmax,double & ymin,double & ymax,double & zmin,double & zmax,double & xscale,double & yscale,double & zscale) const541 bool LammpsInterface::region_bounds(const char * regionName,
542 double & xmin, double & xmax,
543 double & ymin, double & ymax,
544 double & zmin, double & zmax,
545 double & xscale, double & yscale, double & zscale) const
546 {
547 int iRegion = region_id(regionName);
548 xscale = region_xscale(iRegion);
549 yscale = region_yscale(iRegion);
550 zscale = region_zscale(iRegion);
551 xmin = region_xlo(iRegion);
552 xmax = region_xhi(iRegion);
553 ymin = region_ylo(iRegion);
554 ymax = region_yhi(iRegion);
555 zmin = region_zlo(iRegion);
556 zmax = region_zhi(iRegion);
557 if (strcmp(region_style(iRegion),"block")==0) { return true; }
558 else { return false; }
559 }
560
minimum_image(double & dx,double & dy,double & dz) const561 void LammpsInterface::minimum_image(double & dx, double & dy, double & dz) const {
562 lammps_->domain->minimum_image(dx,dy,dz);
563 }
564
closest_image(const double * const xi,const double * const xj,double * const xjImage) const565 void LammpsInterface::closest_image(const double * const xi, const double * const xj, double * const xjImage) const {
566 lammps_->domain->closest_image(xi,xj,xjImage);
567 }
568
569
570 // -----------------------------------------------------------------
571 // update interface methods
572 // -----------------------------------------------------------------
units_style(void) const573 LammpsInterface::UnitsType LammpsInterface::units_style(void) const
574 {
575 if (strcmp(lammps_->update->unit_style,"lj") == 0) return LJ;
576 else if (strcmp(lammps_->update->unit_style,"real") == 0) return REAL;
577 else if (strcmp(lammps_->update->unit_style,"metal") == 0) return METAL;
578 else return UNKNOWN;
579 }
580
convert_units(double value,UnitsType in,UnitsType out,int massExp,int lenExp,int timeExp,int engExp) const581 double LammpsInterface::convert_units(double value, UnitsType in, UnitsType out, int massExp, int lenExp, int timeExp, int engExp) const
582 {
583 double ps2fs = 1.e3;
584 double eV2kcal = 23.069;
585 if (in==REAL) {
586 if (out==METAL) {
587 return value*pow(ps2fs,-timeExp)*pow(eV2kcal,-engExp);
588 }
589 else if (out==ATC) {
590 if (units_style()==REAL) {
591 return value;
592 }
593 else if (units_style()==METAL) {
594 return convert_units(value, METAL, out, massExp, lenExp, timeExp)*1.0;
595 }
596 }
597 else throw ATC_Error("can't convert");
598 }
599 else if (in==METAL) {
600 if (out==REAL) {
601 return value*pow(ps2fs,timeExp)*pow(eV2kcal,engExp);
602 }
603 else if (out==ATC) {
604 if (units_style()==REAL) {
605 return convert_units(value, REAL, out, massExp, lenExp, timeExp)*1.0;
606 }
607 else if (units_style()==METAL) {
608 return value;
609 }
610 }
611 else throw ATC_Error("can't convert");
612 }
613 else throw ATC_Error("can't convert");
614 return value;
615 }
616
617
618 // -----------------------------------------------------------------
619 // lattice interface methods
620 // -----------------------------------------------------------------
621
xlattice() const622 double LammpsInterface::xlattice() const { return lammps_->domain->lattice->xlattice; }
623
ylattice() const624 double LammpsInterface::ylattice() const { return lammps_->domain->lattice->ylattice; }
625
zlattice() const626 double LammpsInterface::zlattice() const { return lammps_->domain->lattice->zlattice; }
627
lattice_style() const628 LammpsInterface::LatticeType LammpsInterface::lattice_style() const
629 {
630 if (lammps_->domain->lattice)
631 return (LammpsInterface::LatticeType)lammps_->domain->lattice->style;
632 else
633 throw ATC_Error("Lattice has not been defined");
634 }
635
636 //* returns the number of basis vectors
n_basis() const637 int LammpsInterface::n_basis() const
638 {
639 return lammps_->domain->lattice->nbasis;
640 }
641
642 //* returns the basis vectors, transformed to the box coords
basis_vectors(double ** basis) const643 void LammpsInterface::basis_vectors(double **basis) const
644 {
645 LAMMPS_NS::Lattice *lattice = lammps_->domain->lattice;
646 int i,j;
647 double origin[3] = {0.0, 0.0, 0.0};
648 lattice->lattice2box(origin[0], origin[1], origin[2]);
649 for (i=0; i<n_basis(); i++)
650 {
651 memcpy(basis[i],lattice->basis[i],3*sizeof(double));
652 lattice->lattice2box(basis[i][0], basis[i][1], basis[i][2]);
653 for (j=0; j<3; j++) basis[i][j] -= origin[j];
654 }
655 }
656
657 //* gets the (max) lattice constant
max_lattice_constant(void) const658 double LammpsInterface::max_lattice_constant(void) const
659 {
660 double a1[3], a2[3], a3[3];
661 unit_cell(a1,a2,a3);
662 double a = norm(a1);
663 a = max(a,norm(a2));
664 a = max(a,norm(a3));
665 return a;
666 }
667
668 //* computes a cutoff distance halfway between 1st and 2nd nearest neighbors
near_neighbor_cutoff(void) const669 double LammpsInterface::near_neighbor_cutoff(void) const
670 {
671 double cutoff;
672 double alat = LammpsInterface::max_lattice_constant();
673 LatticeType type = lattice_style();
674 if (type == LammpsInterface::SC) {
675 cutoff = 0.5*(1.0+sqrt(2.0))*alat;
676 } else if (type == LammpsInterface::BCC) {
677 cutoff = 0.5*(0.5*sqrt(3.0)+1.0)*alat;
678 } else if (type == LammpsInterface::FCC) {
679 cutoff = 0.5*(1.0/sqrt(2.0)+1.0)*alat;
680 } else if (type == LammpsInterface::HCP) {
681 cutoff = 0.5*(1.0/sqrt(2.0)+1.0)*alat;
682 } else if (type == LammpsInterface::DIAMOND) {
683 cutoff = 0.5*(0.25*sqrt(3.0)+1.0/sqrt(2.0))*alat;
684 } else if (type == LammpsInterface::SQ) {
685 cutoff = 0.5*(1.0+sqrt(2.0))*alat;
686 } else if (type == LammpsInterface::SQ2) {
687 cutoff = 0.5*(1.0/sqrt(2.0)+1.0)*alat;
688 } else if (type == LammpsInterface::HEX) {
689 cutoff = 0.5*(1.0/sqrt(3.0)+1.0)*alat;
690 } else {
691 throw ATC_Error("Unknown lattice type");
692 }
693 return cutoff;
694 }
695
696 //* gets the unit cell vectors
unit_cell(double * a1,double * a2,double * a3) const697 void LammpsInterface::unit_cell(double *a1, double *a2, double *a3) const
698 {
699 int i, j;
700 double *a[3] = {a1,a2,a3};
701 double origin[3] = {0.0,0.0,0.0};
702 LAMMPS_NS::Lattice *lattice = lammps_->domain->lattice;
703 // transform origin
704 lattice->lattice2box(origin[0], origin[1], origin[2]);
705
706 // copy reference lattice vectors
707 memcpy(a[0], lattice->a1, 3*sizeof(double));
708 memcpy(a[1], lattice->a2, 3*sizeof(double));
709 memcpy(a[2], lattice->a3, 3*sizeof(double));
710
711 for (i=0; i<3; i++)
712 {
713 lattice->lattice2box(a[i][0], a[i][1], a[i][2]);
714 for (j=0; j<3; j++) a[i][j] -= origin[j];
715 }
716 }
717
718 //* gets number of atoms in a unit cell
num_atoms_per_cell(void) const719 int LammpsInterface::num_atoms_per_cell(void) const
720 {
721 int naCell = 0;
722 LatticeType type = lattice_style();
723 if (type == LammpsInterface::SC) naCell = 1;
724 else if (type == LammpsInterface::BCC) naCell = 2;
725 else if (type == LammpsInterface::FCC) naCell = 4;
726 else if (type == LammpsInterface::DIAMOND) naCell = 8;
727 else if (comm_rank()==0) {
728 //{throw ATC_Error("lattice style not currently supported by ATC");}
729 print_msg_once("WARNING: Cannot get number of atoms per cell from lattice");
730 naCell = 1;
731 }
732 return naCell;
733 }
734
735 //* gets tributary volume for an atom
volume_per_atom(void) const736 double LammpsInterface::volume_per_atom(void) const
737 {
738 double naCell = num_atoms_per_cell();
739 double volPerAtom =
740 xlattice() * ylattice() * zlattice() / naCell;
741 return volPerAtom;
742 }
743
744 //* gets lattice basis
lattice(MATRIX & N,MATRIX & B) const745 void LammpsInterface::lattice(MATRIX &N, MATRIX &B) const
746 {
747 int nbasis = n_basis();
748 double **basis = new double*[nbasis];
749 N.reset(3,3);
750 B.reset(3,nbasis);
751 for (int i=0; i<nbasis; i++) basis[i] = column(B,i).ptr();
752 basis_vectors(basis);
753 unit_cell(column(N,0).ptr(),
754 column(N,1).ptr(),
755 column(N,2).ptr());
756 delete [] basis;
757 }
758
759 // -----------------------------------------------------------------
760 // force interface methods
761 // -----------------------------------------------------------------
762
boltz() const763 double LammpsInterface::boltz() const{ return lammps_->force->boltz; }
764
mvv2e() const765 double LammpsInterface::mvv2e() const{ return lammps_->force->mvv2e; }
766
ftm2v() const767 double LammpsInterface::ftm2v()const { return lammps_->force->ftm2v; }
768
nktv2p() const769 double LammpsInterface::nktv2p()const{ return lammps_->force->nktv2p; }
770
qqr2e() const771 double LammpsInterface::qqr2e() const{ return lammps_->force->qqr2e; }
772
qe2f() const773 double LammpsInterface::qe2f() const{ return lammps_->force->qe2f; }
dielectric() const774 double LammpsInterface::dielectric()const{return lammps_->force->dielectric; }
775
qqrd2e() const776 double LammpsInterface::qqrd2e()const{ return lammps_->force->qqrd2e; }
777
qv2e() const778 double LammpsInterface::qv2e() const{ return qe2f()*ftm2v(); }
779
pair_force(int i,int j,double rsq,double & fmag_over_rmag) const780 double LammpsInterface::pair_force(int i, int j, double rsq,
781 double & fmag_over_rmag) const
782 {
783 int itype = (lammps_->atom->type)[i];
784 int jtype = (lammps_->atom->type)[j];
785 // return value is the energy
786 if (rsq < (lammps_->force->pair->cutsq)[itype][jtype]) {
787 return lammps_->force->pair->single(i,j,itype,jtype,
788 rsq,1.0,1.0,fmag_over_rmag);
789 }
790 return 0.0;
791 }
pair_force(int n,double rsq,double & fmag_over_rmag) const792 double LammpsInterface::pair_force(int n, double rsq,
793 double & fmag_over_rmag) const
794 {
795 int i = bond_list_i(n);
796 int j = bond_list_j(n);
797 int type = bond_list_type(n);
798 // return value is the energy
799 return lammps_->force->bond->single(type,rsq,i,j,fmag_over_rmag);
800 }
pair_force(map<std::pair<int,int>,int>::const_iterator itr,double rsq,double & fmag_over_rmag,int nbonds) const801 double LammpsInterface::pair_force(
802 map< std::pair< int,int >,int >::const_iterator itr, double rsq,
803 double & fmag_over_rmag, int nbonds) const
804 {
805 int n = itr->second;
806 if (n < nbonds) {
807 return pair_force(n, rsq,fmag_over_rmag);
808 }
809 else {
810 std::pair <int,int> ij = itr->first;
811 int i = ij.first;
812 int j = ij.second;
813 return pair_force(i,j, rsq,fmag_over_rmag);
814 }
815 }
pair_force(std::pair<std::pair<int,int>,int> apair,double rsq,double & fmag_over_rmag,int nbonds) const816 double LammpsInterface::pair_force(
817 std::pair< std::pair< int,int >,int > apair, double rsq,
818 double & fmag_over_rmag, int nbonds) const
819 {
820 int n = apair.second;
821 if (n < nbonds) {
822 return pair_force(n, rsq,fmag_over_rmag);
823 }
824 else {
825 std::pair <int,int> ij = apair.first;
826 int i = ij.first;
827 int j = ij.second;
828 return pair_force(i,j, rsq,fmag_over_rmag);
829 }
830 }
bond_stiffness(int i,int j,double rsq0) const831 double LammpsInterface::bond_stiffness(int i, int j, double rsq0) const
832 {
833 const double perturbation = 1.e-8;
834 double rsq1 = sqrt(rsq0)+perturbation;
835 rsq1 *= rsq1;
836 double f0,f1;
837 pair_force(i,j,rsq0,f0);
838 pair_force(i,j,rsq1,f1);
839 double k = (f1-f0)/perturbation;
840 return k;
841 }
842
pair_cutoff() const843 double LammpsInterface::pair_cutoff() const
844 {
845 return lammps_->force->pair->cutforce;
846 }
847
pair_reinit() const848 void LammpsInterface::pair_reinit() const
849 {
850 lammps_->force->pair->reinit();
851 }
852
single_enable() const853 int LammpsInterface::single_enable() const
854 {
855 return lammps_->force->pair->single_enable; // all bonds have a single
856 }
857
858 //* insertion/deletion functions : see FixGCMC
859 //* delete atom
delete_atom(int id) const860 int LammpsInterface::delete_atom(int id) const
861 {
862 LAMMPS_NS::Atom * atom = lammps_->atom;
863 atom->avec->copy(atom->nlocal-1,id,1);
864 atom->nlocal--;
865 return atom->nlocal;
866 }
867
868 //* insert atom
insert_atom(int atype,int amask,double * ax,double * av,double aq) const869 int LammpsInterface::insert_atom(int atype, int amask,
870 double *ax, double *av, double aq) const
871 {
872 LAMMPS_NS::Atom * atom = lammps_->atom;
873 atom->avec->create_atom(atype,ax);
874 int m = atom->nlocal - 1;
875 atom->mask[m] = amask;
876 atom->v[m][0] = av[0];
877 atom->v[m][1] = av[1];
878 atom->v[m][2] = av[2];
879 if (aq != 0) atom->q[m] = aq;
880
881 int nfix = lammps_->modify->nfix;
882 LAMMPS_NS::Fix **fix = lammps_->modify->fix;
883 for (int j = 0; j < nfix; j++) {
884 if (fix[j]->create_attribute) fix[j]->set_arrays(m);
885 }
886 return m;
887 }
888
reset_ghosts(int deln) const889 int LammpsInterface::reset_ghosts(int deln) const
890 {
891 LAMMPS_NS::Atom * atom = lammps_->atom;
892 //ATC::LammpsInterface::instance()->print_msg("reset_ghosts del n "+to_string(deln));
893 if (atom->tag_enable) {
894 atom->natoms += deln;
895 //ATC::LammpsInterface::instance()->print_msg("reset_ghosts natoms "+to_string(atom->natoms));
896 if (deln > 0) { atom->tag_extend(); }
897 if (atom->map_style) atom->map_init();
898 }
899 atom->nghost = 0;
900 lammps_->comm->borders();
901 //ATC::LammpsInterface::instance()->print_msg("reset_ghosts nghosts "+to_string(atom->nghost));
902 return atom->nghost;
903 }
904
905 //* energy for interactions within the shortrange cutoff
shortrange_energy(double * coord,int itype,int id,double) const906 double LammpsInterface::shortrange_energy(double *coord,
907 int itype, int id, double /* max */) const
908 {
909 LAMMPS_NS::Atom * atom = lammps_->atom;
910 double **x = atom->x;
911 int *type = atom->type;
912 int nall = atom->nlocal+ atom->nghost;
913 LAMMPS_NS::Pair *pair = lammps_->force->pair;
914 double **cutsq = lammps_->force->pair->cutsq;
915
916 double fpair = 0.0; // an output of single
917 double factor_coul = 1.0;
918 double factor_lj = 1.0;
919
920 double total_energy = 0.0;
921 for (int j = 0; j < nall; j++) {
922 if (id == j) continue;
923 // factor_lj = special_lj[sbmask(j)];
924 // factor_coul = special_coul[sbmask(j)];
925 //j &= NEIGHMASK;
926 double delx = coord[0] - x[j][0];
927 double dely = coord[1] - x[j][1];
928 double delz = coord[2] - x[j][2];
929 double rsq = delx*delx + dely*dely + delz*delz;
930 int jtype = type[j];
931 double cut2 = cutsq[itype][jtype];
932 if (rsq < cut2) {
933 total_energy += pair->single(id,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); //id is for charge lookup
934 }
935 }
936 return total_energy;
937 }
shortrange_energy(int id,double max) const938 double LammpsInterface::shortrange_energy(int id, double max) const
939 {
940 double *x = (lammps_->atom->x)[id];
941 int type = (lammps_->atom->type)[id];
942 return shortrange_energy(x,type,id,max);
943 }
944
potential() const945 POTENTIAL LammpsInterface::potential() const
946 {
947 // find pair style - FRAGILE
948 const int nStyles = 4;
949 string pairStyles[nStyles] = {"lj/cut",
950 "lj/cut/coul/long",
951 "lj/cut/coul/cut",
952 "lj/charmm/coul/long"};
953 LAMMPS_NS::Pair *pair = nullptr;
954 for (int i = 0; i < nStyles; i++){
955 pair = lammps_->force->pair_match(pairStyles[i].c_str(),1);
956 if (pair != nullptr) break;
957 }
958 return pair;
959 }
960
type_to_groupbit(int itype) const961 int LammpsInterface::type_to_groupbit(int itype) const
962 {
963 LAMMPS_NS::Atom * atom = lammps_->atom;
964 int groupbit = -1;
965 int *type = atom->type;
966 int *mask = atom->mask;
967 for (int i = 0; i < nlocal(); i++) {
968 if (type[i] == itype) {
969 groupbit = mask[i];
970 break;
971 }
972 }
973 return int_allmax(groupbit);
974 }
975
epsilons(int itype,POTENTIAL pair,double * epsilon0) const976 bool LammpsInterface::epsilons(int itype, POTENTIAL pair, double * epsilon0) const
977 {
978 // grab energy parameters
979 char * pair_parameter = new char[8];
980 strcpy(pair_parameter,"epsilon");
981 int dim = 2; // a return value for extract
982 double ** epsilons = (double**) ( pair->extract(pair_parameter,dim) );
983 delete [] pair_parameter;
984 if (epsilons == nullptr) return false;
985 //if (epsilons == nullptr) error->all(FLERR,"Fix concentration adapted pair style parameter not supported");
986 int i1,i2;
987 for (int i=1; i < ntypes()+1; i++) {
988 if (i < itype) { i1 = i; i2 = itype; }
989 else { i2 = i; i1 = itype; }
990 epsilon0[i-1] = epsilons[i1][i2];
991 }
992 return true;
993 }
994
set_epsilons(int itype,POTENTIAL pair,double * epsilon) const995 bool LammpsInterface::set_epsilons(int itype, POTENTIAL pair, double * epsilon) const
996 {
997 // grab energy parameters
998 char * pair_parameter = new char[8];
999 strcpy(pair_parameter,"epsilon");
1000 int dim = 2; // a return value for extract
1001 double ** epsilons = (double**) ( pair->extract(pair_parameter,dim) );
1002 delete [] pair_parameter;
1003 if (epsilons == nullptr) return false;
1004 //if (epsilons == nullptr) error->all(FLERR,"Fix concentration adapted pair style parameter not supported");
1005 // scale interactions
1006 int i1,i2;
1007 for (int i = 1; i < ntypes()+1; i++) {
1008 if (i < itype) { i1 = i; i2 = itype; }
1009 else { i2 = i; i1 = itype; }
1010 epsilons[i1][i2] = epsilon[i-1];
1011 }
1012
1013 return true;
1014 }
1015
set_charge(int itype,double charge) const1016 int LammpsInterface::set_charge(int itype, double charge) const
1017 {
1018 int nlocal = lammps_->atom->nlocal;
1019 int *type = lammps_->atom->type;
1020 double *q = lammps_->atom->q;
1021 int count = 0;
1022 for (int i = 0; i < nlocal; i++) {
1023 if (type[i] == itype) {
1024 q[i] = charge;
1025 count++;
1026 }
1027 }
1028 return count;
1029 }
1030
change_type(int itype,int jtype) const1031 int LammpsInterface::change_type(int itype, int jtype) const
1032 {
1033 int nlocal = lammps_->atom->nlocal;
1034 int *type = lammps_->atom->type;
1035 int count = 0;
1036 for (int i = 0; i < nlocal; i++) {
1037 if (type[i] == itype) {
1038 type[i] = jtype;
1039 count++;
1040 }
1041 }
1042 return count;
1043 }
1044
count_type(int itype) const1045 int LammpsInterface::count_type(int itype) const
1046 {
1047 int nlocal = lammps_->atom->nlocal;
1048 int *type = lammps_->atom->type;
1049 int count = 0;
1050 for (int i = 0; i < nlocal; i++) {
1051 if (type[i] == itype) { count++; }
1052 }
1053 return int_allsum(count);
1054 }
1055
1056 // random number generators
random_number_generator() const1057 RNG_POINTER LammpsInterface::random_number_generator() const {
1058 RNG_POINTER p = new LAMMPS_NS::RanPark(lammps_,seed_);
1059 return p;
1060 }
random_uniform(RNG_POINTER p) const1061 double LammpsInterface::random_uniform(RNG_POINTER p) const {
1062 return p->uniform();
1063 }
random_normal(RNG_POINTER p) const1064 double LammpsInterface::random_normal (RNG_POINTER p) const {
1065 return p->gaussian();
1066 }
random_state(RNG_POINTER p) const1067 int LammpsInterface::random_state (RNG_POINTER p) const {
1068 return p->state();
1069 }
set_random_state(RNG_POINTER p,int seed) const1070 void LammpsInterface::set_random_state (RNG_POINTER p, int seed) const {
1071 return p->reset(seed);
1072 }
advance_random_generator(RNG_POINTER p,int n) const1073 void LammpsInterface::advance_random_generator (RNG_POINTER p, int n) const {
1074 advance_random_uniform(p,n);
1075 }
advance_random_uniform(RNG_POINTER p,int n) const1076 void LammpsInterface::advance_random_uniform (RNG_POINTER p, int n) const {
1077 for (int i = 0; i < n; i++) p->uniform();
1078 }
advance_random_normal(RNG_POINTER p,int n) const1079 void LammpsInterface::advance_random_normal (RNG_POINTER p, int n) const {
1080 for (int i = 0; i < n; i++) p->gaussian();
1081 }
1082
1083 //* Boltzmann's constant in M,L,T,t units
kBoltzmann() const1084 double LammpsInterface::kBoltzmann() const {
1085 return (lammps_->force->boltz)/(lammps_->force->mvv2e);
1086 }
1087
1088 //* Planck's constant
hbar() const1089 double LammpsInterface::hbar() const {
1090 const int UNITS_STYLE = (int) units_style();
1091 double hbar = 1.0; // LJ: Dimensionless
1092 if (UNITS_STYLE == 2) hbar = 15.1685792814; // Real: KCal/mol-fs
1093 else if (UNITS_STYLE == 3) hbar = 0.000658212202469; // Metal: eV-ps
1094 return hbar;
1095 }
1096
1097 //* Dulong-Petit heat capacity
heat_capacity() const1098 double LammpsInterface::heat_capacity() const {
1099 double rhoCp = dimension()*kBoltzmann()/volume_per_atom();
1100 return rhoCp;
1101 }
1102
1103 //* reference mass density for a *unit cell*
1104 // all that is needed is a unit cell: volume, types, mass per type
mass_density(int * numPerType) const1105 double LammpsInterface::mass_density(int* numPerType) const
1106 {
1107 const double *mass = lammps_->atom->mass;
1108 if (!mass) throw ATC_Error("cannot compute a mass density: no mass");
1109 const int ntypes = lammps_->atom->ntypes;
1110 const int *mass_setflag = lammps_->atom->mass_setflag;
1111 const int *type = lammps_->atom->type;
1112 double naCell = num_atoms_per_cell();
1113 double vol = volume_per_atom();
1114 if (numPerType) {
1115 double m = 0.0;
1116 double n = 0;
1117 for (int i = 0; i < ntypes; i++) {
1118 m += numPerType[i]*mass[i+1];
1119 n += numPerType[i];
1120 }
1121 if (n>naCell) throw ATC_Error("cannot compute a mass density: too many basis atoms");
1122 return m/n/vol;
1123 }
1124 // since basis->type map not stored only monatomic lattices are automatic
1125 // if not given a basis try to guess
1126 else {
1127 if (ntypes == 1) {
1128 if ((this->natoms()==0) && mass_setflag[1]) {
1129 return mass[1]/vol;
1130 }
1131 else {
1132 if (type) return mass[type[0]]/vol;
1133 else if (mass_setflag[1]) return mass[1]/vol;
1134 }
1135 }
1136 }
1137 throw ATC_Error("cannot compute a mass density");
1138 return 0.0;
1139 }
1140
1141 //* permittivity of free space
epsilon0() const1142 double LammpsInterface::epsilon0() const
1143 {
1144 return qe2f()/(4.*PI*qqr2e());
1145 }
1146
1147 //* Coulomb's constant
coulomb_constant() const1148 double LammpsInterface::coulomb_constant() const
1149 {
1150 return qqr2e()/qe2f();
1151 }
1152
1153 //* special coulombic interactions
special_coul() const1154 double * LammpsInterface::special_coul() const
1155 {
1156 return lammps_->force->special_coul;
1157 }
1158
1159 //* flag for newton
newton_pair() const1160 int LammpsInterface::newton_pair() const
1161 {
1162 return lammps_->force->newton_pair;
1163 }
1164
1165 // -----------------------------------------------------------------
1166 // group interface methods
1167 // -----------------------------------------------------------------
1168
ngroup() const1169 int LammpsInterface::ngroup() const { return lammps_->group->ngroup; }
1170
group_bit(string name) const1171 int LammpsInterface::group_bit(string name) const
1172 {
1173 return group_bit(group_index(name));
1174 }
1175
group_bit(int iGroup) const1176 int LammpsInterface::group_bit(int iGroup) const
1177 {
1178 int mybit = 0;
1179 mybit |= lammps_->group->bitmask[iGroup];
1180 if (mybit < 0 || mybit > MAX_GROUP_BIT) {
1181 string msg("LammpsInterface::group_bit() lammps group bit "+to_string(mybit)+" is out of range 0:"+to_string(MAX_GROUP_BIT));
1182 throw ATC_Error(msg);
1183 }
1184
1185 return mybit;
1186 }
1187
group_index(string name) const1188 int LammpsInterface::group_index(string name) const
1189 {
1190 int igroup = lammps_->group->find(name.c_str());
1191 if (igroup == -1) {
1192 string msg("LammpsInterface::group_index() lammps group "+name+" does not exist");
1193 throw ATC_Error(msg);
1194 }
1195 return igroup;
1196 }
1197
group_inverse_mask(int iGroup) const1198 int LammpsInterface::group_inverse_mask(int iGroup) const
1199 {
1200 return lammps_->group->inversemask[iGroup];
1201 }
1202
group_name(int iGroup) const1203 char * LammpsInterface::group_name(int iGroup) const
1204 {
1205 return lammps_->group->names[iGroup];
1206 }
1207
group_bounds(int iGroup,double * b) const1208 void LammpsInterface::group_bounds(int iGroup, double * b) const
1209 {
1210 lammps_->group->bounds(iGroup, b);
1211 }
1212
1213 // -----------------------------------------------------------------
1214 // memory interface methods
1215 // -----------------------------------------------------------------
1216
1217
create_1d_double_array(int length,const char * name) const1218 double * LammpsInterface::create_1d_double_array(int length, const char *name) const {
1219 double * myArray;
1220 return lammps_->memory->create(myArray, length, name);
1221 }
1222
grow_1d_double_array(double * array,int length,const char * name) const1223 double *LammpsInterface::grow_1d_double_array(double *array,
1224 int length,
1225 const char *name) const
1226 {
1227 return lammps_->memory->grow(array, length, name);
1228 }
1229
destroy_1d_double_array(double * d) const1230 void LammpsInterface::destroy_1d_double_array(double * d) const {
1231 lammps_->memory->destroy(d);
1232 }
1233
create_2d_double_array(int n1,int n2,const char * name) const1234 double ** LammpsInterface::create_2d_double_array(int n1, int n2, const char *name) const {
1235 double ** myArray;
1236 return lammps_->memory->create(myArray, n1, n2, name);
1237 }
1238
destroy_2d_double_array(double ** d) const1239 void LammpsInterface::destroy_2d_double_array(double **d) const {
1240 lammps_->memory->destroy(d);
1241 }
1242
grow_2d_double_array(double ** array,int n1,int n2,const char * name) const1243 double **LammpsInterface::grow_2d_double_array(double **array,
1244 int n1,
1245 int n2,
1246 const char *name) const
1247 {
1248 return lammps_->memory->grow(array, n1, n2, name);
1249 }
1250
1251
create_1d_int_array(int length,const char * name) const1252 int * LammpsInterface::create_1d_int_array(int length, const char *name) const {
1253 int * myArray;
1254 return lammps_->memory->create(myArray, length, name);
1255 }
1256
grow_1d_int_array(int * array,int length,const char * name) const1257 int *LammpsInterface::grow_1d_int_array(int *array,
1258 int length,
1259 const char *name) const
1260 {
1261 return lammps_->memory->grow(array, length, name);
1262 }
1263
destroy_1d_int_array(int * d) const1264 void LammpsInterface::destroy_1d_int_array(int * d) const {
1265 lammps_->memory->destroy(d);
1266 }
1267
create_2d_int_array(int n1,int n2,const char * name) const1268 int ** LammpsInterface::create_2d_int_array(int n1, int n2, const char *name) const {
1269 int ** myArray;
1270 return lammps_->memory->create(myArray, n1, n2, name);
1271 }
1272
destroy_2d_int_array(int ** i) const1273 void LammpsInterface::destroy_2d_int_array(int **i) const {
1274 lammps_->memory->destroy(i);
1275 }
1276
grow_2d_int_array(int ** array,int n1,int n2,const char * name) const1277 int ** LammpsInterface::grow_2d_int_array(int **array, int n1, int n2, const char *name) const {
1278 return lammps_->memory->grow(array, n1, n2, name);
1279 }
1280
1281 // -----------------------------------------------------------------
1282 // update interface methods
1283 // -----------------------------------------------------------------
1284
dt() const1285 double LammpsInterface::dt() const { return lammps_->update->dt; }
1286
ntimestep() const1287 bigint LammpsInterface::ntimestep() const { return lammps_->update->ntimestep; }
1288
nsteps() const1289 int LammpsInterface::nsteps() const { return lammps_->update->nsteps; }
1290
1291
1292 // -----------------------------------------------------------------
1293 // neighbor list interface methods
1294 // -----------------------------------------------------------------
1295
sbmask(int j) const1296 int LammpsInterface::sbmask(int j) const {return j >> SBBITS & 3; }
1297
set_list(int,LAMMPS_NS::NeighList * ptr)1298 void LammpsInterface::set_list(int /* id */, LAMMPS_NS::NeighList *ptr) { list_ = ptr; }
1299
neighbor_list_inum() const1300 int LammpsInterface::neighbor_list_inum() const { return list_->inum; }
1301
neighbor_list_numneigh() const1302 int * LammpsInterface::neighbor_list_numneigh() const { return list_->numneigh; }
1303
neighbor_list_ilist() const1304 int * LammpsInterface::neighbor_list_ilist() const { return list_->ilist; }
1305
neighbor_list_firstneigh() const1306 int ** LammpsInterface::neighbor_list_firstneigh() const { return list_->firstneigh; }
1307
neighbor_ago() const1308 int LammpsInterface::neighbor_ago() const { return lammps_->neighbor->ago; }
1309
reneighbor_frequency() const1310 int LammpsInterface::reneighbor_frequency() const {return lammps_->neighbor->every; }
1311
1312 // -----------------------------------------------------------------
1313 // bond list interface methods
1314 // -----------------------------------------------------------------
1315
bond_list_length() const1316 int LammpsInterface::bond_list_length() const { return lammps_->neighbor->nbondlist; }
bond_list() const1317 int** LammpsInterface::bond_list() const { return lammps_->neighbor->bondlist; }
1318
1319 // -----------------------------------------------------------------
1320 // region interface methods
1321 // -----------------------------------------------------------------
1322
region_name(int iRegion) const1323 char * LammpsInterface::region_name(int iRegion) const
1324 {
1325 return lammps_->domain->regions[iRegion]->id;
1326 }
1327
region_style(int iRegion) const1328 char * LammpsInterface::region_style(int iRegion) const
1329 {
1330 return lammps_->domain->regions[iRegion]->style;
1331 }
1332
region_xlo(int iRegion) const1333 double LammpsInterface::region_xlo(int iRegion) const
1334 {
1335 return lammps_->domain->regions[iRegion]->extent_xlo;
1336 }
1337
region_xhi(int iRegion) const1338 double LammpsInterface::region_xhi(int iRegion) const
1339 {
1340 return lammps_->domain->regions[iRegion]->extent_xhi;
1341 }
1342
region_ylo(int iRegion) const1343 double LammpsInterface::region_ylo(int iRegion) const
1344 {
1345 return lammps_->domain->regions[iRegion]->extent_ylo;
1346 }
1347
region_yhi(int iRegion) const1348 double LammpsInterface::region_yhi(int iRegion) const
1349 {
1350 return lammps_->domain->regions[iRegion]->extent_yhi;
1351 }
1352
region_zlo(int iRegion) const1353 double LammpsInterface::region_zlo(int iRegion) const
1354 {
1355 return lammps_->domain->regions[iRegion]->extent_zlo;
1356 }
1357
region_zhi(int iRegion) const1358 double LammpsInterface::region_zhi(int iRegion) const
1359 {
1360 return lammps_->domain->regions[iRegion]->extent_zhi;
1361 }
1362
region_xscale(int iRegion) const1363 double LammpsInterface::region_xscale(int iRegion) const
1364 {
1365 return lammps_->domain->regions[iRegion]->xscale;
1366 }
1367
region_yscale(int iRegion) const1368 double LammpsInterface::region_yscale(int iRegion) const
1369 {
1370 return lammps_->domain->regions[iRegion]->yscale;
1371 }
1372
region_zscale(int iRegion) const1373 double LammpsInterface::region_zscale(int iRegion) const
1374 {
1375 return lammps_->domain->regions[iRegion]->zscale;
1376 }
1377
region_match(int iRegion,double x,double y,double z) const1378 int LammpsInterface::region_match(int iRegion, double x, double y, double z) const {
1379 return lammps_->domain->regions[iRegion]->match(x,y,z);
1380 }
1381
1382 // -----------------------------------------------------------------
1383 // compute methods
1384 // -----------------------------------------------------------------
compute_pointer(string tag) const1385 COMPUTE_POINTER LammpsInterface::compute_pointer(string tag) const
1386 {
1387 // get the compute id
1388 char * name = const_cast <char*> (tag.c_str());
1389 int id = lammps_->modify->find_compute(name);
1390 if (id < 0) {
1391 string msg("Could not find compute "+tag);
1392 msg += tag;
1393 throw ATC_Error(msg);
1394 }
1395 // get the compute
1396 LAMMPS_NS::Compute* cmpt = lammps_->modify->compute[id];
1397 // insert it into our set, recall it won't be added if it already exists
1398 computePointers_.insert(cmpt);
1399 return cmpt;
1400 }
1401
computes_addstep(int step) const1402 void LammpsInterface::computes_addstep(int step) const
1403 {
1404 set<LAMMPS_NS::Compute * >::iterator iter;
1405 for (iter = computePointers_.begin(); iter != computePointers_.end(); iter++) {
1406 (*iter)->addstep(step);
1407 }
1408 }
1409
compute_addstep(COMPUTE_POINTER computePointer,int step) const1410 void LammpsInterface::compute_addstep(COMPUTE_POINTER computePointer, int step) const
1411 {
1412 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1413 cmpt->addstep(step);
1414 }
1415
compute_matchstep(COMPUTE_POINTER computePointer,int step) const1416 int LammpsInterface::compute_matchstep(COMPUTE_POINTER computePointer, int step) const
1417 {
1418 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1419 return cmpt->matchstep(step);
1420 }
1421
reset_invoked_flag(COMPUTE_POINTER computePointer) const1422 void LammpsInterface::reset_invoked_flag(COMPUTE_POINTER computePointer) const
1423 {
1424 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1425 cmpt->invoked_flag = 0;
1426 }
1427
compute_ncols_peratom(COMPUTE_POINTER computePointer) const1428 int LammpsInterface::compute_ncols_peratom(COMPUTE_POINTER computePointer) const
1429 {
1430 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1431 int ndof = cmpt->size_peratom_cols;
1432 if (ndof == 0 ) ndof = 1;
1433 return ndof;
1434 }
1435
compute_vector_peratom(COMPUTE_POINTER computePointer) const1436 double* LammpsInterface::compute_vector_peratom(COMPUTE_POINTER computePointer) const
1437 {
1438 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1439 if (!(cmpt->invoked_flag & INVOKED_PERATOM)) {
1440 cmpt->compute_peratom();
1441 cmpt->invoked_flag |= INVOKED_PERATOM;
1442 }
1443 return cmpt->vector_atom;
1444 }
1445
compute_array_peratom(COMPUTE_POINTER computePointer) const1446 double** LammpsInterface::compute_array_peratom(COMPUTE_POINTER computePointer) const
1447 {
1448 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1449 if (!(cmpt->invoked_flag & INVOKED_PERATOM)) {
1450 cmpt->compute_peratom();
1451 cmpt->invoked_flag |= INVOKED_PERATOM;
1452 }
1453 return cmpt->array_atom;
1454 }
1455
const_to_active(COMPUTE_POINTER computePointer) const1456 LAMMPS_NS::Compute * LammpsInterface::const_to_active(COMPUTE_POINTER computePointer) const
1457 {
1458 LAMMPS_NS::Compute* cmpt = const_cast<LAMMPS_NS::Compute* >(computePointer);
1459 set<LAMMPS_NS::Compute * >::iterator cmptPtr;
1460 cmptPtr = computePointers_.find(cmpt);
1461 if (cmptPtr != computePointers_.end())
1462 return *cmptPtr;
1463 else
1464 throw ATC_Error("Requested bad computer pointer");
1465 }
1466
1467 // -----------------------------------------------------------------
1468 // compute pe/atom interface methods
1469 // - the only compute "owned" by ATC
1470 // -----------------------------------------------------------------
create_compute_pe_peratom(void) const1471 int LammpsInterface::create_compute_pe_peratom(void) const
1472 {
1473 char **list = new char*[4];
1474 string atomPeName = compute_pe_name();
1475 list[0] = (char *) atomPeName.c_str();
1476 list[1] = (char *) "all";
1477 list[2] = (char *) "pe/atom";
1478 list[3] = (char *) "pair";
1479
1480 int icompute = lammps_->modify->find_compute(list[0]);
1481 if (icompute < 0) {
1482 lammps_->modify->add_compute(4,list);
1483 icompute = lammps_->modify->find_compute(list[0]);
1484 }
1485 delete [] list;
1486 if (! atomPE_ ) {
1487 atomPE_ = lammps_->modify->compute[icompute];
1488 }
1489 computePointers_.insert(atomPE_);
1490 stringstream ss;
1491 ss << "peratom PE compute created with ID: " << icompute;
1492 ATC::LammpsInterface::instance()->print_msg_once(ss.str());
1493 return icompute;
1494 }
1495
compute_pe_peratom(void) const1496 double * LammpsInterface::compute_pe_peratom(void) const
1497 {
1498 if (atomPE_) {
1499 atomPE_->compute_peratom();
1500 return atomPE_->vector_atom;
1501 }
1502 else {
1503 return nullptr;
1504 }
1505 }
1506
1507 /* ---------------------------------------------------------------------- */
1508
unwrap_coordinates(int iatom,double * xatom) const1509 void LammpsInterface::unwrap_coordinates(int iatom, double* xatom) const
1510 {
1511 double **x = lammps_->atom->x;
1512 int *image = lammps_->atom->image;
1513
1514 double *h = lammps_->domain->h;
1515 double xprd = lammps_->domain->xprd;
1516 double yprd = lammps_->domain->yprd;
1517 double zprd = lammps_->domain->zprd;
1518 int xbox,ybox,zbox;
1519
1520 // for triclinic, need to unwrap current atom coord via h matrix
1521
1522 if (lammps_->domain->triclinic == 0) {
1523 xbox = (image[iatom] & 1023) - 512;
1524 ybox = (image[iatom] >> 10 & 1023) - 512;
1525 zbox = (image[iatom] >> 20) - 512;
1526 xatom[0] = x[iatom][0] + xbox*xprd;
1527 xatom[1] = x[iatom][1] + ybox*yprd;
1528 xatom[2] = x[iatom][2] + zbox*zprd;
1529 } else {
1530 xbox = (image[iatom] & 1023) - 512;
1531 ybox = (image[iatom] >> 10 & 1023) - 512;
1532 zbox = (image[iatom] >> 20) - 512;
1533 xatom[0] = x[iatom][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
1534 xatom[1] = x[iatom][1] + h[1]*ybox + h[3]*zbox;
1535 xatom[2] = x[iatom][2] + h[2]*zbox;
1536 }
1537 }
1538
1539 /* ---------------------------------------------------------------------- */
1540
pair_eam() const1541 LAMMPS_NS::PairEAM* LammpsInterface::pair_eam() const
1542 {
1543 //if (typeid(lammps_->force->pair) == typeid(LAMMPS_NS::PairEAM)) {
1544 // return lammps_->force->pair;
1545 //}
1546 LAMMPS_NS::PairEAM* pair_eam = dynamic_cast<LAMMPS_NS::PairEAM*> (lammps_->force->pair);
1547 if (pair_eam != nullptr) {
1548 return pair_eam;
1549 }
1550 else {
1551 throw ATC_Error("LAMMPS Pair object is not of the derived class type PairEAM");
1552 }
1553 }
1554
1555 // -----------------------------------------------------------------
1556 // other methods
1557 // -----------------------------------------------------------------
1558
1559 /** Return lammps pointer -- only as a last resort! */
lammps_pointer() const1560 LAMMPS_NS::LAMMPS * LammpsInterface::lammps_pointer() const { return lammps_; }
1561
1562 }
1563