Skip to content

GitLab

Commit f5ec3142 authored by ChayaSt's avatar ChayaSt
Browse files

Merge branch 'master' of https://github.com/pandegroup/openmm into nbfix

parents 6b7ef87f 9b5b52f2
Hide whitespace changes
Inline Side-by-side
Showing with 86 additions and 45 deletions
platforms/cpu/src/CpuNeighborList.cpp
View file @ f5ec3142
......@@ -225,38 +225,78 @@ public:
if (usePeriodic)
voxelIndex.y = (y < 0 ? y+ny : (y >= ny ? y-ny : y));
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
// conditions, there may be two separate ranges.
float minx = centerPos[0];
float maxx = centerPos[0];
fvec4 offset(-xoffset, -yoffset+voxelSizeY*y+(usePeriodic ? 0.0f : miny), voxelSizeZ*z+(usePeriodic ? 0.0f : minz), 0);
for (int k = 0; k < (int) blockAtoms.size(); k++) {
const float* atomPos = &sortedPositions[4*(blockSize*blockIndex+k)];
fvec4 posVec(atomPos);
fvec4 delta1 = offset-posVec;
fvec4 delta2 = delta1+fvec4(0, voxelSizeY, voxelSizeZ, 0);
if (usePeriodic) {
delta1 -= round(delta1*invBoxSize)*boxSize;
delta2 -= round(delta2*invBoxSize)*boxSize;
if (usePeriodic && triclinic) {
for (int k = 0; k < (int) blockAtoms.size(); k++) {
const float* atomPos = &sortedPositions[4*(blockSize*blockIndex+k)];
fvec4 delta1(0, voxelSizeY*voxelIndex.y-atomPos[1], voxelSizeZ*voxelIndex.z-atomPos[2], 0);
fvec4 delta2 = delta1+fvec4(0, voxelSizeY, 0, 0);
fvec4 delta3 = delta1+fvec4(0, 0, voxelSizeZ, 0);
fvec4 delta4 = delta1+fvec4(0, voxelSizeY, voxelSizeZ, 0);
delta1 -= periodicBoxVec4[2]*floorf(delta1[2]*recipBoxSize[2]+0.5f);
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]);
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);
}
else {
float xoffset = (float) (usePeriodic ? boxy*periodicBoxVectors[1][0]+boxz*periodicBoxVectors[2][0] : 0);
fvec4 offset(-xoffset, -yoffset+voxelSizeY*y+(usePeriodic ? 0.0f : miny), voxelSizeZ*z+(usePeriodic ? 0.0f : minz), 0);
for (int k = 0; k < (int) blockAtoms.size(); k++) {
const float* atomPos = &sortedPositions[4*(blockSize*blockIndex+k)];
fvec4 posVec(atomPos);
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)
continue;
bool needPeriodic = (centerPos[1]-blockWidth[1] < maxDistance || centerPos[1]+blockWidth[1] > periodicBoxSize[1]-maxDistance ||
centerPos[2]-blockWidth[2] < maxDistance || centerPos[2]+blockWidth[2] > periodicBoxSize[2]-maxDistance ||
minx < 0.0f || maxx > periodicBoxVectors[0][0]);
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 ||
minx < 0.0f || maxx > periodicBoxVectors[0][0]);
int numRanges;
int rangeStart[2];
int rangeEnd[2];
......@@ -294,7 +334,7 @@ public:
continue;
fvec4 atomPos(&sortedPositions[4*sortedIndex]);
fvec4 delta = atomPos-centerPos;
fvec4 delta = atomPos-blockCenter;
if (periodicRectangular)
delta -= round(delta*invBoxSize)*boxSize;
else if (needPeriodic) {
......
platforms/cpu/tests/TestCpuNeighborList.cpp
View file @ f5ec3142
......@@ -53,14 +53,14 @@ void testNeighborList(bool periodic, bool triclinic) {
const float cutoff = 2.0f;
RealVec boxVectors[3];
if (triclinic) {
boxVectors[0] = RealVec(20, 0, 0);
boxVectors[1] = RealVec(5, 15, 0);
boxVectors[2] = RealVec(-3, -7, 22);
boxVectors[0] = RealVec(10, 0, 0);
boxVectors[1] = RealVec(4, 9, 0);
boxVectors[2] = RealVec(-3, -3.5, 11);
}
else {
boxVectors[0] = RealVec(20, 0, 0);
boxVectors[1] = RealVec(0, 15, 0);
boxVectors[2] = RealVec(0, 0, 22);
boxVectors[0] = RealVec(10, 0, 0);
boxVectors[1] = RealVec(0, 9, 0);
boxVectors[2] = RealVec(0, 0, 11);
}
const float boxSize[3] = {(float) boxVectors[0][0], (float) boxVectors[1][1], (float) boxVectors[2][2]};
const int blockSize = 8;
......
platforms/cuda/src/CudaKernels.cpp
View file @ f5ec3142
......@@ -1495,12 +1495,13 @@ public:
addEnergyKernel(addEnergyKernel), pmeEnergyBuffer(pmeEnergyBuffer), forceGroup(forceGroup) {
}
double computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) {
if ((groups&(1<<forceGroup)) != 0)
if ((groups&(1<<forceGroup)) != 0) {
cuStreamWaitEvent(cu.getCurrentStream(), event, 0);
if (includeEnergy) {
int bufferSize = pmeEnergyBuffer.getSize();
void* args[] = {&pmeEnergyBuffer.getDevicePointer(), &cu.getEnergyBuffer().getDevicePointer(), &bufferSize};
cu.executeKernel(addEnergyKernel, args, bufferSize);
if (includeEnergy) {
int bufferSize = pmeEnergyBuffer.getSize();
void* args[] = {&pmeEnergyBuffer.getDevicePointer(), &cu.getEnergyBuffer().getDevicePointer(), &bufferSize};
cu.executeKernel(addEnergyKernel, args, bufferSize);
}
}
return 0.0;
}
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ f5ec3142
......@@ -1502,9 +1502,9 @@ public:
events[0] = event;
event = cl::Event();
cl.getQueue().enqueueWaitForEvents(events);
if (includeEnergy)
cl.executeKernel(addEnergyKernel, pmeEnergyBuffer.getSize());
}
if (includeEnergy)
cl.executeKernel(addEnergyKernel, pmeEnergyBuffer.getSize());
return 0.0;
}
private:
......
plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h
View file @ f5ec3142
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Mark Friedrichs, Peter Eastman *
* Contributors: *
* *
......@@ -87,7 +87,7 @@ public:
* Extrapolated perturbation theory approximation. The dipoles are iterated a few times, and then an analytic
* approximation is used to extrapolate to the fully converged values. Call setExtrapolationCoefficients()
* to set the coefficients used for the extrapolation. The default coefficients used in this release are
* [0, -0.3, 0, 1.3], but be aware that those may change in a future release.
* [-0.154, 0.017, 0.658, 0.474], but be aware that those may change in a future release.
*/
Extrapolated = 2
......@@ -322,7 +322,7 @@ public:
/**
* Get the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the extrapolation
* algorithm for induced dipoles. In this release, the default values for the coefficients are
* [0, -0.3, 0, 1.3], but be aware that those may change in a future release.
* [-0.154, 0.017, 0.658, 0.474], but be aware that those may change in a future release.
*/
const std::vector<double>& getExtrapolationCoefficients() const;
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ f5ec3142
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: *
* Contributors: *
* *
......@@ -43,10 +43,10 @@ AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), polari
mutualInducedTargetEpsilon(1.0e-02), scalingDistanceCutoff(100.0), electricConstant(138.9354558456), aewald(0.0) {
pmeGridDimension.resize(3);
pmeGridDimension[0] = pmeGridDimension[1] = pmeGridDimension[2];
extrapolationCoefficients.push_back(0.0);
extrapolationCoefficients.push_back(-0.3);
extrapolationCoefficients.push_back(0.0);
extrapolationCoefficients.push_back(1.3);
extrapolationCoefficients.push_back(-0.154);
extrapolationCoefficients.push_back(0.017);
extrapolationCoefficients.push_back(0.658);
extrapolationCoefficients.push_back(0.474);
}
AmoebaMultipoleForce::NonbondedMethod AmoebaMultipoleForce::getNonbondedMethod() const {
......
wrappers/python/simtk/openmm/app/forcefield.py
View file @ f5ec3142
......@@ -856,13 +856,13 @@ class ForceField(object):
uniquePropers = set()
for angle in data.angles:
for atom in bondedToAtom[angle[0]]:
if atom != angle[1]:
if atom not in angle:
if atom < angle[2]:
uniquePropers.add((atom, angle[0], angle[1], angle[2]))
else:
uniquePropers.add((angle[2], angle[1], angle[0], atom))
for atom in bondedToAtom[angle[2]]:
if atom != angle[1]:
if atom not in angle:
if atom > angle[0]:
uniquePropers.add((angle[0], angle[1], angle[2], atom))
else:
......
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