Commit 62c4fd53 authored by Peter Eastman's avatar Peter Eastman
Browse files

Renamed "atoms" to "particles"

parent ad75a390
......@@ -53,8 +53,8 @@ const double TOL = 1e-5;
void testSingleBond() {
ReferencePlatform platform;
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
system.setParticleMass(0, 2.0);
system.setParticleMass(1, 2.0);
double dt = 0.01;
BrownianIntegrator integrator(0, 0.1, dt);
HarmonicBondForce* forceField = new HarmonicBondForce(1);
......@@ -85,23 +85,23 @@ void testSingleBond() {
}
void testTemperature() {
const int numAtoms = 8;
const int numBonds = numAtoms-1;
const int numParticles = 8;
const int numBonds = numParticles-1;
const double temp = 100.0;
ReferencePlatform platform;
System system(numAtoms, 0);
System system(numParticles, 0);
BrownianIntegrator integrator(temp, 2.0, 0.01);
HarmonicBondForce* forceField = new HarmonicBondForce(numBonds);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
// forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0);
// forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
for (int i = 0; i < numBonds; ++i)
forceField->setBondParameters(i, i, i+1, 1.0, i);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
for (int i = 0; i < numAtoms; ++i)
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3(i, 0, 0);
context.setPositions(positions);
......@@ -123,24 +123,24 @@ void testTemperature() {
}
void testConstraints() {
const int numAtoms = 8;
const int numParticles = 8;
const double temp = 100.0;
ReferencePlatform platform;
System system(numAtoms, numAtoms-1);
System system(numParticles, numParticles-1);
BrownianIntegrator integrator(temp, 2.0, 0.001);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 10.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numAtoms-1; ++i)
for (int i = 0; i < numParticles-1; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
......@@ -151,7 +151,7 @@ void testConstraints() {
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numAtoms-1; ++j) {
for (int j = 0; j < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
......
......@@ -51,37 +51,37 @@ using namespace std;
Vec3 calcCM(const vector<Vec3>& values, System& system) {
Vec3 cm;
for (int j = 0; j < system.getNumAtoms(); ++j) {
cm[0] += values[j][0]*system.getAtomMass(j);
cm[1] += values[j][1]*system.getAtomMass(j);
cm[2] += values[j][2]*system.getAtomMass(j);
for (int j = 0; j < system.getNumParticles(); ++j) {
cm[0] += values[j][0]*system.getParticleMass(j);
cm[1] += values[j][1]*system.getParticleMass(j);
cm[2] += values[j][2]*system.getParticleMass(j);
}
return cm;
}
void testMotionRemoval() {
const int numAtoms = 8;
const int numParticles = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
ReferencePlatform platform;
System system(numAtoms, 0);
System system(numParticles, 0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* bonds = new HarmonicBondForce(1);
bonds->setBondParameters(0, 2, 3, 2.0, 0.5);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(numAtoms, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, i+1);
nonbonded->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
NonbondedForce* nonbonded = new NonbondedForce(numParticles, 0);
for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, i+1);
nonbonded->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(nonbonded);
CMMotionRemover* remover = new CMMotionRemover();
system.addForce(remover);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
......
......@@ -49,13 +49,13 @@ using namespace std;
const double TOL = 1e-5;
void testSingleAtom() {
void testSingleParticle() {
ReferencePlatform platform;
System system(1, 0);
system.setAtomMass(0, 2.0);
system.setParticleMass(0, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(1);
forceField->setAtomParameters(0, 0.5, 0.15, 1);
forceField->setParticleParameters(0, 0.5, 0.15, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(1);
......@@ -72,20 +72,20 @@ void testSingleAtom() {
void testForce() {
ReferencePlatform platform;
const int numAtoms = 10;
System system(numAtoms, 0);
const int numParticles = 10;
System system(numParticles, 0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(numAtoms);
for (int i = 0; i < numAtoms; ++i)
forceField->setAtomParameters(i, i%2 == 0 ? -1 : 1, 0.15, 1);
GBSAOBCForceField* forceField = new GBSAOBCForceField(numParticles);
for (int i = 0; i < numParticles; ++i)
forceField->setParticleParameters(i, i%2 == 0 ? -1 : 1, 0.15, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
// Set random positions for all the atoms.
// Set random positions for all the particles.
vector<Vec3> positions(numAtoms);
vector<Vec3> positions(numParticles);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i)
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -93,14 +93,14 @@ void testForce() {
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numAtoms; ++i) {
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numAtoms; ++i) {
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
......@@ -115,7 +115,7 @@ void testForce() {
int main() {
try {
testSingleAtom();
testSingleParticle();
testForce();
}
catch(const exception& e) {
......
......@@ -50,7 +50,7 @@ void testCalcKE() {
ReferencePlatform platform;
System system(4, 0);
for (int i = 0; i < 4; ++i)
system.setAtomMass(i, i+1);
system.setParticleMass(i, i+1);
VerletIntegrator integrator(0.01);
OpenMMContext context(system, integrator, platform);
vector<Vec3> velocities(4);
......
......@@ -53,8 +53,8 @@ const double TOL = 1e-5;
void testSingleBond() {
ReferencePlatform platform;
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
system.setParticleMass(0, 2.0);
system.setParticleMass(1, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
HarmonicBondForce* forceField = new HarmonicBondForce(1);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
......@@ -95,20 +95,20 @@ void testSingleBond() {
}
void testTemperature() {
const int numAtoms = 8;
const int numParticles = 8;
const double temp = 100.0;
ReferencePlatform platform;
System system(numAtoms, 0);
System system(numParticles, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
for (int i = 0; i < numAtoms; ++i)
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
......@@ -125,29 +125,29 @@ void testTemperature() {
integrator.step(1);
}
ke /= 1000;
double expected = 0.5*numAtoms*3*BOLTZ*temp;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
}
void testConstraints() {
const int numAtoms = 8;
const int numParticles = 8;
const double temp = 100.0;
ReferencePlatform platform;
System system(numAtoms, numAtoms-1);
System system(numParticles, numParticles-1);
LangevinIntegrator integrator(temp, 2.0, 0.01);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 10.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numAtoms-1; ++i)
for (int i = 0; i < numParticles-1; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
......@@ -158,7 +158,7 @@ void testConstraints() {
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numAtoms-1; ++j) {
for (int j = 0; j < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
......
......@@ -9,35 +9,35 @@ using namespace OpenMM;
void testNeighborList()
{
RealOpenMM* atomList[2];
atomList[0] = new RealOpenMM[3];
atomList[1] = new RealOpenMM[3];
atomList[2] = new RealOpenMM[3];
atomList[0][0] = 13.6f;
atomList[0][1] = 0;
atomList[0][2] = 0;
atomList[1][0] = 0;
atomList[1][1] = 0;
atomList[1][2] = 0;
RealOpenMM* particleList[2];
particleList[0] = new RealOpenMM[3];
particleList[1] = new RealOpenMM[3];
particleList[2] = new RealOpenMM[3];
particleList[0][0] = 13.6f;
particleList[0][1] = 0;
particleList[0][2] = 0;
particleList[1][0] = 0;
particleList[1][1] = 0;
particleList[1][2] = 0;
vector<set<int> > exclusions(2);
NeighborList neighborList;
computeNeighborListNaive(neighborList, 2, atomList, exclusions, NULL, 13.7, 0.01);
computeNeighborListNaive(neighborList, 2, particleList, exclusions, NULL, 13.7, 0.01);
assert(neighborList.size() == 1);
computeNeighborListNaive(neighborList, 2, atomList, exclusions, NULL, 13.5, 0.01);
computeNeighborListNaive(neighborList, 2, particleList, exclusions, NULL, 13.5, 0.01);
assert(neighborList.size() == 0);
computeNeighborListVoxelHash(neighborList, 2, atomList, exclusions, NULL, 13.7, 0.01);
computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, NULL, 13.7, 0.01);
assert(neighborList.size() == 1);
computeNeighborListVoxelHash(neighborList, 2, atomList, exclusions, NULL, 13.5, 0.01);
computeNeighborListVoxelHash(neighborList, 2, particleList, exclusions, NULL, 13.5, 0.01);
assert(neighborList.size() == 0);
delete[] atomList[0];
delete[] atomList[1];
delete[] atomList[2];
delete[] particleList[0];
delete[] particleList[1];
delete[] particleList[2];
}
double periodicDifference(double val1, double val2, double period) {
......@@ -53,40 +53,40 @@ double distance2(RealOpenMM* pos1, RealOpenMM* pos2, const RealOpenMM* periodicB
return dx*dx+dy*dy+dz*dz;
}
void verifyNeighborList(NeighborList& list, int numAtoms, RealOpenMM** positions, const RealOpenMM* periodicBoxSize, double cutoff) {
void verifyNeighborList(NeighborList& list, int numParticles, RealOpenMM** positions, const RealOpenMM* periodicBoxSize, double cutoff) {
for (int i = 0; i < (int) list.size(); i++) {
int atom1 = list[i].first;
int atom2 = list[i].second;
ASSERT(distance2(positions[atom1], positions[atom2], periodicBoxSize) <= cutoff*cutoff);
int particle1 = list[i].first;
int particle2 = list[i].second;
ASSERT(distance2(positions[particle1], positions[particle2], periodicBoxSize) <= cutoff*cutoff);
}
int count = 0;
for (int i = 0; i < numAtoms; i++)
for (int j = i+1; j < numAtoms; j++)
for (int i = 0; i < numParticles; i++)
for (int j = i+1; j < numParticles; j++)
if (distance2(positions[i], positions[j], periodicBoxSize) <= cutoff*cutoff)
count++;
ASSERT(count == list.size());
}
void testPeriodic() {
const int numAtoms = 100;
const int numParticles = 100;
const double cutoff = 3.0;
const RealOpenMM periodicBoxSize[3] = {20.0, 15.0, 22.0};
RealOpenMM* atomList[numAtoms];
RealOpenMM* particleList[numParticles];
init_gen_rand(0);
for (int i = 0; i <numAtoms; i++) {
atomList[i] = new RealOpenMM[3];
atomList[i][0] = (RealOpenMM) (genrand_real2()*periodicBoxSize[0]*3);
atomList[i][1] = (RealOpenMM) (genrand_real2()*periodicBoxSize[1]*3);
atomList[i][2] = (RealOpenMM) (genrand_real2()*periodicBoxSize[2]*3);
for (int i = 0; i <numParticles; i++) {
particleList[i] = new RealOpenMM[3];
particleList[i][0] = (RealOpenMM) (genrand_real2()*periodicBoxSize[0]*3);
particleList[i][1] = (RealOpenMM) (genrand_real2()*periodicBoxSize[1]*3);
particleList[i][2] = (RealOpenMM) (genrand_real2()*periodicBoxSize[2]*3);
}
vector<set<int> > exclusions(numAtoms);
vector<set<int> > exclusions(numParticles);
NeighborList neighborList;
computeNeighborListNaive(neighborList, numAtoms, atomList, exclusions, periodicBoxSize, cutoff);
verifyNeighborList(neighborList, numAtoms, atomList, periodicBoxSize, cutoff);
computeNeighborListVoxelHash(neighborList, numAtoms, atomList, exclusions, periodicBoxSize, cutoff);
verifyNeighborList(neighborList, numAtoms, atomList, periodicBoxSize, cutoff);
for (int i = 0; i <numAtoms; i++)
delete[] atomList[i];
computeNeighborListNaive(neighborList, numParticles, particleList, exclusions, periodicBoxSize, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
computeNeighborListVoxelHash(neighborList, numParticles, particleList, exclusions, periodicBoxSize, cutoff);
verifyNeighborList(neighborList, numParticles, particleList, periodicBoxSize, cutoff);
for (int i = 0; i <numParticles; i++)
delete[] particleList[i];
}
int main()
......
......@@ -54,8 +54,8 @@ void testCoulomb() {
System system(2, 0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce(2, 0);
forceField->setAtomParameters(0, 0.5, 1, 0);
forceField->setAtomParameters(1, -1.5, 1, 0);
forceField->setParticleParameters(0, 0.5, 1, 0);
forceField->setParticleParameters(1, -1.5, 1, 0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
......@@ -75,8 +75,8 @@ void testLJ() {
System system(2, 0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce(2, 0);
forceField->setAtomParameters(0, 0, 1.2, 1);
forceField->setAtomParameters(1, 0, 1.4, 2);
forceField->setParticleParameters(0, 0, 1.2, 1);
forceField->setParticleParameters(1, 0, 1.4, 2);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(2);
......@@ -113,11 +113,11 @@ void testExclusionsAnd14() {
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
nonbonded->setAtomParameters(j, 0, 1.5, 0);
nonbonded->setParticleParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
nonbonded->setAtomParameters(0, 0, 1.5, 1);
nonbonded->setAtomParameters(i, 0, 1.5, 1);
nonbonded->setParticleParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
......@@ -143,8 +143,8 @@ void testExclusionsAnd14() {
// Test Coulomb forces
nonbonded->setAtomParameters(0, 2, 1.5, 0);
nonbonded->setAtomParameters(i, 2, 1.5, 0);
nonbonded->setParticleParameters(0, 2, 1.5, 0);
nonbonded->setParticleParameters(i, 2, 1.5, 0);
nonbonded->setNonbonded14Parameters(0, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
context.reinitialize();
......@@ -172,9 +172,9 @@ void testCutoff() {
System system(3, 0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce(3, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
forceField->setParticleParameters(0, 1.0, 1, 0);
forceField->setParticleParameters(1, 1.0, 1, 0);
forceField->setParticleParameters(2, 1.0, 1, 0);
forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
......@@ -226,10 +226,10 @@ void testCutoff14() {
// Test LJ forces
nonbonded->setAtomParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
nonbonded->setAtomParameters(j, 0, 1.5, 0);
nonbonded->setAtomParameters(i, 0, 1.5, 1);
nonbonded->setParticleParameters(j, 0, 1.5, 0);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
......@@ -256,8 +256,8 @@ void testCutoff14() {
// Test Coulomb forces
const double q = 0.7;
nonbonded->setAtomParameters(0, q, 1.5, 0);
nonbonded->setAtomParameters(i, q, 1.5, 0);
nonbonded->setParticleParameters(0, q, 1.5, 0);
nonbonded->setParticleParameters(i, q, 1.5, 0);
nonbonded->setNonbonded14Parameters(0, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
context.reinitialize();
......@@ -291,9 +291,9 @@ void testPeriodic() {
bonds->setBondParameters(0, 0, 1, 1, 0);
system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(3, 0);
nonbonded->setAtomParameters(0, 1.0, 1, 0);
nonbonded->setAtomParameters(1, 1.0, 1, 0);
nonbonded->setAtomParameters(2, 1.0, 1, 0);
nonbonded->setParticleParameters(0, 1.0, 1, 0);
nonbonded->setParticleParameters(1, 1.0, 1, 0);
nonbonded->setParticleParameters(2, 1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
......
......@@ -53,8 +53,8 @@ const double TOL = 1e-5;
void testSingleBond() {
ReferencePlatform platform;
System system(2, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
system.setParticleMass(0, 2.0);
system.setParticleMass(1, 2.0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce(1);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
......@@ -86,24 +86,24 @@ void testSingleBond() {
}
void testConstraints() {
const int numAtoms = 8;
const int numParticles = 8;
const double temp = 100.0;
ReferencePlatform platform;
System system(numAtoms, numAtoms-1);
System system(numParticles, numParticles-1);
VerletIntegrator integrator(0.002);
NonbondedForce* forceField = new NonbondedForce(numAtoms, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0);
for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 10.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numAtoms-1; ++i)
for (int i = 0; i < numParticles-1; ++i)
system.setConstraintParameters(i, i, i+1, 1.0);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(numAtoms);
vector<Vec3> velocities(numAtoms);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i) {
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
......@@ -115,7 +115,7 @@ void testConstraints() {
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy);
for (int j = 0; j < numAtoms-1; ++j) {
for (int j = 0; j < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
......
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