Skip to content

GitLab

Commit 83ed602e authored by peastman's avatar peastman
Browse files

Merge pull request #797 from peastman/triclinic 

C++ libraries support triclinic boxes
parents b51e05a8 050e1262
Hide whitespace changes
Inline Side-by-side
Showing with 708 additions and 402 deletions
platforms/reference/src/ReferenceKernels.cpp
View file @ 83ed602e
......@@ -134,6 +134,11 @@ static RealVec& extractBoxSize(ContextImpl& context) {
return *(RealVec*) data->periodicBoxSize;
}
static RealVec* extractBoxVectors(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return (RealVec*) data->periodicBoxVectors;
}
static ReferenceConstraints& extractConstraints(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *(ReferenceConstraints*) data->constraints;
......@@ -256,10 +261,10 @@ void ReferenceUpdateStateDataKernel::getForces(ContextImpl& context, std::vector
}
void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
RealVec& box = extractBoxSize(context);
a = Vec3(box[0], 0, 0);
b = Vec3(0, box[1], 0);
c = Vec3(0, 0, box[2]);
RealVec* vectors = extractBoxVectors(context);
a = vectors[0];
b = vectors[1];
c = vectors[2];
}
void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const {
......@@ -267,33 +272,37 @@ void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context,
box[0] = (RealOpenMM) a[0];
box[1] = (RealOpenMM) b[1];
box[2] = (RealOpenMM) c[2];
RealVec* vectors = extractBoxVectors(context);
vectors[0] = a;
vectors[1] = b;
vectors[2] = c;
}
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);
}
......@@ -854,15 +863,15 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
bool ewald = (nonbondedMethod == Ewald);
bool pme = (nonbondedMethod == PME);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
}
if (periodic || ewald || pme) {
RealVec& box = extractBoxSize(context);
RealVec* boxVectors = extractBoxVectors(context);
double minAllowedSize = 1.999999*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.");
clj.setPeriodic(box);
clj.setPeriodic(boxVectors);
}
if (ewald)
clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
......@@ -876,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;
......@@ -1013,19 +1022,19 @@ 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);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
ixn.setUseCutoff(nonbondedCutoff, *neighborList);
}
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);
......@@ -1046,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;
}
......@@ -1110,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);
}
......@@ -1184,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) {
......@@ -1320,9 +1329,9 @@ 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, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
ixn.setUseCutoff(nonbondedCutoff, *neighborList);
}
map<string, double> globalParameters;
......@@ -1499,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++)
......@@ -1652,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;
......@@ -2005,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/ReferencePlatform.cpp
View file @ 83ed602e
......@@ -97,6 +97,7 @@ ReferencePlatform::PlatformData::PlatformData(const System& system) : time(0.0),
velocities = new vector<RealVec>(numParticles);
forces = new vector<RealVec>(numParticles);
periodicBoxSize = new RealVec();
periodicBoxVectors = new RealVec[3];
constraints = new ReferenceConstraints(system);
}
......@@ -105,5 +106,6 @@ ReferencePlatform::PlatformData::~PlatformData() {
delete (vector<RealVec>*) velocities;
delete (vector<RealVec>*) forces;
delete (RealVec*) periodicBoxSize;
delete[] (RealVec*) periodicBoxVectors;
delete (ReferenceConstraints*) constraints;
}
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @ 83ed602e
......@@ -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 @ 83ed602e
......@@ -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 @ 83ed602e
......@@ -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 @ 83ed602e
......@@ -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/ReferenceForce.cpp
View file @ 83ed602e
......@@ -41,13 +41,6 @@ using OpenMM::RealVec;
--------------------------------------------------------------------------------------- */
ReferenceForce::ReferenceForce( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceForce::ReferenceForce";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -57,13 +50,6 @@ ReferenceForce::ReferenceForce( ){
--------------------------------------------------------------------------------------- */
ReferenceForce::~ReferenceForce( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceForce::~ReferenceForce";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -79,27 +65,8 @@ RealOpenMM ReferenceForce::periodicDifference(RealOpenMM val1, RealOpenMM val2,
}
/**---------------------------------------------------------------------------------------
Get deltaR and distance and distance**2 between atomI and atomJ (static method)
deltaR: j - i
@param atomCoordinatesI atom i coordinates
@param atomCoordinatesI atom j coordinates
@param deltaR deltaX, deltaY, deltaZ, R2, R upon return
--------------------------------------------------------------------------------------- */
void ReferenceForce::getDeltaR( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
RealOpenMM* deltaR ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nReferenceForce::getDeltaR";
// ---------------------------------------------------------------------------------------
deltaR[XIndex] = atomCoordinatesJ[0] - atomCoordinatesI[0];
deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1];
deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2];
......@@ -108,27 +75,8 @@ void ReferenceForce::getDeltaR( const RealVec& atomCoordinatesI, const RealVec&
deltaR[RIndex] = (RealOpenMM) SQRT( deltaR[R2Index] );
}
/**---------------------------------------------------------------------------------------
Get deltaR and distance and distance**2 between atomI and atomJ, assuming periodic
boundary conditions (static method); deltaR: j - i
@param atomCoordinatesI atom i coordinates
@param atomCoordinatesI atom j coordinates
@param boxSize X, Y, and Z sizes of the periodic box
@param deltaR deltaX, deltaY, deltaZ, R2, R upon return
--------------------------------------------------------------------------------------- */
void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
const RealOpenMM* boxSize, RealOpenMM* deltaR ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nReferenceForce::getDeltaR";
// ---------------------------------------------------------------------------------------
deltaR[XIndex] = periodicDifference(atomCoordinatesJ[0], atomCoordinatesI[0], boxSize[0]);
deltaR[YIndex] = periodicDifference(atomCoordinatesJ[1], atomCoordinatesI[1], boxSize[1]);
deltaR[ZIndex] = periodicDifference(atomCoordinatesJ[2], atomCoordinatesI[2], boxSize[2]);
......@@ -137,29 +85,17 @@ void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const R
deltaR[RIndex] = (RealOpenMM) SQRT( deltaR[R2Index] );
}
/**---------------------------------------------------------------------------------------
Get deltaR between atomI and atomJ (static method); deltaR: j - i
@param atomCoordinatesI atom i coordinates
@param atomCoordinatesI atom j coordinates
@param deltaR deltaX, deltaY, deltaZ upon return
--------------------------------------------------------------------------------------- */
void ReferenceForce::getDeltaROnly( const RealOpenMM* atomCoordinatesI,
const RealOpenMM* atomCoordinatesJ,
RealOpenMM* deltaR ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nReferenceForce::getDeltaR";
// ---------------------------------------------------------------------------------------
deltaR[XIndex] = atomCoordinatesJ[0] - atomCoordinatesI[0];
deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1];
deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2];
void ReferenceForce::getDeltaRPeriodic( const RealVec& atomCoordinatesI, const RealVec& atomCoordinatesJ,
const RealVec* boxVectors, RealOpenMM* deltaR ){
RealVec diff = atomCoordinatesJ-atomCoordinatesI;
diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5);
diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5);
diff -= boxVectors[0]*floor(diff[0]/boxVectors[0][0]+0.5);
deltaR[XIndex] = diff[0];
deltaR[YIndex] = diff[1];
deltaR[ZIndex] = diff[2];
deltaR[R2Index] = diff.dot(diff);
deltaR[RIndex] = SQRT(deltaR[R2Index]);
}
double* ReferenceForce::getVariablePointer(Lepton::CompiledExpression& expression, const std::string& name) {
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ 83ed602e
......@@ -112,20 +112,20 @@ void ReferenceLJCoulombIxn::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 ReferenceLJCoulombIxn::setPeriodic( RealVec& boxSize ) {
void ReferenceLJCoulombIxn::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];
}
/**---------------------------------------------------------------------------------------
......@@ -197,7 +197,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
RealOpenMM factorEwald = -1 / (4*alphaEwald*alphaEwald);
RealOpenMM SQRT_PI = sqrt(PI_M);
RealOpenMM TWO_PI = 2.0 * PI_M;
RealOpenMM recipCoeff = (RealOpenMM)(ONE_4PI_EPS0*4*PI_M/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM recipCoeff = (RealOpenMM)(ONE_4PI_EPS0*4*PI_M/(periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2]) /epsilon);
RealOpenMM totalSelfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
......@@ -230,14 +230,13 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
if (pme && includeReciprocal) {
pme_t pmedata; /* abstract handle for PME data */
RealOpenMM virial[3][3];
pme_init(&pmedata,alphaEwald,numberOfAtoms,meshDim,5,1);
vector<RealOpenMM> charges(numberOfAtoms);
for (int i = 0; i < numberOfAtoms; i++)
charges[i] = atomParameters[i][QIndex];
pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxSize,&recipEnergy,virial);
pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxVectors,&recipEnergy);
if( totalEnergy )
*totalEnergy += recipEnergy;
......@@ -255,7 +254,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
// setup reciprocal box
RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxVectors[0][0], TWO_PI / periodicBoxVectors[1][1], TWO_PI / periodicBoxVectors[2][2]};
// setup K-vectors
......@@ -370,7 +369,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
int jj = pair.second;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM switchValue = 1, switchDeriv = 0;
......@@ -556,7 +555,7 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
// get deltaR, R2, and R between 2 atoms
if (periodic)
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0] );
else
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
......
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
View file @ 83ed602e
......@@ -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/SimTKReference/ReferenceNeighborList.cpp
View file @ 83ed602e
......@@ -12,26 +12,15 @@ namespace OpenMM {
typedef std::vector<AtomIndex> AtomList;
static double periodicDifference(double val1, double val2, double period) {
double diff = val1-val2;
double base = floor(diff/period+0.5)*period;
return diff-base;
}
// squared distance between two points
static double compPairDistanceSquared(const RealVec& pos1, const RealVec& pos2, const RealVec& periodicBoxSize, bool usePeriodic) {
double dx, dy, dz;
if (!usePeriodic) {
dx = pos2[0] - pos1[0];
dy = pos2[1] - pos1[1];
dz = pos2[2] - pos1[2];
}
else {
dx = periodicDifference(pos2[0], pos1[0], periodicBoxSize[0]);
dy = periodicDifference(pos2[1], pos1[1], periodicBoxSize[1]);
dz = periodicDifference(pos2[2], pos1[2], periodicBoxSize[2]);
static double compPairDistanceSquared(const RealVec& pos1, const RealVec& pos2, const RealVec* periodicBoxVectors, bool usePeriodic) {
RealVec diff = pos2-pos1;
if (usePeriodic) {
diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
}
return dx*dx + dy*dy + dz*dz;
return diff.dot(diff);
}
// Ridiculous O(n^2) version of neighbor list
......@@ -41,7 +30,7 @@ void OPENMM_EXPORT computeNeighborListNaive(
int nAtoms,
const AtomLocationList& atomLocations,
const vector<set<int> >& exclusions,
const RealVec& periodicBoxSize,
const RealVec* periodicBoxVectors,
bool usePeriodic,
double maxDistance,
double minDistance,
......@@ -55,7 +44,7 @@ void OPENMM_EXPORT computeNeighborListNaive(
{
for (AtomIndex atomJ = atomI + 1; atomJ < (AtomIndex) nAtoms; ++atomJ)
{
double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxSize, usePeriodic);
double pairDistanceSquared = compPairDistanceSquared(atomLocations[atomI], atomLocations[atomJ], periodicBoxVectors, usePeriodic);
if ( (pairDistanceSquared <= maxDistanceSquared) && (pairDistanceSquared >= minDistanceSquared))
if (exclusions[atomI].find(atomJ) == exclusions[atomI].end())
{
......@@ -95,15 +84,15 @@ typedef std::vector< VoxelItem > Voxel;
class VoxelHash
{
public:
VoxelHash(double vsx, double vsy, double vsz, const RealVec& periodicBoxSize, bool usePeriodic) :
voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxSize(periodicBoxSize), usePeriodic(usePeriodic) {
VoxelHash(double vsx, double vsy, double vsz, const RealVec* periodicBoxVectors, bool usePeriodic) :
voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxVectors(periodicBoxVectors), usePeriodic(usePeriodic) {
if (usePeriodic) {
nx = (int) floor(periodicBoxSize[0]/voxelSizeX+0.5);
ny = (int) floor(periodicBoxSize[1]/voxelSizeY+0.5);
nz = (int) floor(periodicBoxSize[2]/voxelSizeZ+0.5);
voxelSizeX = periodicBoxSize[0]/nx;
voxelSizeY = periodicBoxSize[1]/ny;
voxelSizeZ = periodicBoxSize[2]/nz;
nx = (int) floor(periodicBoxVectors[0][0]/voxelSizeX+0.5);
ny = (int) floor(periodicBoxVectors[1][1]/voxelSizeY+0.5);
nz = (int) floor(periodicBoxVectors[2][2]/voxelSizeZ+0.5);
voxelSizeX = periodicBoxVectors[0][0]/nx;
voxelSizeY = periodicBoxVectors[1][1]/ny;
voxelSizeZ = periodicBoxVectors[2][2]/nz;
}
}
......@@ -118,21 +107,16 @@ public:
VoxelIndex getVoxelIndex(const RealVec& location) const {
double xperiodic, yperiodic, zperiodic;
if (!usePeriodic) {
xperiodic = location[0];
yperiodic = location[1];
zperiodic = location[2];
}
else {
xperiodic = location[0]-periodicBoxSize[0]*floor(location[0]/periodicBoxSize[0]);
yperiodic = location[1]-periodicBoxSize[1]*floor(location[1]/periodicBoxSize[1]);
zperiodic = location[2]-periodicBoxSize[2]*floor(location[2]/periodicBoxSize[2]);
RealVec r = location;
if (usePeriodic) {
r -= periodicBoxVectors[2]*floor(r[2]/periodicBoxVectors[2][2]);
r -= periodicBoxVectors[1]*floor(r[1]/periodicBoxVectors[1][1]);
r -= periodicBoxVectors[0]*floor(r[0]/periodicBoxVectors[0][0]);
}
int x = int(floor(xperiodic / voxelSizeX));
int y = int(floor(yperiodic / voxelSizeY));
int z = int(floor(zperiodic / voxelSizeZ));
int x = int(floor(r[0]/voxelSizeX));
int y = int(floor(r[1]/voxelSizeY));
int z = int(floor(r[2]/voxelSizeZ));
return VoxelIndex(x, y, z);
}
......@@ -163,19 +147,33 @@ public:
int dIndexY = int(maxDistance / voxelSizeY) + 1;
int dIndexZ = int(maxDistance / voxelSizeZ) + 1;
VoxelIndex centerVoxelIndex = getVoxelIndex(locationI);
int lastx = centerVoxelIndex.x+dIndexX;
int lasty = centerVoxelIndex.y+dIndexY;
int lastz = centerVoxelIndex.z+dIndexZ;
if (usePeriodic) {
lastx = min(lastx, centerVoxelIndex.x-dIndexX+nx-1);
lasty = min(lasty, centerVoxelIndex.y-dIndexY+ny-1);
lastz = min(lastz, centerVoxelIndex.z-dIndexZ+nz-1);
}
for (int x = centerVoxelIndex.x - dIndexX; x <= lastx; ++x)
int minz = centerVoxelIndex.z-dIndexZ;
int maxz = centerVoxelIndex.z+dIndexZ;
if (usePeriodic)
maxz = min(maxz, minz+nz-1);
for (int z = minz; z <= maxz; ++z)
{
for (int y = centerVoxelIndex.y - dIndexY; y <= lasty; ++y)
int boxz = (int) floor((float) z/nz);
int miny = centerVoxelIndex.y-dIndexY;
int maxy = centerVoxelIndex.y+dIndexY;
if (usePeriodic) {
double yoffset = boxz*periodicBoxVectors[2][1]/voxelSizeY;
miny -= (int) ceil(yoffset);
maxy -= (int) floor(yoffset);
maxy = min(maxy, miny+ny-1);
}
for (int y = miny; y <= maxy; ++y)
{
for (int z = centerVoxelIndex.z - dIndexZ; z <= lastz; ++z)
int boxy = (int) floor((float) y/ny);
int minx = centerVoxelIndex.x-dIndexX;
int maxx = centerVoxelIndex.x+dIndexX;
if (usePeriodic) {
double xoffset = (boxy*periodicBoxVectors[1][0]+boxz*periodicBoxVectors[2][0])/voxelSizeX;
minx -= (int) ceil(xoffset);
maxx -= (int) floor(xoffset);
maxx = min(maxx, minx+nx-1);
}
for (int x = minx; x <= maxx; ++x)
{
VoxelIndex voxelIndex(x, y, z);
if (usePeriodic) {
......@@ -194,7 +192,7 @@ public:
// Ignore self hits
if (atomI == atomJ) continue;
double dSquared = compPairDistanceSquared(locationI, locationJ, periodicBoxSize, usePeriodic);
double dSquared = compPairDistanceSquared(locationI, locationJ, periodicBoxVectors, usePeriodic);
if (dSquared > maxDistanceSquared) continue;
if (dSquared < minDistanceSquared) continue;
......@@ -213,7 +211,7 @@ public:
private:
double voxelSizeX, voxelSizeY, voxelSizeZ;
int nx, ny, nz;
const RealVec& periodicBoxSize;
const RealVec* periodicBoxVectors;
const bool usePeriodic;
std::map<VoxelIndex, Voxel> voxelMap;
};
......@@ -223,9 +221,9 @@ private:
void OPENMM_EXPORT computeNeighborListVoxelHash(
NeighborList& neighborList,
int nAtoms,
const AtomLocationList& atomLocations,
const AtomLocationList& atomLocations,
const vector<set<int> >& exclusions,
const RealVec& periodicBoxSize,
const RealVec* periodicBoxVectors,
bool usePeriodic,
double maxDistance,
double minDistance,
......@@ -238,11 +236,11 @@ void OPENMM_EXPORT computeNeighborListVoxelHash(
if (!usePeriodic)
edgeSizeX = edgeSizeY = edgeSizeZ = maxDistance; // TODO - adjust this as needed
else {
edgeSizeX = 0.5*periodicBoxSize[0]/floor(periodicBoxSize[0]/maxDistance);
edgeSizeY = 0.5*periodicBoxSize[1]/floor(periodicBoxSize[1]/maxDistance);
edgeSizeZ = 0.5*periodicBoxSize[2]/floor(periodicBoxSize[2]/maxDistance);
edgeSizeX = 0.5*periodicBoxVectors[0][0]/floor(periodicBoxVectors[0][0]/maxDistance);
edgeSizeY = 0.5*periodicBoxVectors[1][1]/floor(periodicBoxVectors[1][1]/maxDistance);
edgeSizeZ = 0.5*periodicBoxVectors[2][2]/floor(periodicBoxVectors[2][2]/maxDistance);
}
VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxSize, usePeriodic);
VoxelHash voxelHash(edgeSizeX, edgeSizeY, edgeSizeZ, periodicBoxVectors, usePeriodic);
for (AtomIndex atomJ = 0; atomJ < (AtomIndex) nAtoms; ++atomJ) // use "j", because j > i for pairs
{
// 1) Find other atoms that are close to this one
......
platforms/reference/src/SimTKReference/ReferencePME.cpp
View file @ 83ed602e
......@@ -192,11 +192,21 @@ pme_calculate_bsplines_moduli(pme_t pme)
}
static void invert_box_vectors(const RealVec boxVectors[3], RealVec recipBoxVectors[3])
{
RealOpenMM determinant = boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2];
assert(determinant > 0);
RealOpenMM scale = 1.0/determinant;
recipBoxVectors[0] = RealVec(boxVectors[1][1]*boxVectors[2][2], 0, 0)*scale;
recipBoxVectors[1] = RealVec(-boxVectors[1][0]*boxVectors[2][2], boxVectors[0][0]*boxVectors[2][2], 0)*scale;
recipBoxVectors[2] = RealVec(boxVectors[1][0]*boxVectors[2][1]-boxVectors[1][1]*boxVectors[2][0], -boxVectors[0][0]*boxVectors[2][1], boxVectors[0][0]*boxVectors[1][1])*scale;
}
static void
pme_update_grid_index_and_fraction(pme_t pme,
const vector<RealVec>& atomCoordinates,
const RealOpenMM periodicBoxSize[3])
const RealVec periodicBoxVectors[3],
const RealVec recipBoxVectors[3])
{
int i;
int d;
......@@ -205,47 +215,47 @@ pme_update_grid_index_and_fraction(pme_t pme,
for(i=0;i<pme->natoms;i++)
{
/* Index calculation (Look mom, no conditionals!):
*
* Both for Cuda and modern CPUs it is nice to avoid conditionals, but we still need to apply periodic boundary conditions.
* Instead of having loops to add/subtract the box dimension, we do it this way:
*
* 1. First add the box size, to make sure this atom coordinate isnt -0.1 or something.
* After this we assume all fractional box positions are *positive*.
* The reason for this is that we always want to round coordinates _down_ to get
* their grid index, and when taking the integer part of -3.4 we would get -3, not -4 as we want.
* Since we anyway need the grid indices to fall in the central box, it is more convenient
* to first manipulate the coordinates to be positive.
* 2. Convert to integer grid index
* Since we have added a whole box unit in step 1, this index might actually be larger than
* the grid dimension. Examples, assuming 10*10*10nm box and grid dimension 100*100*100 (spacing 0.1 nm):
*
* coordinate is { 0.543 , 6.235 , -0.73 }
*
* x[i][d]/box[d] becomes { 0.0543 , 0.6235 , -0.073 }
* (x[i][d]/box[d] + 1.0) becomes { 1.0543 , 1.6235 , 0.927 }
* (x[i][d]/box[d] + 1.0)*ngrid[d] becomes { 105.43 , 162.35 , 92.7 }
*
* integer part is now { 105 , 162 , 92 }
*
* The fraction is calculates as t-ti, which becomes { 0.43 , 0.35 , 0.7 }
*
* 3. Take the first integer index part (which can be larger than the grid) modulo the grid dimension
*
* Now we get { 5 , 62 , 92 }
*
* Voila, both index and fraction, entirely without conditionals. The one limitation here is that
* we only add one box length, so if the particle had a coordinate <=-10.0, we would be screwed.
* In principle we can of course add 100.0, but that just moves the problem, it doesnt solve it.
* In practice, MD programs will apply PBC to reset particles inside the central box to avoid
* numerical problems, so this shouldnt cause any problems.
* (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!)
*/
RealVec coord = atomCoordinates[i];
for(d=0;d<3;d++)
{
/* Index calculation (Look mom, no conditionals!):
*
* Both for Cuda and modern CPUs it is nice to avoid conditionals, but we still need to apply periodic boundary conditions.
* Instead of having loops to add/subtract the box dimension, we do it this way:
*
* 1. First add the box size, to make sure this atom coordinate isnt -0.1 or something.
* After this we assume all fractional box positions are *positive*.
* The reason for this is that we always want to round coordinates _down_ to get
* their grid index, and when taking the integer part of -3.4 we would get -3, not -4 as we want.
* Since we anyway need the grid indices to fall in the central box, it is more convenient
* to first manipulate the coordinates to be positive.
* 2. Convert to integer grid index
* Since we have added a whole box unit in step 1, this index might actually be larger than
* the grid dimension. Examples, assuming 10*10*10nm box and grid dimension 100*100*100 (spacing 0.1 nm):
*
* coordinate is { 0.543 , 6.235 , -0.73 }
*
* x[i][d]/box[d] becomes { 0.0543 , 0.6235 , -0.073 }
* (x[i][d]/box[d] + 1.0) becomes { 1.0543 , 1.6235 , 0.927 }
* (x[i][d]/box[d] + 1.0)*ngrid[d] becomes { 105.43 , 162.35 , 92.7 }
*
* integer part is now { 105 , 162 , 92 }
*
* The fraction is calculates as t-ti, which becomes { 0.43 , 0.35 , 0.7 }
*
* 3. Take the first integer index part (which can be larger than the grid) modulo the grid dimension
*
* Now we get { 5 , 62 , 92 }
*
* Voila, both index and fraction, entirely without conditionals. The one limitation here is that
* we only add one box length, so if the particle had a coordinate <=-10.0, we would be screwed.
* In principle we can of course add 100.0, but that just moves the problem, it doesnt solve it.
* In practice, MD programs will apply PBC to reset particles inside the central box to avoid
* numerical problems, so this shouldnt cause any problems.
* (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!)
*/
RealOpenMM coord = atomCoordinates[i][d];
coord -= floor(coord/periodicBoxSize[d])*periodicBoxSize[d];
t = (coord / periodicBoxSize[d])*pme->ngrid[d];
t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d];
t = (t-floor(t))*pme->ngrid[d];
ti = (int) t;
pme->particlefraction[i][d] = t - ti;
......@@ -397,9 +407,9 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges)
static void
pme_reciprocal_convolution(pme_t pme,
const RealOpenMM periodicBoxSize[3],
RealOpenMM * energy,
RealOpenMM pme_virial[3][3])
const RealVec periodicBoxVectors[3],
const RealVec recipBoxVectors[3],
RealOpenMM * energy)
{
int kx,ky,kz;
int nx,ny,nz;
......@@ -424,7 +434,7 @@ pme_reciprocal_convolution(pme_t pme,
one_4pi_eps = (RealOpenMM) (ONE_4PI_EPS0/pme->epsilon_r);
factor = (RealOpenMM) (M_PI*M_PI/(pme->ewaldcoeff*pme->ewaldcoeff));
boxfactor = (RealOpenMM) (M_PI*periodicBoxSize[0]*periodicBoxSize[1]*periodicBoxSize[2]);
boxfactor = (RealOpenMM) (M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
esum = 0;
virxx = 0;
......@@ -442,14 +452,14 @@ pme_reciprocal_convolution(pme_t pme,
{
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mx = (RealOpenMM) ((kx<maxkx) ? kx : (kx-nx));
mhx = mx/periodicBoxSize[0];
mhx = mx*recipBoxVectors[0][0];
bx = boxfactor*pme->bsplines_moduli[0][kx];
for(ky=0;ky<ny;ky++)
{
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
my = (RealOpenMM) ((ky<maxky) ? ky : (ky-ny));
mhy = my/periodicBoxSize[1];
mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
by = pme->bsplines_moduli[1][ky];
for(kz=0;kz<nz;kz++)
......@@ -469,7 +479,7 @@ pme_reciprocal_convolution(pme_t pme,
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mz = (RealOpenMM) ((kz<maxkz) ? kz : (kz-nz));
mhz = mz/periodicBoxSize[2];
mhz = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
/* Pointer to the grid cell in question */
ptr = pme->grid + kx*ny*nz + ky*nz + kz;
......@@ -494,24 +504,9 @@ pme_reciprocal_convolution(pme_t pme,
/* Long-range PME contribution to the energy for this frequency */
ets2 = eterm*struct2;
esum += ets2;
/* PME long-range contribution to atomic virial. Since it is symmetric, we only calculate half the matrix inside this loop. */
vfactor = (factor*m2+1)*2/m2;
virxx += ets2*(vfactor*mhx*mhx-1);
virxy += ets2*vfactor*mhx*mhy;
virxz += ets2*vfactor*mhx*mhz;
viryy += ets2*(vfactor*mhy*mhy-1);
viryz += ets2*vfactor*mhy*mhz;
virzz += ets2*(vfactor*mhz*mhz-1);
}
}
}
pme_virial[0][0] = (RealOpenMM) (0.25*virxx);
pme_virial[1][1] = (RealOpenMM) (0.25*viryy);
pme_virial[2][2] = (RealOpenMM) (0.25*virzz);
pme_virial[0][1] = pme_virial[1][0] = (RealOpenMM) (0.25*virxy);
pme_virial[0][2] = pme_virial[2][0] = (RealOpenMM) (0.25*virxz);
pme_virial[1][2] = pme_virial[2][1] = (RealOpenMM) (0.25*viryz);
/* The factor 0.5 is nothing special, but it is better to have it here than inside the loop :-) */
*energy = (RealOpenMM) (0.5*esum);
......@@ -520,7 +515,7 @@ pme_reciprocal_convolution(pme_t pme,
static void
pme_grid_interpolate_force(pme_t pme,
const RealOpenMM periodicBoxSize[3],
const RealVec recipBoxVectors[3],
const vector<RealOpenMM>& charges,
vector<RealVec>& forces)
{
......@@ -610,9 +605,9 @@ pme_grid_interpolate_force(pme_t pme,
}
}
/* Update memory force, note that we multiply by charge and some box stuff */
forces[i][0] -= q*fx*nx/periodicBoxSize[0];
forces[i][1] -= q*fy*ny/periodicBoxSize[1];
forces[i][2] -= q*fz*nz/periodicBoxSize[2];
forces[i][0] -= q*(fx*nx*recipBoxVectors[0][0]);
forces[i][1] -= q*(fx*nx*recipBoxVectors[1][0]+fy*ny*recipBoxVectors[1][1]);
forces[i][2] -= q*(fx*nx*recipBoxVectors[2][0]+fy*ny*recipBoxVectors[2][1]+fz*nz*recipBoxVectors[2][2]);
}
}
......@@ -669,12 +664,14 @@ int pme_exec(pme_t pme,
const vector<RealVec>& atomCoordinates,
vector<RealVec>& forces,
const vector<RealOpenMM>& charges,
const RealOpenMM periodicBoxSize[3],
RealOpenMM* energy,
RealOpenMM pme_virial[3][3])
const RealVec periodicBoxVectors[3],
RealOpenMM* energy)
{
/* Routine is called with coordinates in x, a box, and charges in q */
RealVec recipBoxVectors[3];
invert_box_vectors(periodicBoxVectors, recipBoxVectors);
/* Before we can do the actual interpolation, we need to recalculate and update
* the indices for each particle in the charge grid (initialized in pme_init()),
* and what its fractional offset in this grid cell is.
......@@ -683,7 +680,7 @@ int pme_exec(pme_t pme,
/* Update charge grid indices and fractional offsets for each atom.
* The indices/fractions are stored internally in the pme datatype
*/
pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxSize);
pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxVectors,recipBoxVectors);
/* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
pme_update_bsplines(pme);
......@@ -695,13 +692,13 @@ int pme_exec(pme_t pme,
fftpack_exec_3d(pme->fftplan,FFTPACK_FORWARD,pme->grid,pme->grid);
/* solve in k-space */
pme_reciprocal_convolution(pme,periodicBoxSize,energy,pme_virial);
pme_reciprocal_convolution(pme,periodicBoxVectors,recipBoxVectors,energy);
/* do 3d-invfft */
fftpack_exec_3d(pme->fftplan,FFTPACK_BACKWARD,pme->grid,pme->grid);
/* Get the particle forces from the grid and bsplines in the pme structure */
pme_grid_interpolate_force(pme,periodicBoxSize,charges,forces);
pme_grid_interpolate_force(pme,recipBoxVectors,charges,forces);
return 0;
}
......
platforms/reference/src/gbsa/GBVIParameters.cpp
View file @ 83ed602e
......@@ -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 @ 83ed602e
......@@ -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/TestReferenceCustomManyParticleForce.cpp
View file @ 83ed602e
......@@ -53,7 +53,17 @@ using namespace std;
const double TOL = 1e-5;
void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
Vec3 computeDelta(const Vec3& pos1, const Vec3& pos2, bool periodic, const Vec3* periodicBoxVectors) {
Vec3 diff = pos1-pos2;
if (periodic) {
diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
}
return diff;
}
void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize, bool triclinic) {
// Create a System and Context.
int numParticles = force->getNumParticles();
......@@ -61,7 +71,18 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
Vec3 boxVectors[3];
if (triclinic) {
boxVectors[0] = Vec3(boxSize, 0, 0);
boxVectors[1] = Vec3(0.2*boxSize, boxSize, 0);
boxVectors[2] = Vec3(-0.3*boxSize, -0.1*boxSize, boxSize);
}
else {
boxVectors[0] = Vec3(boxSize, 0, 0);
boxVectors[1] = Vec3(0, boxSize, 0);
boxVectors[2] = Vec3(0, 0, boxSize);
}
system.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
system.addForce(force);
if (force->getNonbondedMethod() == CustomManyParticleForce::CutoffPeriodic) {
ASSERT(force->usesPeriodicBoundaryConditions());
......@@ -81,20 +102,14 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
// See if the energy matches the expected value.
double expectedEnergy = 0;
bool periodic = (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic);
for (int i = 0; i < (int) expectedSets.size(); i++) {
int p1 = expectedSets[i][0];
int p2 = expectedSets[i][1];
int p3 = expectedSets[i][2];
Vec3 d12 = positions[p2]-positions[p1];
Vec3 d13 = positions[p3]-positions[p1];
Vec3 d23 = positions[p3]-positions[p2];
if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
for (int j = 0; j < 3; j++) {
d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
}
}
Vec3 d12 = computeDelta(positions[p2], positions[p1], periodic, boxVectors);
Vec3 d13 = computeDelta(positions[p3], positions[p1], periodic, boxVectors);
Vec3 d23 = computeDelta(positions[p3], positions[p2], periodic, boxVectors);
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
......@@ -218,7 +233,7 @@ void testNoCutoff() {
positions.push_back(Vec3(0.4, 0, -0.8));
int sets[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
vector<const int*> expectedSets(&sets[0], &sets[4]);
validateAxilrodTeller(force, positions, expectedSets, 2.0);
validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
}
void testCutoff() {
......@@ -243,7 +258,7 @@ void testCutoff() {
positions.push_back(Vec3(0.2, 0.5, -0.1));
int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[7]);
validateAxilrodTeller(force, positions, expectedSets, 2.0);
validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
}
void testPeriodic() {
......@@ -269,7 +284,33 @@ void testPeriodic() {
double boxSize = 2.1;
int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[5]);
validateAxilrodTeller(force, positions, expectedSets, boxSize);
validateAxilrodTeller(force, positions, expectedSets, boxSize, false);
}
void testTriclinic() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
force->setCutoffDistance(1.05);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
positions.push_back(Vec3(0.2, 0.5, -0.1));
double boxSize = 2.1;
int sets[4][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {1,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[4]);
validateAxilrodTeller(force, positions, expectedSets, boxSize, true);
}
void testExclusions() {
......@@ -294,7 +335,7 @@ void testExclusions() {
force->addExclusion(0, 3);
int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[5]);
validateAxilrodTeller(force, positions, expectedSets, 2.0);
validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
}
void testAllTerms() {
......@@ -600,6 +641,7 @@ int main() {
testNoCutoff();
testCutoff();
testPeriodic();
testTriclinic();
testExclusions();
testAllTerms();
testParameters();
......
platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp
View file @ 83ed602e
......@@ -231,6 +231,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;
......@@ -869,6 +929,7 @@ int main() {
testExclusions();
testCutoff();
testPeriodic();
testTriclinic();
testContinuous1DFunction();
testContinuous2DFunction();
testContinuous3DFunction();
......
platforms/reference/tests/TestReferenceEwald.cpp
View file @ 83ed602e
......@@ -302,6 +302,58 @@ void testWaterSystem() {
}
void testTriclinic() {
// Create a triclinic box containing eight particles.
ReferencePlatform platform;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(2.5, 0, 0), Vec3(0.5, 3.0, 0), Vec3(0.7, 0.9, 3.5));
for (int i = 0; i < 8; i++)
system.addParticle(1.0);
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
force->setNonbondedMethod(NonbondedForce::PME);
force->setCutoffDistance(1.0);
force->setPMEParameters(3.45891, 32, 40, 48);
for (int i = 0; i < 4; i++)
force->addParticle(-1, 0.440104, 0.4184); // Cl parameters
for (int i = 0; i < 4; i++)
force->addParticle(1, 0.332840, 0.0115897); // Na parameters
vector<Vec3> positions(8);
positions[0] = Vec3(1.744, 2.788, 3.162);
positions[1] = Vec3(1.048, 0.762, 2.340);
positions[2] = Vec3(2.489, 1.570, 2.817);
positions[3] = Vec3(1.027, 1.893, 3.271);
positions[4] = Vec3(0.937, 0.825, 0.009);
positions[5] = Vec3(2.290, 1.887, 3.352);
positions[6] = Vec3(1.266, 1.111, 2.894);
positions[7] = Vec3(0.933, 1.862, 3.490);
// Compute the forces and energy.
VerletIntegrator integ(0.001);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// Compare them to values computed by Gromacs.
double expectedEnergy = -963.370;
vector<Vec3> expectedForce(8);
expectedForce[0] = Vec3(4.25253e+01, -1.23503e+02, 1.22139e+02);
expectedForce[1] = Vec3(9.74752e+01, 1.68213e+02, 1.93169e+02);
expectedForce[2] = Vec3(-1.50348e+02, 1.29165e+02, 3.70435e+02);
expectedForce[3] = Vec3(9.18644e+02, -3.52571e+00, -1.34772e+03);
expectedForce[4] = Vec3(-1.61193e+02, 9.01528e+01, -7.12904e+01);
expectedForce[5] = Vec3(2.82630e+02, 2.78029e+01, -3.72864e+02);
expectedForce[6] = Vec3(-1.47454e+02, -2.14448e+02, -3.55789e+02);
expectedForce[7] = Vec3(-8.82195e+02, -7.39132e+01, 1.46202e+03);
for (int i = 0; i < 8; i++) {
ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-4);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-4);
}
void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
// Create a cloud of random point charges.
......@@ -409,6 +461,7 @@ int main() {
testEwaldPME();
// testEwald2Ions();
// testWaterSystem();
testTriclinic();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
testPMEParameters();
......
platforms/reference/tests/TestReferenceMonteCarloAnisotropicBarostat.cpp
View file @ 83ed602e
......@@ -243,6 +243,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:
......@@ -395,6 +472,7 @@ int main() {
testIdealGasAxis(1);
testIdealGasAxis(2);
testRandomSeed();
testTriclinic();
//testEinsteinCrystal();
}
catch(const exception& e) {
......
platforms/reference/tests/TestReferenceNeighborList.cpp
View file @ 83ed602e
......@@ -48,43 +48,38 @@ void testNeighborList()
NeighborList neighborList;
RealVec boxSize;
computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxSize, false, 13.7, 0.01);
RealVec boxVectors[3];
computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxVectors, false, 13.7, 0.01);
assert(neighborList.size() == 1);
computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxSize, false, 13.5, 0.01);
computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxVectors, false, 13.5, 0.01);
assert(neighborList.size() == 0);
computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxSize, false, 13.7, 0.01);
computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxVectors, false, 13.7, 0.01);
assert(neighborList.size() == 1);
computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxSize, false, 13.5, 0.01);
computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, boxVectors, false, 13.5, 0.01);
assert(neighborList.size() == 0);
}
double periodicDifference(double val1, double val2, double period) {
double diff = val1-val2;
double base = floor(diff/period+0.5)*period;
return diff-base;
double distance2(RealVec& pos1, RealVec& pos2, const RealVec* periodicBoxVectors) {
RealVec diff = pos1-pos2;
diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
return diff.dot(diff);
}
double distance2(RealVec& pos1, RealVec& pos2, const RealVec& periodicBoxSize) {
double dx = periodicDifference(pos1[0], pos2[0], periodicBoxSize[0]);
double dy = periodicDifference(pos1[1], pos2[1], periodicBoxSize[1]);
double dz = periodicDifference(pos1[2], pos2[2], periodicBoxSize[2]);
return dx*dx+dy*dy+dz*dz;
}
void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& positions, const RealVec& periodicBoxSize, double cutoff) {
void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& positions, const RealVec* periodicBoxVectors, double cutoff) {
for (int i = 0; i < (int) list.size(); i++) {
int particle1 = list[i].first;
int particle2 = list[i].second;
ASSERT(distance2(positions[particle1], positions[particle2], periodicBoxSize) <= cutoff*cutoff);
ASSERT(distance2(positions[particle1], positions[particle2], periodicBoxVectors) <= cutoff*cutoff);
}
int count = 0;
for (int i = 0; i < numParticles; i++)
for (int j = i+1; j < numParticles; j++)
if (distance2(positions[i], positions[j], periodicBoxSize) <= cutoff*cutoff)
if (distance2(positions[i], positions[j], periodicBoxVectors) <= cutoff*cutoff)
count++;
ASSERT_EQUAL(count, list.size());
}
......@@ -92,22 +87,49 @@ void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& p
void testPeriodic() {
const int numParticles = 100;
const double cutoff = 3.0;
const RealVec periodicBoxSize(20.0, 15.0, 22.0);
RealVec periodicBoxVectors[3];
periodicBoxVectors[0] = RealVec(20, 0, 0);
periodicBoxVectors[1] = RealVec(0, 15, 0);
periodicBoxVectors[2] = RealVec(0, 0, 22);
vector<RealVec> particleList(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i <numParticles; i++) {
particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[0][0]*3);
particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[1][1]*3);
particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[2][2]*3);
}
vector<set<int> > exclusions(numParticles);
NeighborList neighborList;
computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
}
void testTriclinic() {
const int numParticles = 1000;
const double cutoff = 3.0;
RealVec periodicBoxVectors[3];
periodicBoxVectors[0] = RealVec(20, 0, 0);
periodicBoxVectors[1] = RealVec(5, 15, 0);
periodicBoxVectors[2] = RealVec(-3, -7, 22);
vector<RealVec> particleList(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i <numParticles; i++) {
particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxSize[0]*3);
particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxSize[1]*3);
particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxSize[2]*3);
particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[0][0]*3);
particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[1][1]*3);
particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[2][2]*3);
}
vector<set<int> > exclusions(numParticles);
NeighborList neighborList;
computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxSize, true, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxSize, true, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxVectors, true, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxVectors, cutoff);
}
int main()
......@@ -115,6 +137,7 @@ int main()
try {
testNeighborList();
testPeriodic();
testTriclinic();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/reference/tests/TestReferenceNonbondedForce.cpp
View file @ 83ed602e
......@@ -41,6 +41,7 @@
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "openmm/HarmonicBondForce.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
......@@ -361,6 +362,69 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), 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);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::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);
const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
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*ONE_4PI_EPS0*(-1.0/(distance*distance*distance)+2.0*krf);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(1.0/distance+krf*distance*distance-crf), state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(force, state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(-force, state.getForces()[1], TOL);
}
}
}
void testDispersionCorrection() {
// Create a box full of identical particles.
......@@ -503,6 +567,7 @@ int main() {
testCutoff();
testCutoff14();
testPeriodic();
testTriclinic();
testDispersionCorrection();
testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
testSwitchingFunction(NonbondedForce::PME);
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 83ed602e
......@@ -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-2013 Stanford University and the Authors. *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
......@@ -796,12 +796,12 @@ private:
CudaCalcAmoebaMultipoleForceKernel::CudaCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) :
CalcAmoebaMultipoleForceKernel(name, platform), cu(cu), system(system), hasInitializedScaleFactors(false), hasInitializedFFT(false), multipolesAreValid(false),
multipoleParticles(NULL), molecularDipoles(NULL), molecularQuadrupoles(NULL), labFrameDipoles(NULL), labFrameQuadrupoles(NULL),
field(NULL), fieldPolar(NULL), inducedField(NULL), inducedFieldPolar(NULL), torque(NULL), dampingAndThole(NULL), inducedDipole(NULL),
multipoleParticles(NULL), molecularDipoles(NULL), molecularQuadrupoles(NULL), labFrameDipoles(NULL), labFrameQuadrupoles(NULL), fracDipoles(NULL),
fracQuadrupoles(NULL), field(NULL), fieldPolar(NULL), inducedField(NULL), inducedFieldPolar(NULL), torque(NULL), dampingAndThole(NULL), inducedDipole(NULL),
diisCoefficients(NULL), inducedDipolePolar(NULL), inducedDipoleErrors(NULL), prevDipoles(NULL), prevDipolesPolar(NULL), prevDipolesGk(NULL),
prevDipolesGkPolar(NULL), prevErrors(NULL), diisMatrix(NULL), polarizability(NULL), covalentFlags(NULL), polarizationGroupFlags(NULL),
pmeGrid(NULL), pmeBsplineModuliX(NULL), pmeBsplineModuliY(NULL), pmeBsplineModuliZ(NULL), pmeIgrid(NULL), pmePhi(NULL),
pmePhid(NULL), pmePhip(NULL), pmePhidp(NULL), pmeAtomGridIndex(NULL), lastPositions(NULL), sort(NULL), gkKernel(NULL) {
pmePhid(NULL), pmePhip(NULL), pmePhidp(NULL), pmeCphi(NULL), pmeAtomGridIndex(NULL), lastPositions(NULL), sort(NULL), gkKernel(NULL) {
}
CudaCalcAmoebaMultipoleForceKernel::~CudaCalcAmoebaMultipoleForceKernel() {
......@@ -816,6 +816,10 @@ CudaCalcAmoebaMultipoleForceKernel::~CudaCalcAmoebaMultipoleForceKernel() {
delete labFrameDipoles;
if (labFrameQuadrupoles != NULL)
delete labFrameQuadrupoles;
if (fracDipoles != NULL)
delete fracDipoles;
if (fracQuadrupoles != NULL)
delete fracQuadrupoles;
if (field != NULL)
delete field;
if (fieldPolar != NULL)
......@@ -872,6 +876,8 @@ CudaCalcAmoebaMultipoleForceKernel::~CudaCalcAmoebaMultipoleForceKernel() {
delete pmePhip;
if (pmePhidp != NULL)
delete pmePhidp;
if (pmeCphi != NULL)
delete pmeCphi;
if (pmeAtomGridIndex != NULL)
delete pmeAtomGridIndex;
if (lastPositions != NULL)
......@@ -968,7 +974,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
int elementSize = (cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
labFrameDipoles = new CudaArray(cu, 3*paddedNumAtoms, elementSize, "labFrameDipoles");
labFrameQuadrupoles = new CudaArray(cu, 9*paddedNumAtoms, elementSize, "labFrameQuadrupoles");
labFrameQuadrupoles = new CudaArray(cu, 5*paddedNumAtoms, elementSize, "labFrameQuadrupoles");
fracDipoles = new CudaArray(cu, 3*paddedNumAtoms, elementSize, "fracDipoles");
fracQuadrupoles = new CudaArray(cu, 6*paddedNumAtoms, elementSize, "fracQuadrupoles");
field = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "field");
fieldPolar = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "fieldPolar");
torque = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "torque");
......@@ -1185,6 +1193,8 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
pmeDefines["DIRECT_POLARIZATION"] = "";
CUmodule module = cu.createModule(CudaKernelSources::vectorOps+CudaAmoebaKernelSources::multipolePme, pmeDefines);
pmeGridIndexKernel = cu.getKernel(module, "findAtomGridIndex");
pmeTransformMultipolesKernel = cu.getKernel(module, "transformMultipolesToFractionalCoordinates");
pmeTransformPotentialKernel = cu.getKernel(module, "transformPotentialToCartesianCoordinates");
pmeSpreadFixedMultipolesKernel = cu.getKernel(module, "gridSpreadFixedMultipoles");
pmeSpreadInducedDipolesKernel = cu.getKernel(module, "gridSpreadInducedDipoles");
pmeFinishSpreadChargeKernel = cu.getKernel(module, "finishSpreadCharge");
......@@ -1212,6 +1222,7 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
pmePhid = new CudaArray(cu, 10*numMultipoles, elementSize, "pmePhid");
pmePhip = new CudaArray(cu, 10*numMultipoles, elementSize, "pmePhip");
pmePhidp = new CudaArray(cu, 20*numMultipoles, elementSize, "pmePhidp");
pmeCphi = new CudaArray(cu, 10*numMultipoles, elementSize, "pmeCphi");
pmeAtomRange = CudaArray::create<int>(cu, gridSizeX*gridSizeY*gridSizeZ+1, "pmeAtomRange");
pmeAtomGridIndex = CudaArray::create<int2>(cu, numMultipoles, "pmeAtomGridIndex");
sort = new CudaSort(cu, new SortTrait(), cu.getNumAtoms());
......@@ -1483,15 +1494,46 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
gkKernel->finishComputation(*torque, *labFrameDipoles, *labFrameQuadrupoles, *inducedDipole, *inducedDipolePolar, *dampingAndThole, *covalentFlags, *polarizationGroupFlags);
}
else {
// Compute reciprocal box vectors.
Vec3 boxVectors[3];
cu.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
double determinant = boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2];
double scale = 1.0/determinant;
double3 recipBoxVectors[3];
recipBoxVectors[0] = make_double3(boxVectors[1][1]*boxVectors[2][2]*scale, 0, 0);
recipBoxVectors[1] = make_double3(-boxVectors[1][0]*boxVectors[2][2]*scale, boxVectors[0][0]*boxVectors[2][2]*scale, 0);
recipBoxVectors[2] = make_double3((boxVectors[1][0]*boxVectors[2][1]-boxVectors[1][1]*boxVectors[2][0])*scale, -boxVectors[0][0]*boxVectors[2][1]*scale, boxVectors[0][0]*boxVectors[1][1]*scale);
float3 recipBoxVectorsFloat[3];
void* recipBoxVectorPointer[3];
if (cu.getUseDoublePrecision()) {
recipBoxVectorPointer[0] = &recipBoxVectors[0];
recipBoxVectorPointer[1] = &recipBoxVectors[1];
recipBoxVectorPointer[2] = &recipBoxVectors[2];
}
else {
recipBoxVectorsFloat[0] = make_float3((float) recipBoxVectors[0].x, 0, 0);
recipBoxVectorsFloat[1] = make_float3((float) recipBoxVectors[1].x, (float) recipBoxVectors[1].y, 0);
recipBoxVectorsFloat[2] = make_float3((float) recipBoxVectors[2].x, (float) recipBoxVectors[2].y, (float) recipBoxVectors[2].z);
recipBoxVectorPointer[0] = &recipBoxVectorsFloat[0];
recipBoxVectorPointer[1] = &recipBoxVectorsFloat[1];
recipBoxVectorPointer[2] = &recipBoxVectorsFloat[2];
}
// Reciprocal space calculation.
unsigned int maxTiles = nb.getInteractingTiles().getSize();
void* gridIndexArgs[] = {&cu.getPosq().getDevicePointer(), &pmeAtomGridIndex->getDevicePointer(),
cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer()};
cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeGridIndexKernel, gridIndexArgs, cu.getNumAtoms(), cu.ThreadBlockSize, cu.ThreadBlockSize*PmeOrder*PmeOrder*elementSize);
sort->sort(*pmeAtomGridIndex);
void* pmeSpreadFixedMultipolesArgs[] = {&cu.getPosq().getDevicePointer(), &labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(),
&pmeGrid->getDevicePointer(), &pmeAtomGridIndex->getDevicePointer(), cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer()};
void* pmeTransformMultipolesArgs[] = {&labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(),
&fracDipoles->getDevicePointer(), &fracQuadrupoles->getDevicePointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeTransformMultipolesKernel, pmeTransformMultipolesArgs, cu.getNumAtoms());
void* pmeSpreadFixedMultipolesArgs[] = {&cu.getPosq().getDevicePointer(), &fracDipoles->getDevicePointer(), &fracQuadrupoles->getDevicePointer(),
&pmeGrid->getDevicePointer(), &pmeAtomGridIndex->getDevicePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeSpreadFixedMultipolesKernel, pmeSpreadFixedMultipolesArgs, cu.getNumAtoms());
void* finishSpreadArgs[] = {&pmeGrid->getDevicePointer()};
if (cu.getUseDoublePrecision())
......@@ -1501,7 +1543,7 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
else
cufftExecC2C(fft, (float2*) pmeGrid->getDevicePointer(), (float2*) pmeGrid->getDevicePointer(), CUFFT_FORWARD);
void* pmeConvolutionArgs[] = {&pmeGrid->getDevicePointer(), &pmeBsplineModuliX->getDevicePointer(), &pmeBsplineModuliY->getDevicePointer(),
&pmeBsplineModuliZ->getDevicePointer(), cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer()};
&pmeBsplineModuliZ->getDevicePointer(), cu.getPeriodicBoxSizePointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeConvolutionKernel, pmeConvolutionArgs, cu.getNumAtoms());
if (cu.getUseDoublePrecision())
cufftExecZ2Z(fft, (double2*) pmeGrid->getDevicePointer(), (double2*) pmeGrid->getDevicePointer(), CUFFT_INVERSE);
......@@ -1509,11 +1551,15 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
cufftExecC2C(fft, (float2*) pmeGrid->getDevicePointer(), (float2*) pmeGrid->getDevicePointer(), CUFFT_INVERSE);
void* pmeFixedPotentialArgs[] = {&pmeGrid->getDevicePointer(), &pmePhi->getDevicePointer(), &field->getDevicePointer(),
&fieldPolar ->getDevicePointer(), &cu.getPosq().getDevicePointer(), &labFrameDipoles->getDevicePointer(),
cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(), &pmeAtomGridIndex->getDevicePointer()};
cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2], &pmeAtomGridIndex->getDevicePointer()};
cu.executeKernel(pmeFixedPotentialKernel, pmeFixedPotentialArgs, cu.getNumAtoms());
void* pmeTransformFixedPotentialArgs[] = {&pmePhi->getDevicePointer(), &pmeCphi->getDevicePointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeTransformPotentialKernel, pmeTransformFixedPotentialArgs, cu.getNumAtoms());
void* pmeFixedForceArgs[] = {&cu.getPosq().getDevicePointer(), &cu.getForce().getDevicePointer(), &torque->getDevicePointer(),
&cu.getEnergyBuffer().getDevicePointer(), &labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(),
&pmePhi->getDevicePointer(), cu.getInvPeriodicBoxSizePointer()};
&fracDipoles->getDevicePointer(), &fracQuadrupoles->getDevicePointer(), &pmePhi->getDevicePointer(), &pmeCphi->getDevicePointer(),
recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeFixedForceKernel, pmeFixedForceArgs, cu.getNumAtoms());
// Direct space calculation.
......@@ -1521,7 +1567,8 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
void* computeFixedFieldArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(), &cu.getPosq().getDevicePointer(),
&covalentFlags->getDevicePointer(), &polarizationGroupFlags->getDevicePointer(), &nb.getExclusionTiles().getDevicePointer(), &startTileIndex, &numTileIndices,
&nb.getInteractingTiles().getDevicePointer(), &nb.getInteractionCount().getDevicePointer(), cu.getPeriodicBoxSizePointer(),
cu.getInvPeriodicBoxSizePointer(), &maxTiles, &nb.getBlockCenters().getDevicePointer(), &nb.getInteractingAtoms().getDevicePointer(),
cu.getInvPeriodicBoxSizePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
&maxTiles, &nb.getBlockCenters().getDevicePointer(), &nb.getInteractingAtoms().getDevicePointer(),
&labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(), &dampingAndThole->getDevicePointer()};
cu.executeKernel(computeFixedFieldKernel, computeFixedFieldArgs, numForceThreadBlocks*fixedFieldThreads, fixedFieldThreads);
void* recordInducedDipolesArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(),
......@@ -1532,7 +1579,8 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
cu.clearBuffer(*pmeGrid);
void* pmeSpreadInducedDipolesArgs[] = {&cu.getPosq().getDevicePointer(), &inducedDipole->getDevicePointer(), &inducedDipolePolar->getDevicePointer(),
&pmeGrid->getDevicePointer(), &pmeAtomGridIndex->getDevicePointer(), cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer()};
&pmeGrid->getDevicePointer(), &pmeAtomGridIndex->getDevicePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeSpreadInducedDipolesKernel, pmeSpreadInducedDipolesArgs, cu.getNumAtoms());
if (cu.getUseDoublePrecision())
cu.executeKernel(pmeFinishSpreadChargeKernel, finishSpreadArgs, pmeGrid->getSize());
......@@ -1546,7 +1594,8 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
else
cufftExecC2C(fft, (float2*) pmeGrid->getDevicePointer(), (float2*) pmeGrid->getDevicePointer(), CUFFT_INVERSE);
void* pmeInducedPotentialArgs[] = {&pmeGrid->getDevicePointer(), &pmePhid->getDevicePointer(), &pmePhip->getDevicePointer(),
&pmePhidp->getDevicePointer(), &cu.getPosq().getDevicePointer(), cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(),
&pmePhidp->getDevicePointer(), &cu.getPosq().getDevicePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(),
cu.getPeriodicBoxVecZPointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2],
&pmeAtomGridIndex->getDevicePointer()};
cu.executeKernel(pmeInducedPotentialKernel, pmeInducedPotentialArgs, cu.getNumAtoms());
......@@ -1559,7 +1608,8 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
void* computeInducedFieldArgs[] = {&inducedField->getDevicePointer(), &inducedFieldPolar->getDevicePointer(), &cu.getPosq().getDevicePointer(),
&nb.getExclusionTiles().getDevicePointer(), &inducedDipole->getDevicePointer(), &inducedDipolePolar->getDevicePointer(), &startTileIndex, &numTileIndices,
&nb.getInteractingTiles().getDevicePointer(), &nb.getInteractionCount().getDevicePointer(), cu.getPeriodicBoxSizePointer(),
cu.getInvPeriodicBoxSizePointer(), &maxTiles, &nb.getBlockCenters().getDevicePointer(), &nb.getInteractingAtoms().getDevicePointer(),
cu.getInvPeriodicBoxSizePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
&maxTiles, &nb.getBlockCenters().getDevicePointer(), &nb.getInteractingAtoms().getDevicePointer(),
&dampingAndThole->getDevicePointer()};
cu.executeKernel(computeInducedFieldKernel, computeInducedFieldArgs, numForceThreadBlocks*inducedFieldThreads, inducedFieldThreads);
cu.clearBuffer(*pmeGrid);
......@@ -1577,7 +1627,7 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
cufftExecC2C(fft, (float2*) pmeGrid->getDevicePointer(), (float2*) pmeGrid->getDevicePointer(), CUFFT_INVERSE);
cu.executeKernel(pmeInducedPotentialKernel, pmeInducedPotentialArgs, cu.getNumAtoms());
void* pmeRecordInducedFieldDipolesArgs[] = {&pmePhid->getDevicePointer(), &pmePhip->getDevicePointer(),
&inducedField->getDevicePointer(), &inducedFieldPolar->getDevicePointer(), cu.getInvPeriodicBoxSizePointer()};
&inducedField->getDevicePointer(), &inducedFieldPolar->getDevicePointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeRecordInducedFieldDipolesKernel, pmeRecordInducedFieldDipolesArgs, cu.getNumAtoms());
bool converged = iterateDipolesByDIIS(i);
if (converged)
......@@ -1590,14 +1640,18 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
&cu.getPosq().getDevicePointer(), &covalentFlags->getDevicePointer(), &polarizationGroupFlags->getDevicePointer(),
&nb.getExclusionTiles().getDevicePointer(), &startTileIndex, &numTileIndices,
&nb.getInteractingTiles().getDevicePointer(), &nb.getInteractionCount().getDevicePointer(),
cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(), &maxTiles, &nb.getBlockCenters().getDevicePointer(), &nb.getInteractingAtoms().getDevicePointer(),
cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
&maxTiles, &nb.getBlockCenters().getDevicePointer(), &nb.getInteractingAtoms().getDevicePointer(),
&labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(), &inducedDipole->getDevicePointer(),
&inducedDipolePolar->getDevicePointer(), &dampingAndThole->getDevicePointer()};
cu.executeKernel(electrostaticsKernel, electrostaticsArgs, numForceThreadBlocks*electrostaticsThreads, electrostaticsThreads);
void* pmeTransformInducedPotentialArgs[] = {&pmePhidp->getDevicePointer(), &pmeCphi->getDevicePointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeTransformPotentialKernel, pmeTransformInducedPotentialArgs, cu.getNumAtoms());
void* pmeInducedForceArgs[] = {&cu.getPosq().getDevicePointer(), &cu.getForce().getDevicePointer(), &torque->getDevicePointer(),
&cu.getEnergyBuffer().getDevicePointer(), &labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(),
&fracDipoles->getDevicePointer(), &fracQuadrupoles->getDevicePointer(),
&inducedDipole->getDevicePointer(), &inducedDipolePolar->getDevicePointer(), &pmePhi->getDevicePointer(), &pmePhid->getDevicePointer(),
&pmePhip->getDevicePointer(), &pmePhidp->getDevicePointer(), cu.getInvPeriodicBoxSizePointer()};
&pmePhip->getDevicePointer(), &pmePhidp->getDevicePointer(), &pmeCphi->getDevicePointer(), recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
cu.executeKernel(pmeInducedForceKernel, pmeInducedForceArgs, cu.getNumAtoms());
}
......@@ -1785,7 +1839,8 @@ void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl&
void* computePotentialArgs[] = {&cu.getPosq().getDevicePointer(), &labFrameDipoles->getDevicePointer(),
&labFrameQuadrupoles->getDevicePointer(), &inducedDipole->getDevicePointer(), &points.getDevicePointer(),
&potential.getDevicePointer(), &numPoints, cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer()};
&potential.getDevicePointer(), &numPoints, cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(),
cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer()};
int blockSize = 128;
cu.executeKernel(computePotentialKernel, computePotentialArgs, numPoints, blockSize, blockSize*15*elementSize);
outputElectrostaticPotential.resize(numPoints);
......
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