diff --git a/api/coot-molecule.cc b/api/coot-molecule.cc index e900a4f71..dbba63457 100644 --- a/api/coot-molecule.cc +++ b/api/coot-molecule.cc @@ -4740,3 +4740,105 @@ coot::molecule_t::assign_sequence(const clipper::Xmap &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 +coot::molecule_t::get_residue_CA_position(const std::string &cid) const { + + std::vector 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; iatisTer()) { + 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 +coot::molecule_t::get_residue_average_position(const std::string &cid) const { + + std::vector v; + mmdb::Residue *residue_p = cid_to_residue(cid); + if (residue_p) { + std::vector atom_positions; + mmdb::Atom **residue_atoms = 0; + int n_residue_atoms = 0; + residue_p->GetAtomTable(residue_atoms, n_residue_atoms); + for (int iat=0; iatisTer()) { + 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(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 +coot::molecule_t::get_residue_sidechain_average_position(const std::string &cid) const { + + std::vector v; + mmdb::Residue *residue_p = cid_to_residue(cid); + if (residue_p) { + std::vector side_chain_positions; + std::set 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; iatisTer()) { + 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(side_chain_positions.size()); + v = {sum.x() * is, sum.y() * is, sum.z() * is}; + } + } + return v; +} + diff --git a/api/coot-molecule.hh b/api/coot-molecule.hh index 34a6f4b1b..f6c16f957 100644 --- a/api/coot-molecule.hh +++ b/api/coot-molecule.hh @@ -631,6 +631,21 @@ namespace coot { //! Get the chains that are related by NCS: std::vector > 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 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 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 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, diff --git a/api/molecules-container-nanobind.cc b/api/molecules-container-nanobind.cc index 952dbb79c..c40adfed2 100644 --- a/api/molecules-container-nanobind.cc +++ b/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) diff --git a/api/molecules-container.cc b/api/molecules-container.cc index 9175c4ae4..c5192779e 100644 --- a/api/molecules-container.cc +++ b/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 +molecules_container_t::get_residue_CA_position(int imol, const std::string &cid) const { + + std::vector 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 +molecules_container_t::get_residue_average_position(int imol, const std::string &cid) const { + + std::vector 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 +molecules_container_t::get_residue_sidechain_average_position(int imol, const std::string &cid) const { + + std::vector 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; +} diff --git a/api/molecules-container.hh b/api/molecules-container.hh index 75f041cb7..05c5c8672 100644 --- a/api/molecules-container.hh +++ b/api/molecules-container.hh @@ -1179,6 +1179,21 @@ public: std::pair 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 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 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 get_residue_sidechain_average_position(int imol, const std::string &cid) const; + //! Get number of atoms //! //! @param imol is the model molecule index diff --git a/api/test-molecules-container.cc b/api/test-molecules-container.cc index 5fd1acddc..64e8ad003 100644 --- a/api/test-molecules-container.cc +++ b/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 &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 ap_1 = mc.get_residue_average_position(imol, "//A/15"); + std::vector ap_2 = mc.get_residue_sidechain_average_position(imol, "//A/15"); + std::vector 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(); diff --git a/rel-todo-gtk4 b/rel-todo-gtk4 index 1b4a972f0..8dba878f7 100644 --- a/rel-todo-gtk4 +++ b/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