Re: amber-developers: About calculating the minimum imaged distance

From: Thomas Cheatham III <tec3.utah.edu>
Date: Wed, 19 Sep 2007 14:34:17 -0600 (Mountain Daylight Time)

> The program ptraj may give out the incorrect minimum imaged distance
> between two atoms when the non-imaged distance between the atoms is
> greater than two times the box size. This is caused by the function
> calculateDistance2. I am not sure whether the programs sander and pmemd
> have the same bug. Can someone give me the answer of this question?

This bug is only relevant to ptraj and it is a bug (or misbalance or
misfeature). It is only relevant to non-orthorhombic periodicity. The
suggested fix for the "GROMACS" truncated octahedral is not general and
not appropriate (as the "AMBER" definition is a triclinic unit cell that
requires interconversion between reciprocal and regular coordinates.

The misfeature results from my mistaken assumption that all of the
coordinates would be imaged prior to using this command. Then, I search
each unit cell +/-1 away to see if this results in a smaller (imaged)
distance.

The command as implemented will work fine if you run

  image origin familar

prior to doing the distance or closestwater commands.

I have fixed my code and will soon incorporate this into the amber10
version...

--tec3

The code fix is to first image each set of coordinates (x[j] and x[i])
prior to doing the "search" over grids (the for loops over ix, iy, iz):


    fx = x[j]*recip[0] + y[j]*recip[1] + z[j]*recip[2];
    fy = x[j]*recip[3] + y[j]*recip[4] + z[j]*recip[5];
    fz = x[j]*recip[6] + y[j]*recip[7] + z[j]*recip[8];

    f2x = x[i]*recip[0] + y[i]*recip[1] + z[i]*recip[2];
    f2y = x[i]*recip[3] + y[i]*recip[4] + z[i]*recip[5];
    f2z = x[i]*recip[6] + y[i]*recip[7] + z[i]*recip[8];

    fx = fx - floor(fx);
    fy = fy - floor(fy);
    fz = fz - floor(fz);

    f2x = f2x - floor(f2x);
    f2y = f2y - floor(f2y);
    f2z = f2z - floor(f2z);

    X = (fx*ucell[0] + fy*ucell[3] + fz*ucell[6]) - (f2x*ucell[0] + f2y*ucell[3] + f2z*ucell[6]);
    Y = (fx*ucell[1] + fy*ucell[4] + fz*ucell[7]) - (f2x*ucell[1] + f2y*ucell[4] + f2z*ucell[7]);
    Z = (fx*ucell[2] + fy*ucell[5] + fz*ucell[8]) - (f2x*ucell[2] + f2y*ucell[5] + f2z*ucell[8]);
Received on Sun Sep 23 2007 - 06:07:17 PDT
Custom Search