Commit 8d956093 authored by Lee-Ping's avatar Lee-Ping
Browse files

Replaced some tabs with spaces

parent 665de794
......@@ -53,7 +53,7 @@ Context::Context(const System& system, Integrator& integrator, Platform& platfor
}
Context::~Context() {
delete impl;
delete impl;
}
const System& Context::getSystem() const {
......
......@@ -64,13 +64,13 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
if (++step < owner.getFrequency() || owner.getFrequency() == 0)
return;
step = 0;
// Compute the current potential energy.
double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
// Modify the periodic box size.
Vec3 box[3];
context.getPeriodicBoxVectors(box[0], box[1], box[2]);
double volume = box[0][0]*box[1][1]*box[2][2];
......@@ -79,7 +79,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
double lengthScale = std::pow(newVolume/volume, 1.0/3.0);
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale, lengthScale, lengthScale);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale, box[1]*lengthScale, box[2]*lengthScale);
// Compute the energy of the modified system.
double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
......@@ -88,7 +88,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
// Reject the step.
kernel.getAs<ApplyMonteCarloBarostatKernel>().restoreCoordinates(context);
context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
volume = newVolume;
......@@ -145,33 +145,33 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
if (owner.getScaleX() == 0 && owner.getScaleY() == 0 && owner.getScaleZ() == 0)
return;
step = 0;
// Compute the current potential energy.
double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
double pressure;
// Choose which axis to modify at random.
double rnd = genrand_real2(random)*3.0;
int axis;
while (1) {
if (rnd < 1.0 && owner.getScaleX()) {
axis = 0;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
break;
} else if (rnd < 2.0 && owner.getScaleY()) {
axis = 1;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureY())*(AVOGADRO*1e-25);
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureY())*(AVOGADRO*1e-25);
break;
} else if (owner.getScaleZ()) {
axis = 2;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureZ())*(AVOGADRO*1e-25);
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureZ())*(AVOGADRO*1e-25);
break;
}
}
// Modify the periodic box size.
Vec3 box[3];
context.getPeriodicBoxVectors(box[0], box[1], box[2]);
double volume = box[0][0]*box[1][1]*box[2][2];
......@@ -179,11 +179,11 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
double newVolume = volume+deltaVolume;
Vec3 lengthScale;
for (int i=0; i<3; i++)
lengthScale[i] = 1.0;
lengthScale[i] = 1.0;
lengthScale[axis] = newVolume/volume;
kernel.getAs<ApplyMonteCarloAnisotropicBarostatKernel>().scaleCoordinates(context, lengthScale[0], lengthScale[1], lengthScale[2]);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale[0], box[1]*lengthScale[1], box[2]*lengthScale[2]);
// Compute the energy of the modified system.
double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
......@@ -191,7 +191,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
// Reject the step.
kernel.getAs<ApplyMonteCarloAnisotropicBarostatKernel>().restoreCoordinates(context);
context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
volume = newVolume;
......
......@@ -98,9 +98,9 @@ void testIdealGas(int aniso) {
const double temp[] = {300.0, 600.0, 1000.0};
const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
vector<Vec3> positions(numParticles);
......@@ -112,32 +112,32 @@ void testIdealGas(int aniso) {
}
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
if (aniso)
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
system.addForce(barostat);
// Test it for three different temperatures.
for (int i = 0; i < 3; i++) {
barostat->setTemperature(temp[i]);
LangevinIntegrator integrator(temp[i], 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];
if (!aniso) {
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
}
if (!aniso) {
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
}
integrator.step(frequency);
}
volume /= steps;
......@@ -168,9 +168,9 @@ void testRandomSeed() {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
barostat->setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
......@@ -182,9 +182,9 @@ void testRandomSeed() {
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
barostat->setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
......@@ -196,9 +196,9 @@ void testRandomSeed() {
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
......@@ -219,9 +219,9 @@ void testWater(int aniso) {
const double angle = 109.47*M_PI/180;
const double dOH = 0.1;
const double dHH = dOH*2*std::sin(0.5*angle);
// Create a box of SPC water molecules.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
NonbondedForce* nonbonded = new NonbondedForce();
......@@ -256,11 +256,11 @@ void testWater(int aniso) {
system.addForce(nonbonded);
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
if (aniso)
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, frequency);
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, frequency);
system.addForce(barostat);
// Simulate it and see if the density matches the expected value (1 g/mL).
LangevinIntegrator integrator(temp, 1.0, 0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
......
......@@ -22,7 +22,7 @@ __kernel void scalePositions(float scaleX, float scaleY, float scaleZ, int numMo
int ycell = (int) floor(center.y*invPeriodicBoxSize.y);
int zcell = (int) floor(center.z*invPeriodicBoxSize.z);
real4 delta = (real4) (xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z, 0);
real4 scaleXYZ = (real4) (scaleX, scaleY, scaleZ, 1);
real4 scaleXYZ = (real4) (scaleX, scaleY, scaleZ, 1);
center -= delta;
// Now scale the position of the molecule center.
......
......@@ -97,9 +97,9 @@ void testIdealGas(int aniso) {
const double temp[] = {300.0, 600.0, 1000.0};
const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
ReferencePlatform platform;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
......@@ -112,32 +112,32 @@ void testIdealGas(int aniso) {
}
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
if (aniso)
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
system.addForce(barostat);
// Test it for three different temperatures.
for (int i = 0; i < 3; i++) {
barostat->setTemperature(temp[i]);
LangevinIntegrator integrator(temp[i], 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];
if (!aniso) {
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
}
if (!aniso) {
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
}
integrator.step(frequency);
}
volume /= steps;
......@@ -169,9 +169,9 @@ void testRandomSeed() {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
barostat->setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
......@@ -183,9 +183,9 @@ void testRandomSeed() {
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
barostat->setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
......@@ -197,9 +197,9 @@ void testRandomSeed() {
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
......
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