SSAGES  0.1
A MetaDynamics Package
Public Member Functions | Static Public Member Functions | Private Attributes | List of all members
SSAGES::AlphaRMSDCV Class Reference

Collective variable to measure alpha helix secondary structure. More...

#include <AlphaRMSDCV.h>

Inheritance diagram for SSAGES::AlphaRMSDCV:
Inheritance graph
[legend]

Public Member Functions

 AlphaRMSDCV (std::vector< int > resids, std::string refpdb, double unitconv)
 Constructor. More...
 
void Initialize (const Snapshot &snapshot) override
 Initialize necessary variables. More...
 
void Evaluate (const Snapshot &snapshot) override
 Evaluate the CV. More...
 
- Public Member Functions inherited from SSAGES::CollectiveVariable
 CollectiveVariable ()
 Constructor.
 
virtual ~CollectiveVariable ()
 Destructor.
 
virtual void Initialize (const class Snapshot &)
 Initialize CV. More...
 
virtual void Evaluate (const class Snapshot &)=0
 Evaluate CV. More...
 
double GetValue () const
 Get current value of the CV. More...
 
virtual double GetMinimumImage (double) const
 Returns the minimum image of a CV based on the input location. More...
 
virtual double GetPeriodicValue (double location) const
 Apply periodic boundaries to a given value. More...
 
const std::vector< Vector3 > & GetGradient () const
 Get current gradient of the CV. More...
 
const Matrix3GetBoxGradient () const
 Get gradient contribution to box.
 
const std::array< double, 2 > & GetBoundaries ()
 Get CV boundaries. More...
 
virtual double GetDifference (double location) const
 

Static Public Member Functions

static AlphaRMSDCVBuild (const Json::Value &json, const std::string &path)
 Set up collective variable. More...
 
- Static Public Member Functions inherited from SSAGES::CollectiveVariable
static CollectiveVariableBuildCV (const Json::Value &json, const std::string &path)
 Set up collective variable. More...
 

Private Attributes

std::vector< int > resids_
 Residue IDs for secondary structure calculation.
 
std::vector< std::vector< std::string > > atomids_
 Atom IDs for secondary structure calculation: backbone of resids_.
 
std::vector< int > ids_only_
 Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
 
std::string refpdb_
 Name of pdb reference for system.
 
double unitconv_
 Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10)
 
std::vector< double > refdists_
 Pairwise distance between backbone atoms for reference structure.
 

Additional Inherited Members

- Protected Attributes inherited from SSAGES::CollectiveVariable
std::vector< Vector3grad_
 Gradient vector dCv/dxi.
 
Matrix3 boxgrad_
 Gradient w.r.t box vectors dCv/dHij.
 
double val_
 Current value of CV.
 
std::array< double, 2 > bounds_
 Bounds on CV.
 

Detailed Description

Collective variable to measure alpha helix secondary structure.

Following treatment in Pietrucci and Laio, "A Collective Variable for the Efficient Exploration of Protein Beta-Sheet Structures: Application to SH3 and GB1", JCTC, 2009, 5(9): 2197-2201.

Check blocks of six consecutive protein residues for RMSD from reference "ideal" alpha helix structure.

Definition at line 43 of file AlphaRMSDCV.h.

Constructor & Destructor Documentation

SSAGES::AlphaRMSDCV::AlphaRMSDCV ( std::vector< int >  resids,
std::string  refpdb,
double  unitconv 
)
inline

Constructor.

Parameters
residsIDs of residues for calculating secondary structure
refpdbString of pdb filename with atom and residue indices.
unitconvConversion for internal MD length unit: 1 nm is equal to unitconv internal units

Construct an AlphaRMSD CV – calculates alpha helix character by summing pairwise RMSD to an ideal helix structure for all possible 6 residue segments.

Definition at line 77 of file AlphaRMSDCV.h.

Referenced by Build().

77  :
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(unsigned int i = resids[0]; i <= resids[1]; i++){
95  resids_.push_back(i);
96  }
97  }
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< int > resids_
Residue IDs for secondary structure calculation.
Definition: AlphaRMSDCV.h:48
std::string refpdb_
Name of pdb reference for system.
Definition: AlphaRMSDCV.h:57

Here is the caller graph for this function:

Member Function Documentation

static AlphaRMSDCV* SSAGES::AlphaRMSDCV::Build ( const Json::Value &  json,
const std::string &  path 
)
inlinestatic

Set up collective variable.

Parameters
jsonJSON input value.
pathPath for JSON path specification.
Returns
Pointer to the CV built by this function. nullptr in case of unknown error.

Builds a CV from a JSON node. Returns a pointer to the built cv. If an unknown error is encountered, this function will return a nullptr, but generally it will throw a BuildException on failure.

Warning
Object lifetime is the caller's responsibility.

Definition at line 233 of file AlphaRMSDCV.h.

References AlphaRMSDCV(), Json::Requirement::GetErrors(), Json::Requirement::HasErrors(), Json::ObjectRequirement::Parse(), and Json::ObjectRequirement::Validate().

Referenced by SSAGES::CollectiveVariable::BuildCV().

234  {
235  Json::ObjectRequirement validator;
236  Json::Value schema;
237  Json::Reader reader;
238 
239  reader.parse(JsonSchema::AlphaRMSDCV, schema);
240  validator.Parse(schema, path);
241 
242  //Validate inputs
243  validator.Validate(json, path);
244  if(validator.HasErrors())
245  throw BuildException(validator.GetErrors());
246 
247  std::vector<int> resids;
248  for(auto& s : json["residue_ids"])
249  resids.push_back(s.asInt());
250  auto reference = json["reference"].asString();
251 
252  double unitconv = json.get("length_unit", 1).asDouble();
253 
254  return new AlphaRMSDCV(resids, reference, unitconv);
255  }
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Requirements on an object.
AlphaRMSDCV(std::vector< int > resids, std::string refpdb, double unitconv)
Constructor.
Definition: AlphaRMSDCV.h:77
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.

Here is the call graph for this function:

Here is the caller graph for this function:

void SSAGES::AlphaRMSDCV::Evaluate ( const Snapshot snapshot)
inlineoverride

Evaluate the CV.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 171 of file AlphaRMSDCV.h.

References SSAGES::Snapshot::GetCommunicator(), SSAGES::Snapshot::GetLocalIndex(), SSAGES::Snapshot::GetNumAtoms(), SSAGES::Snapshot::GetPositions(), SSAGES::CollectiveVariable::grad_, and SSAGES::CollectiveVariable::val_.

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  unsigned 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  }
std::vector< int > resids_
Residue IDs for secondary structure calculation.
Definition: AlphaRMSDCV.h:48
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
std::vector< int > ids_only_
Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
Definition: AlphaRMSDCV.h:54
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
double val_
Current value of CV.
std::vector< double > refdists_
Pairwise distance between backbone atoms for reference structure.
Definition: AlphaRMSDCV.h:63

Here is the call graph for this function:

void SSAGES::AlphaRMSDCV::Initialize ( const Snapshot snapshot)
inlineoverride

Initialize necessary variables.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 103 of file AlphaRMSDCV.h.

References SSAGES::Snapshot::GetCommunicator(), SSAGES::Snapshot::GetLocalIndices(), and SSAGES::ReadBackbone::GetPdbBackbone().

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  int 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  }
std::vector< std::vector< std::string > > atomids_
Atom IDs for secondary structure calculation: backbone of resids_.
Definition: AlphaRMSDCV.h:51
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< int > resids_
Residue IDs for secondary structure calculation.
Definition: AlphaRMSDCV.h:48
std::string refpdb_
Name of pdb reference for system.
Definition: AlphaRMSDCV.h:57
std::vector< int > ids_only_
Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
Definition: AlphaRMSDCV.h:54
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
std::vector< double > refdists_
Pairwise distance between backbone atoms for reference structure.
Definition: AlphaRMSDCV.h:63
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

Here is the call graph for this function:


The documentation for this class was generated from the following file: