User:Karsten Theis/Rigid body interpolation

From Proteopedia

Revision as of 18:15, 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.

function linear(progress, sel1, sel2) {
/* linear interpolation of coordinates in the atom sets sel1 and sel2 */
/* progress is a number between 0 and 1, and indicates the weight for atoms in selection sel2 */
  coord1 = {@sel1}.xyz.all
  coord2 = {@sel2}.xyz.all
  len = coord1.length
  ssergorp = 1 - progress
  for (var i FROM [1,len]) {coord1[i] = (coord2[i] * progress + coord1[i] * ssergorp)}
  {@sel1}.xyz = @coord1
/* interpolated coordinates are stored in selection sel1 */
}

function rigid(progress, sel1, sel2, anch1, anch2) {
/* atoms selected by sel1 and sel2 are transformed onto each other */
/* atoms are rotated into an intermediate orientation determined by progress */
/* the atoms are translated such that a set of atoms (defined by anch1 and anch2) is positioned
on a common position anchorpos, the linear interpolation between the two anchors */
  fourbyfour = compare({@sel1}, {@sel2}) /* calculates lowest RMSD superposition */
  total_quat = quaternion(@fourbyfour%1) /* extract rotation sel1 -> sel2 and convert into quaternion */
  theta = total_quat %"theta"
  ax = total_quat %"vector"
  if (theta > 160 and theta < 170) total_quat = quaternion(ax, 360 - theta) /* special case, needs a parameter*/
  ssergorp = 1 - progress
  quat1 = total_quat * progress      /* partial rotation for atoms in selection sel1 */
  quat2 = total_quat * ssergorp * -1 /* partial rotation for atoms in selection sel2 in the other direction */
  anchorpos = {@anch2}.xyz * progress + {@anch1}.xyz * ssergorp /* this is where the anchors should go */
  select sel1
  rotateselected @quat1 molecular /* rotate the coordinates into common orientation */
  transl1 = anchorpos - {@anch1}.xyz /* figure out how much to translate to have anchor in desired position */
  translateselected @transl1 /* translate (shift) the coordinates */
  select sel2
  rotateselected @quat2 molecular /* same steps for second set of atomic coordinates, using the other rotation */
  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