SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
FiniteTempString.cpp
1 
21 #include "FiniteTempString.h"
22 #include "CVs/CVManager.h"
23 #include "Snapshot.h"
24 #include "spline.h"
25 #include <math.h>
26 #include <iostream>
27 #include <algorithm>
28 
29 namespace SSAGES
30 {
31  // Check whether CV values are within their respective Voronoi cell in CV space
32  bool FiniteTempString::InCell(const CVList& cvs) const
33  {
34  std::vector<double> dists (numnodes_, 0);
35 
36  // Record the difference between all cvs and all nodes
37  for (int i = 0; i < numnodes_; i++)
38  for(size_t j = 0; j < cvs.size(); j++)
39  dists[i]+= pow(cvs[j]->GetDifference(worldstring_[i][j]),2);
40 
41  if(std::min_element(dists.begin(), dists.end()) - dists.begin() == mpiid_)
42  return true;
43 
44  return false;
45  }
46 
47  // Post-integration hook.
48  void FiniteTempString::PostIntegration(Snapshot* snapshot, const CVManager& cvmanager)
49  {
50  auto cvs = cvmanager.GetCVs(cvmask_);
51  auto& forces = snapshot->GetForces();
52  bool insidecell;
53 
55  {
56  //If the system was not going to run the umbrella anymore, reset its position
57  for(auto& force : forces)
58  {
59  force.setZero();
60  }
61 
62  auto& Pos = snapshot->GetPositions();
63 
64  for(size_t i = 0; i < prev_IDs_[0].size(); i++)
65  {
66  auto localindex = snapshot->GetLocalIndex(prev_IDs_[0][i]);
67  if(localindex!= -1)
68  {
69  Pos[localindex][0] = prev_positions_[0][i*3];
70  Pos[localindex][1] = prev_positions_[0][i*3 + 1];
71  Pos[localindex][2] = prev_positions_[0][i*3 + 2];
72  }
73  }
74  }
75  for(auto& cv : cvs)
76  {
77  //Trigger a rebuild of the CVs since we reset the positions
78  cv->Evaluate(*snapshot);
79  }
80  insidecell = InCell(cvs);
81 
82  MPI_Allreduce(MPI_IN_PLACE, &run_umbrella_, 1, MPI_INT, MPI_LOR, world_); //Compares all run_umbrellas on all processors - if any is true, all are set to true
83 
84  if(run_umbrella_)
85  {
87  {
88  //After reaching the minimum umbrella sampling checkpoint, evaluate whether to keep going (per individual image)
89  if(insidecell)
90  {
91  //Stop umbrella sampling if inside cell
92  run_umbrella_ = false;
93  reset_for_umbrella = true;
94 
95  //This node is done initializing; so store this snapshot
96  prev_positions_[0] = snapshot->SerializePositions();
97  prev_IDs_[0] = snapshot->SerializeIDs();
98  }
99  umbrella_iter_ = 1;
100  }
101  else
102  {
103  if(!reset_for_umbrella)
104  {
105  //If not at the minimum step checkpoint AND not a system which was already done, add restraining force
106  for(size_t i = 0; i < cvs.size(); i++)
107  {
108  // Get current cv and gradient
109  auto& cv = cvs[i];
110  auto& grad = cv->GetGradient();
111 
112  // Compute dV/dCV
113  auto D = cvspring_[i]*(cv->GetDifference(centers_[i]));
114 
115  // Update forces
116  for(size_t j = 0; j < forces.size(); j++)
117  forces[j] -= D*grad[j];
118  }
119  umbrella_iter_++; //Progress toward checkpoint
120  }
121  else
122  {
123  //End umbrella sampling of this node if it was already initialized
124  run_umbrella_ = false;
125  }
126  }
127  return; //Don't do any sampling after restraining
128  }
129  else
130  {
131  //If run umbrella was false (only possible if every node was done) then set reset to false for next string method iteration (only proceed past umbrella sampling if no images require further umbrella sampling)
132  reset_for_umbrella = false;
133  }
134 
135  if(!insidecell)
136  {
137  //If last step took the system outside the Voronoi cell, reset the position and zero the force to reset the trajectory
138  for(auto& force : forces)
139  force.setZero();
140 
141  auto& Pos = snapshot->GetPositions();
142 
143  for(size_t i = 0; i < prev_IDs_[0].size(); i++)
144  {
145  auto localindex = snapshot->GetLocalIndex(prev_IDs_[0][i]);
146  if(localindex!= -1)
147  {
148  Pos[localindex][0] = prev_positions_[0][i*3];
149  Pos[localindex][1] = prev_positions_[0][i*3 + 1];
150  Pos[localindex][2] = prev_positions_[0][i*3 + 2];
151  }
152  }
153 
154  // Calculate running averages for each CV at each node based on previous CV
155  for(size_t i = 0; i < newcenters_.size(); i++)
156  {
159  }
160  }
161  else
162  {
163  // Calculate running averages for each CV at each node
164  for(size_t i = 0; i < newcenters_.size(); i++)
165  {
166  newcenters_[i] = newcenters_[i] * (iteration_ * blockiterations_ + iterator_ - 1) + cvs[i]->GetMinimumImage(centers_[i]);
168  }
169 
170  //Store info about this step in prev_vectors
171  prev_CVs_.clear();
172  for(size_t i = 0; i < centers_.size(); i++)
173  {
174  prev_CVs_.push_back(cvs[i]->GetMinimumImage(centers_[i]));
175  }
176  prev_positions_[0] = snapshot->SerializePositions();
177  prev_IDs_[0] = snapshot->SerializeIDs();
178  }
179 
180  // Update the string, every blockiterations_ string method iterations
181  if(iterator_ % blockiterations_ == 0)
182  {
183  MPI_Barrier(world_);
184  StringUpdate();
185  CheckEnd(cvs);
186  MPI_Barrier(world_);
187  UpdateWorldString(cvs);
188  PrintString(cvs);
189 
190  iterator_ = 0; //Number of steps within current method iteration
191  iteration_++; //Increment string method iteration
192 
193  if(!InCell(cvs))
194  {
195  run_umbrella_ = true;
196  reset_for_umbrella = false;
197  }
198  MPI_Allreduce(MPI_IN_PLACE, &run_umbrella_, 1, MPI_INT, MPI_LOR, world_); //Compares all run_umbrellas on all processors - if any is true, all are set to true
199  }
200 
201  iterator_++; //Iterate number of steps within current method iteration
202  }
203 
205  {
206  std::vector<double> lcv0, ucv0;
207  lcv0.resize(centers_.size(), 0);
208  ucv0.resize(centers_.size(), 0);
209 
210  GatherNeighbors(&lcv0, &ucv0);
211 
212  // Update node locations toward running averages:
213  for(size_t i = 0; i < centers_.size(); i++)
214  {
215  if(mpiid_ == 0 || mpiid_ == numnodes_ - 1)
216  centers_[i] =centers_[i] - tau_ * (centers_[i] - newcenters_[i]);
217  else
218  centers_[i] = centers_[i] - tau_ * (centers_[i] - newcenters_[i]) +
219  (kappa_ * numnodes_ * tau_ * (ucv0[i] + lcv0[i] - 2 * centers_[i]));
220  }
221 
222  GatherNeighbors(&lcv0, &ucv0);
223  double alphastar = distance(centers_, lcv0);
224  StringReparam(alphastar);
225  }
226 }
std::vector< CollectiveVariable * > GetCVs(const std::vector< unsigned int > &mask=std::vector< unsigned int >()) const
Get CV iterator.
Definition: CVManager.h:80
Collective variable manager.
Definition: CVManager.h:40
int run_umbrella_
Flag to run umbrella or not during post-integration.
void PostIntegration(Snapshot *snapshot, const class CVManager &cvmanager) override
Post-integration hook.
void StringUpdate() override
Updates the string according to the FTS method.
std::vector< CollectiveVariable * > CVList
List of Collective Variables.
Definition: types.h:51
unsigned int min_num_umbrella_steps_
Minimum number of steps to apply umbrella sampling.
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:362
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:43
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:555
void GatherNeighbors(std::vector< double > *lcv0, std::vector< double > *ucv0)
Gather neighbors over MPI.
std::vector< double > centers_
CV starting location values.
Definition: StringMethod.h:43
std::vector< double > newcenters_
CV starting location values.
Definition: StringMethod.h:46
std::vector< double > SerializePositions()
Return the serialized positions across all local cores.
Definition: Snapshot.h:648
mxx::comm world_
Global MPI communicator.
Definition: Method.h:46
unsigned int blockiterations_
Number of steps to block average the CV&#39;s postions over.
bool CheckEnd(const CVList &CV)
Check if method reached one of the exit criteria.
int numnodes_
Number of nodes on a string.
Definition: StringMethod.h:61
std::vector< int > SerializeIDs()
Return the serialized positions across all local cores.
Definition: Snapshot.h:713
std::vector< std::vector< int > > prev_IDs_
Store atom IDs for starting trajectories.
Definition: StringMethod.h:91
int mpiid_
The node this belongs to.
Definition: StringMethod.h:55
void UpdateWorldString(const CVList &cvs)
Update the world string over MPI.
void StringReparam(double alpha_star)
Reparameterize the string.
void PrintString(const CVList &CV)
Prints the CV positions to file.
bool InCell(const CVList &cvs) const
Checks if CV is in voronoi cell.
unsigned int iterator_
The local method iterator.
Definition: StringMethod.h:70
std::vector< std::vector< double > > worldstring_
The world&#39;s strings centers for each CV.
Definition: StringMethod.h:52
std::vector< unsigned int > cvmask_
Mask which identifies which CVs to act on.
Definition: Method.h:50
double tau_
Time step of string change.
unsigned int iteration_
The global method iteration.
Definition: StringMethod.h:73
std::vector< double > cvspring_
Vector of spring constants.
Definition: StringMethod.h:67
std::vector< double > prev_CVs_
Stores the last positions of the CVs.
bool reset_for_umbrella
Flag for whether a system was to run umbrella sampling before checking against other systems...
unsigned int umbrella_iter_
Iterator that keeps track of umbrella iterations.
double distance(const std::vector< double > &x, const std::vector< double > &y) const
Helper function for calculating distances.
Definition: StringMethod.h:102
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
std::vector< std::vector< double > > prev_positions_
Store positions for starting trajectories.
Definition: StringMethod.h:85
double kappa_
String modification parameter.