Skip to content

GitLab

Commit e1c49ebd authored by peastman's avatar peastman
Browse files

Merge pull request #1403 from peastman/neighbors 

Fixed a bug in CPU neighbor list with triclinic boxes
parents 7633a01c a4f7077b
Hide whitespace changes
Inline Side-by-side
Showing with 68 additions and 28 deletions
platforms/cpu/src/CpuNeighborList.cpp
View file @ e1c49ebd
...@@ -225,38 +225,78 @@ public: ...@@ -225,38 +225,78 @@ public:
if (usePeriodic) if (usePeriodic)
voxelIndex.y = (y < 0 ? y+ny : (y >= ny ? y-ny : y)); voxelIndex.y = (y < 0 ? y+ny : (y >= ny ? y-ny : y));
float boxy = floor((float) y/ny); float boxy = floor((float) y/ny);
float xoffset = (float) (usePeriodic ? boxy*periodicBoxVectors[1][0]+boxz*periodicBoxVectors[2][0] : 0);
// Identify the range of atoms within this bin we need to search. When using periodic boundary // Identify the range of atoms within this bin we need to search. When using periodic boundary
// conditions, there may be two separate ranges. // conditions, there may be two separate ranges.
float minx = centerPos[0]; float minx = centerPos[0];
float maxx = centerPos[0]; float maxx = centerPos[0];
fvec4 offset(-xoffset, -yoffset+voxelSizeY*y+(usePeriodic ? 0.0f : miny), voxelSizeZ*z+(usePeriodic ? 0.0f : minz), 0); if (usePeriodic && triclinic) {
for (int k = 0; k < (int) blockAtoms.size(); k++) { for (int k = 0; k < (int) blockAtoms.size(); k++) {
const float* atomPos = &sortedPositions[4*(blockSize*blockIndex+k)]; const float* atomPos = &sortedPositions[4*(blockSize*blockIndex+k)];
fvec4 posVec(atomPos); fvec4 delta1(0, voxelSizeY*voxelIndex.y-atomPos[1], voxelSizeZ*voxelIndex.z-atomPos[2], 0);
fvec4 delta1 = offset-posVec; fvec4 delta2 = delta1+fvec4(0, voxelSizeY, 0, 0);
fvec4 delta2 = delta1+fvec4(0, voxelSizeY, voxelSizeZ, 0); fvec4 delta3 = delta1+fvec4(0, 0, voxelSizeZ, 0);
if (usePeriodic) { fvec4 delta4 = delta1+fvec4(0, voxelSizeY, voxelSizeZ, 0);
delta1 -= round(delta1*invBoxSize)*boxSize; delta1 -= periodicBoxVec4[2]*floorf(delta1[2]*recipBoxSize[2]+0.5f);
delta2 -= round(delta2*invBoxSize)*boxSize; delta1 -= periodicBoxVec4[1]*floorf(delta1[1]*recipBoxSize[1]+0.5f);
delta1 -= periodicBoxVec4[0]*floorf(delta1[0]*recipBoxSize[0]+0.5f);
delta2 -= periodicBoxVec4[2]*floorf(delta2[2]*recipBoxSize[2]+0.5f);
delta2 -= periodicBoxVec4[1]*floorf(delta2[1]*recipBoxSize[1]+0.5f);
delta2 -= periodicBoxVec4[0]*floorf(delta2[0]*recipBoxSize[0]+0.5f);
delta3 -= periodicBoxVec4[2]*floorf(delta3[2]*recipBoxSize[2]+0.5f);
delta3 -= periodicBoxVec4[1]*floorf(delta3[1]*recipBoxSize[1]+0.5f);
delta3 -= periodicBoxVec4[0]*floorf(delta3[0]*recipBoxSize[0]+0.5f);
delta4 -= periodicBoxVec4[2]*floorf(delta4[2]*recipBoxSize[2]+0.5f);
delta4 -= periodicBoxVec4[1]*floorf(delta4[1]*recipBoxSize[1]+0.5f);
delta4 -= periodicBoxVec4[0]*floorf(delta4[0]*recipBoxSize[0]+0.5f);
if (delta1[1] < 0 && delta1[1]+voxelSizeY > 0)
delta1 = fvec4(delta1[0], 0, delta1[2], 0);
if (delta1[2] < 0 && delta1[2]+voxelSizeZ > 0)
delta1 = fvec4(delta1[0], delta1[1], 0, 0);
if (delta3[1] < 0 && delta3[1]+voxelSizeY > 0)
delta3 = fvec4(delta3[0], 0, delta3[2], 0);
if (delta2[2] < 0 && delta2[2]+voxelSizeZ > 0)
delta2 = fvec4(delta2[0], delta2[1], 0, 0);
fvec4 delta = min(min(min(abs(delta1), abs(delta2)), abs(delta3)), abs(delta4));
float dy = (voxelIndex.y == atomVoxelIndex[k].y ? 0.0f : delta[1]);
float dz = (voxelIndex.z == atomVoxelIndex[k].z ? 0.0f : delta[2]);
float dist2 = maxDistanceSquared-dy*dy-dz*dz;
if (dist2 > 0) {
float dist = sqrtf(dist2);
minx = min(minx, atomPos[0]-dist-max(max(max(delta1[0], delta2[0]), delta3[0]), delta4[0]));
maxx = max(maxx, atomPos[0]+dist-min(min(min(delta1[0], delta2[0]), delta3[0]), delta4[0]));
}
} }
fvec4 delta = min(abs(delta1), abs(delta2)); }
float dy = (y == atomVoxelIndex[k].y ? 0.0f : delta[1]); else {
float dz = (z == atomVoxelIndex[k].z ? 0.0f : delta[2]); float xoffset = (float) (usePeriodic ? boxy*periodicBoxVectors[1][0]+boxz*periodicBoxVectors[2][0] : 0);
float dist2 = maxDistanceSquared-dy*dy-dz*dz; fvec4 offset(-xoffset, -yoffset+voxelSizeY*y+(usePeriodic ? 0.0f : miny), voxelSizeZ*z+(usePeriodic ? 0.0f : minz), 0);
if (dist2 > 0) { for (int k = 0; k < (int) blockAtoms.size(); k++) {
float dist = sqrtf(dist2); const float* atomPos = &sortedPositions[4*(blockSize*blockIndex+k)];
minx = min(minx, atomPos[0]-dist-xoffset); fvec4 posVec(atomPos);
maxx = max(maxx, atomPos[0]+dist-xoffset); fvec4 delta1 = offset-posVec;
fvec4 delta2 = delta1+fvec4(0, voxelSizeY, voxelSizeZ, 0);
if (usePeriodic) {
delta1 -= round(delta1*invBoxSize)*boxSize;
delta2 -= round(delta2*invBoxSize)*boxSize;
}
fvec4 delta = min(abs(delta1), abs(delta2));
float dy = (y == atomVoxelIndex[k].y ? 0.0f : delta[1]);
float dz = (z == atomVoxelIndex[k].z ? 0.0f : delta[2]);
float dist2 = maxDistanceSquared-dy*dy-dz*dz;
if (dist2 > 0) {
float dist = sqrtf(dist2);
minx = min(minx, atomPos[0]-dist-xoffset);
maxx = max(maxx, atomPos[0]+dist-xoffset);
}
} }
} }
if (minx == maxx) if (minx == maxx)
continue; continue;
bool needPeriodic = (centerPos[1]-blockWidth[1] < maxDistance || centerPos[1]+blockWidth[1] > periodicBoxSize[1]-maxDistance || bool needPeriodic = usePeriodic && (centerPos[1]-blockWidth[1] < maxDistance || centerPos[1]+blockWidth[1] > periodicBoxSize[1]-maxDistance ||
centerPos[2]-blockWidth[2] < maxDistance || centerPos[2]+blockWidth[2] > periodicBoxSize[2]-maxDistance || centerPos[2]-blockWidth[2] < maxDistance || centerPos[2]+blockWidth[2] > periodicBoxSize[2]-maxDistance ||
minx < 0.0f || maxx > periodicBoxVectors[0][0]); minx < 0.0f || maxx > periodicBoxVectors[0][0]);
int numRanges; int numRanges;
int rangeStart[2]; int rangeStart[2];
int rangeEnd[2]; int rangeEnd[2];
...@@ -294,7 +334,7 @@ public: ...@@ -294,7 +334,7 @@ public:
continue; continue;
fvec4 atomPos(&sortedPositions[4*sortedIndex]); fvec4 atomPos(&sortedPositions[4*sortedIndex]);
fvec4 delta = atomPos-centerPos; fvec4 delta = atomPos-blockCenter;
if (periodicRectangular) if (periodicRectangular)
delta -= round(delta*invBoxSize)*boxSize; delta -= round(delta*invBoxSize)*boxSize;
else if (needPeriodic) { else if (needPeriodic) {
......
platforms/cpu/tests/TestCpuNeighborList.cpp
View file @ e1c49ebd
...@@ -53,14 +53,14 @@ void testNeighborList(bool periodic, bool triclinic) { ...@@ -53,14 +53,14 @@ void testNeighborList(bool periodic, bool triclinic) {
const float cutoff = 2.0f; const float cutoff = 2.0f;
RealVec boxVectors[3]; RealVec boxVectors[3];
if (triclinic) { if (triclinic) {
boxVectors[0] = RealVec(20, 0, 0); boxVectors[0] = RealVec(10, 0, 0);
boxVectors[1] = RealVec(5, 15, 0); boxVectors[1] = RealVec(4, 9, 0);
boxVectors[2] = RealVec(-3, -7, 22); boxVectors[2] = RealVec(-3, -3.5, 11);
} }
else { else {
boxVectors[0] = RealVec(20, 0, 0); boxVectors[0] = RealVec(10, 0, 0);
boxVectors[1] = RealVec(0, 15, 0); boxVectors[1] = RealVec(0, 9, 0);
boxVectors[2] = RealVec(0, 0, 22); boxVectors[2] = RealVec(0, 0, 11);
} }
const float boxSize[3] = {(float) boxVectors[0][0], (float) boxVectors[1][1], (float) boxVectors[2][2]}; const float boxSize[3] = {(float) boxVectors[0][0], (float) boxVectors[1][1], (float) boxVectors[2][2]};
const int blockSize = 8; const int blockSize = 8;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment