Possible Bug in printing Triclinic Box Coordinates for PDB?

[DEC014]

Hey all,

So after many /head-desk and /head-wall moments these past two weeks I 
discovered something unusual.  In the PDB format, it ASSUMES (we all know 
what happens when you assume things) the A-axis is co-linear with the 
X-axis and the B-axis is in the XY-plane.  CP2K does not make this 
conversion before printing out Cell Vectors and Coordinates to PDB Format.  
Made this discovery trying to figure out why CUBE files visualize correctly 
but PDB files didn't in VMD.  VMD does not make that assumption when 
importing CUBE files and performs the Transformation while importing.  
Makes visualizing things in VMD interesting to say the least.  So is this a 
bug with the code or an oversight?  

Dug through the VMD CUBE Plugin source code and found the relevant code.  
The following is in C, but I'm sure someone can convert it to Fortran and 
implement it into CP2K:

Finding the Transformation Matrix:
// calculate and store origin and rotation matrix to realign everything 
later.
static void cube_buildrotmat(cube_t *cube, float *o, float *a, float *b)
{
  // we rotate first around y and z to align a along the x-axis...
  const double len   = sqrt(a[0]*a[0] + a[1]*a[1]);
  const double phi   = atan2((double) a[2], (double) len);
  const double theta = atan2((double) a[1], (double) a[0]);

  const double cph = cos(phi);
  const double cth = cos(theta);
  const double sph = sin(phi);
  const double sth = sin(theta);

  // ...then we rotate around x to put b into the xy-plane.
  const double psi = atan2(-sph*cth*b[0] - sph*sth*b[1] + 
cph*b[2],-sth*b[0] + cth*b[1]);
  const double cps = cos(psi);
  const double sps = sin(psi);

  const double r[3][3] = {
    {               cph*cth,                    cph*sth,      sph},
    {-sth*cps - sph*cth*sps,      cth*cps - sph*sth*sps,  cph*sps},
    { sth*sps - sph*cth*cps,     -cth*sps - sph*sth*cps,  cph*cps}
  };

  for (int i=0; i<3; ++i) {
    cube->origin[i] = o[i];
    for (int j=0; j<3; ++j) {
      cube->rotmat[i][j] = r[i][j];
    }
  }
}

Rotating Vectors:
 // store aligned axes.
  for (i=0; i<3; ++i) {
    voltmpl.xaxis[i] = cube->rotmat[i][0] * a[0]
      + cube->rotmat[i][1] * a[1] + cube->rotmat[i][2] * a[2];

    voltmpl.yaxis[i] = cube->rotmat[i][0] * b[0]
      + cube->rotmat[i][1] * b[1] + cube->rotmat[i][2] * b[2];

    voltmpl.zaxis[i] = cube->rotmat[i][0] * c[0]
      + cube->rotmat[i][1] * c[1] + cube->rotmat[i][2] * c[2];
  }

Rotating Coordinates is just [transformation matrix] * [coordinate]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20130411/05ce89da/attachment.htm>