diff --git a/GrandBrownTown.cu b/GrandBrownTown.cu
index 4addee478a666f489a41b51568a42e4277662b5f..ca06f6a88af9afad839fb5c0ebc1e1df19cfba3b 100644
--- a/GrandBrownTown.cu
+++ b/GrandBrownTown.cu
@@ -17,7 +17,7 @@ cudaEvent_t START, STOP;
GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
const long int randomSeed, bool debug, bool imd_on, unsigned int imd_port, int numReplicas) :
imd_on(imd_on), imd_port(imd_port), numReplicas(numReplicas),
- conf(c), RBC(RigidBodyController(c)) {
+ conf(c), RBC(RigidBodyController(c,outArg)) {
for (int i = 0; i < numReplicas; i++) {
std::stringstream curr_file, restart_file, out_prefix;
@@ -399,6 +399,8 @@ void GrandBrownTown::run() {
int numBlocks = (num * numReplicas) / NUM_THREADS + (num * numReplicas % NUM_THREADS == 0 ? 0 : 1);
int tl = temperatureGrid.length();
+ RBC.updateForces(); /* update RB forces before update particle positions... */
+
// Call the kernel to update the positions of each particle
updateKernel<<< numBlocks, NUM_THREADS >>>(pos_d, forceInternal_d, type_d,
part_d, kT, kTGrid_d,
@@ -406,20 +408,18 @@ void GrandBrownTown::run() {
sys_d, randoGen_d, numReplicas);
//gpuErrchk(cudaPeekAtLastError()); // Does not work on old GPUs (like mine). TODO: write a better wrapper around Peek
-
/* Time position computations.
rt_timer_stop(cputimer);
float dt2 = rt_timer_time(cputimer);
printf("Position Update Time: %f ms\n", dt2 * 1000);
*/
- RBC.updateForces();
+ RBC.integrate(s);
/* for (int j = 0; j < t->num; j++) { */
/* computeGridGridForce<<< numBlocks, NUM_THREADS >>>(grid1_d, grid2_d); */
// int numBlocks = (numRB ) / NUM_THREADS + (num * numReplicas % NUM_THREADS == 0 ? 0 : 1);
-
Vector3 force0(0.0f);
if (imd_on && clientsock && s % outputPeriod == 0) {
diff --git a/RigidBody.cu b/RigidBody.cu
index 448bf099934478886150d3f58bcdc05b5e06e503..01739b8d7594a6efdbc05bcfbcc3d9363fa2794e 100644
--- a/RigidBody.cu
+++ b/RigidBody.cu
@@ -51,7 +51,7 @@ void RigidBody::addTorque(Force torq) {
}
RigidBody::~RigidBody() {}
-/*===========================================================================\
+ /*===========================================================================\
| Following "Algorithm for rigid-body Brownian dynamics" Dan Gordon, Matthew |
| Hoyles, and Shin-Ho Chung |
| http://langevin.anu.edu.au/publications/PhysRevE_80_066703.pdf |
@@ -82,12 +82,13 @@ void RigidBody::addLangevin(Vector3 w1, Vector3 w2) {
addTorque(torq);
}
-/*==========================================================================\
+ /*==========================================================================\
| from: Dullweber, Leimkuhler, Maclachlan. Symplectic splitting methods for |
| rigid body molecular dynamics. JCP 107. (1997) |
| http://jcp.aip.org/resource/1/jcpsa6/v107/i15/p5840_s1 |
\==========================================================================*/
-void RigidBody::integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll) {
+// void RigidBody::integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll) {}
+void RigidBody::integrate(int startFinishAll) {
Vector3 trans; // = *p_trans;
Matrix3 rot = Matrix3(1); // = *p_rot;
@@ -169,18 +170,18 @@ void RigidBody::integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishA
// update actual orientation
Matrix3 newOrientation = orientation*rot; // not 100% sure; rot could be in rb frame
orientation = newOrientation;
- rot = rot.transpose();
+ /* rot = rot.transpose(); */
- DebugM(2, "trans during: " << trans
- << "\n" << endi);
- DebugM(2, "rot during: " << rot
- << "\n" << endi);
+ /* DebugM(2, "trans during: " << trans */
+ /* << "\n" << endi); */
+ /* DebugM(2, "rot during: " << rot */
+ /* << "\n" << endi); */
- clearForce();
- clearTorque();
+ /* clearForce(); */
+ /* clearTorque(); */
- old_trans = trans;
- old_rot = rot;
+ /* old_trans = trans; */
+ /* old_rot = rot; */
}
DebugM(3, " position after: " << position << "\n" << endi);
}
diff --git a/RigidBody.h b/RigidBody.h
index 6e8af6d2a2206c0ddec2186ed49f2bc67122a13d..cebbf0f80b5070ff99d2aedc9270bca03fd02f0f 100644
--- a/RigidBody.h
+++ b/RigidBody.h
@@ -38,9 +38,11 @@ class RigidBody { // host side representation of rigid bodies
HOST DEVICE inline void clearForce() { force = Force(0.0f); }
HOST DEVICE inline void clearTorque() { torque = Force(0.0f); }
- HOST DEVICE void integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll);
+ // HOST DEVICE void integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll);
+ // HOST DEVICE void integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll);
+ HOST DEVICE void integrate(int startFinishAll);
- HOST DEVICE inline String getKey() { return key; }
+ HOST DEVICE inline String getKey() const { return key; }
HOST DEVICE inline Vector3 getPosition() const { return position; }
HOST DEVICE inline Matrix3 getOrientation() const { return orientation; }
diff --git a/RigidBodyController.cu b/RigidBodyController.cu
index 9bbae16e8e3cced45e8593a409760bf1853ba1e6..1fe941d93c322e2721517b01b40f2d54351d72e0 100644
--- a/RigidBodyController.cu
+++ b/RigidBodyController.cu
@@ -11,8 +11,9 @@
#include "ComputeGridGrid.cuh"
// #include <vector>
+#include "Debug.h"
-/* #include "Random.h" */
+#include "Random.h" /* RBTODO: fix this? */
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, String file, int line, bool abort=true) {
@@ -26,8 +27,8 @@ inline void gpuAssert(cudaError_t code, String file, int line, bool abort=true)
/* #include <cuda_runtime.h> */
/* #include <curand_kernel.h> */
-RigidBodyController::RigidBodyController(const Configuration& c) :
-conf(c) {
+RigidBodyController::RigidBodyController(const Configuration& c, const char* outArg) :
+ conf(c), outArg(outArg) {
if (conf.numRigidTypes > 0) {
copyGridsToDevice();
@@ -45,12 +46,15 @@ conf(c) {
rigidBodyByType.push_back(tmp);
}
+ random = new Random(conf.seed + 1); /* +1 to avoid using same seed as RandomCUDA */
+
initializeForcePairs();
}
RigidBodyController::~RigidBodyController() {
for (int i = 0; i < rigidBodyByType.size(); i++)
rigidBodyByType[i].clear();
rigidBodyByType.clear();
+ delete random;
}
void RigidBodyController::initializeForcePairs() {
@@ -140,7 +144,63 @@ void RigidBodyController::updateForces() {
gpuErrchk(cudaDeviceSynchronize());
}
+void RigidBodyController::integrate(int step) {
+ // tell RBs to integrate
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+
+ // thermostat
+ rb.addLangevin( random->gaussian_vector(), random->gaussian_vector() );
+ }
+ }
+ if ( step % conf.outputPeriod == 0 ) { /* PRINT & INTEGRATE */
+ if (step == 0) { // first step so only start this cycle
+ print(step);
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+ rb.integrate(0);
+ }
+ }
+ } else { // finish last cycle
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+ rb.integrate(1);
+ }
+ }
+ print(step);
+
+ // start this cycle
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+ rb.integrate(0);
+ }
+ }
+ }
+ } else { /* INTEGRATE ONLY */
+ if (step == 0) { // first step so only start this cycle
+ print(step);
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+ rb.integrate(0);
+ }
+ }
+ } else { // integrate end of last step and start of this one
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ RigidBody& rb = rigidBodyByType[i][j];
+ rb.integrate(2);
+ }
+ }
+ }
+ }
+}
+
// RBTODO: bundle several rigidbodypair evaluations in single kernel call
void RigidBodyForcePair::updateForces() {
// get the force/torque between a pair of rigid bodies
@@ -155,6 +215,7 @@ void RigidBodyForcePair::updateForces() {
printf(" Calculating grid forces\n");
+ // RBTODO: add energy
computeGridGridForce<<< nb, numThreads >>>
(type1->rawDensityGrids_d[k1], type2->rawPotentialGrids_d[k2],
rb1->getBasis(), rb1->getPosition(), /* RBTODO: include offset from grid */
@@ -254,11 +315,11 @@ void RigidBodyController::copyGridsToDevice() {
for (int i = 0; i < conf.numRigidTypes; i++) {
printf("Working on RB %d\n",i);
- RigidBodyType& rb = conf.rigidBody[i]; /* convenience... may be trouble... */
+ RigidBodyType& rb = conf.rigidBody[i];
int ng = rb.numPotGrids;
- rb.rawPotentialGrids_d = new RigidBodyGrid*[ng]; /* not sure this is needed */
-
+ rb.rawPotentialGrids_d = new RigidBodyGrid*[ng]; /* not 100% sure this is needed, possible memory leak */
+
printf(" RigidBodyType %d: numGrids = %d\n", i, ng);
// copy potential grid data to device
for (int gid = 0; gid < ng; gid++) {
@@ -288,15 +349,284 @@ void RigidBodyController::copyGridsToDevice() {
}
}
- RigidBodyType& rb = conf.rigidBody[0];
gpuErrchk(cudaDeviceSynchronize());
printf("Done copying RBs\n");
- printRigidBodyGrid<<<1,1>>>(conf.rigidBody[0].rawPotentialGrids_d[0]);
+
+ // DEBUG
+ RigidBodyType& rb = conf.rigidBody[0];
+ printRigidBodyGrid<<<1,1>>>( rb.rawPotentialGrids_d[0] );
gpuErrchk(cudaDeviceSynchronize());
- printRigidBodyGrid<<<1,1>>>(conf.rigidBody[0].rawDensityGrids_d[0]);
+ printRigidBodyGrid<<<1,1>>>( rb.rawDensityGrids_d[0] );
gpuErrchk(cudaDeviceSynchronize());
}
+void RigidBodyController::print(int step) {
+ // modeled after outputExtendedData() in Controller.C
+ if ( step >= 0 ) {
+ // Write RIGID BODY trajectory file
+ if ( step % conf.outputPeriod == 0 ) {
+ if ( ! trajFile.rdbuf()->is_open() ) {
+ // open file
+ printf("OPENING RIGID BODY TRAJECTORY FILE\n");
+ // RBTODO: backup_file(simParams->rigidBodyTrajectoryFile);
+
+ char fname[140];
+ strcpy(fname,outArg);
+ strcat(fname, ".rb-traj");
+ trajFile.open(fname);
+
+ while (!trajFile) {
+ /* if ( errno == EINTR ) {
+ printf("Warning: Interrupted system call opening RIGIDBODY trajectory file, retrying.\n");
+ trajFile.clear();
+ trajFile.open(simParams->rigidBodyTrajectoryFile);
+ continue;
+ }
+ */
+ //char err_msg[257];
+ printf("Error opening RigidBody trajectory file %s",fname);
+ exit(1);
+ }
+ trajFile << "# NAMD RigidBody trajectory file" << std::endl;
+ printLegend(trajFile);
+ }
+ printData(step,trajFile);
+ trajFile.flush();
+ }
+
+ // Write restart File
+ /* if ( simParams->restartFrequency && */
+ /* ((step % simParams->restartFrequency) == 0) && */
+ /* (step != simParams->firstTimestep) ) { */
+ if ( step % conf.outputPeriod == 0 && step != 0 ){
+ printf("RIGID BODY: WRITING RESTART FILE AT STEP %d\n", step);
+ char fname[140];
+ strcpy(fname,outArg);
+ strcat(fname, ".rigid");
+ // RBTODO: NAMD_backup_file(fname,".old"); /* */
+ std::ofstream restartFile(fname);
+ while (!restartFile) {
+ /* RBTODO
+ if ( errno == EINTR ) {
+ printf("Warning: Interrupted system call opening rigid body restart file, retrying.\n");
+ restartFile.clear();
+ restartFile.open(fname);
+ continue;
+ }
+ */
+ printf("Error opening rigid body restart file %s",fname);
+ exit(1); // NAMD_err(err_msg);
+ }
+ restartFile << "# NAMD rigid body restart file" << std::endl;
+ printLegend(restartFile);
+ printData(step,restartFile);
+ if (!restartFile) {
+ printf("Error writing rigid body restart file %s",fname);
+ exit(-1); // NAMD_err(err_msg);
+ }
+ }
+ }
+ /*
+ // Output final coordinates
+ if (step == FILE_OUTPUT || step == END_OF_RUN) {
+ int realstep = ( step == FILE_OUTPUT ?
+ simParams->firstTimestep : simParams->N );
+ iout << "WRITING RIGID BODY OUTPUT FILE AT STEP " << realstep << "\n" << endi;
+ static char fname[140];
+ strcpy(fname, simParams->outputFilename);
+ strcat(fname, ".rigid");
+ NAMD_backup_file(fname);
+ std::ofstream outputFile(fname);
+ while (!outputFile) {
+ if ( errno == EINTR ) {
+ CkPrintf("Warning: Interrupted system call opening rigid body output file, retrying.\n");
+ outputFile.clear();
+ outputFile.open(fname);
+ continue;
+ }
+ char err_msg[257];
+ sprintf(err_msg, "Error opening rigid body output file %s",fname);
+ NAMD_err(err_msg);
+ }
+ outputFile << "# NAMD rigid body output file" << std::endl;
+ printLegend(outputFile);
+ printData(realstep,outputFile);
+ if (!outputFile) {
+ char err_msg[257];
+ sprintf(err_msg, "Error writing rigid body output file %s",fname);
+ NAMD_err(err_msg);
+ }
+ }
+
+ // Close trajectory file
+ if (step == END_OF_RUN) {
+ if ( trajFile.rdbuf()->is_open() ) {
+ trajFile.close();
+ iout << "CLOSING RIGID BODY TRAJECTORY FILE\n" << endi;
+ }
+ }
+ */
+}
+void RigidBodyController::printLegend(std::ofstream &file) {
+ file << "#$LABELS step RigidBodyKey"
+ << " posX posY posZ"
+ << " rotXX rotXY rotXZ"
+ << " rotYX rotYY rotYZ"
+ << " rotZX rotZY rotZZ"
+ << " velX velY velZ"
+ << " angVelX angVelY angVelZ" << std::endl;
+}
+void RigidBodyController::printData(int step,std::ofstream &file) {
+ printf("WRITING RIGID BODY COORDINATES AT STEP %d\n",step);
+ // tell RBs to integrate
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ const RigidBody& rb = rigidBodyByType[i][j];
+
+ Vector3 v = rb.getPosition();
+ Matrix3 t = rb.getOrientation();
+ file << step <<" "<< rb.getKey()
+ <<" "<< v.x <<" "<< v.y <<" "<< v.z;
+ file <<" "<< t.exx <<" "<< t.exy <<" "<< t.exz
+ <<" "<< t.eyx <<" "<< t.eyy <<" "<< t.eyz
+ <<" "<< t.ezx <<" "<< t.ezy <<" "<< t.ezz;
+ v = rb.getVelocity();
+ file <<" "<< v.x <<" "<< v.y <<" "<< v.z;
+ v = rb.getAngularVelocity();
+ file <<" "<< v.x <<" "<< v.y <<" "<< v.z
+ << std::endl;
+ }
+ }
+}
+/*
+void RigidBodyController::integrate(int step) {
+ DebugM(3, "RBC::integrate: step "<< step << "\n" << endi);
+
+ DebugM(1, "RBC::integrate: Waiting for grid reduction\n" << endi);
+ gridReduction->require();
+
+ const Molecule * mol = Node::Object()->molecule;
+
+ // pass reduction force and torque to each grid
+ // DebugM(3, "Summing forces on rigid bodies" << "\n" << endi);
+ for (int i=0; i < mol->rbReductionIdToRigidBody.size(); i++) {
+ Force f;
+ Force t;
+ for (int k = 0; k < 3; k++) {
+ f[k] = gridReduction->item(6*i + k);
+ t[k] = gridReduction->item(6*i + k + 3);
+ }
+
+ if (step % 100 == 1)
+ DebugM(4, "RBC::integrate: reduction/rb " << i <<":"
+ << "\n\tposition: "
+ << rigidBodyList[mol->rbReductionIdToRigidBody[i]]->getPosition()
+ <<"\n\torientation: "
+ << rigidBodyList[mol->rbReductionIdToRigidBody[i]]->getOrientation()
+ << "\n" << endi);
+
+ DebugM(4, "RBC::integrate: reduction/rb " << i <<": "
+ << "force " << f <<": "<< "torque: " << t << "\n" << endi);
+ rigidBodyList[mol->rbReductionIdToRigidBody[i]]->addForce(f);
+ rigidBodyList[mol->rbReductionIdToRigidBody[i]]->addTorque(t);
+ }
+
+ // Langevin
+ for (int i=0; i<rigidBodyList.size(); i++) {
+ // continue; // debug
+ if (rigidBodyList[i]->langevin) {
+ DebugM(1, "RBC::integrate: reduction/rb " << i
+ <<": calling langevin" << "\n" << endi);
+ rigidBodyList[i]->addLangevin(
+ random->gaussian_vector(), random->gaussian_vector() );
+ // DebugM(4, "RBC::integrate: reduction/rb " << i
+ // << " after langevin: force " << rigidBodyList[i]f <<": "<< "torque: " << t << "\n" << endi);
+ // <<": calling langevin" << "\n" << endi);
+ }
+ }
+
+ if ( step >= 0 && simParams->rigidBodyOutputFrequency &&
+ (step % simParams->rigidBodyOutputFrequency) == 0 ) {
+ DebugM(1, "RBC::integrate:integrating for before printing output" << "\n" << endi);
+ // PRINT
+ if ( step == simParams->firstTimestep ) {
+ print(step);
+ // first step so only start this cycle
+ for (int i=0; i<rigidBodyList.size(); i++) {
+ DebugM(2, "RBC::integrate: reduction/rb " << i
+ <<": starting integration cycle of step "
+ << step << "\n" << endi);
+ rigidBodyList[i]->integrate(&trans[i],&rot[i],0);
+ }
+ } else {
+ // finish last cycle
+ // DebugM(1, "RBC::integrate: reduction/rb " << i
+ // <<": firststep: calling rb->integrate" << "\n" << endi);
+ for (int i=0; i<rigidBodyList.size(); i++) {
+ DebugM(2, "RBC::integrate: reduction/rb " << i
+ <<": finishing integration cycle of step "
+ << step-1 << "\n" << endi);
+ rigidBodyList[i]->integrate(&trans[i],&rot[i],1);
+ }
+ print(step);
+ // start this cycle
+ for (int i=0; i<rigidBodyList.size(); i++) {
+ DebugM(2, "RBC::integrate: reduction/rb " << i
+ <<": starting integration cycle of step "
+ << step << "\n" << endi);
+ rigidBodyList[i]->integrate(&trans[i],&rot[i],0);
+ }
+ }
+ } else {
+ DebugM(1, "RBC::integrate: trans[0] before: " << trans[0] << "\n" << endi);
+ if ( step == simParams->firstTimestep ) {
+ // integrate the start of this cycle
+ for (int i=0; i<rigidBodyList.size(); i++) {
+ DebugM(2, "RBC::integrate: reduction/rb " << i
+ <<": starting integration cycle of (first) step "
+ << step << "\n" << endi);
+ rigidBodyList[i]->integrate(&trans[i],&rot[i],0);
+ }
+ } else {
+ // integrate end of last ts and start of this one
+ for (int i=0; i<rigidBodyList.size(); i++) {
+ DebugM(2, "RBC::integrate: reduction/rb " << i
+ <<": ending / starting integration cycle of step "
+ << step-1 << "-" << step << "\n" << endi);
+ rigidBodyList[i]->integrate(&trans[i],&rot[i],2);
+ }
+ }
+ DebugM(1, "RBC::integrate: trans[0] after: " << trans[0] << "\n" << endi);
+ }
+
+ DebugM(3, "sendRigidBodyUpdate on step: " << step << "\n" << endi);
+ if (trans.size() != rot.size())
+ NAMD_die("failed sanity check\n");
+ RigidBodyMsg *msg = new RigidBodyMsg;
+ msg->trans.copy(trans); // perhaps .swap() would cause problems
+ msg->rot.copy(rot);
+ computeMgr->sendRigidBodyUpdate(msg);
+}
+
+
+RigidBodyParams* RigidBodyParamsList::find_key(const char* key) {
+ RBElem* cur = head;
+ RBElem* found = NULL;
+ RigidBodyParams* result = NULL;
+
+ while (found == NULL && cur != NULL) {
+ if (!strcasecmp((cur->elem).rigidBodyKey,key)) {
+ found = cur;
+ } else {
+ cur = cur->nxt;
+ }
+ }
+ if (found != NULL) {
+ result = &(found->elem);
+ }
+ return result;
+}
+*/
/* RigidBodyForcePair::RigidBodyForcePair(RigidBodyType* t1, RigidBodyType* t2, */
/* RigidBody* rb1, RigidBody* rb2, */
diff --git a/RigidBodyController.h b/RigidBodyController.h
index 84599ad49d94b0c63708ea7c153e665362baa276..79ed18b4fd8ce175838ef5802a3eb9f00a5c4b08 100644
--- a/RigidBodyController.h
+++ b/RigidBodyController.h
@@ -1,4 +1,5 @@
#include <vector>
+#include <fstream>
#include "RigidBody.h"
/* #include "RandomCUDA.h" */
#include <cuda.h>
@@ -7,6 +8,7 @@
#define NUMTHREADS 256
class Configuration;
+class Random;
class RigidBodyForcePair {
friend class RigidBodyController;
@@ -64,28 +66,28 @@ public:
/* DEVICE RigidBodyController(const NamdState *s, int reductionTag, SimParameters *sp); */
RigidBodyController();
~RigidBodyController();
- RigidBodyController(const Configuration& c);
+ RigidBodyController(const Configuration& c, const char* outArg);
void integrate(int step);
- DEVICE void print(int step);
-
void updateForces();
private:
void copyGridsToDevice();
void initializeForcePairs();
-
- /* void printLegend(std::ofstream &file); */
- /* void printData(int step, std::ofstream &file); */
+
+ void print(int step);
+ void printLegend(std::ofstream &file);
+ void printData(int step, std::ofstream &file);
public:
RigidBodyType** rbType_d;
private:
- /* std::ofstream trajFile; */
-
+ std::ofstream trajFile;
+
const Configuration& conf;
+ const char* outArg;
- /* Random* random; */
+ Random* random;
/* RequireReduction *gridReduction; */
Vector3* trans; // would have made these static, but
@@ -93,4 +95,7 @@ private:
std::vector< std::vector<RigidBody> > rigidBodyByType;
std::vector< RigidBodyForcePair > forcePairs;
+
+
+
};
diff --git a/common.h b/common.h
index 53ce1c62087f9471ab209ebf10e749cdbbecdd5c..ce32520098039a827c789d570397d096f45430a3 100644
--- a/common.h
+++ b/common.h
@@ -112,7 +112,7 @@ void NAMD_die(const char *);
void NAMD_err(const char *); // also prints strerror(errno)
void NAMD_bug(const char *);
void NAMD_backup_file(const char *filename, const char *extension = 0);
-void NAMD_write(int fd, const void *buf, size_t count); // NAMD_die on error
+// void NAMD_write(int fd, const void *buf, size_t count); // NAMD_die on error
char *NAMD_stringdup(const char *);
FILE *Fopen(const char *filename, const char *mode);
int Fclose(FILE *fout);
diff --git a/runBrownTown.cpp b/runBrownTown.cpp
index 5de068f6ab6a3cdcd9c7f1140ad1d2f75bd17d2e..fd1bd4bbfa0cec2d0022119829d5d53892eff204 100644
--- a/runBrownTown.cpp
+++ b/runBrownTown.cpp
@@ -103,7 +103,8 @@ int main(int argc, char* argv[]) {
GPUManager::init();
GPUManager::safe(safe);
Configuration config(configFile, replicas, debug);
- GPUManager::set(0);
+ // GPUManager::set(0);
+ GPUManager::set(1);
config.copyToCUDA();
GrandBrownTown brown(config, outArg, randomSeed,