VMD-L Mailing List

From: Andrew Dalke (dalke_at_mag.com)
Date: Thu Dec 04 1997 - 18:41:03 CST

[Think of this as a mini-VMDTECH :)]

Nithin Rai <nithin_at_mblab.gla.ac.uk> asked:

> has anyone written a tcl script for writing the xyz co-ordinates of
> a molecule after it has been modified ??

  The reason what you did doesn't work is because rotations with the
mouse do not affect the underlying coordinates. The mouse only
changes the matrix transformation used in the display routines. As
with much of VMD, that transformation information is available to the
scripting language, so you can build what you want.

  First, a bit about graphics transformations. VMD, as with most
graphics programs, uses a 4x4 graphics matrix to convert the 3D
coordinates from "molecule space" (the coordinates in Angstroms) to
scene space (which, according to VMD, is a cube from -1 to 1 on each
side). There a four matricies involved in the process.
  1) perform a translation to center the molecule (bring the center of
the molecule to (0,0,0))
  2) scale the system
  3) rotate the molecule around the origin to the desired orientation
  4) perform the final (global) translation transformation

The process is often written as:

   x' = (global * rotate * scale * center) * x
         \ /
          \_______ The view matrix ____/

For more details about this, see one of the many books on 3D graphics.
The most popular is Foley and van Dam, but there are many others.

To get the information in VMD, use the 'molinfo' command. For
example:

vmd > molinfo top get center_matrix

Info) {{1.000000 0.000000 0.000000 -28.416485} {0.000000 1.000000
0.000000 -9.767003} {0.000000 0.000000 1.000000 -0.333005} {0.000000
0.000000 0.000000 1.000000}}

This says there is a translation by (-28.416485, -9.767003, -0.333005)
to get the center of the molecule to the origin.

The keywords for the different matricies are center_matrix,
rotate_matrix, scale_matrix and global_matrix. There is another
variable, view_matrix, which is the product of all four of these.

The question is, what transformation should be applied to the molecule
to get it to the right coordinates? It isn't simple. You don't need
the scaling, since that doesn't make sense (or do you want a C-C bond
that's 0.03A across?).

  The final transformation, a global translation, looks useful, but
the key thing to remember is that that information is in scene space,
and not coordinate space. Thus some additional transformations would
have to occur to use it. It turns out you would only worry about this
one if you want to compare two different molecules (eg, if you are
doing docking). In that case you would have to compute the inverse
view matrix of one system to get you to the other one.

  Instead, you probably want just the molecule coordinates centered
around the same center, but with the new rotation applied to it. If
that is so, the following will work:

proc save_rotated_coordinates {filename {mol top}} {
        # get the selection and save the current coordinates
        set sel [atomselect $mol all]
        set xyz [$sel get {x y z}]

        # get the center and rotation transformation
        lassign [molinfo $mol get {center_matrix rotate_matrix}] \
              center rotate

        # compute the matrix inverse of the center matrix
        set center_inv [measure inverse $center]

        # the transformation matrix for the rotated coordinates is
        set trans [transmult $center_inv $rotate $center]

        # apply it to the coordinates
        $sel move $trans

        # save the coordinates
        $sel writepdb $filename

        # revert to the original coordinates
        $sel set {x y z} $xyz
}

Then use it like:

  save_rotated_coordinates "rotated.pdb"

Here's what I did to test it:

  mol load pdb test.pdb
  rotate y by 30
  save_rotated_coordinates rotated.pdb
  mol load pdb rotated.pdb
  mol fix top
  rotate y by 30

At the end of this, the two molecules should be overlapping.

Let me point out that the trick I used for the procedure is not very
clean. I replaced the existing coordinates instead of doing
everything independently. This can be a problem when you have complex
structures, like surfaces. These are recomputed when any of the
coordinates are changed. On the other hand, I know the 'move' option
for selections is very fast, and I can use the 'writepdb' command
instead of creating a new PDB output routine.

                                                Andrew Dalke
                                                dalke_at_mag.com
P.S.
  If you wanted to move one molecule into the coordinate space of
another you would need to have two molecule ids ("from_mol" and
"to_mol"), and the transformations would be something like:

  lassign [molinfo $from_mol get view_matrix] view_matrix_from
  lassign [molinfo $to_mol get view_matrix] view_matrix_to

  set trans [transmult [measure inverse $view_matrix_to] $view_matrix_from] 

---
Not speaking for the Molecular Applications Group.