User:Karsten Theis/Rigid body interpolation

From Proteopedia

The methods described here has been further developed into the Storymorph Jmol scrips to perform supoerpositions and to create morph, see Jmol/Storymorph.

This page shows a way to superpose two rigid bodies so that they both have the same rotation, and they are translated such that a common anchor (center of mass, given atom) is in the same location.

The entire coordinates are copied first so that you can play around with the parameters and then reset the coordinates. The starting point is a structure with two models containing the same atoms in the same order, but with different coordinates.

/**
 * Linear interpolation of coordinates in the atom sets sel1 and sel2
 * @param {atom set} sel1 First set of atoms, will be modified by the function
 * @param {atom set} sel2 Second related set of atoms, will not be touched
 * @param {float} progress Weight for atoms in selection sel2
**/
function linear(progress, sel1, sel2) {

  // Copy coordinates to list of points and determine number of points
  coord1 = {@sel1}.xyz.all
  coord2 = {@sel2}.xyz.all
  len = coord1.length

  ssergorp = 1 - progress
  for (var i FROM [1,len]) {

    // Weighted average of coordinates for equivalent points
    coord1[i] = (coord2[i] * progress + coord1[i] * ssergorp)

  }

  // Interpolated coordinates are stored in selection sel1
  {@sel1}.xyz = @coord1
}

/**
 * Rigid body interpolation of coordinates in the equivalent atom sets sel1 and sel2
 * @param {float} progress Relative weight of sel2 vs sel1 to choose rotation and translation
 * @param {atom set} sel1 First set of atoms, will be modified by the function
 * @param {atom set} sel2 Second related set of atoms, will be modified by the function
 * @param {atom set} anch1 Anchor to determine translation of first set 
 * @param {atom set} anch2 Anchor to determine translation of second set
**/
function rigid(progress, sel1, sel2, anch1, anch2) {

  // Calculate lowest RMSD superposition of sel1 and sel2 as 4x4 matrix
  fourbyfour = compare({@sel1}, {@sel2})

  // Extract rotation sel1 -> sel2 and convert into quaternion
  total_quat = quaternion(@fourbyfour%1)

  // For special case, rotate in the other direction (should have a parameter to make explicit 
  theta = total_quat %"theta"
  ax = total_quat %"vector"
  if (theta > 160 and theta < 170) total_quat = quaternion(ax, 360 - theta) 

  ssergorp = 1 - progress

  // Partial rotation for atoms in sel1
  quat1 = total_quat * progress

  // Partial rotation for atoms in sel2 in the other direction */
  quat2 = total_quat * ssergorp * -1

  // Desired postion of anchors store in anchorpos
  anchorpos = {@anch2}.xyz * progress + {@anch1}.xyz * ssergorp 

  // Work of first set
  select sel1

  // Rotate the coordinates into common orientation
  rotateselected @quat1 molecular

  // Calculate and apply translation to match anchor in rotated rigid body with anchorpos
  transl1 = anchorpos - {@anch1}.xyz 
  translateselected @transl1 /* translate (shift) the coordinates */

  // Same steps for second (equivalent) set of atoms 
  select sel2
  rotateselected @quat2 molecular
  transl2 = anchorpos - {@anch2}.xyz
  translateselected @transl2
}

// Define rigid bodies
sel = [{912-984},{1125-1162},{not (1125-1162 or 912-984)}]

//Define anchors for each rigid body
anch = [{984},{1125},{not (1125-1162 or 912-984)}]

// Save coordinates to use in every step of the interpolation
original = {all}.xyz.all

for (var i FROM [1,20]) {
  progress = 0.05 * i
  for (var j FROM [1,sel.length]) {
    s = sel[j]
    a = anch[j]

    // Spike protein is a trimer with chains A, B, and C. Each chain is morphed individually
    for (var c in [{chain='A'},{chain='B'},{chain='C'}]) {
      s1 = {1.1 and @s and @c}
      s2 = {1.2 and @s and @c}
      a1 = {1.1 and @a and @c}
      a2 = {1.2 and @a and @c}

      // Atom sets s1 and s2 will be moved into a common frame on the path from s1 to s2.
      // Position on the path (closer to s1, closer to s2) is determined by parameter progress
      rigid(progress, s1, s2, a1, a2)

      // Linear interpolation between the superimposed coordinates, weighted by parameter progress
      linear(progress, s1, s2)
    }
  }
  delay 1.0
  // To save this step in the morph, un-comment the following three lines while running in local Jmol app
  // fname = "morph" + i + ".pdb"
  // select 1.1
  // write @fname

  // Set coordinates back to original values in preparation of calculating next step in the morph
  {all}.xyz = @original
}
 

Example