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

Collective variable to measure parallel beta sheet secondary structure. More...

#include <ParallelBetaRMSDCV.h>

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

Public Member Functions

 ParallelBetaRMSDCV (std::vector< int > resids, std::string refpdb, double unitconv, int mode)
 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 ParallelBetaRMSDCVBuild (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.
 
std::vector< Vector3refalpha_
 Coordinates for reference structure.
 
double unitconv_
 Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10)
 
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< 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 parallel beta sheet 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 2 blocks of 3 consecutive protein residues for RMSD from reference "ideal" parallel beta sheet structure.

Definition at line 43 of file ParallelBetaRMSDCV.h.

Constructor & Destructor Documentation

SSAGES::ParallelBetaRMSDCV::ParallelBetaRMSDCV ( std::vector< int >  resids,
std::string  refpdb,
double  unitconv,
int  mode 
)
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
modeSpecification for inter/intra mode for beta sheet character

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 85 of file ParallelBetaRMSDCV.h.

Referenced by Build().

85  :
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(unsigned int i = resids[0]; i <= resids[1]; i++){
103  resids_.push_back(i);
104  }
105  }
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< 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, 1 for inter, 2 for intra.

Here is the caller graph for this function:

Member Function Documentation

static ParallelBetaRMSDCV* SSAGES::ParallelBetaRMSDCV::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 260 of file ParallelBetaRMSDCV.h.

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

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

261  {
262  Json::ObjectRequirement validator;
263  Json::Value schema;
264  Json::Reader reader;
265 
266  reader.parse(JsonSchema::ParallelBetaRMSDCV, schema);
267  validator.Parse(schema, path);
268 
269  //Validate inputs
270  validator.Validate(json, path);
271  if(validator.HasErrors())
272  throw BuildException(validator.GetErrors());
273 
274  std::vector<int> resids;
275  for(auto& s : json["residue_ids"])
276  resids.push_back(s.asInt());
277  auto reference = json["reference"].asString();
278 
279  double unitconv = json.get("length_unit", 1).asDouble();
280 
281  int mode = json.get("mode", 0).asInt();
282 
283  return new ParallelBetaRMSDCV(resids, reference, unitconv, mode);
284  }
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
ParallelBetaRMSDCV(std::vector< int > resids, std::string refpdb, double unitconv, int mode)
Constructor.
Requirements on an object.
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::ParallelBetaRMSDCV::Evaluate ( const Snapshot snapshot)
inlineoverride

Evaluate the CV.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 179 of file ParallelBetaRMSDCV.h.

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

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  unsigned 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  }
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_.
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
std::vector< double > refdists_
Pairwise distance between backbone atoms for reference structure.
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< 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.

Here is the call graph for this function:

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

Initialize necessary variables.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 111 of file ParallelBetaRMSDCV.h.

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

112  {
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  int 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  }
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< int > resids_
Residue IDs for secondary structure calculation.
std::vector< int > ids_only_
Atom IDs only, for secondary structure calculation: backbone atoms of resids_.
std::vector< double > refdists_
Pairwise distance between backbone atoms 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
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: