diff --git a/Configuration.cpp b/Configuration.cpp
index a52b56b0c9d5d8f82546544e812c12b3a3da88cb..2ed7f53aab9bea28d382e4e740815979184e1f38 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -113,28 +113,24 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
if (readExcludesFromFile) readExcludes();
if (readAnglesFromFile) readAngles();
if (readDihedralsFromFile) readDihedrals();
- /*
- if (kTGridFile.length() != 0) {
- printf("\nFound kT grid file: %s\n", kTGridFile.val());
- kTGrid = new BaseGrid(kTGridFile.val());
- printf("Loaded `%s'.\n", kTGridFile.val());
- printf("System size %s.\n", kTGrid->getExtent().toString().val());
- }
- */
- if (temperatureGrid.length() != 0) {
- printf("\nFound temperature grid file: %s\n", temperatureGrid.val());
- tGrid = new BaseGrid(temperatureGrid.val());
- printf("Loaded `%s'.\n", temperatureGrid.val());
+
+ kT = temperature * 0.0019872065f; // `units "k K" "kcal_mol"`
+ if (temperatureGridFile.length() != 0) {
+ printf("\nFound temperature grid file: %s\n", temperatureGridFile.val());
+ tGrid = new BaseGrid(temperatureGridFile.val());
+ printf("Loaded `%s'.\n", temperatureGridFile.val());
printf("System size %s.\n", tGrid->getExtent().toString().val());
+ // TODO: ask Max Belkin what this is about and how to remove hard-coded temps
float ToSo = 1.0f / (295.0f * 4.634248239f); // 1 / (To * sigma(To))
- sigmaT = new BaseGrid(temperatureGrid.val());
+ sigmaT = new BaseGrid(*tGrid);
sigmaT->shift(-122.8305f);
sigmaT->scale(0.0269167f);
sigmaT->mult(*tGrid);
sigmaT->scale(ToSo);
- kTGrid = new BaseGrid(temperatureGrid.val());
- float factor = kT/295.0f;
+
+ kTGrid = new BaseGrid(*tGrid);
+ float factor = 0.0019872065f; // `units "k K" "kcal_mol"`
kTGrid->scale(factor);
// char outFile[256];
// char comment[256]; sprintf(comment,"KTGrid");
@@ -195,7 +191,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
printf("System size %s.\n", part[i].diffusionGrid->getExtent().toString().val());
}
- if (temperatureGrid.length() != 0) {
+ if (temperatureGridFile.length() != 0) {
if (partDiffusionGridFile[i].length() != 0) {
part[i].diffusionGrid->mult(*sigmaT);
} else {
@@ -479,7 +475,7 @@ void Configuration::copyToCUDA() {
// kTGrid_d
kTGrid_d = NULL;
- if (temperatureGrid.length() > 0) {
+ if (temperatureGridFile.length() > 0) {
gpuErrchk(cudaMalloc(&kTGrid_d, sizeof(BaseGrid)));
gpuErrchk(cudaMemcpyAsync(kTGrid_d, kTGrid, sizeof(BaseGrid), cudaMemcpyHostToDevice));
}
@@ -530,7 +526,7 @@ void Configuration::copyToCUDA() {
}
void Configuration::setDefaults() {
- // System parameters
+ // System parameters
outputName = "out";
timestep = 1e-5f;
steps = 100;
@@ -542,9 +538,9 @@ void Configuration::setDefaults() {
interparticleForce = 1;
tabulatedPotential = 0;
fullLongRange = 1;
- kT = 1.0f;
// kTGridFile = ""; // Commented out for an unknown reason
- temperatureGrid = "";
+ temperature = 295.0f;
+ temperatureGridFile = "";
coulombConst = 566.440698f/92.0f;
electricField = 0.0f;
cutoff = 10.0f;
@@ -646,15 +642,10 @@ int Configuration::readParameters(const char * config_file) {
inputCoordinates = value;
else if (param == String("restartCoordinates"))
restartCoordinates = value;
- else if (param == String("kT"))
- kT = (float) strtod(value.val(), NULL);
- // else if (param == String("kTGridFile")) kTGridFile = value;
- else if (param == String("temperature")) {
+ else if (param == String("temperature"))
temperature = (float) strtod(value.val(),NULL);
- kT = 0.58622522f * temperature / 295.0f; // kcal/mol (units "295 k K" "kcal_mol")
- }
else if (param == String("temperatureGrid"))
- temperatureGrid = value;
+ temperatureGridFile = value;
else if (param == String("numberFluct"))
numberFluct = atoi(value.val());
else if (param == String("numberFluctPeriod"))
diff --git a/Configuration.h b/Configuration.h
index 7fc75f180a618353d39d45ca7c25381ab9cc2fe5..450cc8b9a52cb9645f3a51febd8376bc2a8760b6 100644
--- a/Configuration.h
+++ b/Configuration.h
@@ -33,6 +33,13 @@
#include <cuda.h>
#include <cuda_runtime.h>
+// Units:
+// Energy: kcal/mol (6.947694e-24 kJ)
+// Temperature: Kelvin
+// Time: nanoseconds
+// Length: nanometers
+
+
class Configuration {
struct compare {
bool operator()(const String& lhs, const String& rhs);
@@ -120,7 +127,7 @@ public:
long int steps;
long int seed;
// String kTGridFile;
- String temperatureGrid;
+ String temperatureGridFile;
String inputCoordinates;
String restartCoordinates;
int numberFluct;
diff --git a/GrandBrownTown.cu b/GrandBrownTown.cu
index bc4788cd471cd789ca3e5b420c3152013b89635f..ba82eddecb4b99d6c048e1e31ecae00f51f7b773 100644
--- a/GrandBrownTown.cu
+++ b/GrandBrownTown.cu
@@ -77,7 +77,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
timestep = c.timestep;
steps = c.steps;
seed = c.seed;
- temperatureGrid = c.temperatureGrid;
+ temperatureGridFile = c.temperatureGridFile;
inputCoordinates = c.inputCoordinates;
restartCoordinates = c.restartCoordinates;
numberFluct = c.numberFluct;
@@ -86,6 +86,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
readBondsFromFile = c.readBondsFromFile;
fullLongRange = c.fullLongRange;
kT = c.kT;
+ temperature = c.temperature;
coulombConst = c.coulombConst;
electricField = c.electricField;
cutoff = c.cutoff;
@@ -449,7 +450,7 @@ void GrandBrownTown::run() {
// int numBlocks = (num * numReplicas) / NUM_THREADS + 1;
int numBlocks = (num * numReplicas) / NUM_THREADS + (num * numReplicas % NUM_THREADS == 0 ? 0 : 1);
- int tl = temperatureGrid.length();
+ int tl = temperatureGridFile.length();
RBC.updateForces(internal->getPos_d(), internal->getForceInternal_d(), s); /* update RB forces before update particle positions... */
diff --git a/GrandBrownTown.cuh b/GrandBrownTown.cuh
index 157978d62bdc0e9f9e1f2c4bd7b669c0a472a997..9ceb183132ed8cf686d770959439c31a99bf0722 100644
--- a/GrandBrownTown.cuh
+++ b/GrandBrownTown.cuh
@@ -32,7 +32,6 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
// Calculate this thread's ID
const int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
// TODO: Make this a grid-stride loop to make efficient reuse of RNG states
// Loop over ALL particles in ALL replicas
if (idx < num * numReplicas) {
@@ -53,7 +52,7 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
// Add a force defined via 3D FORCE maps (not 3D potential maps)
if (pt.forceXGrid != NULL) fe.f.x += pt.forceXGrid->interpolatePotentialLinearly(p);
if (pt.forceYGrid != NULL) fe.f.y += pt.forceYGrid->interpolatePotentialLinearly(p);
- if (pt.forceZGrid != NULL) fe.f.z += pt.forceZGrid->interpolatePotentialLinearly(p);
+ if (pt.forceZGrid != NULL) fe.f.z += pt.forceZGrid->interpolatePotentialLinearly(p);
#endif
// Compute total force:
@@ -76,18 +75,16 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
// printf("atom %d: pos: %f %f %f\n", idx, p.x, p.y, p.z);
// p = pt.diffusionGrid->wrap(p); // illegal mem access; no origin/basis?
- Vector3 gridCenter = pt.diffusionGrid->origin +
- pt.diffusionGrid->basis.transform( Vector3(0.5*pt.diffusionGrid->nx,
- 0.5*pt.diffusionGrid->ny,
- 0.5*pt.diffusionGrid->nz));
- Vector3 p2 = p - gridCenter;
- p2 = sys->wrapDiff( p2 ) + gridCenter;
-
- /* p2 = sys->wrap( p2 ); */
- /* p2 = p2 - gridCenter; */
- /* printf("atom %d: ps2: %f %f %f\n", idx, p2.x, p2.y, p2.z); */
-
- ForceEnergy diff = pt.diffusionGrid->interpolateForceDLinearlyPeriodic(p2);
+ Vector3 gridCenter = pt.diffusionGrid->origin +
+ pt.diffusionGrid->basis.transform( Vector3(0.5*pt.diffusionGrid->nx,
+ 0.5*pt.diffusionGrid->ny,
+ 0.5*pt.diffusionGrid->nz));
+ Vector3 p2 = p - gridCenter;
+ p2 = sys->wrapDiff( p2 ) + gridCenter;
+ /* p2 = sys->wrap( p2 ); */
+ /* p2 = p2 - gridCenter; */
+ /* printf("atom %d: ps2: %f %f %f\n", idx, p2.x, p2.y, p2.z); */
+ ForceEnergy diff = pt.diffusionGrid->interpolateForceDLinearlyPeriodic(p2);
diffusion = diff.e;
diffGrad = diff.f;
}
diff --git a/GrandBrownTown.h b/GrandBrownTown.h
index 44a783ee5b1d3782f1828d2af9a841e571d78e0b..006ae8ad1dfde892ec3c1580ff0c5a555a8ab993 100644
--- a/GrandBrownTown.h
+++ b/GrandBrownTown.h
@@ -171,8 +171,7 @@ private:
float timestep;
long int steps;
long int seed;
- // String kTGridFile;
- String temperatureGrid;
+ String temperatureGridFile;
String inputCoordinates;
String restartCoordinates;
int numberFluct;
@@ -180,6 +179,7 @@ private:
int tabulatedPotential;
int fullLongRange;
float kT;
+ float temperature;
float coulombConst;
float electricField;
float cutoff;
diff --git a/RigidBody.cu b/RigidBody.cu
index bd614d16fb1c6c05f253232e46039682e570e11c..67cea83a5f50ca0074957106d3e244bc7f1f7828 100644
--- a/RigidBody.cu
+++ b/RigidBody.cu
@@ -5,14 +5,13 @@
#include "Debug.h"
-// RBTODO: make units same as those in atomic simulation
RigidBody::RigidBody(String name, const Configuration& cref, RigidBodyType& tref)
: name(name), c(&cref), t(&tref), impulse_to_momentum(4.184e8f) {
// units "(kcal_mol/AA) * ns" "dalton AA/ns" * 4.184e+08
timestep = c->timestep;
- // RBTODO: fix this
- Temp = 295;
+ Temp = c->temperature;
+ // RBTODO: use temperature grids
// tempgrid = c->temperatureGrid;
position = t->initPos;