VMD-L Mailing List

From: Justin Gullingsrud (justin_at_ks.uiuc.edu)
Date: Thu Jul 10 2003 - 09:20:41 CDT

Hi,

On Thu, Jul 10, 2003 at 02:13:07PM +0200, Vlad Cojocaru wrote:
> Dear VMD users,
> I am trying to do the following:
> I have an RNA molecule (from a crystal structure) which containd 2
> helices and 1 loop. I am trying to extend the helices with 7-10 base
> pairs. So I created the A-RNA helices by including a GC base pair common
> to the crystal structure..then I was trying to fit the common base pair
> (by fitting the base atoms) but the problem is that the fit gave a
> strange matrix one of the molecule moves in a very strange position.
> Doing it by hand I can see that the common base pair (in tha crystal and
> in the extended helix) can be overlayed without any problem. Now I have
> a question....Is there another way to overlay common parts of two
> molecules in VMD? When selecting the atoms for fit is VMD respecting the
> order that I gave to atomselect comand...for instance if I say:
> atomselect 0 [{(index 20 14 15}]
> atomselect 1 [{index 16 13 19}]
> is VMD trying to fit atom 20 (mol 0) to atom 16 (mol 1); atom 14 (mol 0)
> to atom 13 (mol 1)...........or does it do it in another way??

You're right - VMD does the matching by increasing atom index.

Here's a way you could cheat: create a copy of the coordinates using
the command "animate dup top", then set the coordinates of the atoms
in the new molecule to have the mixed up order you want. Compute the
fitting matrix from the mixed-up coordinates, then use that matrix to move
the original set of coordinates.

You could also just compute the fitting matrix yourself. Here's
the source code for how VMD handles three points:

static Matrix4 myfit3(const float *x1, const float *x2,
                      const float *y1, const float *y2,
                      const float *comx, const float *comy) {

  Matrix4 mx, my, rx1, ry1;
  double dx1[3], dy1[3], angle;
  float dx2[3], dy2[3], x2t[3], y2t[3];

  for (int i=0; i<3; i++) {
    dx1[i] = x1[i] - comx[i];
    dx2[i] = x2[i] - comx[i];
    dy1[i] = y1[i] - comy[i];
    dy2[i] = y2[i] - comy[i];
  }

  transvecinv(dx1, rx1);
  rx1.multpoint3d(dx2, x2t);
  angle = atan2(x2t[2], x2t[1]);
  mx.translate(-comx[0], -comx[1], -comy[2]);
  mx.multmatrix(rx1);
  mx.rot((float) RADTODEG(angle), 'x');

  transvecinv(dy1, ry1);
  ry1.multpoint3d(dy2, y2t);
  angle = atan2(y2t[2], y2t[1]);
  my.translate(-comy[0], -comy[1], -comy[2]);
  my.multmatrix(ry1);
  my.rot((float) RADTODEG(angle), 'x');

  my.inverse();
  mx.multmatrix(my);
  return mx;
}

It should be pretty easy to convert this into a Tcl script that does the
same thing using the built-in VMD vector/matrix commands (probably because
it used to be a Tcl script ;).

Cheers,
Justin

> Thanks a lot for any advices,
> Best regards,
> Vlad
>
> --
> Vlad Cojocaru
> Max Planck Institute for Biophysical Chemistry
> Department: 060
> Am Fassberg 11, 37077 Goettingen, Germany
> tel: ++49-551-201.1327
> e-mail: Vlad.Cojocaru_at_mpi-bpc.mpg.de
> home tel: ++49-551-9963204
>
> 

-- 
  Justin Gullingsrud        3111 Beckman Institute        217-244-8946
  I been dropping the new science, and I be kicking the new knowledge,
  and I'm seeing to a degree that you can't get in college.  -- b.boys