# 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

## 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 #include #include #include #include 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(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> v; // ---------- Load coordinates to memory ---------- for (int i_start = 0; i_start < n_frmes; ++i_start) { std::vector vi(n_atoms); reader.fill_structure(i_start, vi); v.push_back(vi); } // ---------- Calculate displacement ---------- std::vector avg(t_max); for (size_t i_start = 0; i_start < n_frmes - t_max - 1; ++i_start) { const std::vector &v0 = v[i_start]; for (size_t i_t = 1; i_t <= t_max; ++i_t) { const std::vector &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())); } }