SSAGES  0.8.5
Software Suite for Advanced General Ensemble Simulations
Swarm.cpp
1 
20 #include "Swarm.h"
21 #include "CVs/CVManager.h"
22 #include "Snapshot.h"
23 #include "spline.h"
24 #include <cmath>
25 #include <iostream>
26 #include <iomanip>
27 #include <string>
28 #include <sstream>
29 #include <algorithm>
30 
31 namespace SSAGES
32 {
33 
34  //Helper function to check if CVs are initialized correctly
35  bool Swarm::CVInitialized(const CVList& cvs)
36  {
37  double threshold = 0.2;
38  const double eps = 0.0000000001;
39  double diff;
40 
41  //On the first iteration, check that the CVs are within (threshold*100)% of the center value they're associated to
42  for(size_t i = 0; i < cvs.size(); i++)
43  {
44  if(std::abs(centers_[i]) <= eps)
45  {//e.g. if centers_[i] = 0.0
46  if(std::abs(cvs[i]->GetValue()) <= threshold)
47  {
48  diff = 0; //If current value is approximately zero then go ahead
49  }
50  else
51  {
52  diff = 2*threshold; //Will always be greater than threshold
53  }
54  //diff = std::abs((cvs[i]->GetValue() - (centers_[i]+eps))) / ((cvs[i]->GetValue() + (eps + centers_[i]))/2.0);
55  //std::cout << "MPIID = " << mpiid_ << " and center[" << i << "] = 0.0 and diff = " << diff << std::endl;
56  }
57  else
58  {
59  diff = std::abs((cvs[i]->GetValue() - (centers_[i]))) / ((cvs[i]->GetValue() + (centers_[i]))/2.0);
60  //std::cout << "MPIID = " << mpiid_ << " and center[" << i << "] != 0.0 and diff = " << diff << std::endl;
61  }
62  if(std::abs(diff) >= threshold)
63  {
64  return false; //proceed to initialize again
65  }
66  }
67  return true; //OK to move on to regular sampling
68  }
69 
70  void Swarm::PostIntegration(Snapshot* snapshot, const CVManager& cvmanager)
71  {
72  auto cvs = cvmanager.GetCVs(cvmask_);
73  auto& forces = snapshot->GetForces();
74  auto& positions = snapshot->GetPositions();
75  auto& velocities = snapshot->GetVelocities();
76 
77  //First check if a previous snapshot was stored; the system is definitely initialized if so
78  if(snapshot_stored)
79  {
80  initialized = true;
81  }
82  else
83  {
84  initialized = CVInitialized(cvs); //If no snapshot stored, evalute whether the system finish initializing
85  }
86 
87  original_initialized = initialized; //Will remember whether this system was initialized before the MPI call overrides it
88 
89  //If system finished initializing and hadn't entered the sampling blocks, check if a snapshot was stored. If not, store one and don't store any future snapshots.
91  {
92  if(!snapshot_stored)
93  {
96  prev_IDs_[index_] = snapshot->SerializeIDs();
97  snapshot_stored = true;
98  }
99  }
100  MPI_Allreduce(MPI_IN_PLACE, &initialized, 1, MPI_INT, MPI_LAND, world_); //Check if all systems are initialized; initialize is only true if everything node finished initializing
101  //If any system wasn't initialized and the sampling block were not entered, begin restraining each node
103  {
104  //Do restrained sampling; NO trajectory harvesting
106  {
107  //If the system was not initialized originally, restrain it
108  for(size_t i = 0; i < cvs.size(); i++)
109  {
110  //Get current CV and gradient
111  auto& cv = cvs[i];
112  auto& grad = cv->GetGradient();
113 
114  //Compute dV/dCV
115  auto D = cvspring_[i]*(cv->GetDifference(centers_[i]));
116 
117  //Update forces
118  for(size_t j = 0; j < forces.size(); j++)
119  {
120  for(int k = 0; k < forces[j].size(); k++)
121  {
122  forces[j][k] -= (double)D*grad[j][k];
123  }
124  }
125  }
126  }
127  else
128  {
129  //If the system was already initialized, simply reset the positions and velocities; this keeps them in their initialized states
130  index_ = 0;
131  for(auto& force: forces)
132  force.setZero();
133 
134  for(size_t i = 0; i < prev_IDs_[index_].size(); i++)
135  {
136  auto localindex = snapshot->GetLocalIndex(prev_IDs_[index_][i]);
137  if(localindex!= -1)
138  {
139  positions[localindex][0] = prev_positions_[index_][i*3];
140  positions[localindex][1] = prev_positions_[index_][i*3 + 1];
141  positions[localindex][2] = prev_positions_[index_][i*3 + 2];
142 
143  velocities[localindex][0] = prev_velocities_[index_][i*3];
144  velocities[localindex][1] = prev_velocities_[index_][i*3 + 1];
145  velocities[localindex][2] = prev_velocities_[index_][i*3 + 2];
146 
147  }
148  }
149  }
150  }
151  else
152  {
153  //Enter this block is every node was done initializing, or the sampling had already begun
154  if(!sampling_started)
155  {
156  sampling_started = true; //If the node got here right after initializing, set the flag that sampling has begun (preventing unneeded initializing in the future)
157  }
159  {
160  //Do restrained sampling, and do harvest trajectories
161  for(size_t i = 0; i < cvs.size(); i++)
162  {
163  if(iterator_ == 0)
164  {
165  index_ = 0; //Reset index when starting
166  }
167  //Get current CV and gradient
168  auto& cv = cvs[i];
169  auto& grad = cv->GetGradient();
170 
171  //Compute dV/dCV
172  auto D = cvspring_[i]*(cv->GetDifference(centers_[i]));
173 
174  //Update forces
175  for(size_t j = 0; j < forces.size(); j++)
176  {
177  for(int k = 0; k < forces[j].size(); k++)
178  {
179  forces[j][k] -= (double)D*grad[j][k];
180  }
181  }
182  }
184  {
185  //Harvest a trajectory every harvest_length_ steps
186  if(iterator_ % harvest_length_ == 0)
187  {
190  prev_IDs_[index_] = snapshot->SerializeIDs();
191  index_++;
192  }
193  }
195  {
196  //Reset positions and forces before first call to unrestrained sampling
197  index_ = 0;
198  for(auto& force: forces)
199  force.setZero();
200 
201  for(size_t i = 0; i < prev_IDs_[index_].size(); i++)
202  {
203  auto localindex = snapshot->GetLocalIndex(prev_IDs_[index_][i]);
204  if(localindex!= -1)
205  {
206  positions[localindex][0] = prev_positions_[index_][i*3];
207  positions[localindex][1] = prev_positions_[index_][i*3 + 1];
208  positions[localindex][2] = prev_positions_[index_][i*3 + 2];
209 
210  velocities[localindex][0] = prev_velocities_[index_][i*3];
211  velocities[localindex][1] = prev_velocities_[index_][i*3 + 1];
212  velocities[localindex][2] = prev_velocities_[index_][i*3 + 2];
213 
214  }
215  }
216  }
217  iterator_++;
218  }
220  {
221  //Launch unrestrained trajectories
223  {
224  //End of trajectory, harvest drift
225  for(size_t i = 0; i < cv_drift_.size(); i++)
226  {
227  cv_drift_[i] = (cv_drift_[i]*index_ + cvs[i]->GetMinimumImage(centers_[i]) - centers_[i]) / (index_+1); //Calculate running average of drifts
228  }
229  //Set up for next trajectory
230  index_++;
231 
233  {
234  for(auto& force: forces)
235  force.setZero();
236 
237  for(size_t i = 0; i < prev_IDs_[index_].size(); i++)
238  {
239  auto localindex = snapshot->GetLocalIndex(prev_IDs_[index_][i]);
240  if(localindex!= -1)
241  {
242  positions[localindex][0] = prev_positions_[index_][i*3];
243  positions[localindex][1] = prev_positions_[index_][i*3 + 1];
244  positions[localindex][2] = prev_positions_[index_][i*3 + 2];
245 
246  velocities[localindex][0] = prev_velocities_[index_][i*3];
247  velocities[localindex][1] = prev_velocities_[index_][i*3 + 1];
248  velocities[localindex][2] = prev_velocities_[index_][i*3 + 2];
249 
250  }
251  }
252  }
253  }
254  iterator_++;
255  }
256  else
257  {
258  //Evolve CVs, reparametrize, and reset vectors
259  iteration_++;
260 
261  MPI_Barrier(world_);
262  StringUpdate();
263  CheckEnd(cvs);
264  MPI_Barrier(world_);
265  UpdateWorldString(cvs);
266  PrintString(cvs);
267 
268  iterator_ = 0;
269  index_ = 0;
270  snapshot_stored = false;
271 
272  for(size_t i = 0; i < cv_drift_.size(); i++)
273  {
274  cv_drift_[i] = 0;
275  }
276  }
277  }
278  }
279 
281  {
282  // Update node locations toward running averages:
283  for(size_t i = 0; i < centers_.size(); i++)
284  {
285 
286  centers_[i] = centers_[i] + cv_drift_[i];
287  }
288 
289  std::vector<double> lcv0, ucv0;
290  lcv0.resize(centers_.size(), 0);
291  ucv0.resize(centers_.size(), 0);
292 
293  GatherNeighbors(&lcv0, &ucv0);
294 
295  double alphastar = distance(centers_, lcv0);
296 
297  StringReparam(alphastar);
298 
299  }
300 }
301 
double distance(const std::vector< double > &x, const std::vector< double > &y) const
Helper function for calculating distances.
Definition: StringMethod.h:102
Collective variable manager.
Definition: CVManager.h:42
unsigned int initialize_steps_
Total number of MD steps for initialization for one iteration.
Definition: Swarm.h:40
bool original_initialized
Flag for whether a system was initialized before it checked whether other systems were...
Definition: Swarm.h:84
bool CVInitialized(const CVList &cvs)
Helper function check if CVs are initialized correctly.
Definition: Swarm.cpp:35
std::vector< CollectiveVariable * > CVList
List of Collective Variables.
Definition: types.h:51
std::vector< std::vector< double > > prev_velocities_
Store velocities for starting trajectories.
Definition: StringMethod.h:88
int initialized
Flag for whether a system is initialized at a given iteration.
Definition: Swarm.h:81
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:47
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 > cv_drift_
Drift of CVs for one iteration.
Definition: Swarm.h:37
unsigned int swarm_length_
Length of unrestrained trajectories.
Definition: Swarm.h:52
bool sampling_started
Flag for determing whether to perform initialization or not.
Definition: Swarm.h:75
std::vector< double > SerializePositions()
Return the serialized positions across all local cores.
Definition: Snapshot.h:712
void PostIntegration(Snapshot *snapshot, const class CVManager &cvmanager) override
Method call post integration.
Definition: Swarm.cpp:70
mxx::comm world_
Global MPI communicator.
Definition: Method.h:46
bool CheckEnd(const CVList &CV)
Check if method reached one of the exit criteria.
bool snapshot_stored
Flag for whether a snapshot was stored in the umbrella sampling.
Definition: Swarm.h:78
std::vector< int > SerializeIDs()
Return the serialized IDs across all local cores.
Definition: Snapshot.h:783
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:616
std::vector< std::vector< int > > prev_IDs_
Store atom IDs for starting trajectories.
Definition: StringMethod.h:91
unsigned int number_trajectories_
Number of trajectories per swarm.
Definition: Swarm.h:49
void UpdateWorldString(const CVList &cvs)
Update the world string over MPI.
void StringReparam(double alpha_star)
Reparameterize the string.
unsigned int harvest_length_
Length to run before harvesting a trajectory for unrestrained sampling.
Definition: Swarm.h:43
unsigned int restrained_steps_
Total number of restrained MD steps for one iteration.
Definition: Swarm.h:46
void PrintString(const CVList &CV)
Prints the CV positions to file.
unsigned int iterator_
The local method iterator.
Definition: StringMethod.h:70
std::vector< unsigned int > cvmask_
Mask which identifies which CVs to act on.
Definition: Method.h:50
unsigned int iteration_
The global method iteration.
Definition: StringMethod.h:73
std::vector< double > cvspring_
Vector of spring constants.
Definition: StringMethod.h:67
unsigned int unrestrained_steps_
Total number of unrestrained MD steps for one iteration.
Definition: Swarm.h:55
unsigned int index_
For indexing trajectory vectors.
Definition: Swarm.h:58
void StringUpdate() override
Updates the positions of the string.
Definition: Swarm.cpp:280
std::vector< double > SerializeVelocities()
Return the serialized velocities across all local cores.
Definition: Snapshot.h:748
const std::vector< Vector3 > & GetVelocities() const
Access the particle velocities.
Definition: Snapshot.h:353
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:366
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:327
std::vector< std::vector< double > > prev_positions_
Store positions for starting trajectories.
Definition: StringMethod.h:85
CVList GetCVs(const std::vector< unsigned int > &mask=std::vector< unsigned int >()) const
Get CV iterator.
Definition: CVManager.h:81