> "(Complete) Periodic Imaging of Truncated Octahedra for Dummies" :)
I adapted what is in sander() for use in ptraj(). sander() [PME] is set
up for generic triclinic imaging as opposed to the algorithmic version of
the truncated octahedron in Allen & Tildesley. Both algorithms will
produce a truncated octehedron however the "distance" along a box edge
will be different (I think) as well as the orientation. To image like in
AMBER, you want the triclinic imaging. This has the benefit of working
with other non-orthorhombic shapes (like rhombic dodecahedral, etc).
In ptraj() two versions of imaging are implemented, i.e. the standard
rectangular (since it is faster to imaging by truncating box length) and
the general triclinic.
In ptraj(), imaging is handled by transformImage() and in any routine that
measures a simple distance (transformDistance()) in actions.c
(1) Create reciprocal vectors
box[0->2] is X, Y, Z
box[3->5] is alpha, beta, gamma
ucell[0] = box[0];
ucell[1] = 0.0;
ucell[2] = 0.0;
ucell[3] = box[1]*cos(DEGRAD*box[5]);
ucell[4] = box[1]*sin(DEGRAD*box[5]);
ucell[5] = 0.0;
ucell[6] = box[2]*cos(DEGRAD*box[4]);
ucell[7] = (box[1]*box[2]*cos(DEGRAD*box[3]) - ucell[6]*ucell[3]) / ucell[4];
ucell[8] = sqrt(box[2]*box[2] - ucell[6]*ucell[6] - ucell[7]*ucell[7]);
u23x = ucell[4]*ucell[8] - ucell[5]*ucell[7];
u23y = ucell[5]*ucell[6] - ucell[3]*ucell[8];
u23z = ucell[3]*ucell[7] - ucell[4]*ucell[6];
u31x = ucell[7]*ucell[2] - ucell[8]*ucell[1];
u31y = ucell[8]*ucell[0] - ucell[6]*ucell[2];
u31z = ucell[6]*ucell[1] - ucell[7]*ucell[0];
u12x = ucell[1]*ucell[5] - ucell[2]*ucell[4];
u12y = ucell[2]*ucell[3] - ucell[0]*ucell[5];
u12z = ucell[0]*ucell[4] - ucell[1]*ucell[3];
volume=ucell[0]*u23x + ucell[1]*u23y + ucell[2]*u23z;
recip[0] = u23x/volume;
recip[1] = u23y/volume;
recip[2] = u23z/volume;
recip[3] = u31x/volume;
recip[4] = u31y/volume;
recip[5] = u31z/volume;
recip[6] = u12x/volume;
recip[7] = u12y/volume;
recip[8] = u12z/volume;
(2) Convert to reciprocal coordinates
fcx=cx*recip[0]+cy*recip[1]+cz*recip[2];
fcy=cx*recip[3]+cy*recip[4]+cz*recip[5];
fcz=cx*recip[6]+cy*recip[7]+cz*recip[8];
(3) Offset to origin rather than box center if desired
if( origin ) {
fcx += .5;
fcy += .5;
fcz += .5;
}
(4) Truncate the reciproal coordinates to [0->1] to bring it back to the
primary unit cell
ffcx = floor(fcx);
ffcy = floor(fcy);
ffcz = floor(fcz);
(4) Calculate the offset to add to the coordinates in each dimension to
image...
boxXtrans -= (ffcx*ucell[0] + ffcy*ucell[3] + ffcz*ucell[6]);
boxYtrans -= (ffcx*ucell[1] + ffcy*ucell[4] + ffcz*ucell[7]);
boxZtrans -= (ffcx*ucell[2] + ffcy*ucell[5] + ffcz*ucell[8]);
(5) Do the imaging,
for ( j = firstAtom; j < lastAtom; j++ ) {
if (mask[j]) {
x[j] += boxXtrans;
y[j] += boxYtrans;
z[j] += boxZtrans;
}
}
Note that this will create the skewed rectangle look rather than the more
familiar (and equivalent) "spherical" rendition of the truncated
octahedron...
To get the more familiar version shape, in ptraj(), I do a distance check
to find out what is the closest cell to the central cell to a given
imaging and adjust accordingly, i.e. if you have
_________________
/ x/
/ /
/________________/
imaging point x to:
_________________
/ /
/ /
/________________/
x
may bring it closer to the center. To do this, I go through each of the
adjacent periodic images and calculate the minimum distance (this gets
expensive)... ie. for the calculated boxXtrans, etc. values +/- one imaged
distance, which is the closest and use this. In transformImage() search
for the code above the text "IMAGING, FAMILIAR OFFSETS ARE"...
> though all atoms are certainly in the truncated octahedron. I've tried to
> sift through the general imaging code in ptraj.. but it is not entirely clear
> to me how that is working.
Let me know if you are still not completely understanding, but I think the
Allen&Tildesley is simply in a different orientation and may refer to a
different box size (i.e. original box length before the algorithm rather
than the side of the triclinic unit cell).
Good luck! I and likely the larger community would love to see MDDisplay
handle generic imaging.
--tom
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Mon Apr 05 2004 - 22:53:01 PDT
