SSAGES  0.8.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  inline double roundf(double x)
32  {
33  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
34  }
35 
37 
43  class Snapshot
44  {
45  private:
47  mxx::comm comm_;
48 
49  unsigned int wid_;
50  unsigned int nlocal_;
51 
52  std::string ID_;
53 
56 
58 
60 
62 
63  std::vector<Vector3> positions_;
64  std::vector<Integer3> images_;
65  std::vector<Vector3> velocities_;
66  std::vector<Vector3> forces_;
67  std::vector<double> masses_;
68  std::vector<double> charges_;
71  std::vector<std::vector<double> > sigma_;
72 
73  size_t iteration_;
74  size_t targetiter_;
75  double temperature_;
76  double energy_;
77  double kb_;
78 
79  double dielectric_;
80  double qqrd2e_;
81  bool changed_;
82 
83  public:
85 
92  Snapshot(const MPI_Comm& comm, unsigned int wid) :
93  comm_(comm), wid_(wid), H_(), Hinv_(), virial_(Matrix3::Zero()),
94  origin_({0,0,0}), isperiodic_({true, true, true}), positions_(0),
95  images_(0), velocities_(0), forces_(0), masses_(0), atomids_(0),
96  types_(0), iteration_(0), temperature_(0), energy_(0), kb_(0)
97  {}
98 
100 
103  size_t GetIteration() const {return iteration_; }
104 
106 
109  size_t GetTargetIterations() const {return targetiter_; }
110 
112 
115  double GetTemperature() const {return temperature_; }
116 
118 
121  double GetEnergy() const { return energy_; }
122 
124 
127  const Matrix3& GetHMatrix() const { return H_; }
128 
130 
133  const Matrix3& GetVirial() const { return virial_; }
134 
136 
139  Matrix3& GetVirial() { return virial_; }
140 
142 
145  const Vector3& GetOrigin() const { return origin_; }
146 
148 
151  const Bool3& IsPeriodic() const {return isperiodic_; }
152 
154 
157  double GetVolume() const { return H_.determinant(); }
158 
160 
163  double GetKb() const { return kb_; }
164 
166 
169  double GetDielectric() const { return dielectric_; }
170 
172 
175  double Getqqrd2e() const { return qqrd2e_; }
176 
178 
184  const mxx::comm& GetCommunicator() const
185  {
186  return comm_;
187  }
188 
190 
193  unsigned GetWalkerID() const { return wid_; }
194 
195 
197 
200  unsigned GetNumAtoms() const { return nlocal_; }
201 
203 
206  void SetIteration(size_t iteration)
207  {
208  iteration_ = iteration;
209  changed_ = true;
210  }
211 
213 
216  void SetTargetIterations(int target)
217  {
218  targetiter_ = target;
219  changed_ = true;
220  }
221 
223 
226  void SetTemperature(double temperature)
227  {
228  temperature_ = temperature;
229  changed_ = true;
230  }
231 
233 
236  void SetEnergy(double energy)
237  {
238  energy_ = energy;
239  changed_ = true;
240  }
241 
243 
246  void SetHMatrix(const Matrix3& hmat)
247  {
248  H_ = hmat;
249  Hinv_ = hmat.inverse();
250  changed_ = true;
251  }
252 
254 
257  void SetVirial(const Matrix3& virial)
258  {
259  virial_ = virial;
260  changed_ = true;
261  }
262 
264 
267  void SetOrigin(const Vector3& origin)
268  {
269  origin_ = origin;
270  changed_ = true;
271  }
272 
274 
277  void SetPeriodicity(const Bool3& isperiodic)
278  {
279  isperiodic_ = isperiodic;
280  changed_ = true;
281  }
282 
284 
287  void SetKb(double kb)
288  {
289  kb_ = kb;
290  changed_ = true;
291  }
292 
294 
297  void SetDielectric(double dielectric)
298  {
299  dielectric_ = dielectric;
300  changed_ = true;
301  }
302 
304 
307  void Setqqrd2e(double qqrd2e)
308  {
309  qqrd2e_ = qqrd2e;
310  changed_ = true;
311  }
312 
314 
317  void SetNumAtoms(unsigned int natoms) { nlocal_ = natoms; }
318 
320 
323  const std::vector<Vector3>& GetPositions() const { return positions_; }
324 
326  std::vector<Vector3>& GetPositions()
327  {
328  changed_ = true;
329  return positions_;
330  }
331 
333 
336  const std::vector<Integer3>& GetImageFlags() const { return images_; }
337 
339  std::vector<Integer3>& GetImageFlags()
340  {
341  changed_ = true;
342  return images_;
343  }
344 
346 
349  const std::vector<Vector3>& GetVelocities() const { return velocities_; }
350 
352  std::vector<Vector3>& GetVelocities()
353  {
354  changed_ = true;
355  return velocities_;
356  }
357 
359 
362  const std::vector<Vector3>& GetForces() const { return forces_; }
363 
365  std::vector<Vector3>& GetForces()
366  {
367  changed_ = true;
368  return forces_;
369  }
370 
372 
378  const std::vector<double>& GetMasses() const { return masses_; }
379 
381  std::vector<double>& GetMasses()
382  {
383  changed_ = true;
384  return masses_;
385  }
386 
388 
392  Vector3 ScaleVector(const Vector3& v) const
393  {
394  return Hinv_*(v-origin_);
395  }
396 
398 
409  Vector3 UnwrapVector(const Vector3& v, const Integer3& image) const
410  {
411  return H_*image.cast<double>()+v;
412  }
413 
415 
418  void ApplyMinimumImage(Vector3* v) const
419  {
420  Vector3 scaled = Hinv_*(*v);
421 
422  for(int i = 0; i < 3; ++i)
423  scaled[i] -= isperiodic_[i]*roundf(scaled[i]);
424 
425  *v = H_*scaled;
426  }
427 
429 
434  {
435  Vector3 scaled = Hinv_*v;
436 
437  for(int i = 0; i < 3; ++i)
438  scaled[i] -= isperiodic_[i]*roundf(scaled[i]);
439 
440  return H_*scaled;
441  }
442 
444 
451  double TotalMass(const Label& indices) const
452  {
453  auto mtot = 0.;
454  for(auto& i : indices)
455  mtot += masses_[i];
456  MPI_Allreduce(MPI_IN_PLACE, &mtot, 1, MPI_DOUBLE, MPI_SUM, comm_);
457  return mtot;
458  }
459 
462 
466  Vector3 CenterOfMass(const Label& indices) const
467  {
468  // Get total mass.
469  auto mtot = TotalMass(indices);
470 
471  return CenterOfMass(indices, mtot);
472  }
473 
476 
483  Vector3 CenterOfMass(const Label& indices, double mtot) const
484  {
485  // Store coorinates and masses in vectors to gather.
486  std::vector<double> pos, mass, gpos, gmass;
487  std::vector<int> pcounts(comm_.size(), 0), mcounts(comm_.size(), 0);
488  std::vector<int> pdispls(comm_.size()+1, 0), mdispls(comm_.size()+1, 0);
489 
490  pcounts[comm_.rank()] = 3*indices.size();
491  mcounts[comm_.rank()] = indices.size();
492 
493  // Reduce counts.
494  MPI_Allreduce(MPI_IN_PLACE, pcounts.data(), pcounts.size(), MPI_INT, MPI_SUM, comm_);
495  MPI_Allreduce(MPI_IN_PLACE, mcounts.data(), mcounts.size(), MPI_INT, MPI_SUM, comm_);
496 
497  // Compute displacements.
498  std::partial_sum(pcounts.begin(), pcounts.end(), pdispls.begin() + 1);
499  std::partial_sum(mcounts.begin(), mcounts.end(), mdispls.begin() + 1);
500 
501  // Fill up mass and position vectors.
502  for(auto& idx : indices)
503  {
504  auto& p = positions_[idx];
505  pos.push_back(p[0]);
506  pos.push_back(p[1]);
507  pos.push_back(p[2]);
508  mass.push_back(masses_[idx]);
509  }
510 
511  // Re-size receiving vectors.
512  gpos.resize(pdispls.back(), 0);
513  gmass.resize(mdispls.back(), 0);
514 
515  // All-gather data.
516  MPI_Allgatherv(pos.data(), pos.size(), MPI_DOUBLE, gpos.data(), pcounts.data(), pdispls.data(), MPI_DOUBLE, comm_);
517  MPI_Allgatherv(mass.data(), mass.size(), MPI_DOUBLE, gmass.data(), mcounts.data(), mdispls.data(), MPI_DOUBLE, comm_);
518 
519  // Loop through atoms and compute mass weighted sum.
520  // We march linearly through list and find nearest image
521  // to each successive particle to properly unwrap object.
522  Vector3 ppos = {gpos[0], gpos[1], gpos[2]}; // Previous unwrapped position.
523  Vector3 cpos = ppos;
524  Vector3 xcm = gmass[0]*cpos;
525 
526  for(size_t i = 1, j = 3; i < gmass.size(); ++i, j += 3)
527  {
528  cpos = {gpos[j], gpos[j+1], gpos[j+2]};
529  cpos = ApplyMinimumImage(cpos - ppos) + ppos;
530  xcm += gmass[i]*cpos;
531  ppos = cpos;
532  }
533 
534  return xcm/mtot;
535  }
536 
538 
541  const Label& GetAtomIDs() const { return atomids_; }
542 
545  {
546  changed_ = true;
547  return atomids_;
548  }
549 
551 
555  int GetLocalIndex(int id) const
556  {
557 
558  auto s = std::find(atomids_.begin(), atomids_.end(), id);
559  if(s == atomids_.end())
560  return -1;
561  else
562  return s - atomids_.begin();
563  }
564 
567 
573  void GetLocalIndices(const Label& ids, Label* indices) const
574  {
575  for(auto& id : ids)
576  {
577  auto idx = GetLocalIndex(id);
578  if(idx != -1)
579  indices->push_back(idx);
580  }
581  }
582 
584 
587  const std::vector<double>& GetCharges() const { return charges_; }
588 
590  std::vector<double>& GetCharges()
591  {
592  changed_ = true;
593  return charges_;
594  }
595 
597 
600  const Label& GetAtomTypes() const { return types_; }
601 
604  {
605  changed_ = true;
606  return types_;
607  }
608 
610 
613  const std::vector<std::vector<double>>& GetSigmas() const { return sigma_; }
614 
616  std::vector<std::vector<double>>& GetSigmas()
617  {
618  changed_ = true;
619  return sigma_;
620  }
621 
623 
626  const std::string& GetSnapshotID() const { return ID_; }
627 
629  std::string& GetSnapshotID()
630  {
631  changed_ = true;
632  return ID_;
633  }
634 
636 
639  bool HasChanged() const { return changed_; }
640 
642 
645  void Changed(bool state) { changed_ = state; }
646 
648  std::vector<double> SerializePositions()
649  {
650 
651  std::vector<int> pcounts(comm_.size(), 0);
652  std::vector<int> pdispls(comm_.size()+1, 0);
653 
654  pcounts[comm_.rank()] = 3*nlocal_;
655 
656  // Reduce counts.
657  MPI_Allreduce(MPI_IN_PLACE, pcounts.data(), pcounts.size(), MPI_INT, MPI_SUM, comm_);
658 
659  // Compute displacements.
660  std::partial_sum(pcounts.begin(), pcounts.end(), pdispls.begin() + 1);
661 
662  // Re-size receiving vectors.
663  std::vector<double> positions;
664  positions.resize(pdispls.back(), 0);
665 
666  std::vector<double> ptemp;
667 
668  for(auto& p : positions_)
669  {
670  ptemp.push_back(p[0]);
671  ptemp.push_back(p[1]);
672  ptemp.push_back(p[2]);
673  }
674 
675  // All-gather data.
676  MPI_Allgatherv(ptemp.data(), ptemp.size(), MPI_DOUBLE, positions.data(), pcounts.data(), pdispls.data(), MPI_DOUBLE, comm_);
677  return positions;
678  }
679 
681  std::vector<double> SerializeVelocities()
682  {
683  std::vector<int> vcounts(comm_.size(), 0);
684  std::vector<int> vdispls(comm_.size()+1, 0);
685 
686  vcounts[comm_.rank()] = 3*nlocal_;
687 
688  // Reduce counts.
689  MPI_Allreduce(MPI_IN_PLACE, vcounts.data(), vcounts.size(), MPI_INT, MPI_SUM, comm_);
690 
691  // Compute displacements.
692  std::partial_sum(vcounts.begin(), vcounts.end(), vdispls.begin() + 1);
693 
694  // Re-size receiving vectors.
695  std::vector<double> velocities;
696  velocities.resize(vdispls.back(), 0);
697 
698  std::vector<double> vtemp;
699 
700  for(auto& v : velocities_)
701  {
702  vtemp.push_back(v[0]);
703  vtemp.push_back(v[1]);
704  vtemp.push_back(v[2]);
705  }
706 
707  // All-gather data.
708  MPI_Allgatherv(vtemp.data(), vtemp.size(), MPI_DOUBLE, velocities.data(), vcounts.data(), vdispls.data(), MPI_DOUBLE, comm_);
709  return velocities;
710  }
711 
713  std::vector<int> SerializeIDs()
714  {
715  std::vector<int> mcounts(comm_.size(), 0);
716  std::vector<int> mdispls(comm_.size()+1, 0);
717 
718  mcounts[comm_.rank()] = nlocal_;
719 
720  // Reduce counts.
721  MPI_Allreduce(MPI_IN_PLACE, mcounts.data(), mcounts.size(), MPI_INT, MPI_SUM, comm_);
722 
723  // Compute displacements.
724  std::partial_sum(mcounts.begin(), mcounts.end(), mdispls.begin() + 1);
725 
726  // Re-size receiving vectors.
727  std::vector<int> IDs;
728  IDs.resize(mdispls.back(), 0);
729 
730  // All-gather data.
731  MPI_Allgatherv(atomids_.data(), atomids_.size(), MPI_INT, IDs.data(), mcounts.data(), mdispls.data(), MPI_INT, comm_);
732  return IDs;
733  }
734 
737  };
738 }
739 
bool changed_
TRUE is Simulation state changed
Definition: Snapshot.h:81
std::vector< Integer3 > & GetImageFlags()
Access the particles image flags.
Definition: Snapshot.h:339
void SetOrigin(const Vector3 &origin)
Change the box origin.
Definition: Snapshot.h:267
std::vector< int > Label
List of integers.
Definition: types.h:48
Eigen::Matrix3d Matrix3
3x3 matrix.
Definition: types.h:42
Eigen::Vector3i Integer3
Three-dimensional integer vector.
Definition: types.h:39
std::vector< Integer3 > images_
Unwrapped positions.
Definition: Snapshot.h:64
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:200
void Changed(bool state)
Set the "changed" flag of the Snapshot.
Definition: Snapshot.h:645
void ApplyMinimumImage(Vector3 *v) const
Apply minimum image to a vector.
Definition: Snapshot.h:418
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:184
std::vector< Vector3 > forces_
Forces.
Definition: Snapshot.h:66
void SetNumAtoms(unsigned int natoms)
Set number of atoms in this snapshot.
Definition: Snapshot.h:317
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:362
size_t iteration_
Iteration of Simulation.
Definition: Snapshot.h:73
Matrix3 H_
Parrinello-Rahman box H-matrix.
Definition: Snapshot.h:54
double kb_
Kb from the MD driver.
Definition: Snapshot.h:77
Label & GetAtomIDs()
Access the atom IDs.
Definition: Snapshot.h:544
const std::string & GetSnapshotID() const
Access the snapshot ID.
Definition: Snapshot.h:626
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:43
std::vector< Vector3 > positions_
Positions.
Definition: Snapshot.h:63
std::vector< Vector3 > & GetForces()
Access the per-particle forces.
Definition: Snapshot.h:365
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:555
Matrix3 Hinv_
Parinello-Rahman box inverse.
Definition: Snapshot.h:55
void SetVirial(const Matrix3 &virial)
Change the box virial.
Definition: Snapshot.h:257
const Label & GetAtomTypes() const
Access the atom types.
Definition: Snapshot.h:600
Vector3 ScaleVector(const Vector3 &v) const
Scale a vector into fractional coordinates.
Definition: Snapshot.h:392
void SetTargetIterations(int target)
Set target iterations.
Definition: Snapshot.h:216
double energy_
Average per-particle energy.
Definition: Snapshot.h:76
std::vector< double > SerializePositions()
Return the serialized positions across all local cores.
Definition: Snapshot.h:648
double TotalMass(const Label &indices) const
Compute the total mass of a group of particles based on index.
Definition: Snapshot.h:451
std::vector< double > charges_
Charges.
Definition: Snapshot.h:68
const std::vector< std::vector< double > > & GetSigmas() const
Access the atom sigmas.
Definition: Snapshot.h:613
Vector3 ApplyMinimumImage(const Vector3 &v) const
Apply minimum image to a vector.
Definition: Snapshot.h:433
double temperature_
Current temperature.
Definition: Snapshot.h:75
double GetKb() const
Get system Kb.
Definition: Snapshot.h:163
void SetEnergy(double energy)
Change the energy.
Definition: Snapshot.h:236
size_t GetIteration() const
Get the current iteration.
Definition: Snapshot.h:103
std::vector< int > SerializeIDs()
Return the serialized positions across all local cores.
Definition: Snapshot.h:713
std::vector< double > & GetCharges()
Access the atom charges.
Definition: Snapshot.h:590
Matrix3 virial_
Virial tensor.
Definition: Snapshot.h:57
std::vector< std::vector< double > > & GetSigmas()
Access the atom sigmas.
Definition: Snapshot.h:616
Bool3 isperiodic_
Periodicity of box.
Definition: Snapshot.h:61
Vector3 CenterOfMass(const Label &indices, double mtot) const
Definition: Snapshot.h:483
std::vector< std::vector< double > > sigma_
Sigma.
Definition: Snapshot.h:71
Snapshot(const MPI_Comm &comm, unsigned int wid)
Constructor.
Definition: Snapshot.h:92
double GetVolume() const
Get system volume.
Definition: Snapshot.h:157
void SetPeriodicity(const Bool3 &isperiodic)
Change the periodicity of the system.
Definition: Snapshot.h:277
mxx::comm comm_
Local snapshot (walker) communicator.
Definition: Snapshot.h:47
std::vector< double > masses_
Masses.
Definition: Snapshot.h:67
double GetTemperature() const
Get current temperature.
Definition: Snapshot.h:115
Label types_
List of Atom types.
Definition: Snapshot.h:70
size_t targetiter_
Iteration target of simulation.
Definition: Snapshot.h:74
~Snapshot()
Destructor.
Definition: Snapshot.h:736
const std::vector< Integer3 > & GetImageFlags() const
Access the particles image flags.
Definition: Snapshot.h:336
double dielectric_
Dielectric.
Definition: Snapshot.h:79
size_t GetTargetIterations() const
Get target iterations.
Definition: Snapshot.h:109
Label & GetAtomTypes()
Access the atom types.
Definition: Snapshot.h:603
double GetDielectric() const
Get dielectric constant.
Definition: Snapshot.h:169
double Getqqrd2e() const
Get qqrd2e value.
Definition: Snapshot.h:175
bool HasChanged() const
Query if Snapshot was modified.
Definition: Snapshot.h:639
std::vector< Vector3 > velocities_
Velocities.
Definition: Snapshot.h:65
double GetEnergy() const
Get per-particle energy.
Definition: Snapshot.h:121
void SetDielectric(double dielectric)
Set the dielectric constant.
Definition: Snapshot.h:297
void SetKb(double kb)
Change the kb.
Definition: Snapshot.h:287
void SetTemperature(double temperature)
Change the temperature.
Definition: Snapshot.h:226
Vector3 CenterOfMass(const Label &indices) const
Definition: Snapshot.h:466
double roundf(double x)
Quick helper function to round a double.
Definition: Snapshot.h:31
const std::vector< Vector3 > & GetVelocities() const
Access the particle velocities.
Definition: Snapshot.h:349
unsigned int wid_
Walker ID.
Definition: Snapshot.h:49
const Bool3 & IsPeriodic() const
Get periodicity of three dimensions.
Definition: Snapshot.h:151
std::vector< Vector3 > & GetVelocities()
Access the particle velocities.
Definition: Snapshot.h:352
Vector3 origin_
Box origin.
Definition: Snapshot.h:59
const Matrix3 & GetVirial() const
Get box virial.
Definition: Snapshot.h:133
std::string ID_
ID string.
Definition: Snapshot.h:52
const std::vector< double > & GetMasses() const
Const access to the particle masses.
Definition: Snapshot.h:378
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
void SetHMatrix(const Matrix3 &hmat)
Change the Box H-matrix.
Definition: Snapshot.h:246
std::string & GetSnapshotID()
Access the snapshot ID.
Definition: Snapshot.h:629
unsigned GetWalkerID() const
Get walker ID.
Definition: Snapshot.h:193
std::vector< double > SerializeVelocities()
Return the serialized velocities across all local cores.
Definition: Snapshot.h:681
Label atomids_
List of Atom IDs.
Definition: Snapshot.h:69
const Label & GetAtomIDs() const
Access the atom IDs.
Definition: Snapshot.h:541
Vector3 UnwrapVector(const Vector3 &v, const Integer3 &image) const
Unwrap a vector&#39;s real coordinates according to its image replica count.
Definition: Snapshot.h:409
double qqrd2e_
qqrd2e
Definition: Snapshot.h:80
void Setqqrd2e(double qqrd2e)
Set the value for qqrd2e.
Definition: Snapshot.h:307
const Vector3 & GetOrigin() const
Get origin of the system.
Definition: Snapshot.h:145
void SetIteration(size_t iteration)
Set the iteration.
Definition: Snapshot.h:206
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
const Matrix3 & GetHMatrix() const
Get system H-matrix.
Definition: Snapshot.h:127
std::vector< double > & GetMasses()
Const access to the particle masses.
Definition: Snapshot.h:381
unsigned int nlocal_
Number of atoms located on this snapshot.
Definition: Snapshot.h:50
std::vector< Vector3 > & GetPositions()
Access the particle positions.
Definition: Snapshot.h:326
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:573
Matrix3 & GetVirial()
Get box virial.
Definition: Snapshot.h:139
const std::vector< double > & GetCharges() const
Access the atom charges.
Definition: Snapshot.h:587
Eigen::Matrix< bool, 3, 1 > Bool3
Three-dimensional boolean.
Definition: types.h:36