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
30 // ATC includes
31 #include "ATC_Error.h"
32 #include "MatrixLibrary.h"
33 #include "Utility.h"
34 using ATC_Utility::to_string;
35
36 // Other include files
37 #include "mpi.h"
38 #include <cstring>
39 #include <map>
40 #include <typeinfo>
41
42 using std::max;
43 using std::stringstream;
44 using std::copy;
45 using std::map;
46 using std::pair;
47 using std::string;
48 using std::set;
49 using LAMMPS_NS::bigint;
50
51 namespace ATC {
52
53 const static double PI = 3.141592653589793238;
54
55 const static int seed_ = 3141592;
56
57 const static int MAX_GROUP_BIT = 2147483647; //4294967295; // pow(2,31)-1;
58
norm(double * v)59 double norm(double * v) {return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); }
60
61 LammpsInterface * LammpsInterface::myInstance_ = NULL;
62
63 // -----------------------------------------------------------------
64 // instance()
65 // -----------------------------------------------------------------
instance()66 LammpsInterface * LammpsInterface::instance()
67 {
68 if (myInstance_ == NULL) {
69 myInstance_ = new LammpsInterface();
70 }
71 return myInstance_;
72 }
73
74 // -----------------------------------------------------------------
75 // Destroy()
76 // -----------------------------------------------------------------
Destroy()77 void LammpsInterface::Destroy()
78 {
79 if (myInstance_) delete myInstance_;
80 myInstance_ = NULL;
81 }
82
83
84 // -----------------------------------------------------------------
85 // constructor
86 // -----------------------------------------------------------------
LammpsInterface()87 LammpsInterface::LammpsInterface()
88 : lammps_(NULL),
89 fixPointer_(NULL),
90 commRank_(0),
91 atomPE_(NULL),
92 refBoxIsSet_(false),
93 random_(NULL),
94 globalrandom_(NULL)
95 {
96 }
97
98 // -----------------------------------------------------------------
99 // general interface methods
100 // -----------------------------------------------------------------
101
world() const102 MPI_Comm LammpsInterface::world() const { return lammps_->world; }
set_fix_pointer(LAMMPS_NS::Fix * thisFix)103 void LammpsInterface::set_fix_pointer(LAMMPS_NS::Fix * thisFix) { fixPointer_ = thisFix; }
forward_comm_fix() const104 void LammpsInterface::forward_comm_fix() const { lammps_->comm->forward_comm_fix(fixPointer_); }
comm_borders() const105 void LammpsInterface::comm_borders() const { lammps_->comm->borders(); }
106
107 #ifndef ISOLATE_FE
sparse_allsum(SparseMatrix<double> & toShare) const108 void LammpsInterface::sparse_allsum(SparseMatrix<double> &toShare) const
109 {
110 toShare.compress();
111
112 // initialize MPI information
113 int nProcs;
114 int myRank;
115 MPI_Comm_size(MPI_COMM_WORLD, &nProcs);
116 MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
117
118 int error;
119
120 // get numbers of rows, columns, rowsCRS, and
121 // sizes (number of nonzero elements in matrix)
122 SparseMatInfo *recInfo = new SparseMatInfo[nProcs];
123 SparseMatInfo myInfo;
124 myInfo.rows = toShare.nRows();
125 myInfo.cols = toShare.nCols();
126 myInfo.rowsCRS = toShare.nRowsCRS();
127 myInfo.size = toShare.size();
128
129 error = MPI_Allgather(&myInfo, 4, MPI_INT,
130 recInfo, 4, MPI_INT, lammps_->world);
131 if (error != MPI_SUCCESS) throw ATC_Error("error in sparse_allsum_numrows "+to_string(error));
132
133 // adjust row sendcounts because recRowsCRS is off by one
134 int rowCounts[nProcs];
135 int sizeCounts[nProcs];
136 // set up total size of receive buffers for Allgatherv calls
137 int totalRowsCRS = 0;
138 int totalSize = 0;
139 // set up array of displacements for Allgatherv calls
140 int rowOffsets[nProcs];
141 rowOffsets[0] = 0;
142 int sizeOffsets[nProcs];
143 sizeOffsets[0] = 0;
144 for (int i = 0; i < nProcs; i++) {
145 // find the total number of entries to share in the mpi calls below
146 rowCounts[i] = recInfo[i].rowsCRS + 1;
147 sizeCounts[i] = recInfo[i].size;
148 totalRowsCRS += rowCounts[i];
149 totalSize += recInfo[i].size;
150 // these already have their 0th slot filled in
151 if (i == 0) continue;
152 rowOffsets[i] = rowOffsets[i-1] + rowCounts[i-1];
153 sizeOffsets[i] = sizeOffsets[i-1] + sizeCounts[i-1];
154 }
155
156 // get actual rows
157 INDEX *rec_ia = new INDEX[totalRowsCRS];
158 if (toShare.size() == 0) {
159 double dummy[0];
160 error = MPI_Allgatherv(dummy, 0, MPI_INT,
161 rec_ia, rowCounts, rowOffsets, MPI_INT, lammps_->world);
162 }
163 else
164 error = MPI_Allgatherv(toShare.rows(), rowCounts[myRank], MPI_INT,
165 rec_ia, rowCounts, rowOffsets, MPI_INT, lammps_->world);
166 if (error != MPI_SUCCESS)
167 throw ATC_Error("error in sparse_allsum_rowarray "+to_string(error));
168
169 // get actual cols
170 INDEX *rec_ja = new INDEX[totalSize];
171 error = MPI_Allgatherv(toShare.cols(), sizeCounts[myRank], MPI_INT,
172 rec_ja, sizeCounts, sizeOffsets, MPI_INT, lammps_->world);
173 if (error != MPI_SUCCESS)
174 throw ATC_Error("error in sparse_allsum_colarray "+to_string(error));
175
176 // get the array of values
177 double *rec_vals = new double[totalSize];
178 error = MPI_Allgatherv(toShare.ptr(), sizeCounts[myRank], MPI_DOUBLE,
179 rec_vals, sizeCounts, sizeOffsets, MPI_DOUBLE, lammps_->world);
180 if (error != MPI_SUCCESS)
181 throw ATC_Error("error in sparse_allsum_valarray "+to_string(error));
182
183 INDEX *rec_ia_proc;
184 INDEX *rec_ja_proc;
185 double *rec_vals_proc;
186 for (int i = 0; i < nProcs; i++) {
187 if (myRank != i) {
188 // deallocated when tempMat is deleted since it wraps them
189 rec_ia_proc = new INDEX[rowCounts[i]];
190 rec_ja_proc = new INDEX[sizeCounts[i]];
191 rec_vals_proc = new double[sizeCounts[i]];
192
193 // copy the data passed with MPI into the new spots
194 copy(rec_ia + rowOffsets[i],
195 rec_ia + rowOffsets[i] + rowCounts[i],
196 rec_ia_proc);
197 copy(rec_ja + sizeOffsets[i],
198 rec_ja + sizeOffsets[i] + sizeCounts[i],
199 rec_ja_proc);
200 copy(rec_vals + sizeOffsets[i],
201 rec_vals + sizeOffsets[i] + sizeCounts[i],
202 rec_vals_proc);
203
204 // Does anyone know why we have to declare tempMat here (as well as set it equal to
205 // something) to avoid segfaults? there are still segfaults, but they happen at a much
206 // later stage of the game now (and for less benchmarks overall).
207 SparseMatrix<double> tempMat =
208 SparseMatrix<double>(rec_ia_proc, rec_ja_proc, rec_vals_proc,
209 recInfo[i].size, recInfo[i].rows,
210 recInfo[i].cols, recInfo[i].rowsCRS);
211 toShare += tempMat;
212 }
213 }
214
215 delete[] recInfo;
216 delete[] rec_ia;
217 delete[] rec_ja;
218 delete[] rec_vals;
219 }
220 #endif
221
222 // -----------------------------------------------------------------
223 // atom interface methods
224 // -----------------------------------------------------------------
fix_id() const225 string LammpsInterface::fix_id() const { return string(fixPointer_->id); }
226
nlocal() const227 int LammpsInterface::nlocal() const { return lammps_->atom->nlocal; }
228
nghost() const229 int LammpsInterface::nghost() const { return lammps_->atom->nghost; }
230
atoms_sorted() const231 bool LammpsInterface::atoms_sorted() const
232 {
233 int sortfreq = lammps_->atom->sortfreq;
234 if (sortfreq > 0) { return true; }
235 else { return false;
236 }
237 }
238
239
natoms() const240 bigint LammpsInterface::natoms() const { return lammps_->atom->natoms; }
241
nmax() const242 int LammpsInterface::nmax() const { return lammps_->atom->nmax; }
243
ntypes() const244 int LammpsInterface::ntypes() const { return lammps_->atom->ntypes; }
245
xatom() const246 double ** LammpsInterface::xatom() const { return lammps_->atom->x; }
247
type_to_charge(int atype) const248 int LammpsInterface::type_to_charge(int atype) const {
249 double *q = lammps_->atom->q;
250 if (! q) return 0;
251 int nlocal = lammps_->atom->nlocal;
252 int *type = lammps_->atom->type;
253 double aq = 0.0;
254 for (int i = 0; i < nlocal; i++) {
255 if (type[i] == atype) {
256 aq = q[i];
257 break;
258 }
259 }
260 double pcharge;
261 MPI_Allreduce(&aq,&pcharge,1,MPI_DOUBLE,MPI_MAX,world());
262 double ncharge;
263 MPI_Allreduce(&aq,&ncharge,1,MPI_DOUBLE,MPI_MIN,world());
264 double charge = (pcharge == 0.0) ? ncharge : pcharge;
265 return charge;
266 }
267
268 //const double ** LammpsInterface::xatom() const { return (const double**)(lammps_->atom->x); }
269
vatom() const270 double ** LammpsInterface::vatom() const { return lammps_->atom->v; }
271
fatom() const272 double ** LammpsInterface::fatom() const { return lammps_->atom->f; }
273
atom_mask() const274 const int * LammpsInterface::atom_mask() const { return (const int*)lammps_->atom->mask; }
275
atom_mask()276 int * LammpsInterface::atom_mask() { return lammps_->atom->mask; }
277
atom_type() const278 int * LammpsInterface::atom_type() const { return lammps_->atom->type; }
279
atom_tag() const280 int * LammpsInterface::atom_tag() const { return lammps_->atom->tag; }
281
atom_to_molecule() const282 int * LammpsInterface::atom_to_molecule() const { return lammps_->atom->molecule; }
283
284
num_bond() const285 int * LammpsInterface::num_bond() const { return lammps_->atom->num_bond; }
286
bond_atom() const287 int ** LammpsInterface::bond_atom() const { return lammps_->atom->bond_atom; }
288
image() const289 int * LammpsInterface::image() const { return lammps_->atom->image; }
290
bond_per_atom() const291 int LammpsInterface::bond_per_atom() const { return lammps_->atom->bond_per_atom; }
292
newton_bond() const293 int LammpsInterface::newton_bond() const { return lammps_->force->newton_bond; }
294
local_to_global_map(int global) const295 int LammpsInterface::local_to_global_map(int global) const { return lammps_->atom->map(global); }
296
atom_mass() const297 double * LammpsInterface::atom_mass() const { return lammps_->atom->mass; }
298
atom_mass(int iType) const299 double LammpsInterface::atom_mass(int iType) const { return lammps_->atom->mass[iType]; }
300
atom_rmass() const301 double * LammpsInterface::atom_rmass() const { return lammps_->atom->rmass; }
302
atom_charge() const303 double * LammpsInterface::atom_charge() const { return lammps_->atom->q; }
304
atom_scalar(FundamentalAtomQuantity quantityType) const305 double * LammpsInterface::atom_scalar(FundamentalAtomQuantity quantityType) const
306 {
307 if (quantityType==ATOM_MASS) {
308 if (atom_mass())
309 throw ATC_Error("Atom mass array requested but not defined");
310 return atom_rmass();
311 }
312 else if (quantityType==ATOM_CHARGE) {
313 double * atomCharge = atom_charge();
314 if (!atomCharge)
315 throw ATC_Error("Atom charge array requested but not defined");
316 return atomCharge;
317 }
318 else
319 throw ATC_Error("BAD type requested in atom_scalar");
320 return NULL;
321 }
322
atom_vector(FundamentalAtomQuantity quantityType) const323 double ** LammpsInterface::atom_vector(FundamentalAtomQuantity quantityType) const
324 {
325 if (quantityType==ATOM_POSITION)
326 return xatom();
327 else if (quantityType==ATOM_VELOCITY)
328 return vatom();
329 else if (quantityType==ATOM_FORCE)
330 return fatom();
331 else
332 throw ATC_Error("BAD type requested in atom_vector");
333 return NULL;
334 }
335
atom_quantity_ndof(FundamentalAtomQuantity quantityType) const336 int LammpsInterface::atom_quantity_ndof(FundamentalAtomQuantity quantityType) const
337 {
338 if (quantityType==ATOM_MASS || quantityType==ATOM_CHARGE)
339 return 1;
340 else if (quantityType==ATOM_POSITION || quantityType==ATOM_VELOCITY || quantityType==ATOM_FORCE)
341 return dimension();
342 else
343 throw ATC_Error("BAD type requested in atom_quantity_ndof");
344 }
345
atom_quantity_conversion(FundamentalAtomQuantity quantityType) const346 double LammpsInterface::atom_quantity_conversion(FundamentalAtomQuantity quantityType) const
347 {
348 if (quantityType==ATOM_MASS || quantityType==ATOM_CHARGE || quantityType==ATOM_POSITION || quantityType==ATOM_VELOCITY)
349 return 1;
350 else if ( quantityType==ATOM_FORCE)
351 return ftm2v();
352 else
353 throw ATC_Error("BAD type requested in atom_quantity_conversion");
354 }
355
356 // -----------------------------------------------------------------
357 // domain interface methods
358 // -----------------------------------------------------------------
359
dimension() const360 int LammpsInterface::dimension() const { return lammps_->domain->dimension; }
361
nregion() const362 int LammpsInterface::nregion() const { return lammps_->domain->nregion; }
363
box_bounds(double & boxxlo,double & boxxhi,double & boxylo,double & boxyhi,double & boxzlo,double & boxzhi) const364 void LammpsInterface::box_bounds(double & boxxlo, double & boxxhi,
365 double & boxylo, double & boxyhi,
366 double & boxzlo, double &boxzhi) const
367 {
368 if (lammps_->domain->triclinic == 0) {
369 boxxlo = lammps_->domain->boxlo[0];
370 boxxhi = lammps_->domain->boxhi[0];
371 boxylo = lammps_->domain->boxlo[1];
372 boxyhi = lammps_->domain->boxhi[1];
373 boxzlo = lammps_->domain->boxlo[2];
374 boxzhi = lammps_->domain->boxhi[2];
375 }
376 else {
377 boxxlo = lammps_->domain->boxlo_bound[0];
378 boxxhi = lammps_->domain->boxhi_bound[0];
379 boxylo = lammps_->domain->boxlo_bound[1];
380 boxyhi = lammps_->domain->boxhi_bound[1];
381 boxzlo = lammps_->domain->boxlo_bound[2];
382 boxzhi = lammps_->domain->boxhi_bound[2];
383 }
384 }
385
in_box(double * x) const386 bool LammpsInterface::in_box(double * x) const
387 {
388 double xlo,xhi,ylo,yhi,zlo,zhi;
389 box_bounds(xlo,xhi,ylo,yhi,zlo,zhi);
390 if (x[0] >= xlo && x[0] < xhi &&
391 x[1] >= ylo && x[1] < yhi &&
392 x[2] >= zlo && x[2] < zhi)
393 return true;
394 return false;
395 }
396
in_my_processor_box(double * x) const397 bool LammpsInterface::in_my_processor_box(double * x) const
398 {
399 if (x[0] >= lammps_->domain->sublo[0] && x[0] < lammps_->domain->subhi[0] &&
400 x[1] >= lammps_->domain->sublo[1] && x[1] < lammps_->domain->subhi[1] &&
401 x[2] >= lammps_->domain->sublo[2] && x[2] < lammps_->domain->subhi[2])
402 return true;
403 if (! in_box(x))
404 throw ATC_Error("point is in no processors box");
405 return false;
406 }
407
408
sub_bounds(double & subxlo,double & subxhi,double & subylo,double & subyhi,double & subzlo,double & subzhi) const409 void LammpsInterface::sub_bounds(double & subxlo, double & subxhi,
410 double & subylo, double & subyhi,
411 double & subzlo, double & subzhi) const
412 {
413 if (lammps_->domain->triclinic == 0) {
414 subxlo = lammps_->domain->sublo[0];
415 subxhi = lammps_->domain->subhi[0];
416 subylo = lammps_->domain->sublo[1];
417 subyhi = lammps_->domain->subhi[1];
418 subzlo = lammps_->domain->sublo[2];
419 subzhi = lammps_->domain->subhi[2];
420 }
421 else {
422 ATC_Error("Subboxes not accurate when triclinic != 0.");
423 }
424 }
425
xperiodic() const426 int LammpsInterface::xperiodic() const { return lammps_->domain->xperiodic; }
427
yperiodic() const428 int LammpsInterface::yperiodic() const { return lammps_->domain->yperiodic; }
429
zperiodic() const430 int LammpsInterface::zperiodic() const { return lammps_->domain->zperiodic; }
431
nperiodic() const432 int LammpsInterface::nperiodic() const
433 {
434 int nprd = 0;
435 if ( lammps_->domain->xperiodic > 0 ) { nprd++ ; }
436 if ( lammps_->domain->yperiodic > 0 ) { nprd++ ; }
437 if ( lammps_->domain->zperiodic > 0 ) { nprd++ ; }
438 return nprd;
439 }
440
441 // correct posistions for periodic box
periodicity_correction(double * x) const442 void LammpsInterface::periodicity_correction(double * x) const
443 {
444 int* periodicity = lammps_->domain->periodicity;
445 if (!refBoxIsSet_) set_reference_box();
446 for (int m = 0; m < 3; m++) {
447 if ((bool) periodicity[m]) {
448 if (x[m] < lower_[m] || x[m] > upper_[m]) {
449 x[m] -= length_[m]*floor((x[m]-lower_[m])/length_[m]);
450 }
451 if (x[m] < lower_[m] || x[m] > upper_[m]) {
452 throw ATC_Error("periodicity_correction: still out of box bounds");
453 }
454 }
455 }
456 }
457
set_reference_box(void) const458 void LammpsInterface::set_reference_box(void) const
459 {
460 double * hi = lammps_->domain->boxhi;
461 double * lo = lammps_->domain->boxlo;
462 double * len = lammps_->domain->prd;
463 for (int i = 0; i < 3; i++) {
464 upper_[i] = hi[i];
465 lower_[i] = lo[i];
466 length_[i] = len[i];
467 }
468 refBoxIsSet_ = true;
469 }
470
471
domain_xprd() const472 double LammpsInterface::domain_xprd() const { return lammps_->domain->xprd; }
473
domain_yprd() const474 double LammpsInterface::domain_yprd() const { return lammps_->domain->yprd; }
475
domain_zprd() const476 double LammpsInterface::domain_zprd() const { return lammps_->domain->zprd; }
477
domain_volume() const478 double LammpsInterface::domain_volume() const
479 {
480 return (lammps_->domain->xprd)*
481 (lammps_->domain->yprd)*
482 (lammps_->domain->zprd);
483 }
484
domain_xy() const485 double LammpsInterface::domain_xy() const { return lammps_->domain->xy; }
486
domain_xz() const487 double LammpsInterface::domain_xz() const { return lammps_->domain->xz; }
488
domain_yz() const489 double LammpsInterface::domain_yz() const { return lammps_->domain->yz; }
490
domain_triclinic() const491 int LammpsInterface::domain_triclinic() const { return lammps_->domain->triclinic; }
492
box_periodicity(int & xperiodic,int & yperiodic,int & zperiodic) const493 void LammpsInterface::box_periodicity(int & xperiodic,
494 int & yperiodic,
495 int & zperiodic) const
496 {
497 xperiodic = lammps_->domain->xperiodic;
498 yperiodic = lammps_->domain->yperiodic;
499 zperiodic = lammps_->domain->zperiodic;
500 }
501
region_id(const char * regionName) const502 int LammpsInterface::region_id(const char * regionName) const {
503 int nregion = this->nregion();
504 for (int iregion = 0; iregion < nregion; iregion++) {
505 if (strcmp(regionName, region_name(iregion)) == 0) {
506 return iregion;
507 }
508 }
509 throw ATC_Error("Region has not been defined");
510 return -1;
511 }
512
region_bounds(const char * regionName,double & xmin,double & xmax,double & ymin,double & ymax,double & zmin,double & zmax,double & xscale,double & yscale,double & zscale) const513 bool LammpsInterface::region_bounds(const char * regionName,
514 double & xmin, double & xmax,
515 double & ymin, double & ymax,
516 double & zmin, double & zmax,
517 double & xscale, double & yscale, double & zscale) const
518 {
519 int iRegion = region_id(regionName);
520 xscale = region_xscale(iRegion);
521 yscale = region_yscale(iRegion);
522 zscale = region_zscale(iRegion);
523 xmin = region_xlo(iRegion);
524 xmax = region_xhi(iRegion);
525 ymin = region_ylo(iRegion);
526 ymax = region_yhi(iRegion);
527 zmin = region_zlo(iRegion);
528 zmax = region_zhi(iRegion);
529 if (strcmp(region_style(iRegion),"block")==0) { return true; }
530 else { return false; }
531 }
532
minimum_image(double & dx,double & dy,double & dz) const533 void LammpsInterface::minimum_image(double & dx, double & dy, double & dz) const {
534 lammps_->domain->minimum_image(dx,dy,dz);
535 }
536
closest_image(const double * const xi,const double * const xj,double * const xjImage) const537 void LammpsInterface::closest_image(const double * const xi, const double * const xj, double * const xjImage) const {
538 lammps_->domain->closest_image(xi,xj,xjImage);
539 }
540
541
542 // -----------------------------------------------------------------
543 // update interface methods
544 // -----------------------------------------------------------------
units_style(void) const545 LammpsInterface::UnitsType LammpsInterface::units_style(void) const
546 {
547 if (strcmp(lammps_->update->unit_style,"lj") == 0) return LJ;
548 else if (strcmp(lammps_->update->unit_style,"real") == 0) return REAL;
549 else if (strcmp(lammps_->update->unit_style,"metal") == 0) return METAL;
550 else return UNKNOWN;
551 }
552
convert_units(double value,UnitsType in,UnitsType out,int massExp,int lenExp,int timeExp,int engExp) const553 double LammpsInterface::convert_units(double value, UnitsType in, UnitsType out, int massExp, int lenExp, int timeExp, int engExp) const
554 {
555 double ps2fs = 1.e3;
556 double eV2kcal = 23.069;
557 if (in==REAL) {
558 if (out==METAL) {
559 return value*pow(ps2fs,-timeExp)*pow(eV2kcal,-engExp);
560 }
561 else if (out==ATC) {
562 if (units_style()==REAL) {
563 return value;
564 }
565 else if (units_style()==METAL) {
566 return convert_units(value, METAL, out, massExp, lenExp, timeExp)*1.0;
567 }
568 }
569 else throw ATC_Error("can't convert");
570 }
571 else if (in==METAL) {
572 if (out==REAL) {
573 return value*pow(ps2fs,timeExp)*pow(eV2kcal,engExp);
574 }
575 else if (out==ATC) {
576 if (units_style()==REAL) {
577 return convert_units(value, REAL, out, massExp, lenExp, timeExp)*1.0;
578 }
579 else if (units_style()==METAL) {
580 return value;
581 }
582 }
583 else throw ATC_Error("can't convert");
584 }
585 else throw ATC_Error("can't convert");
586 return value;
587 }
588
589
590 // -----------------------------------------------------------------
591 // lattice interface methods
592 // -----------------------------------------------------------------
593
xlattice() const594 double LammpsInterface::xlattice() const { return lammps_->domain->lattice->xlattice; }
595
ylattice() const596 double LammpsInterface::ylattice() const { return lammps_->domain->lattice->ylattice; }
597
zlattice() const598 double LammpsInterface::zlattice() const { return lammps_->domain->lattice->zlattice; }
599
lattice_style() const600 LammpsInterface::LatticeType LammpsInterface::lattice_style() const
601 {
602 if (lammps_->domain->lattice)
603 return (LammpsInterface::LatticeType)lammps_->domain->lattice->style;
604 else
605 throw ATC_Error("Lattice has not been defined");
606 }
607
608 //* retuns the number of basis vectors
n_basis() const609 int LammpsInterface::n_basis() const
610 {
611 return lammps_->domain->lattice->nbasis;
612 }
613
614 //* returns the basis vectors, transformed to the box coords
basis_vectors(double ** basis) const615 void LammpsInterface::basis_vectors(double **basis) const
616 {
617 LAMMPS_NS::Lattice *lattice = lammps_->domain->lattice;
618 int i,j;
619 double origin[3] = {0.0, 0.0, 0.0};
620 lattice->lattice2box(origin[0], origin[1], origin[2]);
621 for (i=0; i<n_basis(); i++)
622 {
623 memcpy(basis[i],lattice->basis[i],3*sizeof(double));
624 lattice->lattice2box(basis[i][0], basis[i][1], basis[i][2]);
625 for (j=0; j<3; j++) basis[i][j] -= origin[j];
626 }
627 }
628
629 //* gets the (max) lattice constant
max_lattice_constant(void) const630 double LammpsInterface::max_lattice_constant(void) const
631 {
632 double a1[3], a2[3], a3[3];
633 unit_cell(a1,a2,a3);
634 double a = norm(a1);
635 a = max(a,norm(a2));
636 a = max(a,norm(a3));
637 return a;
638 }
639
640 //* computes a cutoff distance halfway between 1st and 2nd nearest neighbors
near_neighbor_cutoff(void) const641 double LammpsInterface::near_neighbor_cutoff(void) const
642 {
643 double cutoff;
644 double alat = LammpsInterface::max_lattice_constant();
645 LatticeType type = lattice_style();
646 if (type == LammpsInterface::SC) {
647 cutoff = 0.5*(1.0+sqrt(2.0))*alat;
648 } else if (type == LammpsInterface::BCC) {
649 cutoff = 0.5*(0.5*sqrt(3.0)+1.0)*alat;
650 } else if (type == LammpsInterface::FCC) {
651 cutoff = 0.5*(1.0/sqrt(2.0)+1.0)*alat;
652 } else if (type == LammpsInterface::HCP) {
653 cutoff = 0.5*(1.0/sqrt(2.0)+1.0)*alat;
654 } else if (type == LammpsInterface::DIAMOND) {
655 cutoff = 0.5*(0.25*sqrt(3.0)+1.0/sqrt(2.0))*alat;
656 } else if (type == LammpsInterface::SQ) {
657 cutoff = 0.5*(1.0+sqrt(2.0))*alat;
658 } else if (type == LammpsInterface::SQ2) {
659 cutoff = 0.5*(1.0/sqrt(2.0)+1.0)*alat;
660 } else if (type == LammpsInterface::HEX) {
661 cutoff = 0.5*(1.0/sqrt(3.0)+1.0)*alat;
662 } else {
663 throw ATC_Error("Unknown lattice type");
664 }
665 return cutoff;
666 }
667
668 //* gets the unit cell vectors
unit_cell(double * a1,double * a2,double * a3) const669 void LammpsInterface::unit_cell(double *a1, double *a2, double *a3) const
670 {
671 int i, j;
672 double *a[3] = {a1,a2,a3};
673 double origin[3] = {0.0,0.0,0.0};
674 LAMMPS_NS::Lattice *lattice = lammps_->domain->lattice;
675 // transform origin
676 lattice->lattice2box(origin[0], origin[1], origin[2]);
677
678 // copy reference lattice vectors
679 memcpy(a[0], lattice->a1, 3*sizeof(double));
680 memcpy(a[1], lattice->a2, 3*sizeof(double));
681 memcpy(a[2], lattice->a3, 3*sizeof(double));
682
683 for (i=0; i<3; i++)
684 {
685 lattice->lattice2box(a[i][0], a[i][1], a[i][2]);
686 for (j=0; j<3; j++) a[i][j] -= origin[j];
687 }
688 }
689
690 //* gets number of atoms in a unit cell
num_atoms_per_cell(void) const691 int LammpsInterface::num_atoms_per_cell(void) const
692 {
693 int naCell = 0;
694 LatticeType type = lattice_style();
695 if (type == LammpsInterface::SC) naCell = 1;
696 else if (type == LammpsInterface::BCC) naCell = 2;
697 else if (type == LammpsInterface::FCC) naCell = 4;
698 else if (type == LammpsInterface::DIAMOND) naCell = 8;
699 else if (comm_rank()==0) {
700 //{throw ATC_Error("lattice style not currently supported by ATC");}
701 print_msg_once("WARNING: Cannot get number of atoms per cell from lattice");
702 naCell = 1;
703 }
704 return naCell;
705 }
706
707 //* gets tributary volume for an atom
volume_per_atom(void) const708 double LammpsInterface::volume_per_atom(void) const
709 {
710 double naCell = num_atoms_per_cell();
711 double volPerAtom =
712 xlattice() * ylattice() * zlattice() / naCell;
713 return volPerAtom;
714 }
715
716 //* gets lattice basis
lattice(MATRIX & N,MATRIX & B) const717 void LammpsInterface::lattice(MATRIX &N, MATRIX &B) const
718 {
719 int nbasis = n_basis();
720 double **basis = new double*[nbasis];
721 N.reset(3,3);
722 B.reset(3,nbasis);
723 for (int i=0; i<nbasis; i++) basis[i] = column(B,i).ptr();
724 basis_vectors(basis);
725 unit_cell(column(N,0).ptr(),
726 column(N,1).ptr(),
727 column(N,2).ptr());
728 delete [] basis;
729 }
730
731 // -----------------------------------------------------------------
732 // force interface methods
733 // -----------------------------------------------------------------
734
boltz() const735 double LammpsInterface::boltz() const{ return lammps_->force->boltz; }
736
mvv2e() const737 double LammpsInterface::mvv2e() const{ return lammps_->force->mvv2e; }
738
ftm2v() const739 double LammpsInterface::ftm2v()const { return lammps_->force->ftm2v; }
740
nktv2p() const741 double LammpsInterface::nktv2p()const{ return lammps_->force->nktv2p; }
742
qqr2e() const743 double LammpsInterface::qqr2e() const{ return lammps_->force->qqr2e; }
744
qe2f() const745 double LammpsInterface::qe2f() const{ return lammps_->force->qe2f; }
dielectric() const746 double LammpsInterface::dielectric()const{return lammps_->force->dielectric; }
747
qqrd2e() const748 double LammpsInterface::qqrd2e()const{ return lammps_->force->qqrd2e; }
749
qv2e() const750 double LammpsInterface::qv2e() const{ return qe2f()*ftm2v(); }
751
pair_force(int i,int j,double rsq,double & fmag_over_rmag) const752 double LammpsInterface::pair_force(int i, int j, double rsq,
753 double & fmag_over_rmag) const
754 {
755 int itype = (lammps_->atom->type)[i];
756 int jtype = (lammps_->atom->type)[j];
757 // return value is the energy
758 if (rsq < (lammps_->force->pair->cutsq)[itype][jtype]) {
759 return lammps_->force->pair->single(i,j,itype,jtype,
760 rsq,1.0,1.0,fmag_over_rmag);
761 }
762 return 0.0;
763 }
pair_force(int n,double rsq,double & fmag_over_rmag) const764 double LammpsInterface::pair_force(int n, double rsq,
765 double & fmag_over_rmag) const
766 {
767 int i = bond_list_i(n);
768 int j = bond_list_j(n);
769 int type = bond_list_type(n);
770 // return value is the energy
771 return lammps_->force->bond->single(type,rsq,i,j,fmag_over_rmag);
772 }
pair_force(map<std::pair<int,int>,int>::const_iterator itr,double rsq,double & fmag_over_rmag,int nbonds) const773 double LammpsInterface::pair_force(
774 map< std::pair< int,int >,int >::const_iterator itr, double rsq,
775 double & fmag_over_rmag, int nbonds) const
776 {
777 int n = itr->second;
778 if (n < nbonds) {
779 return pair_force(n, rsq,fmag_over_rmag);
780 }
781 else {
782 std::pair <int,int> ij = itr->first;
783 int i = ij.first;
784 int j = ij.second;
785 return pair_force(i,j, rsq,fmag_over_rmag);
786 }
787 }
pair_force(std::pair<std::pair<int,int>,int> apair,double rsq,double & fmag_over_rmag,int nbonds) const788 double LammpsInterface::pair_force(
789 std::pair< std::pair< int,int >,int > apair, double rsq,
790 double & fmag_over_rmag, int nbonds) const
791 {
792 int n = apair.second;
793 if (n < nbonds) {
794 return pair_force(n, rsq,fmag_over_rmag);
795 }
796 else {
797 std::pair <int,int> ij = apair.first;
798 int i = ij.first;
799 int j = ij.second;
800 return pair_force(i,j, rsq,fmag_over_rmag);
801 }
802 }
bond_stiffness(int i,int j,double rsq0) const803 double LammpsInterface::bond_stiffness(int i, int j, double rsq0) const
804 {
805 const double perturbation = 1.e-8;
806 double rsq1 = sqrt(rsq0)+perturbation;
807 rsq1 *= rsq1;
808 double f0,f1;
809 pair_force(i,j,rsq0,f0);
810 pair_force(i,j,rsq1,f1);
811 double k = (f1-f0)/perturbation;
812 return k;
813 }
814
pair_cutoff() const815 double LammpsInterface::pair_cutoff() const
816 {
817 return lammps_->force->pair->cutforce;
818 }
819
pair_reinit() const820 void LammpsInterface::pair_reinit() const
821 {
822 lammps_->force->pair->reinit();
823 }
824
single_enable() const825 int LammpsInterface::single_enable() const
826 {
827 return lammps_->force->pair->single_enable; // all bonds have a single
828 }
829
830 //* insertion/deletion functions : see FixGCMC
831 //* delete atom
delete_atom(int id) const832 int LammpsInterface::delete_atom(int id) const
833 {
834 LAMMPS_NS::Atom * atom = lammps_->atom;
835 atom->avec->copy(atom->nlocal-1,id,1);
836 atom->nlocal--;
837 return atom->nlocal;
838 }
839
840 //* insert atom
insert_atom(int atype,int amask,double * ax,double * av,double aq) const841 int LammpsInterface::insert_atom(int atype, int amask,
842 double *ax, double *av, double aq) const
843 {
844 LAMMPS_NS::Atom * atom = lammps_->atom;
845 atom->avec->create_atom(atype,ax);
846 int m = atom->nlocal - 1;
847 atom->mask[m] = amask;
848 atom->v[m][0] = av[0];
849 atom->v[m][1] = av[1];
850 atom->v[m][2] = av[2];
851 if (aq != 0) atom->q[m] = aq;
852
853 int nfix = lammps_->modify->nfix;
854 LAMMPS_NS::Fix **fix = lammps_->modify->fix;
855 for (int j = 0; j < nfix; j++) {
856 if (fix[j]->create_attribute) fix[j]->set_arrays(m);
857 }
858 return m;
859 }
860
reset_ghosts(int deln) const861 int LammpsInterface::reset_ghosts(int deln) const
862 {
863 LAMMPS_NS::Atom * atom = lammps_->atom;
864 //ATC::LammpsInterface::instance()->print_msg("reset_ghosts del n "+to_string(deln));
865 if (atom->tag_enable) {
866 atom->natoms += deln;
867 //ATC::LammpsInterface::instance()->print_msg("reset_ghosts natoms "+to_string(atom->natoms));
868 if (deln > 0) { atom->tag_extend(); }
869 if (atom->map_style) atom->map_init();
870 }
871 atom->nghost = 0;
872 lammps_->comm->borders();
873 //ATC::LammpsInterface::instance()->print_msg("reset_ghosts nghosts "+to_string(atom->nghost));
874 return atom->nghost;
875 }
876
877 //* energy for interactions within the shortrange cutoff
shortrange_energy(double * coord,int itype,int id,double max) const878 double LammpsInterface::shortrange_energy(double *coord,
879 int itype, int id, double max) const
880 {
881 LAMMPS_NS::Atom * atom = lammps_->atom;
882 double **x = atom->x;
883 int *type = atom->type;
884 int nall = atom->nlocal+ atom->nghost;
885 LAMMPS_NS::Pair *pair = lammps_->force->pair;
886 double **cutsq = lammps_->force->pair->cutsq;
887
888 double fpair = 0.0; // an output of single
889 double factor_coul = 1.0;
890 double factor_lj = 1.0;
891
892 double total_energy = 0.0;
893 for (int j = 0; j < nall; j++) {
894 if (id == j) continue;
895 // factor_lj = special_lj[sbmask(j)];
896 // factor_coul = special_coul[sbmask(j)];
897 //j &= NEIGHMASK;
898 double delx = coord[0] - x[j][0];
899 double dely = coord[1] - x[j][1];
900 double delz = coord[2] - x[j][2];
901 double rsq = delx*delx + dely*dely + delz*delz;
902 int jtype = type[j];
903 double cut2 = cutsq[itype][jtype];
904 if (rsq < cut2) {
905 total_energy += pair->single(id,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); //id is for charge lookup
906 }
907 }
908 return total_energy;
909 }
shortrange_energy(int id,double max) const910 double LammpsInterface::shortrange_energy(int id, double max) const
911 {
912 double *x = (lammps_->atom->x)[id];
913 int type = (lammps_->atom->type)[id];
914 return shortrange_energy(x,type,id,max);
915 }
916
potential() const917 POTENTIAL LammpsInterface::potential() const
918 {
919 // find pair style - FRAGILE
920 const int nStyles = 4;
921 string pairStyles[nStyles] = {"lj/cut",
922 "lj/cut/coul/long",
923 "lj/cut/coul/cut",
924 "lj/charmm/coul/long"};
925 LAMMPS_NS::Pair *pair = NULL;
926 for (int i = 0; i < nStyles; i++){
927 pair = lammps_->force->pair_match(pairStyles[i].c_str(),1);
928 if (pair != NULL) break;
929 }
930 return pair;
931 }
932
type_to_groupbit(int itype) const933 int LammpsInterface::type_to_groupbit(int itype) const
934 {
935 LAMMPS_NS::Atom * atom = lammps_->atom;
936 int groupbit = -1;
937 int *type = atom->type;
938 int *mask = atom->mask;
939 for (int i = 0; i < nlocal(); i++) {
940 if (type[i] == itype) {
941 groupbit = mask[i];
942 break;
943 }
944 }
945 return int_allmax(groupbit);
946 }
947
epsilons(int itype,POTENTIAL pair,double * epsilon0) const948 bool LammpsInterface::epsilons(int itype, POTENTIAL pair, double * epsilon0) const
949 {
950 // grab energy parameters
951 char * pair_parameter = new char[8];
952 strcpy(pair_parameter,"epsilon");
953 int dim = 2; // a return value for extract
954 double ** epsilons = (double**) ( pair->extract(pair_parameter,dim) );
955 delete [] pair_parameter;
956 if (epsilons == NULL) return false;
957 //if (epsilons == NULL) error->all(FLERR,"Fix concentration adapted pair style parameter not supported");
958 int i1,i2;
959 for (int i=1; i < ntypes()+1; i++) {
960 if (i < itype) { i1 = i; i2 = itype; }
961 else { i2 = i; i1 = itype; }
962 epsilon0[i-1] = epsilons[i1][i2];
963 }
964 return true;
965 }
966
set_epsilons(int itype,POTENTIAL pair,double * epsilon) const967 bool LammpsInterface::set_epsilons(int itype, POTENTIAL pair, double * epsilon) const
968 {
969 // grab energy parameters
970 char * pair_parameter = new char[8];
971 strcpy(pair_parameter,"epsilon");
972 int dim = 2; // a return value for extract
973 double ** epsilons = (double**) ( pair->extract(pair_parameter,dim) );
974 delete [] pair_parameter;
975 if (epsilons == NULL) return false;
976 //if (epsilons == NULL) error->all(FLERR,"Fix concentration adapted pair style parameter not supported");
977 // scale interactions
978 int i1,i2;
979 for (int i = 1; i < ntypes()+1; i++) {
980 if (i < itype) { i1 = i; i2 = itype; }
981 else { i2 = i; i1 = itype; }
982 epsilons[i1][i2] = epsilon[i-1];
983 }
984
985 return true;
986 }
987
set_charge(int itype,double charge) const988 int LammpsInterface::set_charge(int itype, double charge) const
989 {
990 int nlocal = lammps_->atom->nlocal;
991 int *type = lammps_->atom->type;
992 double *q = lammps_->atom->q;
993 int count = 0;
994 for (int i = 0; i < nlocal; i++) {
995 if (type[i] == itype) {
996 q[i] = charge;
997 count++;
998 }
999 }
1000 return count;
1001 }
1002
change_type(int itype,int jtype) const1003 int LammpsInterface::change_type(int itype, int jtype) const
1004 {
1005 int nlocal = lammps_->atom->nlocal;
1006 int *type = lammps_->atom->type;
1007 int count = 0;
1008 for (int i = 0; i < nlocal; i++) {
1009 if (type[i] == itype) {
1010 type[i] = jtype;
1011 count++;
1012 }
1013 }
1014 return count;
1015 }
1016
count_type(int itype) const1017 int LammpsInterface::count_type(int itype) const
1018 {
1019 int nlocal = lammps_->atom->nlocal;
1020 int *type = lammps_->atom->type;
1021 int count = 0;
1022 for (int i = 0; i < nlocal; i++) {
1023 if (type[i] == itype) { count++; }
1024 }
1025 return int_allsum(count);
1026 }
1027
1028 // random number generators
random_number_generator() const1029 RNG_POINTER LammpsInterface::random_number_generator() const {
1030 RNG_POINTER p = new LAMMPS_NS::RanPark(lammps_,seed_);
1031 return p;
1032 }
random_uniform(RNG_POINTER p) const1033 double LammpsInterface::random_uniform(RNG_POINTER p) const {
1034 return p->uniform();
1035 }
random_normal(RNG_POINTER p) const1036 double LammpsInterface::random_normal (RNG_POINTER p) const {
1037 return p->gaussian();
1038 }
random_state(RNG_POINTER p) const1039 int LammpsInterface::random_state (RNG_POINTER p) const {
1040 return p->state();
1041 }
set_random_state(RNG_POINTER p,int seed) const1042 void LammpsInterface::set_random_state (RNG_POINTER p, int seed) const {
1043 return p->reset(seed);
1044 }
advance_random_generator(RNG_POINTER p,int n) const1045 void LammpsInterface::advance_random_generator (RNG_POINTER p, int n) const {
1046 advance_random_uniform(p,n);
1047 }
advance_random_uniform(RNG_POINTER p,int n) const1048 void LammpsInterface::advance_random_uniform (RNG_POINTER p, int n) const {
1049 for (int i = 0; i < n; i++) p->uniform();
1050 }
advance_random_normal(RNG_POINTER p,int n) const1051 void LammpsInterface::advance_random_normal (RNG_POINTER p, int n) const {
1052 for (int i = 0; i < n; i++) p->gaussian();
1053 }
1054
1055 //* Boltzmann's constant in M,L,T,t units
kBoltzmann() const1056 double LammpsInterface::kBoltzmann() const {
1057 return (lammps_->force->boltz)/(lammps_->force->mvv2e);
1058 }
1059
1060 //* Planck's constant
hbar() const1061 double LammpsInterface::hbar() const {
1062 const int UNITS_STYLE = (int) units_style();
1063 double hbar = 1.0; // LJ: Dimensionless
1064 if (UNITS_STYLE == 2) hbar = 15.1685792814; // Real: KCal/mol-fs
1065 else if (UNITS_STYLE == 3) hbar = 0.000658212202469; // Metal: eV-ps
1066 return hbar;
1067 }
1068
1069 //* Dulong-Petit heat capacity
heat_capacity() const1070 double LammpsInterface::heat_capacity() const {
1071 double rhoCp = dimension()*kBoltzmann()/volume_per_atom();
1072 return rhoCp;
1073 }
1074
1075 //* reference mass density for a *unit cell*
1076 // all that is needed is a unit cell: volume, types, mass per type
mass_density(int * numPerType) const1077 double LammpsInterface::mass_density(int* numPerType) const
1078 {
1079 const double *mass = lammps_->atom->mass;
1080 if (!mass) throw ATC_Error("cannot compute a mass density: no mass");
1081 const int ntypes = lammps_->atom->ntypes;
1082 const int *mass_setflag = lammps_->atom->mass_setflag;
1083 const int *type = lammps_->atom->type;
1084 double naCell = num_atoms_per_cell();
1085 double vol = volume_per_atom();
1086 if (numPerType) {
1087 double m = 0.0;
1088 double n = 0;
1089 for (int i = 0; i < ntypes; i++) {
1090 m += numPerType[i]*mass[i+1];
1091 n += numPerType[i];
1092 }
1093 if (n>naCell) throw ATC_Error("cannot compute a mass density: too many basis atoms");
1094 return m/n/vol;
1095 }
1096 // since basis->type map not stored only monatomic lattices are automatic
1097 // if not given a basis try to guess
1098 else {
1099 if (ntypes == 1) {
1100 if ((this->natoms()==0) && mass_setflag[1]) {
1101 return mass[1]/vol;
1102 }
1103 else {
1104 if (type) return mass[type[0]]/vol;
1105 else if (mass_setflag[1]) return mass[1]/vol;
1106 }
1107 }
1108 }
1109 throw ATC_Error("cannot compute a mass density");
1110 return 0.0;
1111 }
1112
1113 //* permittivity of free space
epsilon0() const1114 double LammpsInterface::epsilon0() const
1115 {
1116 return qe2f()/(4.*PI*qqr2e());
1117 }
1118
1119 //* Coulomb's constant
coulomb_constant() const1120 double LammpsInterface::coulomb_constant() const
1121 {
1122 return qqr2e()/qe2f();
1123 }
1124
1125 //* special coulombic interactions
special_coul() const1126 double * LammpsInterface::special_coul() const
1127 {
1128 return lammps_->force->special_coul;
1129 }
1130
1131 //* flag for newton
newton_pair() const1132 int LammpsInterface::newton_pair() const
1133 {
1134 return lammps_->force->newton_pair;
1135 }
1136
1137 // -----------------------------------------------------------------
1138 // group interface methods
1139 // -----------------------------------------------------------------
1140
ngroup() const1141 int LammpsInterface::ngroup() const { return lammps_->group->ngroup; }
1142
group_bit(string name) const1143 int LammpsInterface::group_bit(string name) const
1144 {
1145 return group_bit(group_index(name));
1146 }
1147
group_bit(int iGroup) const1148 int LammpsInterface::group_bit(int iGroup) const
1149 {
1150 int mybit = 0;
1151 mybit |= lammps_->group->bitmask[iGroup];
1152 if (mybit < 0 || mybit > MAX_GROUP_BIT) {
1153 string msg("LammpsInterface::group_bit() lammps group bit "+to_string(mybit)+" is out of range 0:"+to_string(MAX_GROUP_BIT));
1154 throw ATC_Error(msg);
1155 }
1156
1157 return mybit;
1158 }
1159
group_index(string name) const1160 int LammpsInterface::group_index(string name) const
1161 {
1162 int igroup = lammps_->group->find(name.c_str());
1163 if (igroup == -1) {
1164 string msg("LammpsInterface::group_index() lammps group "+name+" does not exist");
1165 throw ATC_Error(msg);
1166 }
1167 return igroup;
1168 }
1169
group_inverse_mask(int iGroup) const1170 int LammpsInterface::group_inverse_mask(int iGroup) const
1171 {
1172 return lammps_->group->inversemask[iGroup];
1173 }
1174
group_name(int iGroup) const1175 char * LammpsInterface::group_name(int iGroup) const
1176 {
1177 return lammps_->group->names[iGroup];
1178 }
1179
group_bounds(int iGroup,double * b) const1180 void LammpsInterface::group_bounds(int iGroup, double * b) const
1181 {
1182 lammps_->group->bounds(iGroup, b);
1183 }
1184
1185 // -----------------------------------------------------------------
1186 // memory interface methods
1187 // -----------------------------------------------------------------
1188
1189
create_1d_double_array(int length,const char * name) const1190 double * LammpsInterface::create_1d_double_array(int length, const char *name) const {
1191 double * myArray;
1192 return lammps_->memory->create(myArray, length, name);
1193 }
1194
grow_1d_double_array(double * array,int length,const char * name) const1195 double *LammpsInterface::grow_1d_double_array(double *array,
1196 int length,
1197 const char *name) const
1198 {
1199 return lammps_->memory->grow(array, length, name);
1200 }
1201
destroy_1d_double_array(double * d) const1202 void LammpsInterface::destroy_1d_double_array(double * d) const {
1203 lammps_->memory->destroy(d);
1204 }
1205
create_2d_double_array(int n1,int n2,const char * name) const1206 double ** LammpsInterface::create_2d_double_array(int n1, int n2, const char *name) const {
1207 double ** myArray;
1208 return lammps_->memory->create(myArray, n1, n2, name);
1209 }
1210
destroy_2d_double_array(double ** d) const1211 void LammpsInterface::destroy_2d_double_array(double **d) const {
1212 lammps_->memory->destroy(d);
1213 }
1214
grow_2d_double_array(double ** array,int n1,int n2,const char * name) const1215 double **LammpsInterface::grow_2d_double_array(double **array,
1216 int n1,
1217 int n2,
1218 const char *name) const
1219 {
1220 return lammps_->memory->grow(array, n1, n2, name);
1221 }
1222
1223
create_1d_int_array(int length,const char * name) const1224 int * LammpsInterface::create_1d_int_array(int length, const char *name) const {
1225 int * myArray;
1226 return lammps_->memory->create(myArray, length, name);
1227 }
1228
grow_1d_int_array(int * array,int length,const char * name) const1229 int *LammpsInterface::grow_1d_int_array(int *array,
1230 int length,
1231 const char *name) const
1232 {
1233 return lammps_->memory->grow(array, length, name);
1234 }
1235
destroy_1d_int_array(int * d) const1236 void LammpsInterface::destroy_1d_int_array(int * d) const {
1237 lammps_->memory->destroy(d);
1238 }
1239
create_2d_int_array(int n1,int n2,const char * name) const1240 int ** LammpsInterface::create_2d_int_array(int n1, int n2, const char *name) const {
1241 int ** myArray;
1242 return lammps_->memory->create(myArray, n1, n2, name);
1243 }
1244
destroy_2d_int_array(int ** i) const1245 void LammpsInterface::destroy_2d_int_array(int **i) const {
1246 lammps_->memory->destroy(i);
1247 }
1248
grow_2d_int_array(int ** array,int n1,int n2,const char * name) const1249 int ** LammpsInterface::grow_2d_int_array(int **array, int n1, int n2, const char *name) const {
1250 return lammps_->memory->grow(array, n1, n2, name);
1251 }
1252
1253 // -----------------------------------------------------------------
1254 // update interface methods
1255 // -----------------------------------------------------------------
1256
dt() const1257 double LammpsInterface::dt() const { return lammps_->update->dt; }
1258
ntimestep() const1259 bigint LammpsInterface::ntimestep() const { return lammps_->update->ntimestep; }
1260
nsteps() const1261 int LammpsInterface::nsteps() const { return lammps_->update->nsteps; }
1262
1263
1264 // -----------------------------------------------------------------
1265 // neighbor list interface methods
1266 // -----------------------------------------------------------------
1267
sbmask(int j) const1268 int LammpsInterface::sbmask(int j) const {return j >> SBBITS & 3; }
1269
set_list(int id,LAMMPS_NS::NeighList * ptr)1270 void LammpsInterface::set_list(int id, LAMMPS_NS::NeighList *ptr) { list_ = ptr; }
1271
neighbor_list_inum() const1272 int LammpsInterface::neighbor_list_inum() const { return list_->inum; }
1273
neighbor_list_numneigh() const1274 int * LammpsInterface::neighbor_list_numneigh() const { return list_->numneigh; }
1275
neighbor_list_ilist() const1276 int * LammpsInterface::neighbor_list_ilist() const { return list_->ilist; }
1277
neighbor_list_firstneigh() const1278 int ** LammpsInterface::neighbor_list_firstneigh() const { return list_->firstneigh; }
1279
neighbor_ago() const1280 int LammpsInterface::neighbor_ago() const { return lammps_->neighbor->ago; }
1281
reneighbor_frequency() const1282 int LammpsInterface::reneighbor_frequency() const {return lammps_->neighbor->every; }
1283
1284 // -----------------------------------------------------------------
1285 // bond list interface methods
1286 // -----------------------------------------------------------------
1287
bond_list_length() const1288 int LammpsInterface::bond_list_length() const { return lammps_->neighbor->nbondlist; }
bond_list() const1289 int** LammpsInterface::bond_list() const { return lammps_->neighbor->bondlist; }
1290
1291 // -----------------------------------------------------------------
1292 // region interface methods
1293 // -----------------------------------------------------------------
1294
region_name(int iRegion) const1295 char * LammpsInterface::region_name(int iRegion) const
1296 {
1297 return lammps_->domain->regions[iRegion]->id;
1298 }
1299
region_style(int iRegion) const1300 char * LammpsInterface::region_style(int iRegion) const
1301 {
1302 return lammps_->domain->regions[iRegion]->style;
1303 }
1304
region_xlo(int iRegion) const1305 double LammpsInterface::region_xlo(int iRegion) const
1306 {
1307 return lammps_->domain->regions[iRegion]->extent_xlo;
1308 }
1309
region_xhi(int iRegion) const1310 double LammpsInterface::region_xhi(int iRegion) const
1311 {
1312 return lammps_->domain->regions[iRegion]->extent_xhi;
1313 }
1314
region_ylo(int iRegion) const1315 double LammpsInterface::region_ylo(int iRegion) const
1316 {
1317 return lammps_->domain->regions[iRegion]->extent_ylo;
1318 }
1319
region_yhi(int iRegion) const1320 double LammpsInterface::region_yhi(int iRegion) const
1321 {
1322 return lammps_->domain->regions[iRegion]->extent_yhi;
1323 }
1324
region_zlo(int iRegion) const1325 double LammpsInterface::region_zlo(int iRegion) const
1326 {
1327 return lammps_->domain->regions[iRegion]->extent_zlo;
1328 }
1329
region_zhi(int iRegion) const1330 double LammpsInterface::region_zhi(int iRegion) const
1331 {
1332 return lammps_->domain->regions[iRegion]->extent_zhi;
1333 }
1334
region_xscale(int iRegion) const1335 double LammpsInterface::region_xscale(int iRegion) const
1336 {
1337 return lammps_->domain->regions[iRegion]->xscale;
1338 }
1339
region_yscale(int iRegion) const1340 double LammpsInterface::region_yscale(int iRegion) const
1341 {
1342 return lammps_->domain->regions[iRegion]->yscale;
1343 }
1344
region_zscale(int iRegion) const1345 double LammpsInterface::region_zscale(int iRegion) const
1346 {
1347 return lammps_->domain->regions[iRegion]->zscale;
1348 }
1349
region_match(int iRegion,double x,double y,double z) const1350 int LammpsInterface::region_match(int iRegion, double x, double y, double z) const {
1351 return lammps_->domain->regions[iRegion]->match(x,y,z);
1352 }
1353
1354 // -----------------------------------------------------------------
1355 // compute methods
1356 // -----------------------------------------------------------------
compute_pointer(string tag) const1357 COMPUTE_POINTER LammpsInterface::compute_pointer(string tag) const
1358 {
1359 // get the compute id
1360 char * name = const_cast <char*> (tag.c_str());
1361 int id = lammps_->modify->find_compute(name);
1362 if (id < 0) {
1363 string msg("Could not find compute "+tag);
1364 msg += tag;
1365 throw ATC_Error(msg);
1366 }
1367 // get the compute
1368 LAMMPS_NS::Compute* cmpt = lammps_->modify->compute[id];
1369 // insert it into our set, recall it won't be added if it already exists
1370 computePointers_.insert(cmpt);
1371 return cmpt;
1372 }
1373
computes_addstep(int step) const1374 void LammpsInterface::computes_addstep(int step) const
1375 {
1376 set<LAMMPS_NS::Compute * >::iterator iter;
1377 for (iter = computePointers_.begin(); iter != computePointers_.end(); iter++) {
1378 (*iter)->addstep(step);
1379 }
1380 }
1381
compute_addstep(COMPUTE_POINTER computePointer,int step) const1382 void LammpsInterface::compute_addstep(COMPUTE_POINTER computePointer, int step) const
1383 {
1384 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1385 cmpt->addstep(step);
1386 }
1387
compute_matchstep(COMPUTE_POINTER computePointer,int step) const1388 int LammpsInterface::compute_matchstep(COMPUTE_POINTER computePointer, int step) const
1389 {
1390 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1391 return cmpt->matchstep(step);
1392 }
1393
reset_invoked_flag(COMPUTE_POINTER computePointer) const1394 void LammpsInterface::reset_invoked_flag(COMPUTE_POINTER computePointer) const
1395 {
1396 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1397 cmpt->invoked_flag = 0;
1398 }
1399
compute_ncols_peratom(COMPUTE_POINTER computePointer) const1400 int LammpsInterface::compute_ncols_peratom(COMPUTE_POINTER computePointer) const
1401 {
1402 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1403 int ndof = cmpt->size_peratom_cols;
1404 if (ndof == 0 ) ndof = 1;
1405 return ndof;
1406 }
1407
compute_vector_peratom(COMPUTE_POINTER computePointer) const1408 double* LammpsInterface::compute_vector_peratom(COMPUTE_POINTER computePointer) const
1409 {
1410 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1411 if (!(cmpt->invoked_flag & INVOKED_PERATOM)) {
1412 cmpt->compute_peratom();
1413 cmpt->invoked_flag |= INVOKED_PERATOM;
1414 }
1415 return cmpt->vector_atom;
1416 }
1417
compute_array_peratom(COMPUTE_POINTER computePointer) const1418 double** LammpsInterface::compute_array_peratom(COMPUTE_POINTER computePointer) const
1419 {
1420 LAMMPS_NS::Compute* cmpt = const_to_active(computePointer);
1421 if (!(cmpt->invoked_flag & INVOKED_PERATOM)) {
1422 cmpt->compute_peratom();
1423 cmpt->invoked_flag |= INVOKED_PERATOM;
1424 }
1425 return cmpt->array_atom;
1426 }
1427
const_to_active(COMPUTE_POINTER computePointer) const1428 LAMMPS_NS::Compute * LammpsInterface::const_to_active(COMPUTE_POINTER computePointer) const
1429 {
1430 LAMMPS_NS::Compute* cmpt = const_cast<LAMMPS_NS::Compute* >(computePointer);
1431 set<LAMMPS_NS::Compute * >::iterator cmptPtr;
1432 cmptPtr = computePointers_.find(cmpt);
1433 if (cmptPtr != computePointers_.end())
1434 return *cmptPtr;
1435 else
1436 throw ATC_Error("Requested bad computer pointer");
1437 }
1438
1439 // -----------------------------------------------------------------
1440 // compute pe/atom interface methods
1441 // - the only compute "owned" by ATC
1442 // -----------------------------------------------------------------
create_compute_pe_peratom(void) const1443 int LammpsInterface::create_compute_pe_peratom(void) const
1444 {
1445 char **list = new char*[4];
1446 string atomPeName = compute_pe_name();
1447 list[0] = (char *) atomPeName.c_str();
1448 list[1] = (char *) "all";
1449 list[2] = (char *) "pe/atom";
1450 list[3] = (char *) "pair";
1451
1452 int icompute = lammps_->modify->find_compute(list[0]);
1453 if (icompute < 0) {
1454 lammps_->modify->add_compute(4,list);
1455 icompute = lammps_->modify->find_compute(list[0]);
1456 }
1457 delete [] list;
1458 if (! atomPE_ ) {
1459 atomPE_ = lammps_->modify->compute[icompute];
1460 }
1461 computePointers_.insert(atomPE_);
1462 stringstream ss;
1463 ss << "peratom PE compute created with ID: " << icompute;
1464 ATC::LammpsInterface::instance()->print_msg_once(ss.str());
1465 return icompute;
1466 }
1467
compute_pe_peratom(void) const1468 double * LammpsInterface::compute_pe_peratom(void) const
1469 {
1470 if (atomPE_) {
1471 atomPE_->compute_peratom();
1472 return atomPE_->vector_atom;
1473 }
1474 else {
1475 return NULL;
1476 }
1477 }
1478
1479 /* ---------------------------------------------------------------------- */
1480
unwrap_coordinates(int iatom,double * xatom) const1481 void LammpsInterface::unwrap_coordinates(int iatom, double* xatom) const
1482 {
1483 double **x = lammps_->atom->x;
1484 int *image = lammps_->atom->image;
1485
1486 double *h = lammps_->domain->h;
1487 double xprd = lammps_->domain->xprd;
1488 double yprd = lammps_->domain->yprd;
1489 double zprd = lammps_->domain->zprd;
1490 int xbox,ybox,zbox;
1491
1492 // for triclinic, need to unwrap current atom coord via h matrix
1493
1494 if (lammps_->domain->triclinic == 0) {
1495 xbox = (image[iatom] & 1023) - 512;
1496 ybox = (image[iatom] >> 10 & 1023) - 512;
1497 zbox = (image[iatom] >> 20) - 512;
1498 xatom[0] = x[iatom][0] + xbox*xprd;
1499 xatom[1] = x[iatom][1] + ybox*yprd;
1500 xatom[2] = x[iatom][2] + zbox*zprd;
1501 } else {
1502 xbox = (image[iatom] & 1023) - 512;
1503 ybox = (image[iatom] >> 10 & 1023) - 512;
1504 zbox = (image[iatom] >> 20) - 512;
1505 xatom[0] = x[iatom][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
1506 xatom[1] = x[iatom][1] + h[1]*ybox + h[3]*zbox;
1507 xatom[2] = x[iatom][2] + h[2]*zbox;
1508 }
1509 }
1510
1511 /* ---------------------------------------------------------------------- */
1512
pair_eam() const1513 LAMMPS_NS::PairEAM* LammpsInterface::pair_eam() const
1514 {
1515 //if (typeid(lammps_->force->pair) == typeid(LAMMPS_NS::PairEAM)) {
1516 // return lammps_->force->pair;
1517 //}
1518 LAMMPS_NS::PairEAM* pair_eam = dynamic_cast<LAMMPS_NS::PairEAM*> (lammps_->force->pair);
1519 if (pair_eam != NULL) {
1520 return pair_eam;
1521 }
1522 else {
1523 throw ATC_Error("LAMMPS Pair object is not of the derived class type PairEAM");
1524 }
1525 }
1526
1527 // -----------------------------------------------------------------
1528 // other methods
1529 // -----------------------------------------------------------------
1530
1531 /** Return lammps pointer -- only as a last resort! */
lammps_pointer() const1532 LAMMPS_NS::LAMMPS * LammpsInterface::lammps_pointer() const { return lammps_; }
1533
1534 }
1535