ap_molecule_diffusion

ap_molecule_diffusion calculates average displacement of a small molecule as a function of time over a trajectory

If a multi-model PDB file was given, the program prints contact count observed in all models

USAGE:
ap_molecule_diffusion  trajectory.pdb HOH box_side

where trajectory.pdb is the input file multimodel-PDB file HOH is the PDB-id of molecules for which the displacement will be evaluated

Categories:

  • core::data::basic::Vec3Cubic

Input files:

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
#include <iostream>

#include <core/index.hh>
#include <core/data/io/Pdb.hh>
#include <core/calc/statistics/OnlineStatistics.hh>
#include <utils/exit.hh>

std::string program_info = R"(

ap_molecule_diffusion calculates average displacement of a small molecule as a function of time over a trajectory

If a multi-model PDB file was given, the program prints contact count observed in all models

USAGE:
    ap_molecule_diffusion  trajectory.pdb HOH box_side

where trajectory.pdb is the input file multimodel-PDB file HOH is the PDB-id of molecules for which the displacement
will be evaluated

)";

/** @brief Calculates  average displacement of a small molecule as a function of time over a trajectory
 *
 * CATEGORIES: core::data::basic::Vec3Cubic
 * KEYWORDS: PDB input; simulation
 */
int main(const int argc, const char *argv[]) {

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

  double L = utils::from_string<double>(argv[3]); // The third parameter is the box width (in Angstroms)
  core::data::basic::Vec3Cubic::set_box_len(L);
  core::data::io::Pdb reader(argv[1]); // --- file name (PDB format, may be gzip-ped)

  core::index4 n_atoms = reader.count_atoms(0);
  core::index4 n_frmes = reader.count_models();
  core::index4 t_max = n_frmes / 5;
  std::vector<std::vector<core::data::basic::Vec3Cubic>> v;

  // ---------- Load coordinates to memory ----------
  for (int i_start = 0; i_start < n_frmes; ++i_start) {
    std::vector<core::data::basic::Vec3Cubic> vi(n_atoms);
    reader.fill_structure(i_start, vi);
    v.push_back(vi);
  }

  // ---------- Calculate displacement ----------
  std::vector<core::calc::statistics::OnlineStatistics> avg(t_max);
  for (size_t i_start = 0; i_start < n_frmes - t_max - 1; ++i_start) {
    const std::vector<core::data::basic::Vec3Cubic> &v0 = v[i_start];
    for (size_t i_t = 1; i_t <= t_max; ++i_t) {
      const std::vector<core::data::basic::Vec3Cubic> &vi = v[i_start + i_t];
      for (size_t i_atom = 0; i_atom < n_atoms; ++i_atom)
        avg[i_t - 1](sqrt(v0[i_atom].closest_distance_square_to(vi[i_atom])));
    }
  }

  for (size_t i_t = 1; i_t <= t_max; ++i_t) {
    std::cout << utils::string_format("%5d %f %f\n", i_t, avg[i_t - 1].avg(), sqrt(avg[i_t - 1].var()));
  }
}
../_images/file_icon.png