PYMOL-5418: Add USAlign implementation for structural alignment

Author: JarrettSJohnson

layer2/USalign.cpp

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include <glm/glm.hpp>
+
+namespace pymol::usalign {
+
+struct Superposition {
+  glm::dvec3 translation{0.0};
+  glm::dmat3 rotation{1.0};
+};
+
+struct TMAlignResult {
+  double tm_score_mobile = 0.0;  // normalized by mobile length
+  double tm_score_target = 0.0;  // normalized by target length
+  double d0_target = 0.0;        // d0 used for target normalization
+  double d0_mobile = 0.0;        // d0 used for mobile normalization
+  double rmsd = 0.0;
+  int aligned_length = 0;        // aligned pairs within distance cutoff
+  double seq_identity = 0.0;
+  Superposition transform;
+  std::vector<int> mobile_indices;  // paired residue indices into mobile CA array
+  std::vector<int> target_indices;  // paired residue indices into target CA array
+  std::string seq_mobile;           // alignment string for mobile
+  std::string seq_target;           // alignment string for target
+  std::string seq_match;            // ':' close, '.' far, ' ' gap
+};
+
+// DP workspace — single allocation reused across all seeds
+struct DPWorkspace {
+  std::vector<double> score_flat;
+  std::vector<double> val_flat;
+  std::vector<char> path_flat;
+  int rows = 0;
+  int cols = 0;
+
+  // Scratch buffers for TMscore8_search and scoring.
+  // xtm, ytm, r1, r2 are sized to min(xlen, ylen); xt is sized to xlen.
+  // All score_fun8 n_cut values are bounded by min(xlen, ylen).
+  std::vector<glm::dvec3> xtm, ytm, xt, r1, r2;
+
+  void resize(int xlen, int ylen);
+
+  double& score(int i, int j) { return score_flat[i * cols + j]; }
+  double& val(int i, int j) { return val_flat[i * cols + j]; }
+  bool path(int i, int j) const { return path_flat[i * cols + j] != 0; }
+  void set_path(int i, int j, bool v) { path_flat[i * cols + j] = v ? 1 : 0; }
+};
+
+/**
+ * Perform TM-score structural alignment between two protein structures.
+ * 
+ * @param target_ca Target structure CA coordinates (remains fixed)
+ * @param mobile_ca Mobile structure CA coordinates (will be aligned to target)
+ * @param target_seq Target sequence (single-letter amino acid codes)
+ * @param mobile_seq Mobile sequence (single-letter amino acid codes)
+ * @param fast Use fast mode with fewer iterations (default: false)
+ * @return TMAlignResult containing TM-scores, RMSD, alignment, and transform
+ * 
+ * @note TM-score ranges from 0 to 1; score > 0.5 indicates same fold
+ * @note Complexity: O(n²) where n = min(target_len, mobile_len)
+ * @note may use a lot of memory for structures above 10K residues.
+ */
+TMAlignResult TMalign(
+    const std::vector<glm::dvec3>& target_ca,
+    const std::vector<glm::dvec3>& mobile_ca,
+    const std::string& target_seq,
+    const std::string& mobile_seq,
+    bool fast = false);
+
+} // namespace pymol::usalign

layer2/USalign.h

@@ -77,6 +77,7 @@
 #include "SceneRay.h"
 #include "ScrollBar.h"
 #include "SculptCache.h"
+#include "Seeker.h"
 #include "Selector.h"
 #include "Seq.h"
 #include "Setting.h"
@@ -114,6 +115,7 @@
 #include "ce_types.h"
 #endif
 
+#include <glm/gtc/quaternion.hpp>
 #include <glm/gtc/type_ptr.hpp>
 #include <glm/vec3.hpp>
 
@@ -17560,3 +17562,197 @@ pymol::Result<std::unordered_set<const pymol::CObject*>> ExecutiveGetObjectDeps(
   obj_set.erase(&obj);
   return obj_set;
 }
+
+/**
+ * Run TM-align on two selections and return results.
+ *
+ * @param mobile_sele mobile selection (will be transformed)
+ * @param target_sele target selection (stays fixed)
+ * @param mobile_state state of mobile selection (0-based)
+ * @param target_state state of target selection (0-based)
+ * @param quiet suppress output
+ * @param transform apply superposition transform
+ * @param oname name for alignment object (empty = don't create)
+ * @param fast use fast mode (fewer iterations)
+ */
+pymol::Result<pymol::usalign::TMAlignResult> ExecutiveUSalign(PyMOLGlobals* G,
+    const char* mobile_sele, const char* target_sele, int mobile_state,
+    int target_state, int quiet, int transform, const char* oname, int fast)
+{
+  // Resolve selections
+  auto sele_mobile = SelectorIndexByName(G, mobile_sele);
+  if (sele_mobile < 0)
+    return pymol::make_error("Invalid mobile selection: ", mobile_sele);
+
+  auto sele_target = SelectorIndexByName(G, target_sele);
+  if (sele_target < 0)
+    return pymol::make_error("Invalid target selection: ", target_sele);
+
+  // Extract CA coordinates and sequences
+  struct ResidueInfo {
+    glm::dvec3 coord;
+    char seq_char;
+    AtomInfoType* ai;
+  };
+
+  auto extract_ca = [&](SelectorID_t sele,
+                        int state) -> std::vector<ResidueInfo> {
+    std::vector<ResidueInfo> residues;
+    SeleCoordIterator iter(G, sele, state);
+    while (iter.next()) {
+      auto* ai = iter.getAtomInfo();
+      if (ai->flags & cAtomFlag_guide) {
+        float* c = iter.getCoord();
+        ResidueInfo ri;
+        ri.coord = glm::dvec3(c[0], c[1], c[2]);
+        ri.seq_char = SeekerGetAbbr(G, LexStr(G, ai->resn), 'O', 'X');
+        ri.ai = ai;
+        residues.push_back(ri);
+      }
+    }
+    return residues;
+  };
+
+  auto mobile_res = extract_ca(sele_mobile, mobile_state);
+  auto target_res = extract_ca(sele_target, target_state);
+
+  if (mobile_res.size() < 3) {
+    return pymol::make_error("Mobile selection has fewer than 3 guide atoms (",
+        mobile_res.size(), ")");
+  }
+  if (target_res.size() < 3) {
+    return pymol::make_error("Target selection has fewer than 3 guide atoms (",
+        target_res.size(), ")");
+  }
+
+  // Build coordinate vectors and sequences
+  std::vector<glm::dvec3> mobile_ca, target_ca;
+  std::string mobile_seq, target_seq;
+  mobile_ca.reserve(mobile_res.size());
+  target_ca.reserve(target_res.size());
+
+  for (const auto& r : mobile_res) {
+    mobile_ca.push_back(r.coord);
+    mobile_seq.push_back(r.seq_char);
+  }
+  for (const auto& r : target_res) {
+    target_ca.push_back(r.coord);
+    target_seq.push_back(r.seq_char);
+  }
+
+  // Run TM-align
+  auto result = pymol::usalign::TMalign(
+      target_ca, mobile_ca, target_seq, mobile_seq, fast != 0);
+
+  if (result.aligned_length < 1) {
+    return pymol::make_error("TM-align failed to find any alignment");
+  }
+
+  // Print results
+  if (!quiet) {
+    PRINTFB(G, FB_Executive, FB_Results)
+    " USalign: TM-score= %6.4f (normalized by target, N=%d, d0=%.2f)\n",
+        result.tm_score_target, static_cast<int>(target_ca.size()),
+        result.d0_target ENDFB(G);
+    PRINTFB(G, FB_Executive, FB_Results)
+    " USalign: TM-score= %6.4f (normalized by mobile, N=%d, d0=%.2f)\n",
+        result.tm_score_mobile, static_cast<int>(mobile_ca.size()),
+        result.d0_mobile ENDFB(G);
+    PRINTFB(G, FB_Executive, FB_Results)
+    " USalign: Aligned length= %d, RMSD= %5.2f, Seq_ID=n_identical/n_aligned= "
+    "%4.3f\n",
+        result.aligned_length, result.rmsd, result.seq_identity ENDFB(G);
+  }
+
+  // Apply transform to mobile object
+  if (transform) {
+    // Convert double-precision Superposition to float TTT
+    // USalign convention: y_aligned = R * x + t
+    // where x = mobile coords, y = target coords
+    // The rotation R and translation t transform mobile -> target space
+
+    const auto& sup = result.transform;
+
+    // Build a legacy-style 16-float TTT matrix
+    // TTT format: [R00 R01 R02 pre_x] [R10 R11 R12 pre_y]
+    //             [R20 R21 R22 pre_z] [tx  ty  tz  1]
+    // where pre is the pre-translation (origin), and t is post-translation
+    // For a simple rotation+translation (no origin): pre=0, R=rotation,
+    // t=translation
+
+    glm::mat3 rot_f(sup.rotation);
+    glm::quat q = glm::quat_cast(rot_f);
+    glm::vec3 t(sup.translation);
+
+    // Create TTT: pretranslate=0, rotate=q, posttranslate=t
+    pymol::TTT ttt(glm::vec3(0.0f), q, t);
+
+    // Convert to legacy float[16] format for ExecuteCombineObjectTTT
+    auto legacy = pymol::TTT::as_pymol_2_legacy(ttt);
+    float tttf[16];
+    std::memcpy(tttf, glm::value_ptr(legacy), 16 * sizeof(float));
+
+    // Follow the same pattern as ExecutiveAlign:
+    // 1. Copy target's TTT and state matrix to mobile (reset to same frame)
+    // 2. Combine the alignment transform (reverse_order=true)
+    // Note: Only the first object in the mobile selection is transformed,
+    // matching ExecutiveAlign behavior for multi-object selections.
+    ObjectMolecule* mobile_obj = SelectorGetFirstObjectMolecule(G, sele_mobile);
+    ObjectMolecule* target_obj =
+        SelectorGetSingleObjectMolecule(G, sele_target);
+    if (mobile_obj && target_obj) {
+      ExecutiveMatrixCopy(G, target_obj->Name, mobile_obj->Name, 1, 1,
+          target_state, mobile_state, false, 0, quiet);
+      ExecutiveMatrixCopy(G, target_obj->Name, mobile_obj->Name, 2, 2,
+          target_state, mobile_state, false, 0, quiet);
+      ExecutiveCombineObjectTTT(G, mobile_obj->Name, tttf, true, -1);
+    }
+  }
+
+  // Create alignment object
+  if (oname && oname[0]) {
+    int align_state = target_state;
+    if (align_state < 0) {
+      align_state = SceneGetState(G);
+    }
+
+    ObjectMolecule* trg_obj = SelectorGetSingleObjectMolecule(G, sele_target);
+    ObjectMolecule* mob_obj = SelectorGetFirstObjectMolecule(G, sele_mobile);
+
+    if (trg_obj && mob_obj) {
+      int n_pair = result.aligned_length;
+      pymol::vla<int> align_vla(n_pair * 3);
+      int* id_p = align_vla.data();
+
+      for (int k = 0; k < n_pair; k++) {
+        int mi = result.mobile_indices[k];
+        int ti = result.target_indices[k];
+        if (mi < static_cast<int>(mobile_res.size()) &&
+            ti < static_cast<int>(target_res.size())) {
+          id_p[0] = AtomInfoCheckUniqueID(G, target_res[ti].ai);
+          id_p[1] = AtomInfoCheckUniqueID(G, mobile_res[mi].ai);
+          id_p[2] = 0;
+          id_p += 3;
+        }
+      }
+
+      ObjectAlignment* obj = nullptr;
+      {
+        pymol::CObject* execObj = ExecutiveFindObjectByName(G, oname);
+        if (execObj && execObj->type != cObjectAlignment) {
+          ExecutiveDelete(G, oname);
+        } else {
+          obj = dynamic_cast<ObjectAlignment*>(execObj);
+        }
+      }
+      obj = ObjectAlignmentDefine(
+          G, obj, align_vla, align_state, true, trg_obj, mob_obj);
+      obj->Color = ColorGetIndex(G, "yellow");
+      ObjectSetName(obj, oname);
+      ExecutiveManageObject(G, obj, 0, quiet);
+      SceneInvalidate(G);
+    }
+  }
+
+  return result;
+}

layer3/Executive.cpp

@@ -38,6 +38,7 @@ Z* -------------------------------------------------------------------
 #include "SpecRecSpecial.h"
 #include "Tracker.h"
 #include "TrackerList.h"
+#include "USalign.h"
 #include "Word.h"
 #include "vla.h"
 
@@ -255,6 +256,11 @@ int ExecutiveAlign(PyMOLGlobals* G, const char* s1, const char* s2,
     float seq_wt, float radius, float scale, float base, float coord_wt,
     float expect, int window, float ante);
 
+pymol::Result<pymol::usalign::TMAlignResult> ExecutiveUSalign(
+    PyMOLGlobals* G, const char* mobile_sele, const char* target_sele,
+    int mobile_state, int target_state, int quiet, int transform,
+    const char* oname, int fast);
+
 void ExecutiveUpdateColorDepends(PyMOLGlobals* G, ObjectMolecule* mol);
 void ExecutiveUpdateCoordDepends(PyMOLGlobals* G, ObjectMolecule* mol);
 pymol::Result<float> ExecutiveDistance(PyMOLGlobals* G, const char* nam,

layer3/Executive.h

@@ -1951,6 +1951,61 @@ static PyObject *CmdAlign(PyObject * self, PyObject * args)
   }
 }
 
+static PyObject *CmdUSalign(PyObject * self, PyObject * args)
+{
+  PyMOLGlobals *G = nullptr;
+  const char *mobile, *target, *oname;
+  int mobile_state, target_state, quiet, transform, fast;
+  API_SETUP_ARGS(G, self, args, "Ossiiiisi", &self,
+      &mobile, &target, &mobile_state, &target_state,
+      &quiet, &transform, &oname, &fast);
+  API_ASSERT(APIEnterNotModal(G));
+
+  OrthoLineType s1, s2;
+  int ok = (SelectorGetTmp(G, mobile, s1) >= 0) &&
+           (SelectorGetTmp(G, target, s2) >= 0);
+
+  // Store results in locals — can't call Py_BuildValue until after APIExit
+  double tm_target = 0, tm_mobile = 0, rmsd = 0, seq_id = 0;
+  int ali_len = 0;
+  bool have_result = false;
+  std::string err_msg;
+
+  if (ok) {
+    auto res = ExecutiveUSalign(G, s1, s2,
+        mobile_state, target_state, quiet, transform, oname, fast);
+    if (res) {
+      auto& r = res.result();
+      tm_target = r.tm_score_target;
+      tm_mobile = r.tm_score_mobile;
+      rmsd = r.rmsd;
+      ali_len = r.aligned_length;
+      seq_id = r.seq_identity;
+      have_result = true;
+    } else {
+      err_msg = res.error().what();
+    }
+  }
+
+  SelectorFreeTmp(G, s1);
+  SelectorFreeTmp(G, s2);
+  APIExit(G);
+
+  if (have_result) {
+    return Py_BuildValue("{s:d,s:d,s:d,s:i,s:d}",
+        "tm_score_target", tm_target,
+        "tm_score_mobile", tm_mobile,
+        "RMSD", rmsd,
+        "alignment_length", ali_len,
+        "seq_identity", seq_id);
+  }
+  if (!err_msg.empty()) {
+    PyErr_SetString(PyExc_RuntimeError, err_msg.c_str());
+    return nullptr;
+  }
+  return APIFailure();
+}
+
 static PyObject *CmdGetCoordsAsNumPy(PyObject * self, PyObject * args)
 {
   PyMOLGlobals *G = nullptr;
@@ -6613,6 +6668,7 @@ static PyMethodDef Cmd_methods[] = {
   {"unset", CmdUnset, METH_VARARGS},
   {"unset_bond", CmdUnsetBond, METH_VARARGS},
   {"update", CmdUpdate, METH_VARARGS},
+  {"usalign", CmdUSalign, METH_VARARGS},
   {"window", CmdWindow, METH_VARARGS},
   {"zoom", CmdZoom, METH_VARARGS},
   {NULL, nullptr}                  /* sentinel */

layer4/Cmd.cpp

@@ -302,7 +302,8 @@
       intra_rms,         \
       intra_rms_cur,     \
       cealign,          \
-      pair_fit
+      pair_fit,          \
+      usalign
 
 #--------------------------------------------------------------------
 # ARE ALL OF THESE UNUSED AND/OR DEPRECATED (?)