Possible Bug in printing Triclinic Box Coordinates for PDB?
DEC014
dcoss... at gmail.com
Thu Apr 11 15:07:14 UTC 2013
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>
More information about the CP2K-user
mailing list