SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
ParallelBetaRMSDCV.h
1 
21 #pragma once
22 
23 #include "CollectiveVariable.h"
24 #include "Validator/ObjectRequirement.h"
25 #include "Drivers/DriverException.h"
26 #include "Snapshot.h"
27 #include "schema.h"
28 #include "Utility/ReadBackbone.h"
29 #include <stdexcept>
30 
31 namespace SSAGES
32 {
34 
44  {
45  private:
46 
48  std::vector<int> resids_;
49 
51  //std::vector<int> atomids_;
52  std::vector< std::vector<std::string> > atomids_;
53 
55  std::vector<int> ids_only_;
56 
58  std::string refpdb_;
59 
61  std::vector<Vector3> refalpha_;
62 
64  double unitconv_;
65 
67  int mode_;
68 
70  std::vector<double> refdists_;
71 
72  public:
74 
85  ParallelBetaRMSDCV(std::vector<int> resids, std::string refpdb, double unitconv, int mode) :
86  resids_(resids), refpdb_(refpdb), unitconv_(unitconv), mode_(mode)
87  {
88  if(resids_.size() != 2 ){
89  throw std::invalid_argument("ParallelBetaRMSDCV: Input must designate range of residues with 2 residue numbers.");
90  }
91 
92  resids_.clear();
93 
94  if(resids[0] >= resids[1]){
95  throw std::invalid_argument("ParallelBetaRMSDCV: Input must list lower residue index first: please reverse residue range.");
96  } else if(resids[1] - resids[0] < 5) {
97  throw std::invalid_argument("ParallelBetaRMSDCV: Residue range must span at least 6 residues for secondary structure calculation.");
98  }
99 
100  std::cout << "ParallelBetaRMSDCV: Calculating parallel beta sheet character from residue " << resids[0] << " to " << resids[1] << "." << std::endl;
101 
102  for(int i = resids[0]; i <= resids[1]; i++){
103  resids_.push_back(i);
104  }
105  }
106 
108 
111  void Initialize(const Snapshot& snapshot) override
112  {
113  atomids_ = ReadBackbone::GetPdbBackbone(refpdb_, resids_);
114  ids_only_.resize(atomids_[0].size());
115 
116  // check to find all necessary atoms
117  std::transform(atomids_[0].begin(), atomids_[0].end(), ids_only_.begin(), [](std::string val) {return std::stoi(val);});
118  using std::to_string;
119  std::vector<int> found;
120  snapshot.GetLocalIndices(ids_only_, &found);
121  size_t nfound = found.size();
122  MPI_Allreduce(MPI_IN_PLACE, &nfound, 1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
123  if(nfound != resids_.size() * 5)
124  throw BuildException({
125  "ParallelBetaRMSDCV: Expected to find " +
126  to_string(resids_.size() * 5) +
127  " atoms, but only found " +
128  to_string(nfound) + "."
129  });
130 
131  // reference structure for parallel beta sheet, values in nanometers
132  // Reference structure taken from values used in Plumed 2.2.0
133  std::vector<Vector3> refalpha;
134  refalpha.push_back(unitconv_ * Vector3{-.1439, -.5122, -.1144}); // N residue i
135  refalpha.push_back(unitconv_ * Vector3{-.0816, -.3803, -.1013}); // CA
136  refalpha.push_back(unitconv_ * Vector3{ .0099, -.3509, -.2206}); // CB
137  refalpha.push_back(unitconv_ * Vector3{-.1928, -.2770, -.0952}); // C
138  refalpha.push_back(unitconv_ * Vector3{-.2991, -.2970, -.1551}); // O
139  refalpha.push_back(unitconv_ * Vector3{-.1698, -.1687, -.0215}); // N residue i+1
140  refalpha.push_back(unitconv_ * Vector3{-.2681, -.0613, -.0143}); // CA
141  refalpha.push_back(unitconv_ * Vector3{-.3323, -.0477, .1267}); // CB
142  refalpha.push_back(unitconv_ * Vector3{-.1984, .0681, -.0574}); // C
143  refalpha.push_back(unitconv_ * Vector3{-.0807, .0921, -.0273}); // O
144  refalpha.push_back(unitconv_ * Vector3{-.2716, .1492, -.1329}); // N residue i+2
145  refalpha.push_back(unitconv_ * Vector3{-.2196, .2731, -.1883}); // CA
146  refalpha.push_back(unitconv_ * Vector3{-.2263, .2692, -.3418}); // CB
147  refalpha.push_back(unitconv_ * Vector3{-.2989, .3949, -.1433}); // C
148  refalpha.push_back(unitconv_ * Vector3{-.4214, .3989, -.1583}); // O
149  refalpha.push_back(unitconv_ * Vector3{ .2464, -.4352, .2149}); // N residue h
150  refalpha.push_back(unitconv_ * Vector3{ .3078, -.3170, .1541}); // CA
151  refalpha.push_back(unitconv_ * Vector3{ .3398, -.3415, .0060}); // CB
152  refalpha.push_back(unitconv_ * Vector3{ .2080, -.2021, .1639}); // C
153  refalpha.push_back(unitconv_ * Vector3{ .0938, -.2178, .1225}); // O
154  refalpha.push_back(unitconv_ * Vector3{ .2525, -.0886, .2183}); // N residue h+1
155  refalpha.push_back(unitconv_ * Vector3{ .1692, .0303, .2346}); // CA
156  refalpha.push_back(unitconv_ * Vector3{ .1541, .0665, .3842}); // CB
157  refalpha.push_back(unitconv_ * Vector3{ .2420, .1410, .1608}); // C
158  refalpha.push_back(unitconv_ * Vector3{ .3567, .1733, .1937}); // O
159  refalpha.push_back(unitconv_ * Vector3{ .1758, .1976, .0600}); // N residue h+2
160  refalpha.push_back(unitconv_ * Vector3{ .2373, .2987, -.0238}); // CA
161  refalpha.push_back(unitconv_ * Vector3{ .2367, .2527, -.1720}); // CB
162  refalpha.push_back(unitconv_ * Vector3{ .1684, .4331, -.0148}); // C
163  refalpha.push_back(unitconv_ * Vector3{ .0486, .4430, -.0415}); // O
164 
165  // calculate all necessary pairwise distances for reference structure
166  std::fill(refdists_.begin(), refdists_.end(), 0.0);
167  refdists_.resize(900, 0.0);
168  for(size_t k = 0; k < 29; k++){
169  for(size_t h = k + 1; h < 30; h++){
170  refdists_[29 * k + h] = (refalpha[k] - refalpha[h]).norm();
171  }
172  }
173  }
174 
176 
179  void Evaluate(const Snapshot& snapshot) override
180  {
181  // need atom positions for all atoms
182  const auto& pos = snapshot.GetPositions();
183  auto& comm = snapshot.GetCommunicator();
184 
185  unsigned int resgroups = resids_.size() - 2;
186  double rmsd, dist_norm, dxgrouprmsd;
187  int localidx;
188  Vector3 dist_xyz, dist_ref;
189  std::vector<Vector3> refxyz;
190  std::vector< std::vector< Vector3 > > deriv(30, std::vector<Vector3>(30, Vector3{0,0,0}));
191  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
192  grad_.resize(snapshot.GetNumAtoms(), Vector3{0,0,0});
193  val_ = 0.0;
194 
195  for(size_t i = 0; i < resgroups - 3; i++){
196  for(size_t j = i + 3; j < resgroups; j++){
197  // mode: 0 for all, 1 for inter, 2 for intra
198  if((mode_ == 0) || (mode_ == 1 && atomids_[1][5 * j] != atomids_[1][5 * i]) || (mode_ == 2 && atomids_[1][5 * j] == atomids_[1][5 * i])){
199  rmsd = 0.0;
200  std::fill(refxyz.begin(), refxyz.end(), Vector3{0,0,0});
201  refxyz.resize(30, Vector3{0,0,0});
202 
203  for(size_t k = 0; k < 15; k++){
204  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
205  if(localidx != -1) refxyz[k] = pos[localidx];
206  }
207  for(size_t k = 0; k < 15; k++){
208  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + k]);
209  if(localidx != -1) refxyz[k + 15] = pos[localidx];
210  }
211 
212  MPI_Allreduce(MPI_IN_PLACE, refxyz.data(), 90, MPI_DOUBLE, MPI_SUM, comm);
213 
214  // sum over all pairs to calculate CV
215  for(size_t k = 0; k < 29; k++){
216  for(size_t h = k + 1; h < 30; h++){
217  dist_xyz = refxyz[k] - refxyz[h];
218  dist_norm = dist_xyz.norm() - refdists_[29 * k + h];
219  rmsd += dist_norm * dist_norm;
220  deriv[k][h] = dist_xyz * dist_norm / dist_xyz.norm();
221  }
222  }
223 
224  rmsd = pow(rmsd / 435, 0.5) / 0.1;
225  val_ += (1 - pow(rmsd, 8.0)) / (1 - pow(rmsd, 12.0));
226 
227  dxgrouprmsd = pow(rmsd, 11.0) + pow(rmsd, 7.0);
228  dxgrouprmsd /= pow(rmsd, 8.0) + pow(rmsd, 4.0) + 1;
229  dxgrouprmsd /= pow(rmsd, 8.0) + pow(rmsd, 4.0) + 1;
230  dxgrouprmsd *= -40. / 435;
231 
232  for(size_t k = 0; k < 15; k++){
233  for(size_t h = k + 1; h < 15; h++){
234  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
235  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h];
236  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + h]);
237  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h];
238  }
239  for(size_t h = 0; h < 15; h++){
240  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
241  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h + 15];
242  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + h]);
243  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h + 15];
244  }
245  }
246  for(size_t k = 0; k < 14; k++){
247  for(size_t h = k + 1; h < 15; h++){
248  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + k]);
249  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k + 15][h + 15];
250  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + h]);
251  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k + 15][h + 15];
252  }
253  }
254  }
255  }
256  }
257  }
258 
260  static ParallelBetaRMSDCV* Build(const Json::Value& json, const std::string& path)
261  {
262  Json::ObjectRequirement validator;
263  Json::Value schema;
264  Json::CharReaderBuilder rbuilder;
265  Json::CharReader* reader = rbuilder.newCharReader();
266 
267  reader->parse(JsonSchema::ParallelBetaRMSDCV.c_str(),
268  JsonSchema::ParallelBetaRMSDCV.c_str() + JsonSchema::ParallelBetaRMSDCV.size(),
269  &schema, NULL);
270  validator.Parse(schema, path);
271 
272  //Validate inputs
273  validator.Validate(json, path);
274  if(validator.HasErrors())
275  throw BuildException(validator.GetErrors());
276 
277  std::vector<int> resids;
278  for(auto& s : json["residue_ids"])
279  resids.push_back(s.asInt());
280  auto reference = json["reference"].asString();
281 
282  double unitconv = json.get("length_unit", 1).asDouble();
283 
284  int mode = json.get("mode", 0).asInt();
285 
286  return new ParallelBetaRMSDCV(resids, reference, unitconv, mode);
287  }
288  };
289 }
double unitconv_
Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10) ...
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:200
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:184
std::string refpdb_
Name of pdb reference for system.
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
std::vector< int > resids_
Residue IDs for secondary structure calculation.
std::vector< int > ids_only_
Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
Exception to be thrown when building the Driver fails.
Collective variable to measure parallel beta sheet secondary structure.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
std::vector< double > refdists_
Pairwise distance between backbone atoms for reference structure.
ParallelBetaRMSDCV(std::vector< int > resids, std::string refpdb, double unitconv, int mode)
Constructor.
Requirements on an object.
int mode_
Specify whether to calculate beta sheet character in intra or inter mode: 0 for either, 1 for inter, 2 for intra.
std::vector< Vector3 > refalpha_
Coordinates for reference structure.
std::vector< std::vector< std::string > > atomids_
Atom IDs for secondary structure calculation: backbone of resids_.
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
double val_
Current value of CV.
Abstract class for a collective variable.
static std::vector< std::vector< std::string > > GetPdbBackbone(std::string refpdb, std::vector< int > resids)
Read protein backbone atoms from reference PDB file.
Definition: ReadBackbone.h:57
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
static ParallelBetaRMSDCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:573
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.