Commit 18501459 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed various bugs that caused energy to come out wrong after per-particle parameters were changed

parent 233c6a70
...@@ -796,9 +796,11 @@ void OpenCLContext::validateMolecules() { ...@@ -796,9 +796,11 @@ void OpenCLContext::validateMolecules() {
velm->upload(); velm->upload();
atomIndex->upload(); atomIndex->upload();
findMoleculeGroups(); findMoleculeGroups();
for (int i = 0; i < (int) reorderListeners.size(); i++)
reorderListeners[i]->execute();
} }
void OpenCLContext::reorderAtoms() { void OpenCLContext::reorderAtoms(bool enforcePeriodic) {
if (numAtoms == 0 || nonbonded == NULL || !nonbonded->getUseCutoff()) if (numAtoms == 0 || nonbonded == NULL || !nonbonded->getUseCutoff())
return; return;
if (moleculesInvalid) if (moleculesInvalid)
...@@ -873,16 +875,18 @@ void OpenCLContext::reorderAtoms() { ...@@ -873,16 +875,18 @@ void OpenCLContext::reorderAtoms() {
molPos[i].x -= dx; molPos[i].x -= dx;
molPos[i].y -= dy; molPos[i].y -= dy;
molPos[i].z -= dz; molPos[i].z -= dz;
for (int j = 0; j < (int) atoms.size(); j++) { if (enforcePeriodic) {
int atom = atoms[j]+mol.offsets[i]; for (int j = 0; j < (int) atoms.size(); j++) {
mm_float4 p = posq->get(atom); int atom = atoms[j]+mol.offsets[i];
p.x -= dx; mm_float4 p = posq->get(atom);
p.y -= dy; p.x -= dx;
p.z -= dz; p.y -= dy;
posq->set(atom, p); p.z -= dz;
posCellOffsets[atom].x -= xcell; posq->set(atom, p);
posCellOffsets[atom].y -= ycell; posCellOffsets[atom].x -= xcell;
posCellOffsets[atom].z -= zcell; posCellOffsets[atom].y -= ycell;
posCellOffsets[atom].z -= zcell;
}
} }
} }
} }
......
...@@ -466,8 +466,10 @@ public: ...@@ -466,8 +466,10 @@ public:
/** /**
* Reorder the internal arrays of atoms to try to keep spatially contiguous atoms close * Reorder the internal arrays of atoms to try to keep spatially contiguous atoms close
* together in the arrays. * together in the arrays.
*
* @param enforcePeriodic if true, the atom positions may be altered to enforce periodic boundary conditions
*/ */
void reorderAtoms(); void reorderAtoms(bool enforcePeriodic);
/** /**
* Add a listener that should be called whenever atoms get reordered. The OpenCLContext * Add a listener that should be called whenever atoms get reordered. The OpenCLContext
* assumes ownership of the object, and deletes it when the context itself is deleted. * assumes ownership of the object, and deletes it when the context itself is deleted.
......
...@@ -98,8 +98,8 @@ void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, boo ...@@ -98,8 +98,8 @@ void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, boo
OpenCLNonbondedUtilities& nb = cl.getNonbondedUtilities(); OpenCLNonbondedUtilities& nb = cl.getNonbondedUtilities();
bool includeNonbonded = ((groups&(1<<nb.getForceGroup())) != 0); bool includeNonbonded = ((groups&(1<<nb.getForceGroup())) != 0);
cl.setAtomsWereReordered(false); cl.setAtomsWereReordered(false);
if (nb.getUseCutoff() && includeNonbonded && (cl.getComputeForceCount()%100 == 0 || cl.getMoleculesAreInvalid())) { if (cl.getMoleculesAreInvalid() || (nb.getUseCutoff() && includeNonbonded && cl.getComputeForceCount()%100 == 0)) {
cl.reorderAtoms(); cl.reorderAtoms(!cl.getMoleculesAreInvalid());
nb.updateNeighborListSize(); nb.updateNeighborListSize();
cl.setComputeForceCount(cl.getComputeForceCount()+1); cl.setComputeForceCount(cl.getComputeForceCount()+1);
} }
...@@ -4645,9 +4645,10 @@ double OpenCLCalcKineticEnergyKernel::execute(ContextImpl& context) { ...@@ -4645,9 +4645,10 @@ double OpenCLCalcKineticEnergyKernel::execute(ContextImpl& context) {
OpenCLArray<mm_float4>& velm = cl.getVelm(); OpenCLArray<mm_float4>& velm = cl.getVelm();
velm.download(); velm.download();
double energy = 0.0; double energy = 0.0;
OpenCLArray<cl_int>& order = cl.getAtomIndex();
for (size_t i = 0; i < masses.size(); ++i) { for (size_t i = 0; i < masses.size(); ++i) {
mm_float4 v = velm[i]; mm_float4 v = velm[i];
energy += masses[i]*(v.x*v.x+v.y*v.y+v.z*v.z); energy += masses[order[i]]*(v.x*v.x+v.y*v.y+v.z*v.z);
} }
return 0.5*energy; return 0.5*energy;
} }
......
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