Add residue and side-chain average position functions to the api

api/coot-molecule.cc

@@ -4740,3 +4740,105 @@ coot::molecule_t::assign_sequence(const clipper::Xmap<float> &xmap, const coot::
 }
 
 
+
+
+//! get the residue CA position
+//!
+//! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+std::vector<double>
+coot::molecule_t::get_residue_CA_position(const std::string &cid) const {
+
+   std::vector<double> v;
+   mmdb::Residue *residue_p = cid_to_residue(cid);
+   if (residue_p) {
+      mmdb::Atom **residue_atoms = 0;
+      int n_residue_atoms = 0;
+      residue_p->GetAtomTable(residue_atoms, n_residue_atoms);
+      for (int iat=0; iat<n_residue_atoms; iat++) {
+         mmdb::Atom *at = residue_atoms[iat];
+         if (! at->isTer()) {
+            std::string name = at->GetAtomName();
+            if (name == " CA ") {
+               v.push_back(at->x);
+               v.push_back(at->y);
+               v.push_back(at->z);
+               break;
+            }
+         }
+      }
+   }
+   return v;
+}
+
+//! get the avarge residue position
+//!
+//! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+std::vector<double>
+coot::molecule_t::get_residue_average_position(const std::string &cid) const {
+
+   std::vector<double> v;
+   mmdb::Residue *residue_p = cid_to_residue(cid);
+   if (residue_p) {
+      std::vector<clipper::Coord_orth> atom_positions;
+      mmdb::Atom **residue_atoms = 0;
+      int n_residue_atoms = 0;
+      residue_p->GetAtomTable(residue_atoms, n_residue_atoms);
+      for (int iat=0; iat<n_residue_atoms; iat++) {
+         mmdb::Atom *at = residue_atoms[iat];
+         if (! at->isTer()) {
+            clipper::Coord_orth p = co(at);
+            atom_positions.push_back(p);
+         }
+      }
+      if (! atom_positions.empty()) {
+         clipper::Coord_orth sum(0,0,0);
+         for (const auto &pos : atom_positions)
+            sum += pos;
+         double is = 1.0/static_cast<double>(atom_positions.size());
+         v = {sum.x() * is, sum.y() * is, sum.z() * is};
+      }
+   }
+   return v;
+}
+
+//! get the avarge residue side-chain position
+//!
+//! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+std::vector<double>
+coot::molecule_t::get_residue_sidechain_average_position(const std::string &cid) const {
+
+   std::vector<double> v;
+   mmdb::Residue *residue_p = cid_to_residue(cid);
+   if (residue_p) {
+      std::vector<clipper::Coord_orth> side_chain_positions;
+      std::set<std::string> main_chain_atoms;
+      main_chain_atoms.insert(" CA ");
+      main_chain_atoms.insert(" C  ");
+      main_chain_atoms.insert(" N  ");
+      main_chain_atoms.insert(" O  ");
+      main_chain_atoms.insert(" H  ");
+      main_chain_atoms.insert(" HA ");
+      mmdb::Atom **residue_atoms = 0;
+      int n_residue_atoms = 0;
+      residue_p->GetAtomTable(residue_atoms, n_residue_atoms);
+      for (int iat=0; iat<n_residue_atoms; iat++) {
+         mmdb::Atom *at = residue_atoms[iat];
+         if (! at->isTer()) {
+            std::string atom_name(at->GetAtomName());
+            if (main_chain_atoms.find(atom_name) == main_chain_atoms.end()) {
+               clipper::Coord_orth p = co(at);
+               side_chain_positions.push_back(p);
+            }
+         }
+      }
+      if (! side_chain_positions.empty()) {
+         clipper::Coord_orth sum(0,0,0);
+         for (const auto &pos : side_chain_positions)
+            sum += pos;
+         double is = 1.0/static_cast<double>(side_chain_positions.size());
+         v = {sum.x() * is, sum.y() * is, sum.z() * is};
+      }
+   }
+   return v;
+}
+

api/coot-molecule.hh

@@ -631,6 +631,21 @@ namespace coot {
       //! Get the chains that are related by NCS:
       std::vector<std::vector<std::string> > get_ncs_related_chains() const;
 
+      //! get the residue CA position
+      //!
+      //! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+      std::vector<double> get_residue_CA_position(const std::string &cid) const;
+
+      //! get the avarge residue position
+      //!
+      //! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+      std::vector<double> get_residue_average_position(const std::string &cid) const;
+
+      //! get the avarge residue side-chain position
+      //!
+      //! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+      std::vector<double> get_residue_sidechain_average_position(const std::string &cid) const;
+
       // ----------------------- model bonds
 
       simple_mesh_t get_bonds_mesh(const std::string &mode, protein_geometry *geom,

api/molecules-container-nanobind.cc

@@ -314,8 +314,11 @@ NB_MODULE(chapi, m) {
     .def("get_ramachandran_validation_markup_mesh",&molecules_container_t::get_ramachandran_validation_markup_mesh)
     //Using allow_raw_pointers(). Perhaps suggests we need to do something different from exposing mmdb pointers to JS.
     .def("get_residue",&molecules_container_t::get_residue, nb::rv_policy::reference)
+    .def("get_residue_average_position",&molecules_container_t::get_residue_average_position)
+    .def("get_residue_CA_position",&molecules_container_t::get_residue_CA_position)
     .def("get_residue_name",&molecules_container_t::get_residue_name)
     .def("get_residue_names_with_no_dictionary",&molecules_container_t::get_residue_names_with_no_dictionary)
+    .def("get_residue_sidechain_average_position",&molecules_container_t::get_residue_sidechain_average_position)
     .def("get_residue_using_cid",&molecules_container_t::get_residue_using_cid)
     .def("get_residues_near_residue",&molecules_container_t::get_residues_near_residue)
     .def("get_rotamer_dodecs",&molecules_container_t::get_rotamer_dodecs)

api/molecules-container.cc

@@ -5900,3 +5900,50 @@ molecules_container_t::dictionary_atom_name_map(const std::string &comp_id_1, in
    return m;
 }
 
+
+
+//! get the residue CA position
+//!
+//! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+std::vector<double>
+molecules_container_t::get_residue_CA_position(int imol, const std::string &cid) const {
+
+   std::vector<double> v;
+   if (is_valid_model_molecule(imol)) {
+      v = molecules[imol].get_residue_CA_position(cid);
+   } else {
+      std::cout << "WARNING:: " << __FUNCTION__ << "(): not a valid model molecule " << imol << std::endl;
+   }
+   return v;
+
+}
+
+//! get the avarge residue position
+//!
+//! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+std::vector<double>
+molecules_container_t::get_residue_average_position(int imol, const std::string &cid) const {
+
+   std::vector<double> v;
+   if (is_valid_model_molecule(imol)) {
+      v = molecules[imol].get_residue_average_position(cid);
+   } else {
+      std::cout << "WARNING:: " << __FUNCTION__ << "(): not a valid model molecule " << imol << std::endl;
+   }
+   return v;
+}
+
+//! get the avarge residue side-chain position
+//!
+//! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+std::vector<double>
+molecules_container_t::get_residue_sidechain_average_position(int imol, const std::string &cid) const {
+
+   std::vector<double> v;
+   if (is_valid_model_molecule(imol)) {
+      v = molecules[imol].get_residue_sidechain_average_position(cid);
+   } else {
+      std::cout << "WARNING:: " << __FUNCTION__ << "(): not a valid model molecule " << imol << std::endl;
+   }
+   return v;
+}

api/molecules-container.hh

@@ -1179,6 +1179,21 @@ public:
    std::pair<bool, coot::Cartesian> get_atom_position(int imol, coot::atom_spec_t &atom_spec);
 #endif
 
+   //! get the residue CA position
+   //!
+   //! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+   std::vector<double> get_residue_CA_position(int imol, const std::string &cid) const;
+
+   //! get the avarge residue position
+   //!
+   //! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+   std::vector<double> get_residue_average_position(int imol, const std::string &cid) const;
+
+   //! get the avarge residue side-chain position
+   //!
+   //! @return a vector. The length of the vector is 0 on failure, otherwise it is the x,y,z values
+   std::vector<double> get_residue_sidechain_average_position(int imol, const std::string &cid) const;
+
    //! Get number of atoms
    //!
    //! @param imol is the model molecule index

api/test-molecules-container.cc

@@ -6145,6 +6145,39 @@ int test_dictionary_atom_name_match(molecules_container_t &mc) {
    return status;
 }
 
+int test_average_position_functions(molecules_container_t &mc) {
+
+   auto close_position = [] (const std::vector<double> &p,
+                             const clipper::Coord_orth &r) {
+      double t = 0.9;
+      if (abs(p[0] - r.x()) < t) {
+         if (abs(p[1] - r.y()) < t) {
+            if (abs(p[2] - r.z()) < t) {
+               return true;
+            }
+         }
+      }
+      return false;
+   };
+
+   starting_test(__FUNCTION__);
+   int status = 0;
+
+   int imol = mc.read_pdb(reference_data("moorhen-tutorial-structure-number-1.pdb"));
+   std::vector<double> ap_1 = mc.get_residue_average_position(imol, "//A/15");
+   std::vector<double> ap_2 = mc.get_residue_sidechain_average_position(imol, "//A/15");
+   std::vector<double> ap_3 = mc.get_residue_CA_position(imol, "//A/15");
+
+   int n_pass = 0;
+   if (close_position(ap_1, clipper::Coord_orth(16.4, 10.1, 57.5))) n_pass += 1;
+   if (close_position(ap_2, clipper::Coord_orth(16.2,  8.9, 56.0))) n_pass += 2;
+   if (close_position(ap_3, clipper::Coord_orth(16.4, 11.5, 58.7))) n_pass += 4;
+
+   std::cout << "here with n_pass " << n_pass << std::endl;
+   if (n_pass == 7) status = 1;
+   return status;
+}
+
 
 int test_template(molecules_container_t &mc) {
 
@@ -6463,7 +6496,8 @@ int main(int argc, char **argv) {
          // status += run_test(test_long_name_ligand_cif_merge, "test long name ligand cif merge", mc);
          // status += run_test(test_merge_ligand_and_gemmi_parse_mmcif, "test_merge_ligand_and_gemmi_parse_mmcif", mc);
          // status += run_test(test_delete_two_add_one_using_gemmi, "test_delete_two_add_one_using_gemmi", mc);
-         status += run_test(test_dictionary_atom_name_match, "dictionary atom names match", mc);
+         // status += run_test(test_dictionary_atom_name_match, "dictionary atom names match", mc);
+         status += run_test(test_average_position_functions, "average position functions", mc);
          if (status == n_tests) all_tests_status = 0;
 
          print_results_summary();

rel-todo-gtk4

@@ -1,5 +1,11 @@
 ----
 
+add to api: centre_of_residue() centre_of_sidechain(), CA_of_residue()
+
+Easy function: for Coot API pentakis dodec vertices and convert them to
+cylinders and make a generic display object out of them given centre and radius.
+Allow the user to change the centre and radius using some (python?) gui.
+
 Use alignment text-anchor="middle" for text in molecule svg
 
 For tubes repo, we need loops - reuse the Stuart loops