ap_Crmsd

ap_Crmsd calculates crmsd value on C-alpha coordinates. The program prints just the crmsd value.

USAGE:
ap_Crmsd structureA.pdb [structureB.pdb]

If two structures are provided, the program calculates crmsd between the first model of structure A and the first model of structure B. If only one input PDB file is given, crsmd is computed for every pair of models found in the input file (each-vs-each)

Keywords:

Categories:

  • core/calc/structural/transformations/Crmsd

Output files:

Program source:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
#include <iostream>

#include <core/data/io/Pdb.hh>
#include <core/calc/structural/transformations/Crmsd.hh>

#include <utils/exit.hh>

std::string program_info = R"(

ap_Crmsd calculates crmsd value on C-alpha coordinates. The program prints just the crmsd value.
USAGE:
    ap_Crmsd structureA.pdb [structureB.pdb]

If two structures are provided, the program calculates crmsd between the first model of structure A
and the first model of structure B. If only one input PDB file is given, crsmd is computed for every pair
of models found in the input file (each-vs-each)

)";

/** @brief Calculates crmsd value on C-alpha coordinates. The program prints just the crmsd value.
 *
 * CATEGORIES: core/calc/structural/transformations/Crmsd
 * KEYWORDS:   PDB input; crmsd
 * IMG: ap_Crmsd_deepteal_brown_1.png
 * IMG_ALT: 2GB1 model structure superimposed on the native, crmsd = 4.93952
 */
int main(const int argc, const char *argv[]) {

  if(argc < 2) utils::exit_OK_with_message(program_info); // --- complain about missing program parameter

  using core::data::basic::Vec3;
  using namespace core::calc::structural::transformations;

  Crmsd<std::vector<Vec3>,std::vector<Vec3>> rms;

  if(argc==2) { // --- The case of each-vs-each calculations beteen models of a single PDB file
    core::data::io::Pdb q_reader(argv[1],core::data::io::is_ca, false);
    core::data::structural::Structure_SP q_strctr = q_reader.create_structure(0); // --- create a structure object

    std::vector<std::vector<core::data::basic::Vec3>> models(q_reader.count_models());
    for(int i=0;i<q_reader.count_models();++i) {
      models[i].resize(q_strctr->count_atoms());
      q_reader.fill_structure(i,models[i]);
      for (int j = 0; j < i; ++j)
        std::cout << i<<" "<<j << " "<<rms.crmsd(models[i], models[j],models[j].size()) << "\n";
    }
  } else { // --- The case when two PDB files are given
    core::data::io::Pdb q_reader(argv[1], core::data::io::is_ca, false);
    core::data::structural::Structure_SP q_strctr = q_reader.create_structure(0); // --- create a structure object

    core::data::io::Pdb t_reader(argv[2], core::data::io::is_ca, false);
    core::data::structural::Structure_SP t_strctr = t_reader.create_structure(0); // --- create a structure object

    if (q_strctr->count_atoms() != t_strctr->count_atoms())
      utils::exit_OK_with_message("The two structures have different number of CA atoms!\n");

    std::vector<Vec3> q, t;
    for (auto atom_it = q_strctr->first_atom(); atom_it != q_strctr->last_atom(); ++atom_it) q.push_back(**atom_it);
    for (auto atom_it = t_strctr->first_atom(); atom_it != t_strctr->last_atom(); ++atom_it) t.push_back(**atom_it);
    std::cout << "crmsd: " << rms.crmsd(q, t, q_strctr->count_atoms()) << "\n";
  }
}
../_images/ap_Crmsd_deepteal_brown_1.png