- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: Thomas E. Cheatham, III <cheatham.chpc.utah.edu>

Date: Mon, 5 Apr 2004 15:21:36 -0600 (Mountain Daylight Time)

*> "(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

Date: Mon, 5 Apr 2004 15:21:36 -0600 (Mountain Daylight Time)

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"...

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

Custom Search