SSAGES  0.9.3
Software Suite for Advanced General Ensemble Simulations
AntiBetaRMSDCV.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  int mode_;
64 
66  std::vector<double> refdists_;
67 
68  public:
70 
80  AntiBetaRMSDCV(std::vector<int> resids, std::string refpdb, double unitconv, int mode) :
81  resids_(resids), refpdb_(refpdb), unitconv_(unitconv), mode_(mode)
82  {
83  if(resids_.size() != 2 ){
84  throw std::invalid_argument("AntiBetaRMSDCV: Input must designate range of residues with 2 residue numbers.");
85  }
86 
87  resids_.clear();
88 
89  if(resids[0] >= resids[1]){
90  throw std::invalid_argument("AntiBetaRMSDCV: Input must list lower residue index first: please reverse residue range.");
91  } else if(resids[1] - resids[0] < 5){
92  throw std::invalid_argument("AntiBetaRMSDCV: Residue range must span at least 6 residues for secondary structure calculation.");
93  }
94 
95  std::cout << "AntiBetaRMSDCV: Calculating antiparallel beta sheet character from residue " << resids[0] << " to " << resids[1] << "." << std::endl;
96 
97  for(int i = resids[0]; i <= resids[1]; i++){
98  resids_.push_back(i);
99  }
100  }
101 
103 
106  void Initialize(const Snapshot& snapshot) override
107  {
109  ids_only_.resize(atomids_[0].size());
110 
111  // check to find all necessary atoms
112  std::transform(atomids_[0].begin(), atomids_[0].end(), ids_only_.begin(), [](std::string val) {return std::stoi(val);});
113  using std::to_string;
114  std::vector<int> found;
115  snapshot.GetLocalIndices(ids_only_, &found);
116  size_t nfound = found.size();
117  MPI_Allreduce(MPI_IN_PLACE, &nfound, 1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
118  if(nfound != resids_.size() * 5)
119  throw BuildException({
120  "AntiBetaRMSDCV: Expected to find " +
121  to_string(resids_.size() * 5) +
122  " atoms, but only found " +
123  to_string(nfound) + "."
124  });
125 
126  // reference 'ideal' antiparallel beta sheet structure in nanometers
127  // Reference structure taken from values used in Plumed 2.2.0
128  std::vector<Vector3> refalpha;
129  refalpha.push_back(unitconv_ * Vector3{ .2263, -.3795, .1722}); // N i
130  refalpha.push_back(unitconv_ * Vector3{ .2493, -.2426, .2263}); // CA
131  refalpha.push_back(unitconv_ * Vector3{ .3847, -.1838, .1761}); // CB
132  refalpha.push_back(unitconv_ * Vector3{ .1301, -.1517, .1921}); // C
133  refalpha.push_back(unitconv_ * Vector3{ .0852, -.1504, .0739}); // O
134  refalpha.push_back(unitconv_ * Vector3{ .0818, -.0738, .2917}); // N i+1
135  refalpha.push_back(unitconv_ * Vector3{-.0299, .0243, .2748}); // CA
136  refalpha.push_back(unitconv_ * Vector3{-.1421, -.0076, .3757}); // CB
137  refalpha.push_back(unitconv_ * Vector3{ .0273, .1680, .2854}); // C
138  refalpha.push_back(unitconv_ * Vector3{ .0902, .1993, .3888}); // O
139  refalpha.push_back(unitconv_ * Vector3{ .0119, .2532, .1813}); // N i+2
140  refalpha.push_back(unitconv_ * Vector3{ .0683, .3916, .1680}); // CA
141  refalpha.push_back(unitconv_ * Vector3{ .1580, .3940, .0395}); // CB
142  refalpha.push_back(unitconv_ * Vector3{-.0394, .5011, .1630}); // C
143  refalpha.push_back(unitconv_ * Vector3{-.1459, .4814, .0982}); // O
144  refalpha.push_back(unitconv_ * Vector3{-.2962, .3559, -.1359}); // N j - 2
145  refalpha.push_back(unitconv_ * Vector3{-.2439, .2526, -.2287}); // CA
146  refalpha.push_back(unitconv_ * Vector3{-.1189, .3006, -.3087}); // CB
147  refalpha.push_back(unitconv_ * Vector3{-.2081, .1231, -.1520}); // C
148  refalpha.push_back(unitconv_ * Vector3{-.1524, .1324, -.0409}); // O
149  refalpha.push_back(unitconv_ * Vector3{-.2326, .0037, -.2095}); // N j - 1
150  refalpha.push_back(unitconv_ * Vector3{-.1858, -.1269, -.1554}); // CA
151  refalpha.push_back(unitconv_ * Vector3{-.3053, -.2199, -.1291}); // CB
152  refalpha.push_back(unitconv_ * Vector3{-.0869, -.1949, -.2512}); // C
153  refalpha.push_back(unitconv_ * Vector3{-.1255, -.2070, -.3710}); // O
154  refalpha.push_back(unitconv_ * Vector3{ .0326, -.2363, -.2072}); // N j
155  refalpha.push_back(unitconv_ * Vector3{ .1405, -.2992, -.2872}); // CA
156  refalpha.push_back(unitconv_ * Vector3{ .2699, -.2129, -.2917}); // CB
157  refalpha.push_back(unitconv_ * Vector3{ .1745, -.4399, -.2330}); // C
158  refalpha.push_back(unitconv_ * Vector3{ .1899, -.4545, -.1102}); // O
159 
160  // calculate all necessary pairwise distances for reference structure
161  std::fill(refdists_.begin(), refdists_.end(), 0.0);
162  refdists_.resize(900, 0.0);
163  for(size_t k = 0; k < 29; k++){
164  for(size_t h = k + 1; h < 30; h++){
165  refdists_[29 * k + h] = (refalpha[k] - refalpha[h]).norm();
166  }
167  }
168  }
169 
171 
174  void Evaluate(const Snapshot& snapshot) override
175  {
176  // need atom positions for all atoms
177  const auto& pos = snapshot.GetPositions();
178  auto& comm = snapshot.GetCommunicator();
179 
180  unsigned int resgroups = resids_.size() - 2;
181  double rmsd, dist_norm, dxgrouprmsd;
182  int localidx;
183  Vector3 dist_xyz, dist_ref;
184  std::vector<Vector3> refxyz;
185  std::vector< std::vector< Vector3 > > deriv(30, std::vector<Vector3>(30, Vector3{0,0,0}));
186  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
187  grad_.resize(snapshot.GetNumAtoms(), Vector3{0,0,0});
188  val_ = 0.0;
189 
190  for(size_t i = 0; i < resgroups - 3; i++){
191  for(size_t j = i + 3; j < resgroups; j++){
192  // mode: 0 for all, 1 for inter, 2 for intra
193  if((mode_ == 0) || (mode_ == 1 && atomids_[1][5 * j] != atomids_[1][5 * i]) || (mode_ == 2 && atomids_[1][5 * j] == atomids_[1][5 * i])){
194  rmsd = 0.0;
195  std::fill(refxyz.begin(), refxyz.end(), Vector3{0,0,0});
196  refxyz.resize(30, Vector3{0,0,0});
197 
198  for(size_t k = 0; k < 15; k++){
199  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
200  if(localidx != -1) refxyz[k] = pos[localidx];
201  }
202 
203  for(size_t k = 0; k < 15; k++){
204  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + k]);
205  if(localidx != -1) refxyz[k + 15] = pos[localidx];
206  }
207 
208  MPI_Allreduce(MPI_IN_PLACE, refxyz.data(), 90, MPI_DOUBLE, MPI_SUM, comm);
209 
210  // sum over all pairs to calculate CV
211  for(size_t k = 0; k < 29; k++){
212  for(size_t h = k + 1; h < 30; h++){
213  dist_xyz = refxyz[k] - refxyz[h];
214  dist_norm = dist_xyz.norm() - refdists_[29 * k + h];
215  rmsd += dist_norm * dist_norm;
216  deriv[k][h] = dist_xyz * dist_norm / dist_xyz.norm();
217  }
218  }
219 
220  rmsd = pow(rmsd / 435, 0.5) / 0.1;
221  val_ += (1 - pow(rmsd, 8.0)) / (1 - pow(rmsd, 12.0));
222 
223  dxgrouprmsd = pow(rmsd, 11.0) + pow(rmsd, 7.0);
224  dxgrouprmsd /= pow(rmsd, 8.0) + pow(rmsd, 4.0) + 1;
225  dxgrouprmsd /= pow(rmsd, 8.0) + pow(rmsd, 4.0) + 1;
226  dxgrouprmsd *= -40. / 435;
227 
228  for(size_t k = 0; k < 15; k++){
229  for(size_t h = k + 1; h < 15; h++){
230  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
231  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h];
232  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + h]);
233  if (localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h];
234  }
235  for(size_t h = 0; h < 15; h++){
236  localidx = snapshot.GetLocalIndex(ids_only_[5 * i + k]);
237  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h + 15];
238  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + h]);
239  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k][h + 15];
240  }
241  }
242  for(size_t k = 0; k < 14; k++){
243  for(size_t h = k + 1; h < 15; h++){
244  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + k]);
245  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k + 15][h + 15];
246  localidx = snapshot.GetLocalIndex(ids_only_[5 * j + h]);
247  if(localidx != -1) grad_[localidx] += dxgrouprmsd * deriv[k + 15][h + 15];
248  }
249  }
250  }
251  }
252  }
253  }
254 
256  static AntiBetaRMSDCV* Build(const Json::Value& json, const std::string& path)
257  {
258  Json::ObjectRequirement validator;
259  Json::Value schema;
260  Json::CharReaderBuilder rbuilder;
261  Json::CharReader* reader = rbuilder.newCharReader();
262 
263  reader->parse(JsonSchema::AntiBetaRMSDCV.c_str(),
264  JsonSchema::AntiBetaRMSDCV.c_str() + JsonSchema::AntiBetaRMSDCV.size(),
265  &schema, nullptr);
266  validator.Parse(schema, path);
267 
268  //Validate inputs
269  validator.Validate(json, path);
270  if(validator.HasErrors())
271  throw BuildException(validator.GetErrors());
272 
273  std::vector<int> resids;
274  for(auto& s : json["residue_ids"])
275  resids.push_back(s.asInt());
276  auto reference = json["reference"].asString();
277 
278  double unitconv = json.get("length_unit", 1).asDouble();
279 
280  int mode = json.get("mode", 0).asInt();
281 
282  return new AntiBetaRMSDCV(resids, reference, unitconv, mode);
283  }
284  };
285 }
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 antiparallel beta secondary structure.
std::vector< double > refdists_
Pairwise distance between backbone atoms for reference structure.
std::vector< int > resids_
Residue IDs for secondary structure calculation.
int mode_
Specify whether to calculate beta sheet character in intra or inter mode: 0 for either,...
static AntiBetaRMSDCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
double unitconv_
Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10)
std::string refpdb_
Name of pdb reference for system.
std::vector< std::vector< std::string > > atomids_
Atom IDs and chain IDs for secondary structure calculation: backbone of resids_.
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
std::vector< int > ids_only_
Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
AntiBetaRMSDCV(std::vector< int > resids, std::string refpdb, double unitconv, int mode)
Constructor.
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