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

Collective variable to measure antiparallel beta secondary structure. More...

#include <AntiBetaRMSDCV.h>

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

Public Member Functions

 AntiBetaRMSDCV (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 AntiBetaRMSDCVBuild (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 and chain 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)
 
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 antiparallel beta 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" antiparallel beta sheet structure.

Definition at line 43 of file AntiBetaRMSDCV.h.

Constructor & Destructor Documentation

SSAGES::AntiBetaRMSDCV::AntiBetaRMSDCV ( 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 AntibetaRMSD CV – calculates anti-beta sheet character by summing pairwise RMSD to an ideal anti-beta sheet structure for all possible 6 residue segments.

Definition at line 80 of file AntiBetaRMSDCV.h.

Referenced by Build().

80  :
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(unsigned int i = resids[0]; i <= resids[1]; i++){
98  resids_.push_back(i);
99  }
100  }
double unitconv_
Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10) ...
std::vector< int > resids_
Residue IDs for secondary structure calculation.
std::string refpdb_
Name of pdb reference for system.
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 AntiBetaRMSDCV* SSAGES::AntiBetaRMSDCV::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 256 of file AntiBetaRMSDCV.h.

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

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

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

Evaluate the CV.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 174 of file AntiBetaRMSDCV.h.

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

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  unsigned 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  }
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.
std::vector< int > resids_
Residue IDs for secondary structure calculation.
std::vector< std::vector< std::string > > atomids_
Atom IDs and chain IDs for secondary structure calculation: backbone of resids_.
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
double val_
Current value of CV.
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 call graph for this function:

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

Initialize necessary variables.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 106 of file AntiBetaRMSDCV.h.

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

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  int 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  }
double unitconv_
Length unit conversion: convert 1 nm to your internal MD units (ex. if using angstroms use 10) ...
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< int > resids_
Residue IDs for secondary structure calculation.
std::vector< std::vector< std::string > > atomids_
Atom IDs and chain IDs for secondary structure calculation: backbone of resids_.
std::string refpdb_
Name of pdb reference for system.
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: