User:Karsten Theis/Rigid body interpolation
From Proteopedia
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
| |||||||||||
