SSAGES  0.9.3
Software Suite for Advanced General Ensemble Simulations
AlphaRMSDCV.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< std::vector<std::string> > atomids_;
52 
54  std::vector<int> ids_only_;
55 
57  std::string refpdb_;
58 
60  double unitconv_;
61 
63  std::vector<double> refdists_;
64 
65  public:
67 
77  AlphaRMSDCV(std::vector<int> resids, std::string refpdb, double unitconv) :
78  resids_(resids), refpdb_(refpdb), unitconv_(unitconv)
79  {
80  if(resids_.size() != 2 ){
81  throw std::invalid_argument("AlphaRMSDCV: Input must designate range of residues with 2 residue numbers.");
82  }
83 
84  resids_.clear();
85 
86  if(resids[0] >= resids[1]){
87  throw std::invalid_argument("AlphaRMSDCV: Input must list lower residue index first: please reverse residue range.");
88  } else if(resids[1] - resids[0] < 5) {
89  throw std::invalid_argument("AlphaRMSDCV: Residue range must span at least 6 residues for alpha helix calculation.");
90  }
91 
92  std::cout << "AlphaRMSDCV: Calculating alpha helix character from residue " << resids[0] << " to " << resids[1] << "." << std::endl;
93 
94  for(int i = resids[0]; i <= resids[1]; i++){
95  resids_.push_back(i);
96  }
97  }
98 
100 
103  void Initialize(const Snapshot& snapshot) override
104  {
106  ids_only_.resize(atomids_[0].size());
107 
108  // check to find all necessary atoms
109  std::transform(atomids_[0].begin(), atomids_[0].end(), ids_only_.begin(), [](std::string val) {return std::stoi(val);});
110  using std::to_string;
111  std::vector<int> found;
112  snapshot.GetLocalIndices(ids_only_, &found);
113  size_t nfound = found.size();
114  MPI_Allreduce(MPI_IN_PLACE, &nfound, 1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
115  if(nfound != resids_.size() * 5)
116  throw BuildException({
117  "AlphaRMSDCV: Expected to find " +
118  to_string(resids_.size() * 5) +
119  " atoms, but only found " +
120  to_string(nfound) + "."
121  });
122 
123  // reference 'ideal' alpha helix structure, in nanometers
124  // Reference structure taken from values used in Plumed 2.2.0
125  std::vector<Vector3> refalpha;
126  refalpha.push_back(unitconv_ * Vector3{ .0733, .0519, .5298 }); // N
127  refalpha.push_back(unitconv_ * Vector3{ .1763, .0810, .4301 }); // CA
128  refalpha.push_back(unitconv_ * Vector3{ .3166, .0543, .4881 }); // CB
129  refalpha.push_back(unitconv_ * Vector3{ .1527, -.0045, .3053 }); // C
130  refalpha.push_back(unitconv_ * Vector3{ .1646, .0436, .1928 }); // O
131  refalpha.push_back(unitconv_ * Vector3{ .1180, -.1312, .3254 }); // N
132  refalpha.push_back(unitconv_ * Vector3{ .0924, -.2203, .2126 }); // CA
133  refalpha.push_back(unitconv_ * Vector3{ .0650, -.3626, .2626 }); // CB
134  refalpha.push_back(unitconv_ * Vector3{-.0239, -.1711, .1261 }); // C
135  refalpha.push_back(unitconv_ * Vector3{-.0190, -.1815, .0032 }); // O
136  refalpha.push_back(unitconv_ * Vector3{-.1280, -.1172, .1891 }); // N
137  refalpha.push_back(unitconv_ * Vector3{-.2416, -.0661, .1127 }); // CA
138  refalpha.push_back(unitconv_ * Vector3{-.3548, -.0217, .2056 }); // CB
139  refalpha.push_back(unitconv_ * Vector3{-.1964, .0529, .0276 }); // C
140  refalpha.push_back(unitconv_ * Vector3{-.2364, .0659, -.0880 }); // O
141  refalpha.push_back(unitconv_ * Vector3{-.1130, .1391, .0856 }); // N
142  refalpha.push_back(unitconv_ * Vector3{-.0620, .2565, .0148 }); // CA
143  refalpha.push_back(unitconv_ * Vector3{ .0228, .3439, .1077 }); // CB
144  refalpha.push_back(unitconv_ * Vector3{ .0231, .2129, -.1032 }); // C
145  refalpha.push_back(unitconv_ * Vector3{ .0179, .2733, -.2099 }); // O
146  refalpha.push_back(unitconv_ * Vector3{ .1028, .1084, -.0833 }); // N
147  refalpha.push_back(unitconv_ * Vector3{ .1872, .0593, -.1919 }); // CA
148  refalpha.push_back(unitconv_ * Vector3{ .2850, -.0462, -.1397 }); // CB
149  refalpha.push_back(unitconv_ * Vector3{ .1020, .0020, -.3049 }); // C
150  refalpha.push_back(unitconv_ * Vector3{ .1317, .0227, -.4224 }); // O
151  refalpha.push_back(unitconv_ * Vector3{-.0051, -.0684, -.2696 }); // N
152  refalpha.push_back(unitconv_ * Vector3{-.0927, -.1261, -.3713 }); // CA
153  refalpha.push_back(unitconv_ * Vector3{-.1933, -.2219, -.3074 }); // CB
154  refalpha.push_back(unitconv_ * Vector3{-.1663, -.0171, -.4475 }); // C
155  refalpha.push_back(unitconv_ * Vector3{-.1916, -.0296, -.5673 }); // O
156 
157  // calculate all necessary pairwise distances for reference structure
158  std::fill(refdists_.begin(), refdists_.end(), 0.0);
159  refdists_.resize(900, 0.0);
160  for(size_t k = 0; k < 29; k++){
161  for(size_t h = k + 1; h < 30; h++){
162  refdists_[29 * k + h] = (refalpha[k] - refalpha[h]).norm();
163  }
164  }
165  }
166 
168 
171  void Evaluate(const Snapshot& snapshot) override
172  {
173  // need atom positions for all atoms
174  const auto& pos = snapshot.GetPositions();
175  auto& comm = snapshot.GetCommunicator();
176 
177  double rmsd, dist_norm, dxgrouprmsd;
178  int localidx;
179  Vector3 dist_xyz, dist_ref;
180  std::vector<Vector3> refxyz;
181  std::vector< std::vector< Vector3 > > deriv(30, std::vector<Vector3>(30, Vector3{0,0,0}));
182  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
183  grad_.resize(snapshot.GetNumAtoms(), Vector3{0,0,0});
184  val_ = 0.0;
185  unsigned int resgroups = resids_.size() - 5;
186 
187  // for each set of 6 residues
188  for(size_t i = 0; i < resgroups; i++){
189 
190  // clear temp rmsd calculation
191  rmsd = 0.0;
192 
193  // load refxyz with the correct 30 reference atoms
194  std::fill(refxyz.begin(), refxyz.end(), Vector3{0,0,0});
195  refxyz.resize(30, Vector3{0,0,0});
196  for(size_t j = 0; j < 30; j++){
197  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + j]);
198  if(localidx != -1) refxyz[j] = pos[localidx];
199  }
200 
201  MPI_Allreduce(MPI_IN_PLACE, refxyz.data(), 90, MPI_DOUBLE, MPI_SUM, comm);
202 
203  // sum over 435 pairs in refxyz and refalpha_ to calculate CV
204  for(size_t j = 0; j < 29; j++){
205  for(size_t k = j + 1; k < 30; k++){
206  dist_xyz = refxyz[j] - refxyz[k];
207  dist_norm = dist_xyz.norm() - refdists_[29 * j + k];
208  rmsd += dist_norm * dist_norm;
209  deriv[j][k] = dist_xyz * dist_norm / dist_xyz.norm();
210  }
211  }
212 
213  rmsd = pow(rmsd / 435, 0.5) / 0.1;
214  val_ += (1 - pow(rmsd, 8.0)) / (1 - pow(rmsd, 12.0));
215 
216  dxgrouprmsd = pow(rmsd, 11.0) + pow(rmsd, 7.0);
217  dxgrouprmsd /= pow(rmsd, 8.0) + pow(rmsd, 4.0) + 1;
218  dxgrouprmsd /= pow(rmsd, 8.0) + pow(rmsd, 4.0) + 1;
219  dxgrouprmsd *= -40. / 435;
220 
221  for(size_t j = 0; j < 29; j++){
222  for(size_t k = j + 1; k < 30; k++){
223  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + j]);
224  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[j][k];
225  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
226  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[j][k];
227  }
228  }
229  }
230  }
231 
233  static AlphaRMSDCV* Build(const Json::Value& json, const std::string& path)
234  {
235  Json::ObjectRequirement validator;
236  Json::Value schema;
237  Json::CharReaderBuilder rbuilder;
238  Json::CharReader* reader = rbuilder.newCharReader();
239 
240  reader->parse(JsonSchema::AlphaRMSDCV.c_str(),
241  JsonSchema::AlphaRMSDCV.c_str() + JsonSchema::AlphaRMSDCV.size(),
242  &schema, nullptr);
243  validator.Parse(schema, path);
244 
245  //Validate inputs
246  validator.Validate(json, path);
247  if(validator.HasErrors())
248  throw BuildException(validator.GetErrors());
249 
250  std::vector<int> resids;
251  for(auto& s : json["residue_ids"])
252  resids.push_back(s.asInt());
253  auto reference = json["reference"].asString();
254 
255  double unitconv = json.get("length_unit", 1).asDouble();
256 
257  return new AlphaRMSDCV(resids, reference, unitconv);
258  }
259  };
260 }
Requirements on an object.
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
Collective variable to measure alpha helix secondary structure.
Definition: AlphaRMSDCV.h:44
std::vector< int > ids_only_
Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
Definition: AlphaRMSDCV.h:54
std::vector< std::vector< std::string > > atomids_
Atom IDs for secondary structure calculation: backbone of resids_.
Definition: AlphaRMSDCV.h:51
std::string refpdb_
Name of pdb reference for system.
Definition: AlphaRMSDCV.h:57
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
Definition: AlphaRMSDCV.h:171
std::vector< int > resids_
Residue IDs for secondary structure calculation.
Definition: AlphaRMSDCV.h:48
AlphaRMSDCV(std::vector< int > resids, std::string refpdb, double unitconv)
Constructor.
Definition: AlphaRMSDCV.h:77
double unitconv_
Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10)
Definition: AlphaRMSDCV.h:60
std::vector< double > refdists_
Pairwise distance between backbone atoms for reference structure.
Definition: AlphaRMSDCV.h:63
static AlphaRMSDCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
Definition: AlphaRMSDCV.h:233
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Definition: AlphaRMSDCV.h:103
Exception to be thrown when building the Driver fails.
Abstract class for a collective variable.
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
double val_
Current value of CV.
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
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:48
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:202
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:537
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:186
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:519
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:325
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33