SSAGES  0.9.3
Software Suite for Advanced General Ensemble Simulations
Snapshot.h
1 
21 #pragma once
22 #include <numeric>
23 #include <vector>
24 #include <memory>
25 #include <mxx/comm.hpp>
26 #include "types.h"
27 
28 namespace SSAGES
29 {
31 
35  inline double roundf(double x)
36  {
37  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
38  }
39 
41 
47  class Snapshot
48  {
49  private:
51  mxx::comm comm_;
52 
53  unsigned int wid_;
54  unsigned int nlocal_;
55 
56  std::string ID_;
57 
60 
62 
64 
66 
67  std::vector<Vector3> positions_;
68  std::vector<Vector3> velocities_;
69  std::vector<Vector3> forces_;
70  std::vector<double> masses_;
71  std::vector<double> charges_;
74 
75  size_t iteration_;
76  size_t targetiter_;
77  double temperature_;
78  double energy_;
79  double kb_;
80 
81  double dielectric_;
82  double qqrd2e_;
83  bool changed_;
84 
85  public:
87 
94  Snapshot(const MPI_Comm& comm, unsigned int wid) :
95  comm_(comm), wid_(wid), H_(), Hinv_(), virial_(Matrix3::Zero()),
96  origin_({0,0,0}), isperiodic_({true, true, true}), positions_(0),
97  velocities_(0), forces_(0), masses_(0), atomids_(0), types_(0),
98  iteration_(0), temperature_(0), energy_(0), kb_(0)
99  {}
100 
102 
105  size_t GetIteration() const {return iteration_; }
106 
108 
111  size_t GetTargetIterations() const {return targetiter_; }
112 
114 
117  double GetTemperature() const {return temperature_; }
118 
120 
123  double GetEnergy() const { return energy_; }
124 
126 
129  const Matrix3& GetHMatrix() const { return H_; }
130 
132 
135  const Matrix3& GetVirial() const { return virial_; }
136 
138 
141  Matrix3& GetVirial() { return virial_; }
142 
144 
147  const Vector3& GetOrigin() const { return origin_; }
148 
150 
153  const Bool3& IsPeriodic() const {return isperiodic_; }
154 
156 
159  double GetVolume() const { return H_.determinant(); }
160 
162 
165  double GetKb() const { return kb_; }
166 
168 
171  double GetDielectric() const { return dielectric_; }
172 
174 
177  double Getqqrd2e() const { return qqrd2e_; }
178 
180 
186  const mxx::comm& GetCommunicator() const
187  {
188  return comm_;
189  }
190 
192 
195  unsigned GetWalkerID() const { return wid_; }
196 
197 
199 
202  unsigned GetNumAtoms() const { return nlocal_; }
203 
205 
208  void SetIteration(size_t iteration)
209  {
210  iteration_ = iteration;
211  changed_ = true;
212  }
213 
215 
218  void SetTargetIterations(int target)
219  {
220  targetiter_ = target;
221  changed_ = true;
222  }
223 
225 
228  void SetTemperature(double temperature)
229  {
230  temperature_ = temperature;
231  changed_ = true;
232  }
233 
235 
238  void SetEnergy(double energy)
239  {
240  energy_ = energy;
241  changed_ = true;
242  }
243 
245 
248  void SetHMatrix(const Matrix3& hmat)
249  {
250  H_ = hmat;
251  Hinv_ = hmat.inverse();
252  changed_ = true;
253  }
254 
256 
259  void SetVirial(const Matrix3& virial)
260  {
261  virial_ = virial;
262  changed_ = true;
263  }
264 
266 
269  void SetOrigin(const Vector3& origin)
270  {
271  origin_ = origin;
272  changed_ = true;
273  }
274 
276 
279  void SetPeriodicity(const Bool3& isperiodic)
280  {
281  isperiodic_ = isperiodic;
282  changed_ = true;
283  }
284 
286 
289  void SetKb(double kb)
290  {
291  kb_ = kb;
292  changed_ = true;
293  }
294 
296 
299  void SetDielectric(double dielectric)
300  {
301  dielectric_ = dielectric;
302  changed_ = true;
303  }
304 
306 
309  void Setqqrd2e(double qqrd2e)
310  {
311  qqrd2e_ = qqrd2e;
312  changed_ = true;
313  }
314 
316 
319  void SetNumAtoms(unsigned int natoms) { nlocal_ = natoms; }
320 
322 
325  const std::vector<Vector3>& GetPositions() const { return positions_; }
326 
328  std::vector<Vector3>& GetPositions()
329  {
330  changed_ = true;
331  return positions_;
332  }
333 
335 
338  const std::vector<Vector3>& GetVelocities() const { return velocities_; }
339 
341  std::vector<Vector3>& GetVelocities()
342  {
343  changed_ = true;
344  return velocities_;
345  }
346 
348 
351  const std::vector<Vector3>& GetForces() const { return forces_; }
352 
354  std::vector<Vector3>& GetForces()
355  {
356  changed_ = true;
357  return forces_;
358  }
359 
361 
367  const std::vector<double>& GetMasses() const { return masses_; }
368 
370  std::vector<double>& GetMasses()
371  {
372  changed_ = true;
373  return masses_;
374  }
375 
377 
381  Vector3 ScaleVector(const Vector3& v) const
382  {
383  return Hinv_*(v-origin_);
384  }
385 
387 
390  void ApplyMinimumImage(Vector3* v) const
391  {
392  Vector3 scaled = Hinv_*(*v);
393 
394  for(int i = 0; i < 3; ++i)
395  scaled[i] -= isperiodic_[i]*roundf(scaled[i]);
396 
397  *v = H_*scaled;
398  }
399 
401 
406  {
407  Vector3 scaled = Hinv_*v;
408 
409  for(int i = 0; i < 3; ++i)
410  scaled[i] -= isperiodic_[i]*roundf(scaled[i]);
411 
412  return H_*scaled;
413  }
414 
416 
423  double TotalMass(const Label& indices) const
424  {
425  auto mtot = 0.;
426  for(auto& i : indices)
427  mtot += masses_[i];
428  MPI_Allreduce(MPI_IN_PLACE, &mtot, 1, MPI_DOUBLE, MPI_SUM, comm_);
429  return mtot;
430  }
431 
433 
442  Vector3 CenterOfMass(const Label& indices, bool mass_weight = true) const
443  {
444  // Store coordinates and masses in vectors to gather.
445  std::vector<double> pos, mass, gpos, gmass;
446  std::vector<int> pcounts(comm_.size(), 0), mcounts(comm_.size(), 0);
447  std::vector<int> pdispls(comm_.size()+1, 0), mdispls(comm_.size()+1, 0);
448 
449  // Calculate total mass
450  double mtot = mass_weight ? TotalMass(indices) : indices.size();
451 
452  pcounts[comm_.rank()] = 3*indices.size();
453  mcounts[comm_.rank()] = indices.size();
454 
455  // Reduce counts.
456  MPI_Allreduce(MPI_IN_PLACE, pcounts.data(), pcounts.size(), MPI_INT, MPI_SUM, comm_);
457  MPI_Allreduce(MPI_IN_PLACE, mcounts.data(), mcounts.size(), MPI_INT, MPI_SUM, comm_);
458 
459  // Compute displacements.
460  std::partial_sum(pcounts.begin(), pcounts.end(), pdispls.begin() + 1);
461  std::partial_sum(mcounts.begin(), mcounts.end(), mdispls.begin() + 1);
462 
463  // Fill up mass and position vectors.
464  for(auto& idx : indices)
465  {
466  auto& p = positions_[idx];
467  pos.push_back(p[0]);
468  pos.push_back(p[1]);
469  pos.push_back(p[2]);
470  mass.push_back(masses_[idx]);
471  }
472 
473  // Re-size receiving vectors.
474  gpos.resize(pdispls.back(), 0);
475  gmass.resize(mdispls.back(), 1.0); // fictitious mass of 1
476 
477  // All-gather data.
478  MPI_Allgatherv(pos.data(), pos.size(), MPI_DOUBLE, gpos.data(), pcounts.data(), pdispls.data(), MPI_DOUBLE, comm_);
479  if(mass_weight) { // overwrite gmass with true mass values
480  MPI_Allgatherv(mass.data(), mass.size(), MPI_DOUBLE, gmass.data(), mcounts.data(), mdispls.data(), MPI_DOUBLE, comm_);
481  }
482 
483  // Loop through atoms and compute mass weighted sum.
484  // We march linearly through list and find nearest image
485  // to each successive particle to properly unwrap object.
486  Vector3 ppos = {gpos[0], gpos[1], gpos[2]}; // Previous unwrapped position.
487  Vector3 cpos = ppos;
488  Vector3 xcm = gmass[0]*cpos;
489 
490  for(size_t i = 1, j = 3; i < gmass.size(); ++i, j += 3)
491  {
492  cpos = {gpos[j], gpos[j+1], gpos[j+2]};
493  cpos = ApplyMinimumImage(cpos - ppos) + ppos;
494  xcm += gmass[i]*cpos;
495  ppos = cpos;
496  }
497 
498  return xcm/mtot;
499  }
500 
502 
505  const Label& GetAtomIDs() const { return atomids_; }
506 
509  {
510  changed_ = true;
511  return atomids_;
512  }
513 
515 
519  int GetLocalIndex(int id) const
520  {
521 
522  auto s = std::find(atomids_.begin(), atomids_.end(), id);
523  if(s == atomids_.end())
524  return -1;
525  else
526  return s - atomids_.begin();
527  }
528 
531 
537  void GetLocalIndices(const Label& ids, Label* indices) const
538  {
539  for(auto& id : ids)
540  {
541  auto idx = GetLocalIndex(id);
542  if(idx != -1)
543  indices->push_back(idx);
544  }
545  }
546 
548 
551  const std::vector<double>& GetCharges() const { return charges_; }
552 
554  std::vector<double>& GetCharges()
555  {
556  changed_ = true;
557  return charges_;
558  }
559 
561 
564  const Label& GetAtomTypes() const { return types_; }
565 
568  {
569  changed_ = true;
570  return types_;
571  }
572 
574 
577  const std::string& GetSnapshotID() const { return ID_; }
578 
580  std::string& GetSnapshotID()
581  {
582  changed_ = true;
583  return ID_;
584  }
585 
587 
590  bool HasChanged() const { return changed_; }
591 
593 
596  void Changed(bool state) { changed_ = state; }
597 
599 
602  std::vector<double> SerializePositions()
603  {
604 
605  std::vector<int> pcounts(comm_.size(), 0);
606  std::vector<int> pdispls(comm_.size()+1, 0);
607 
608  pcounts[comm_.rank()] = 3*nlocal_;
609 
610  // Reduce counts.
611  MPI_Allreduce(MPI_IN_PLACE, pcounts.data(), pcounts.size(), MPI_INT, MPI_SUM, comm_);
612 
613  // Compute displacements.
614  std::partial_sum(pcounts.begin(), pcounts.end(), pdispls.begin() + 1);
615 
616  // Re-size receiving vectors.
617  std::vector<double> positions;
618  positions.resize(pdispls.back(), 0);
619 
620  std::vector<double> ptemp;
621 
622  for(auto& p : positions_)
623  {
624  ptemp.push_back(p[0]);
625  ptemp.push_back(p[1]);
626  ptemp.push_back(p[2]);
627  }
628 
629  // All-gather data.
630  MPI_Allgatherv(ptemp.data(), ptemp.size(), MPI_DOUBLE, positions.data(), pcounts.data(), pdispls.data(), MPI_DOUBLE, comm_);
631  return positions;
632  }
633 
635 
638  std::vector<double> SerializeVelocities()
639  {
640  std::vector<int> vcounts(comm_.size(), 0);
641  std::vector<int> vdispls(comm_.size()+1, 0);
642 
643  vcounts[comm_.rank()] = 3*nlocal_;
644 
645  // Reduce counts.
646  MPI_Allreduce(MPI_IN_PLACE, vcounts.data(), vcounts.size(), MPI_INT, MPI_SUM, comm_);
647 
648  // Compute displacements.
649  std::partial_sum(vcounts.begin(), vcounts.end(), vdispls.begin() + 1);
650 
651  // Re-size receiving vectors.
652  std::vector<double> velocities;
653  velocities.resize(vdispls.back(), 0);
654 
655  std::vector<double> vtemp;
656 
657  for(auto& v : velocities_)
658  {
659  vtemp.push_back(v[0]);
660  vtemp.push_back(v[1]);
661  vtemp.push_back(v[2]);
662  }
663 
664  // All-gather data.
665  MPI_Allgatherv(vtemp.data(), vtemp.size(), MPI_DOUBLE, velocities.data(), vcounts.data(), vdispls.data(), MPI_DOUBLE, comm_);
666  return velocities;
667  }
668 
670 
673  std::vector<int> SerializeIDs()
674  {
675  std::vector<int> mcounts(comm_.size(), 0);
676  std::vector<int> mdispls(comm_.size()+1, 0);
677 
678  mcounts[comm_.rank()] = nlocal_;
679 
680  // Reduce counts.
681  MPI_Allreduce(MPI_IN_PLACE, mcounts.data(), mcounts.size(), MPI_INT, MPI_SUM, comm_);
682 
683  // Compute displacements.
684  std::partial_sum(mcounts.begin(), mcounts.end(), mdispls.begin() + 1);
685 
686  // Re-size receiving vectors.
687  std::vector<int> IDs;
688  IDs.resize(mdispls.back(), 0);
689 
690  // All-gather data.
691  MPI_Allgatherv(atomids_.data(), atomids_.size(), MPI_INT, IDs.data(), mcounts.data(), mdispls.data(), MPI_INT, comm_);
692  return IDs;
693  }
694 
697  };
698 }
699 
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:48
double GetKb() const
Get system Kb.
Definition: Snapshot.h:165
double energy_
Average per-particle energy.
Definition: Snapshot.h:78
const Matrix3 & GetVirial() const
Get box virial.
Definition: Snapshot.h:135
std::vector< double > masses_
Masses.
Definition: Snapshot.h:70
double Getqqrd2e() const
Get qqrd2e value.
Definition: Snapshot.h:177
std::string & GetSnapshotID()
Access the snapshot ID.
Definition: Snapshot.h:580
Vector3 origin_
Box origin.
Definition: Snapshot.h:63
void Setqqrd2e(double qqrd2e)
Set the value for qqrd2e.
Definition: Snapshot.h:309
const std::string & GetSnapshotID() const
Access the snapshot ID.
Definition: Snapshot.h:577
unsigned int nlocal_
Number of atoms located on this snapshot.
Definition: Snapshot.h:54
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:202
bool HasChanged() const
Query if Snapshot was modified.
Definition: Snapshot.h:590
std::vector< double > & GetCharges()
Access the atom charges.
Definition: Snapshot.h:554
double dielectric_
Dielectric.
Definition: Snapshot.h:81
std::vector< double > SerializeVelocities()
Return the serialized velocities across all local cores.
Definition: Snapshot.h:638
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:537
Label atomids_
List of Atom IDs.
Definition: Snapshot.h:72
void Changed(bool state)
Set the "changed" flag of the Snapshot.
Definition: Snapshot.h:596
std::vector< Vector3 > & GetPositions()
Access the particle positions.
Definition: Snapshot.h:328
double GetTemperature() const
Get current temperature.
Definition: Snapshot.h:117
double temperature_
Current temperature.
Definition: Snapshot.h:77
Bool3 isperiodic_
Periodicity of box.
Definition: Snapshot.h:65
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:351
std::vector< int > SerializeIDs()
Return the serialized IDs across all local cores.
Definition: Snapshot.h:673
mxx::comm comm_
Local snapshot (walker) communicator.
Definition: Snapshot.h:51
void SetPeriodicity(const Bool3 &isperiodic)
Change the periodicity of the system.
Definition: Snapshot.h:279
const Bool3 & IsPeriodic() const
Get periodicity of three dimensions.
Definition: Snapshot.h:153
void SetHMatrix(const Matrix3 &hmat)
Change the Box H-matrix.
Definition: Snapshot.h:248
std::vector< double > & GetMasses()
Const access to the particle masses.
Definition: Snapshot.h:370
const Matrix3 & GetHMatrix() const
Get system H-matrix.
Definition: Snapshot.h:129
double GetEnergy() const
Get per-particle energy.
Definition: Snapshot.h:123
void SetNumAtoms(unsigned int natoms)
Set number of atoms in this snapshot.
Definition: Snapshot.h:319
double GetVolume() const
Get system volume.
Definition: Snapshot.h:159
std::vector< Vector3 > positions_
Positions.
Definition: Snapshot.h:67
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:186
std::vector< Vector3 > velocities_
Velocities.
Definition: Snapshot.h:68
void SetVirial(const Matrix3 &virial)
Change the box virial.
Definition: Snapshot.h:259
double qqrd2e_
qqrd2e
Definition: Snapshot.h:82
void SetTargetIterations(int target)
Set target iterations.
Definition: Snapshot.h:218
size_t GetTargetIterations() const
Get target iterations.
Definition: Snapshot.h:111
Label & GetAtomTypes()
Access the atom types.
Definition: Snapshot.h:567
Label types_
List of Atom types.
Definition: Snapshot.h:73
void SetIteration(size_t iteration)
Set the iteration.
Definition: Snapshot.h:208
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:519
std::vector< double > charges_
Charges.
Definition: Snapshot.h:71
Vector3 ApplyMinimumImage(const Vector3 &v) const
Apply minimum image to a vector.
Definition: Snapshot.h:405
double TotalMass(const Label &indices) const
Compute the total mass of a group of particles based on index.
Definition: Snapshot.h:423
Vector3 ScaleVector(const Vector3 &v) const
Scale a vector into fractional coordinates.
Definition: Snapshot.h:381
void SetEnergy(double energy)
Change the energy.
Definition: Snapshot.h:238
size_t targetiter_
Iteration target of simulation.
Definition: Snapshot.h:76
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:325
bool changed_
TRUE is Simulation state changed
Definition: Snapshot.h:83
Matrix3 H_
Parrinello-Rahman box H-matrix.
Definition: Snapshot.h:58
const Label & GetAtomTypes() const
Access the atom types.
Definition: Snapshot.h:564
Matrix3 virial_
Virial tensor.
Definition: Snapshot.h:61
Matrix3 & GetVirial()
Get box virial.
Definition: Snapshot.h:141
unsigned int wid_
Walker ID.
Definition: Snapshot.h:53
Matrix3 Hinv_
Parinello-Rahman box inverse.
Definition: Snapshot.h:59
Label & GetAtomIDs()
Access the atom IDs.
Definition: Snapshot.h:508
double kb_
Kb from the MD driver.
Definition: Snapshot.h:79
size_t GetIteration() const
Get the current iteration.
Definition: Snapshot.h:105
std::string ID_
ID string.
Definition: Snapshot.h:56
void SetTemperature(double temperature)
Change the temperature.
Definition: Snapshot.h:228
std::vector< Vector3 > forces_
Forces.
Definition: Snapshot.h:69
const std::vector< double > & GetMasses() const
Const access to the particle masses.
Definition: Snapshot.h:367
const std::vector< double > & GetCharges() const
Access the atom charges.
Definition: Snapshot.h:551
Vector3 CenterOfMass(const Label &indices, bool mass_weight=true) const
Compute center of mass of a group of atoms based on index.
Definition: Snapshot.h:442
std::vector< Vector3 > & GetForces()
Access the per-particle forces.
Definition: Snapshot.h:354
~Snapshot()
Destructor.
Definition: Snapshot.h:696
double GetDielectric() const
Get dielectric constant.
Definition: Snapshot.h:171
void SetOrigin(const Vector3 &origin)
Change the box origin.
Definition: Snapshot.h:269
const Label & GetAtomIDs() const
Access the atom IDs.
Definition: Snapshot.h:505
Snapshot(const MPI_Comm &comm, unsigned int wid)
Constructor.
Definition: Snapshot.h:94
void ApplyMinimumImage(Vector3 *v) const
Apply minimum image to a vector.
Definition: Snapshot.h:390
std::vector< Vector3 > & GetVelocities()
Access the particle velocities.
Definition: Snapshot.h:341
const std::vector< Vector3 > & GetVelocities() const
Access the particle velocities.
Definition: Snapshot.h:338
void SetDielectric(double dielectric)
Set the dielectric constant.
Definition: Snapshot.h:299
const Vector3 & GetOrigin() const
Get origin of the system.
Definition: Snapshot.h:147
std::vector< double > SerializePositions()
Return the serialized positions across all local cores.
Definition: Snapshot.h:602
void SetKb(double kb)
Change the kb.
Definition: Snapshot.h:289
unsigned GetWalkerID() const
Get walker ID.
Definition: Snapshot.h:195
size_t iteration_
Iteration of Simulation.
Definition: Snapshot.h:75
Eigen::Matrix3d Matrix3
3x3 matrix.
Definition: types.h:42
double roundf(double x)
Quick helper function to round a double.
Definition: Snapshot.h:35
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
Eigen::Matrix< bool, 3, 1 > Bool3
Three-dimensional boolean.
Definition: types.h:36
std::vector< int > Label
List of integers.
Definition: types.h:48