Skip to content

GitLab

Commit 9e2b5a12 authored by peastman's avatar peastman
Browse files

Continuing to implement triclinic boxes in reference platform

parent e2977b91
Hide whitespace changes
Inline Side-by-side
Showing with 267 additions and 136 deletions
openmmapi/src/MonteCarloAnisotropicBarostatImpl.cpp
View file @ 9e2b5a12
......@@ -108,7 +108,9 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
Vec3 lengthScale(1.0, 1.0, 1.0);
lengthScale[axis] = newVolume/volume;
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale[0], lengthScale[1], lengthScale[2]);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale[0], box[1]*lengthScale[1], box[2]*lengthScale[2]);
context.getOwner().setPeriodicBoxVectors(Vec3(box[0][0]*lengthScale[0], box[0][1]*lengthScale[1], box[0][2]*lengthScale[2]),
Vec3(box[1][0]*lengthScale[0], box[1][1]*lengthScale[1], box[1][2]*lengthScale[2]),
Vec3(box[2][0]*lengthScale[0], box[2][1]*lengthScale[1], box[2][2]*lengthScale[2]));
// Compute the energy of the modified system.
......
openmmapi/src/MonteCarloMembraneBarostatImpl.cpp
View file @ 9e2b5a12
......@@ -109,7 +109,9 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
}
double deltaArea = box[0][0]*lengthScale[0]*box[1][1]*lengthScale[1] - box[0][0]*box[1][1];
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale[0], lengthScale[1], lengthScale[2]);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale[0], box[1]*lengthScale[1], box[2]*lengthScale[2]);
context.getOwner().setPeriodicBoxVectors(Vec3(box[0][0]*lengthScale[0], box[0][1]*lengthScale[1], box[0][2]*lengthScale[2]),
Vec3(box[1][0]*lengthScale[0], box[1][1]*lengthScale[1], box[1][2]*lengthScale[2]),
Vec3(box[2][0]*lengthScale[0], box[2][1]*lengthScale[1], box[2][2]*lengthScale[2]));
// Compute the energy of the modified system.
......
platforms/reference/include/GBVIParameters.h
View file @ 9e2b5a12
......@@ -65,7 +65,7 @@ class GBVIParameters {
bool _cutoff;
bool _periodic;
RealOpenMM _periodicBoxSize[3];
OpenMM::RealVec _periodicBoxVectors[3];
RealOpenMM _cutoffDistance;
int _bornRadiusScalingMethod;
......@@ -244,11 +244,11 @@ class GBVIParameters {
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic( OpenMM::RealVec& boxSize );
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......@@ -264,7 +264,7 @@ class GBVIParameters {
--------------------------------------------------------------------------------------- */
const RealOpenMM* getPeriodicBox();
const OpenMM::RealVec* getPeriodicBox();
/**---------------------------------------------------------------------------------------
......
platforms/reference/include/ObcParameters.h
View file @ 9e2b5a12
......@@ -64,7 +64,7 @@ class ObcParameters {
bool _cutoff;
bool _periodic;
RealOpenMM _periodicBoxSize[3];
OpenMM::RealVec _periodicBoxVectors[3];
RealOpenMM _cutoffDistance;
/**---------------------------------------------------------------------------------------
......@@ -329,11 +329,11 @@ class ObcParameters {
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic( const OpenMM::RealVec& boxSize );
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......@@ -345,11 +345,11 @@ class ObcParameters {
/**---------------------------------------------------------------------------------------
Get the periodic box dimension
Get the periodic box vectors
--------------------------------------------------------------------------------------- */
const RealOpenMM* getPeriodicBox();
const OpenMM::RealVec* getPeriodicBox();
};
......
platforms/reference/include/ReferenceCustomGBIxn.h
View file @ 9e2b5a12
......@@ -41,7 +41,7 @@ class ReferenceCustomGBIxn {
bool cutoff;
bool periodic;
const OpenMM::NeighborList* neighborList;
RealOpenMM periodicBoxSize[3];
OpenMM::RealVec periodicBoxVectors[3];
RealOpenMM cutoffDistance;
std::vector<Lepton::ExpressionProgram> valueExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
......@@ -263,11 +263,11 @@ class ReferenceCustomGBIxn {
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic( OpenMM::RealVec& boxSize );
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomHbondIxn.h
View file @ 9e2b5a12
......@@ -43,7 +43,7 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
class DihedralTermInfo;
bool cutoff;
bool periodic;
RealOpenMM periodicBoxSize[3];
OpenMM::RealVec periodicBoxVectors[3];
RealOpenMM cutoffDistance;
std::vector<std::vector<int> > donorAtoms, acceptorAtoms;
Lepton::ExpressionProgram energyExpression;
......@@ -111,11 +111,11 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec& boxSize);
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomManyParticleIxn.h
View file @ 9e2b5a12
......@@ -46,7 +46,7 @@ class ReferenceCustomManyParticleIxn {
int numParticlesPerSet, numPerParticleParameters, numTypes;
bool useCutoff, usePeriodic, centralParticleMode;
RealOpenMM cutoffDistance;
RealOpenMM periodicBoxSize[3];
OpenMM::RealVec periodicBoxVectors[3];
Lepton::ExpressionProgram energyExpression;
std::vector<std::vector<std::string> > particleParamNames;
std::vector<std::set<int> > exclusions;
......@@ -118,11 +118,11 @@ class ReferenceCustomManyParticleIxn {
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec& boxSize);
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomNonbondedIxn.h
View file @ 9e2b5a12
......@@ -43,7 +43,7 @@ class ReferenceCustomNonbondedIxn {
bool useSwitch;
bool periodic;
const OpenMM::NeighborList* neighborList;
RealOpenMM periodicBoxSize[3];
OpenMM::RealVec periodicBoxVectors[3];
RealOpenMM cutoffDistance, switchingDistance;
Lepton::CompiledExpression energyExpression;
Lepton::CompiledExpression forceExpression;
......@@ -129,11 +129,11 @@ class ReferenceCustomNonbondedIxn {
already been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic( OpenMM::RealVec& boxSize );
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceMonteCarloBarostat.h
View file @ 9e2b5a12
......@@ -61,14 +61,14 @@ class ReferenceMonteCarloBarostat {
Apply the barostat at the start of a time step, scaling x, y, and z coordinates independently.
@param atomPositions atom positions
@param boxSize the periodic box dimensions
@param boxVectors the periodic box vectors
@param scaleX the factor by which to scale atomic x coordinates
@param scaleY the factor by which to scale atomic y coordinates
@param scaleZ the factor by which to scale atomic z coordinates
--------------------------------------------------------------------------------------- */
void applyBarostat(std::vector<OpenMM::RealVec>& atomPositions, const OpenMM::RealVec& boxSize, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ);
void applyBarostat(std::vector<OpenMM::RealVec>& atomPositions, const OpenMM::RealVec* boxVectors, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ);
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 9e2b5a12
......@@ -279,30 +279,30 @@ void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context,
}
void ReferenceUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
int version = 1;
int version = 2;
stream.write((char*) &version, sizeof(int));
stream.write((char*) &data.time, sizeof(data.time));
vector<RealVec>& posData = extractPositions(context);
stream.write((char*) &posData[0], sizeof(RealVec)*posData.size());
vector<RealVec>& velData = extractVelocities(context);
stream.write((char*) &velData[0], sizeof(RealVec)*velData.size());
RealVec& box = extractBoxSize(context);
stream.write((char*) &box, sizeof(RealVec));
RealVec* vectors = extractBoxVectors(context);
stream.write((char*) vectors, 3*sizeof(RealVec));
SimTKOpenMMUtilities::createCheckpoint(stream);
}
void ReferenceUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& stream) {
int version;
stream.read((char*) &version, sizeof(int));
if (version != 1)
if (version != 2)
throw OpenMMException("Checkpoint was created with a different version of OpenMM");
stream.read((char*) &data.time, sizeof(data.time));
vector<RealVec>& posData = extractPositions(context);
stream.read((char*) &posData[0], sizeof(RealVec)*posData.size());
vector<RealVec>& velData = extractVelocities(context);
stream.read((char*) &velData[0], sizeof(RealVec)*velData.size());
RealVec& box = extractBoxSize(context);
stream.read((char*) &box, sizeof(RealVec));
RealVec* vectors = extractBoxVectors(context);
stream.read((char*) vectors, 3*sizeof(RealVec));
SimTKOpenMMUtilities::loadCheckpoint(stream);
}
......@@ -867,11 +867,11 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
}
if (periodic || ewald || pme) {
RealVec* vectors = extractBoxVectors(context);
RealVec* boxVectors = extractBoxVectors(context);
double minAllowedSize = 1.999999*nonbondedCutoff;
if (vectors[0][0] < minAllowedSize || vectors[1][1] < minAllowedSize || vectors[2][2] < minAllowedSize)
if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
clj.setPeriodic(vectors);
clj.setPeriodic(boxVectors);
}
if (ewald)
clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
......@@ -885,8 +885,8 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
ReferenceLJCoulomb14 nonbonded14;
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
if (periodic || ewald || pme) {
RealVec& boxSize = extractBoxSize(context);
energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
RealVec* boxVectors = extractBoxVectors(context);
energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
}
}
return energy;
......@@ -1022,7 +1022,7 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealVec& box = extractBoxSize(context);
RealVec* boxVectors = extractBoxVectors(context);
RealOpenMM energy = 0;
ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
......@@ -1032,9 +1032,9 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
}
if (periodic) {
double minAllowedSize = 2*nonbondedCutoff;
if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
ixn.setPeriodic(box);
ixn.setPeriodic(boxVectors);
}
if (interactionGroups.size() > 0)
ixn.setInteractionGroups(interactionGroups);
......@@ -1055,7 +1055,7 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
hasInitializedLongRangeCorrection = true;
}
energy += longRangeCoefficient/(box[0]*box[1]*box[2]);
energy += longRangeCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
return energy;
}
......@@ -1119,7 +1119,7 @@ double ReferenceCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool inclu
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
if (isPeriodic)
obc->getObcParameters()->setPeriodic(extractBoxSize(context));
obc->getObcParameters()->setPeriodic(extractBoxVectors(context));
return obc->computeBornEnergyForces(posData, charges, forceData);
}
......@@ -1193,7 +1193,7 @@ double ReferenceCalcGBVIForceKernel::execute(ContextImpl& context, bool includeF
vector<RealVec>& posData = extractPositions(context);
if (isPeriodic)
gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
gbvi->getGBVIParameters()->setPeriodic(extractBoxVectors(context));
RealOpenMM energy;
if (includeForces) {
......@@ -1329,7 +1329,7 @@ double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool incl
energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (periodic)
ixn.setPeriodic(extractBoxSize(context));
ixn.setPeriodic(extractBoxVectors(context));
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
ixn.setUseCutoff(nonbondedCutoff, *neighborList);
......@@ -1508,7 +1508,7 @@ double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool i
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
if (isPeriodic)
ixn->setPeriodic(extractBoxSize(context));
ixn->setPeriodic(extractBoxVectors(context));
RealOpenMM energy = 0;
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
......@@ -1661,11 +1661,11 @@ double ReferenceCalcCustomManyParticleForceKernel::execute(ContextImpl& context,
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
if (nonbondedMethod == CutoffPeriodic) {
RealVec& box = extractBoxSize(context);
RealVec* boxVectors = extractBoxVectors(context);
double minAllowedSize = 2*cutoffDistance;
if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
ixn->setPeriodic(box);
ixn->setPeriodic(boxVectors);
}
ixn->calculateIxn(posData, particleParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
return energy;
......@@ -2014,8 +2014,8 @@ void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& conte
if (barostat == NULL)
barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
vector<RealVec>& posData = extractPositions(context);
RealVec& boxSize = extractBoxSize(context);
barostat->applyBarostat(posData, boxSize, scaleX, scaleY, scaleZ);
RealVec* boxVectors = extractBoxVectors(context);
barostat->applyBarostat(posData, boxVectors, scaleX, scaleY, scaleZ);
}
void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
......
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @ 9e2b5a12
......@@ -118,21 +118,21 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void ReferenceCustomGBIxn::setPeriodic( RealVec& boxSize ) {
void ReferenceCustomGBIxn::setPeriodic(RealVec* vectors) {
if (cutoff) {
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
assert(vectors[0][0] >= 2.0*cutoffDistance);
assert(vectors[1][1] >= 2.0*cutoffDistance);
assert(vectors[2][2] >= 2.0*cutoffDistance);
}
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
periodicBoxVectors[0] = vectors[0];
periodicBoxVectors[1] = vectors[1];
periodicBoxVectors[2] = vectors[2];
}
void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
......@@ -218,7 +218,7 @@ void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2
const map<string, double>& globalParameters, vector<vector<RealOpenMM> >& values) const {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxVectors, deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
......@@ -292,7 +292,7 @@ void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxVectors, deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
......@@ -390,7 +390,7 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vecto
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxVectors, deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
......
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
View file @ 9e2b5a12
......@@ -87,19 +87,19 @@ void ReferenceCustomHbondIxn::setUseCutoff(RealOpenMM distance) {
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void ReferenceCustomHbondIxn::setPeriodic(RealVec& boxSize) {
void ReferenceCustomHbondIxn::setPeriodic(RealVec* vectors) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
assert(vectors[0][0] >= 2.0*cutoffDistance);
assert(vectors[1][1] >= 2.0*cutoffDistance);
assert(vectors[2][2] >= 2.0*cutoffDistance);
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
periodicBoxVectors[0] = vectors[0];
periodicBoxVectors[1] = vectors[1];
periodicBoxVectors[2] = vectors[2];
}
......@@ -288,7 +288,7 @@ void ReferenceCustomHbondIxn::calculateOneIxn(int donor, int acceptor, vector<Re
void ReferenceCustomHbondIxn::computeDelta(int atom1, int atom2, RealOpenMM* delta, vector<RealVec>& atomCoordinates) const {
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], periodicBoxSize, delta);
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], periodicBoxVectors, delta);
else
ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
}
......
platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
View file @ 9e2b5a12
......@@ -120,15 +120,15 @@ void ReferenceCustomManyParticleIxn::setUseCutoff(RealOpenMM distance) {
cutoffDistance = distance;
}
void ReferenceCustomManyParticleIxn::setPeriodic(RealVec& boxSize) {
void ReferenceCustomManyParticleIxn::setPeriodic(RealVec* vectors) {
assert(useCutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
assert(vectors[0][0] >= 2.0*cutoffDistance);
assert(vectors[1][1] >= 2.0*cutoffDistance);
assert(vectors[2][2] >= 2.0*cutoffDistance);
usePeriodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
periodicBoxVectors[0] = vectors[0];
periodicBoxVectors[1] = vectors[1];
periodicBoxVectors[2] = vectors[2];
}
void ReferenceCustomManyParticleIxn::loopOverInteractions(vector<int>& particles, int loopIndex, vector<OpenMM::RealVec>& atomCoordinates,
......@@ -304,7 +304,7 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
void ReferenceCustomManyParticleIxn::computeDelta(int atom1, int atom2, RealOpenMM* delta, vector<RealVec>& atomCoordinates) const {
if (usePeriodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], periodicBoxSize, delta);
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], periodicBoxVectors, delta);
else
ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
}
......
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @ 9e2b5a12
......@@ -131,20 +131,20 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::setPeriodic( RealVec& boxSize ) {
void ReferenceCustomNonbondedIxn::setPeriodic(OpenMM::RealVec* vectors) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
assert(vectors[0][0] >= 2.0*cutoffDistance);
assert(vectors[1][1] >= 2.0*cutoffDistance);
assert(vectors[2][2] >= 2.0*cutoffDistance);
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
periodicBoxVectors[0] = vectors[0];
periodicBoxVectors[1] = vectors[1];
periodicBoxVectors[2] = vectors[2];
}
......@@ -268,7 +268,7 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR );
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
......
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
View file @ 9e2b5a12
......@@ -57,14 +57,14 @@ ReferenceMonteCarloBarostat::~ReferenceMonteCarloBarostat( ) {
Apply the barostat at the start of a time step.
@param atomPositions atom positions
@param boxSize the periodic box dimensions
@param boxVectors the periodic box vectors
@param scaleX the factor by which to scale atom x-coordinates
@param scaleY the factor by which to scale atom y-coordinates
@param scaleZ the factor by which to scale atom z-coordinates
--------------------------------------------------------------------------------------- */
void ReferenceMonteCarloBarostat::applyBarostat(vector<RealVec>& atomPositions, const RealVec& boxSize, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ) {
void ReferenceMonteCarloBarostat::applyBarostat(vector<RealVec>& atomPositions, const RealVec* boxVectors, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ) {
int numAtoms = savedAtomPositions[0].size();
for (int i = 0; i < numAtoms; i++)
for (int j = 0; j < 3; j++)
......@@ -75,39 +75,29 @@ void ReferenceMonteCarloBarostat::applyBarostat(vector<RealVec>& atomPositions,
for (int i = 0; i < (int) molecules.size(); i++) {
// Find the molecule center.
RealOpenMM pos[3] = {0, 0, 0};
RealVec pos(0, 0, 0);
for (int j = 0; j < (int) molecules[i].size(); j++) {
RealVec& atomPos = atomPositions[molecules[i][j]];
pos[0] += atomPos[0];
pos[1] += atomPos[1];
pos[2] += atomPos[2];
pos += atomPos;
}
pos[0] /= molecules[i].size();
pos[1] /= molecules[i].size();
pos[2] /= molecules[i].size();
pos /= molecules[i].size();
// Move it into the first periodic box.
int xcell = (int) floor(pos[0]/boxSize[0]);
int ycell = (int) floor(pos[1]/boxSize[1]);
int zcell = (int) floor(pos[2]/boxSize[2]);
RealOpenMM dx = xcell*boxSize[0];
RealOpenMM dy = ycell*boxSize[1];
RealOpenMM dz = zcell*boxSize[2];
pos[0] -= dx;
pos[1] -= dy;
pos[2] -= dz;
RealVec newPos = pos;
newPos -= boxVectors[2]*floor(newPos[2]/boxVectors[2][2]);
newPos -= boxVectors[1]*floor(newPos[1]/boxVectors[1][1]);
newPos -= boxVectors[0]*floor(newPos[0]/boxVectors[0][0]);
// Now scale the position of the molecule center.
dx = pos[0]*(scaleX-1)-dx;
dy = pos[1]*(scaleY-1)-dy;
dz = pos[2]*(scaleZ-1)-dz;
newPos[0] *= scaleX;
newPos[1] *= scaleY;
newPos[2] *= scaleZ;
RealVec offset = newPos-pos;
for (int j = 0; j < (int) molecules[i].size(); j++) {
RealVec& atomPos = atomPositions[molecules[i][j]];
atomPos[0] += dx;
atomPos[1] += dy;
atomPos[2] += dz;
atomPos += offset;
}
}
}
......
platforms/reference/src/gbsa/GBVIParameters.cpp
View file @ 9e2b5a12
......@@ -286,22 +286,20 @@ RealOpenMM GBVIParameters::getCutoffDistance() {
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void GBVIParameters::setPeriodic( RealVec& boxSize ) {
void GBVIParameters::setPeriodic(RealVec* vectors) {
assert(_cutoff);
assert(boxSize[0] >= 2.0*_cutoffDistance);
assert(boxSize[1] >= 2.0*_cutoffDistance);
assert(boxSize[2] >= 2.0*_cutoffDistance);
_periodic = true;
_periodicBoxSize[0] = boxSize[0];
_periodicBoxSize[1] = boxSize[1];
_periodicBoxSize[2] = boxSize[2];
assert(_cutoff);
assert(vectors[0][0] >= 2.0*_cutoffDistance);
assert(vectors[1][1] >= 2.0*_cutoffDistance);
assert(vectors[2][2] >= 2.0*_cutoffDistance);
_periodic = true;
_periodicBoxVectors[0] = vectors[0];
_periodicBoxVectors[1] = vectors[1];
_periodicBoxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
......@@ -320,8 +318,8 @@ bool GBVIParameters::getPeriodic() {
--------------------------------------------------------------------------------------- */
const RealOpenMM* GBVIParameters::getPeriodicBox() {
return _periodicBoxSize;
const OpenMM::RealVec* GBVIParameters::getPeriodicBox() {
return _periodicBoxVectors;
}
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/gbsa/ObcParameters.cpp
View file @ 9e2b5a12
......@@ -374,22 +374,22 @@ RealOpenMM ObcParameters::getCutoffDistance() const {
also been set, and the smallest side of the periodic box is at least twice the cutoff
distance.
@param boxSize the X, Y, and Z widths of the periodic box
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void ObcParameters::setPeriodic( const OpenMM::RealVec& boxSize ) {
void ObcParameters::setPeriodic(OpenMM::RealVec* vectors) {
assert(_cutoff);
assert(_cutoff);
assert(boxSize[0] >= 2.0*_cutoffDistance);
assert(boxSize[1] >= 2.0*_cutoffDistance);
assert(boxSize[2] >= 2.0*_cutoffDistance);
assert(boxSize[0][0] >= 2.0*_cutoffDistance);
assert(boxSize[1][1] >= 2.0*_cutoffDistance);
assert(boxSize[2][2] >= 2.0*_cutoffDistance);
_periodic = true;
_periodicBoxSize[0] = boxSize[0];
_periodicBoxSize[1] = boxSize[1];
_periodicBoxSize[2] = boxSize[2];
_periodic = true;
_periodicBoxVectors[0] = vectors[0];
_periodicBoxVectors[1] = vectors[1];
_periodicBoxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
......@@ -408,6 +408,6 @@ bool ObcParameters::getPeriodic() {
--------------------------------------------------------------------------------------- */
const RealOpenMM* ObcParameters::getPeriodicBox() {
return _periodicBoxSize;
const OpenMM::RealVec* ObcParameters::getPeriodicBox() {
return _periodicBoxVectors;
}
platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp
View file @ 9e2b5a12
......@@ -227,6 +227,66 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(1.9+1+0.9, state.getPotentialEnergy(), TOL);
}
void testTriclinic() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
Vec3 a(3.1, 0, 0);
Vec3 b(0.4, 3.5, 0);
Vec3 c(-0.1, -0.5, 4.0);
system.setDefaultPeriodicBoxVectors(a, b, c);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("r");
nonbonded->addParticle(vector<double>());
nonbonded->addParticle(vector<double>());
nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
const double cutoff = 1.5;
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int iteration = 0; iteration < 50; iteration++) {
// Generate random positions for the two particles.
positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
context.setPositions(positions);
// Loop over all possible periodic copies and find the nearest one.
Vec3 delta;
double distance2 = 100.0;
for (int i = -1; i < 2; i++)
for (int j = -1; j < 2; j++)
for (int k = -1; k < 2; k++) {
Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
if (d.dot(d) < distance2) {
delta = d;
distance2 = d.dot(d);
}
}
double distance = sqrt(distance2);
// See if the force and energy are correct.
State state = context.getState(State::Forces | State::Energy);
if (distance >= cutoff) {
ASSERT_EQUAL(0.0, state.getPotentialEnergy());
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[0], 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[1], 0);
}
else {
const Vec3 force = delta/sqrt(delta.dot(delta));
ASSERT_EQUAL_TOL(distance, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(force, state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(-force, state.getForces()[1], TOL);
}
}
}
void testContinuous1DFunction() {
ReferencePlatform platform;
System system;
......@@ -863,6 +923,7 @@ int main() {
testExclusions();
testCutoff();
testPeriodic();
testTriclinic();
testContinuous1DFunction();
testContinuous2DFunction();
testContinuous3DFunction();
......
platforms/reference/tests/TestReferenceMonteCarloAnisotropicBarostat.cpp
View file @ 9e2b5a12
......@@ -237,6 +237,83 @@ void testRandomSeed() {
}
}
void testTriclinic() {
const int numParticles = 64;
const int frequency = 10;
const int steps = 1000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
const double temperature = 300.0;
const double initialVolume = numParticles*BOLTZ*temperature/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
ReferencePlatform platform;
System system;
Vec3 initialBox[3];
initialBox[0] = Vec3(initialLength, 0, 0);
initialBox[1] = Vec3(0.2*initialLength, initialLength, 0);
initialBox[2] = Vec3(0.1*initialLength, 0.3*initialLength, initialLength);
system.setDefaultPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(initialLength*genrand_real2(sfmt), initialLength*genrand_real2(sfmt), initialLength*genrand_real2(sfmt));
}
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temperature, true, true, true, frequency);
system.addForce(barostat);
// Run a simulation
LangevinIntegrator integrator(temperature, 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
double expected = (numParticles+1)*BOLTZ*temperature/pressureInMD;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
// Make sure the box vectors have been scaled consistently.
State state = context.getState(State::Positions);
Vec3 box[3];
state.getPeriodicBoxVectors(box[0], box[1], box[2]);
double xscale = box[2][0]/(0.1*initialLength);
double yscale = box[2][1]/(0.3*initialLength);
double zscale = box[2][2]/(1.0*initialLength);
for (int i = 0; i < 3; i++) {
ASSERT_EQUAL_VEC(Vec3(xscale*initialBox[i][0], yscale*initialBox[i][1], zscale*initialBox[i][2]), box[i], 1e-5);
}
// The barostat should have put all particles inside the first periodic box. One integration step
// has happened since then, so they may have moved slightly outside it.
for (int i = 0; i < numParticles; i++) {
Vec3 pos = state.getPositions()[i];
ASSERT(pos[2]/box[2][2] > -1 && pos[2]/box[2][2] < 2);
pos -= box[2]*floor(pos[2]/box[2][2]);
ASSERT(pos[1]/box[1][1] > -1 && pos[1]/box[1][1] < 2);
pos -= box[1]*floor(pos[1]/box[1][1]);
ASSERT(pos[0]/box[0][0] > -1 && pos[0]/box[0][0] < 2);
}
}
/**
* Run a constant pressure simulation on an anisotropic Einstein crystal
* using isotropic and anisotropic barostats. There are a total of 15 simulations:
......@@ -389,6 +466,7 @@ int main() {
testIdealGasAxis(1);
testIdealGasAxis(2);
testRandomSeed();
testTriclinic();
//testEinsteinCrystal();
}
catch(const exception& e) {
......
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