User:Karsten Theis/Rigid body interpolation

From Proteopedia

Revision as of 18:38, 7 March 2021

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
  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
 * @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 in first set 
 * @param {atom set} anch2 Anchor to determine translation in 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
}

sel = [{912-984},{1125-1162},{not (1125-1162 or 912-984)}]
anch = [{984},{1125},{not (1125-1162 or 912-984)}]

original = {all}.xyz.all
for (var i FROM [1,20]) {
  {all}.xyz = @original
  progress = 0.05 * i
  for (var j FROM [1,sel.length]) {
    s = sel[j]
    a = anch[j]
    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}
      rigid(progress, s1, s2, a1, a2)
      linear(progress, s1, s2)
    }
  }
  delay 1.0
}
{all}.xyz = @original
  
 
  # fname = "morph" + i + ".pdb"
  # select 1.1
  # write @fname
 

Example