diff --git a/.gitignore b/.gitignore
index a149351741048ea77b23bb55a661a7661abd4b6b..942163a195a5b09f28795a09f23e0ac58c6e3fbb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,9 @@
+*.org
+*.o
.svn
*~
.backup
TAGS
.dir-locals.el
+BD_example*
+runBrownCUDA
\ No newline at end of file
diff --git a/BaseGrid.cu b/BaseGrid.cu
index 441f83e7038d0320e26f714d18dcdc9b50aa2585..93155ab340c25f466a1e5ef169211557bbb56afb 100644
--- a/BaseGrid.cu
+++ b/BaseGrid.cu
@@ -17,6 +17,10 @@ void BaseGrid::init() {
size = nx*ny*nz;
val = new float[size];
}
+BaseGrid::BaseGrid() {
+ BaseGrid tmp(Matrix3(1),Vector3(0,0,0),1,1,1);
+ *this = tmp; // TODO: verify that this is OK
+}
// The most obvious of constructors.
BaseGrid::BaseGrid(Matrix3 basis0, Vector3 origin0, int nx0, int ny0, int nz0) {
diff --git a/BaseGrid.h b/BaseGrid.h
index 23d7abacf045b348b1e77dfa60f003a3338d6a24..065b71635266ce632151532e7db47a955cb07075 100644
--- a/BaseGrid.h
+++ b/BaseGrid.h
@@ -40,6 +40,8 @@ private:
void init();
public:
+ BaseGrid(); // cmaffeo2 (2015) moved this out of protected, cause I wanted BaseGrid in a struct
+
// The most obvious of constructors.
BaseGrid(Matrix3 basis0, Vector3 origin0, int nx0, int ny0, int nz0);
@@ -547,7 +549,5 @@ public:
Matrix3 basisInv;
public:
float* val;
-protected:
- BaseGrid();
};
#endif
diff --git a/Configuration.cpp b/Configuration.cpp
index 5dc7cf24d2bda80d089ab57df2a5effc76c9e140..d523a06fd748573e8ba6371b8c45ecb6e3909c36 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -1,5 +1,6 @@
#include "Configuration.h"
+
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true) {
if (code != cudaSuccess) {
@@ -173,6 +174,19 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
name = new String[num * simNum];
currSerial = 0;
+ // RigidBodies...
+ if (numRBs > 0) {
+ printf("\nCounting rigid bodies specified in the configuration file.\n");
+ numRB = 0;
+ for (int i = 0; i < numRBs; i++) num += rigidBody[i].num;
+
+ // // state data
+ // rbPos = new Vector3[numRB * simNum];
+ // type = new int[numRB * simNum];
+
+ }
+
+
// Create exclusions from the exclude rule, if it was specified in the config file
if (excludeRule != String("")) {
int oldNumExcludes = numExcludes;
@@ -386,11 +400,14 @@ void Configuration::copyToCUDA() {
printf("Copying to GPU %d\n", GPUManager::current());
BrownianParticleType **part_addr = new BrownianParticleType*[numParts];
-
+ RigidBodyType **rb_addr = new RigidBodyType*[numRBs];
+
// Copy the BaseGrid objects and their member variables/objects
gpuErrchk(cudaMalloc(&part_d, sizeof(BrownianParticleType*) * numParts));
+ gpuErrchk(cudaMalloc(&rbType_d, sizeof(RigidBodyType*) * numRBs));
// TODO: The above line fails when there is not enough memory. If it fails, stop.
+ // TODO: what's going on here?
for (int i = 0; i < numParts; i++) {
BaseGrid *pmf = NULL, *diffusionGrid = NULL;
BrownianParticleType *b = new BrownianParticleType(part[i]);
@@ -431,6 +448,15 @@ void Configuration::copyToCUDA() {
cudaMemcpyHostToDevice));
}
+ // TODO ? utilize Thrust more extensively
+ for (int i = 0; i < numRBs; i++) {
+ gpuErrchk(cudaMalloc(&rbType_d, sizeof(RigidBodyType*) * numRBs));
+ rigidBody[i].potentialGrids_D = rigidBody[i].potentialGrids;
+ rigidBody[i].densityGrids_D = rigidBody[i].densityGrids;
+
+ }
+
+
// kTGrid_d
kTGrid_d = NULL;
if (temperatureGrid.length() > 0) {
@@ -545,6 +571,7 @@ int Configuration::readParameters(const char * config_file) {
// Get the number of particles.
const int numParams = config.length();
numParts = config.countParameter("particle");
+ numRBs = config.countParameter("rigidBody");
// Allocate the particle variables.
part = new BrownianParticleType[numParts];
@@ -559,6 +586,9 @@ int Configuration::readParameters(const char * config_file) {
partTableFile = new String[numParts*numParts];
partTableIndex0 = new int[numParts*numParts];
partTableIndex1 = new int[numParts*numParts];
+
+ // Allocate rigid body types
+ rigidBody = new RigidBodyType[numRBs];
int btfcap = 10;
bondTableFile = new String[btfcap];
@@ -574,9 +604,16 @@ int Configuration::readParameters(const char * config_file) {
int currBond = -1;
int currAngle = -1;
int currDihedral = -1;
+ int currRB = -1;
+
+ int partClassPart = 0;
+ int partClassRB = 1;
+ int currPartClass = -1; // 0 => particle, 1 => rigidBody
+
for (int i = 0; i < numParams; i++) {
String param = config.getParameter(i);
String value = config.getValue(i);
+ // GLOBAL
if (param == String("outputName"))
outputName = value;
else if (param == String("timestep"))
@@ -622,10 +659,11 @@ int Configuration::readParameters(const char * config_file) {
numCap = atoi(value.val());
else if (param == String("decompPeriod"))
decompPeriod = atoi(value.val());
- else if (param == String("particle"))
+ // PARTICLES
+ else if (param == String("particle")) {
part[++currPart] = BrownianParticleType(value);
- else if (param == String("num"))
- part[currPart].num = atoi(value.val());
+ currPartClass = partClassPart;
+ }
else if (param == String("gridFile"))
partGridFile[currPart] = value;
else if (param == String("forceXGridFile"))
@@ -721,12 +759,55 @@ int Configuration::readParameters(const char * config_file) {
}
if (readDihedralFile(value, ++currDihedral))
numTabDihedralFiles++;
- } else {
+ }
+ // RIGID BODY
+ else if (param == String("rigidBody")) {
+ part[++currPart] = BrownianParticleType(value);
+ rigidBody[++currRB] = RigidBodyType(value);
+ currPartClass = partClassRB;
+ }
+ else if (param == String("mass"))
+ rigidBody[currRB].mass = (float) strtod(value.val(), NULL);
+ else if (param == String("inertia"))
+ rigidBody[currRB].inertia = stringToVector3( value );
+ else if (param == String("transDamping"))
+ rigidBody[currRB].transDamping = stringToVector3( value );
+ else if (param == String("rotDamping"))
+ rigidBody[currRB].rotDamping = stringToVector3( value );
+
+ else if (param == String("potentialGrid"))
+ rigidBody[currRB].addPotentialGrid(value);
+ else if (param == String("potentialGrid"))
+ rigidBody[currRB].addPotentialGrid(value);
+
+ // COMMON
+ else if (param == String("num")) {
+ if (currPartClass == partClassPart)
+ part[currPart].num = atoi(value.val());
+ else if (currPartClass == partClassRB)
+ rigidBody[currRB].num = atoi(value.val());
+ }
+ // UNKNOWN
+ else {
printf("WARNING: Unrecognized keyword `%s'.\n", param.val());
}
}
return numParams;
}
+Vector3 Configuration::stringToVector3(String s) {
+ // tokenize and return
+ int numTokens = s.tokenCount();
+ if (numTokens != 3) {
+ printf("ERROR: could not convert input to Vector3.\n"); // TODO improve this message
+ exit(1);
+ }
+ String* token = new String[numTokens];
+ s.tokenize(token);
+ Vector3 v( (float) strtod(token[0], NULL),
+ (float) strtod(token[1], NULL),
+ (float) strtod(token[2], NULL) );
+ return v;
+}
void Configuration::readAtoms() {
// Open the file
diff --git a/Configuration.h b/Configuration.h
index 9e0d18afd1f1ee7ae172c567e4903397e115cf26..ecefd163760334784c2a02901fadd0e11380c8c3 100644
--- a/Configuration.h
+++ b/Configuration.h
@@ -26,6 +26,7 @@
#include "TabulatedAngle.h"
#include "GPUManager.h"
#include "Dihedral.h"
+#include "RigidBodyType.h"
#include <cuda.h>
#include <cuda_runtime.h>
@@ -40,6 +41,7 @@ class Configuration {
};
void setDefaults();
+ Vector3 stringToVector3(String s);
int readParameters(const char* config_file);
@@ -78,6 +80,7 @@ public:
// Device Variables
int *type_d;
BrownianParticleType **part_d;
+ RigidBodyType **rbType_d;
BaseGrid *sys_d, *kTGrid_d;
Bond* bonds_d;
int2* bondMap_d;
@@ -105,6 +108,9 @@ public:
float timeLast; // used with posLast
float minimumSep; // minimum separation allowed with placing new particles
+ // RigidBody variables
+ int numRB;
+
// System parameters
String outputName;
float timestep;
@@ -183,6 +189,12 @@ public:
Dihedral* dihedrals;
String* dihedralTableFile;
int numTabDihedralFiles;
+
+ // RigidBody parameters.
+ RigidBodyType* rigidBody;
+ int numRBs;
+
+
};
#endif
diff --git a/Debug.h b/Debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..c3bee20ae8423c8894d376c99e1d726387378ebb
--- /dev/null
+++ b/Debug.h
@@ -0,0 +1,61 @@
+/**
+*** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
+*** The Board of Trustees of the University of Illinois.
+*** All rights reserved.
+**/
+
+#pragma once
+
+#ifndef MIN_DEBUG_LEVEL
+ #define MIN_DEBUG_LEVEL 0
+#endif
+#ifndef MAX_DEBUG_LEVEL
+ #define MAX_DEBUG_LEVEL 10
+#endif
+#ifndef STDERR_LEVEL
+ /* anything >= this error level goes to stderr */
+ #define STDERR_LEVEL 5
+#endif
+
+
+/*****************************************************************
+ * DebugM(): function to display a debug message.
+ * Messages have different levels. The low numbers are low severity
+ * while the high numbers are really important. Very high numbers
+ * are sent to stderr rather than stdout.
+ * The default severity scale is from 0 to 10.
+ * 0 = plain message
+ * 4 = important message
+ * 5 = warning (stderr)
+ * 10 = CRASH BANG BOOM error (stderr)
+ * The remaining args are like printf: a format string and some args.
+ * This function can be turned off by compiling without the DEBUGM flag
+ * No parameters to this function should have a side effect!
+ * No functions should be passed as parameters! (including inline)
+ *****************************************************************/
+
+#ifdef DEBUGM
+
+#include "InfoStream.h"
+
+ #define Debug(x) (x)
+ #define DebugM(level,format) \
+ { \
+ if ((level >= MIN_DEBUG_LEVEL) && (level <= MAX_DEBUG_LEVEL)) \
+ { \
+ infostream Dout; \
+ if (level >= STDERR_LEVEL) Dout << "[ERROR " << level << "] "; \
+ else if (level > 0) Dout << "[Debug " << level << "] "; \
+ Dout << iPE << ' ' << iFILE; \
+ Dout << format << endi; \
+ } \
+ }
+
+ #else
+ /* make a void function. */
+ /* parameters with side effects will be removed! */
+ #define Debug(x) ;
+ #define DebugM(x,y) ;
+
+ #endif /* DEBUGM */
+
diff --git a/RandomCUDA.cu b/RandomCUDA.cu
index e5b8260a7b945af47bd08a3c0c3031c62c167b29..2d43249149c690ae3b47c58d96dbb720189386c8 100644
--- a/RandomCUDA.cu
+++ b/RandomCUDA.cu
@@ -3,7 +3,7 @@
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true) {
if (code != cudaSuccess) {
- fprintf(stderr,"CUDA Error: %s %s %d\n", cudaGetErrorString(code), file, line);
+ fprintf(stderr,"CUDA Error: %s (%s:%d)\n", cudaGetErrorString(code), file, line);
if (abort) exit(code);
}
}
@@ -11,7 +11,7 @@ inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true) {
#define cuRandchk(ans) { cuRandAssert((ans), __FILE__, __LINE__); }
inline void cuRandAssert(curandStatus code, char *file, int line, bool abort=true) {
if (code != CURAND_STATUS_SUCCESS) {
- fprintf(stderr, "CURAND Error: %d %s %d\n", code, file, line);
+ fprintf(stderr, "CURAND Error: %d (%s:%d)\n", code, file, line);
if (abort) exit(code);
}
}
diff --git a/RigidBody.C b/RigidBody.C
deleted file mode 100644
index f631ceeb0c618e0994d3508e7b11f56fa36be07b..0000000000000000000000000000000000000000
--- a/RigidBody.C
+++ /dev/null
@@ -1,252 +0,0 @@
-/**
-*** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
-*** The Board of Trustees of the University of Illinois.
-*** All rights reserved.
-**/
-
-#ifndef MIN_DEBUG_LEVEL
-#define MIN_DEBUG_LEVEL 5
-#endif
-#define DEBUGM
-#include "Debug.h"
-
-#include <iostream>
-#include <typeinfo>
-
-#include "RigidBody.h"
-#include "Vector.h"
-#include "RigidBodyParams.h"
-
-
-RigidBody::RigidBody(SimParameters *simParams, RigidBodyParams *rbParams)
- : impulse_to_momentum(0.0004184) {
-
- int len = strlen(rbParams->rigidBodyKey);
- key = new char[len+1];
- strncpy(key,rbParams->rigidBodyKey,len+1);
-
- DebugM(4, "Initializing Rigid Body " << key << "\n" << endi);
-
- timestep = simParams->dt;
- Temp = simParams->langevinTemp;
-
- mass = rbParams->mass;
- inertia = rbParams->inertia;
-
- position = rbParams->position;
- orientation = rbParams->orientation; // orientation a matrix that brings a vector from the RB frame to the lab frame
-
- momentum = rbParams->velocity * mass; // lab frame
- DebugM(4, "velocity " << rbParams->velocity << "\n" << endi);
- DebugM(4, "momentum " << momentum << "\n" << endi);
-
- angularMomentum = rbParams->orientationalVelocity; // rigid body frame
- angularMomentum.x *= rbParams->inertia.x;
- angularMomentum.y *= rbParams->inertia.y;
- angularMomentum.z *= rbParams->inertia.z;
-
- isFirstStep = true; // this might not work flawlessly...
-
- clearForce();
- clearTorque();
-
- DebugM(4, "RigidBody initial Force: " << force << "\n" << endi);
-
- // setup for langevin
- langevin = rbParams->langevin;
- if (langevin) {
- // DiffCoeff = kT / dampingCoeff mass
- // rbParams->DampingCoeff has units of (1/ps)
- // f = - dampingCoeff * momentum
- // units "(1/ps) * (amu AA/fs)" "kcal_mol/AA" * 2.3900574
- transDampingCoeff = 2.3900574 * rbParams->transDampingCoeff;
-
- // < f(t) f(t') > = 2 kT dampingCoeff mass delta(t-t')
- // units "sqrt( k K (1/ps) amu / fs )" "kcal_mol/AA" * 0.068916889
- transForceCoeff = 0.068916889 * element_sqrt( 2*Temp*mass*rbParams->transDampingCoeff/timestep );
-
- // T = - dampingCoeff * angularMomentum
- rotDampingCoeff = 2.3900574 * rbParams->rotDampingCoeff;
- // < f(t) f(t') > = 2 kT dampingCoeff inertia delta(t-t')
- rotTorqueCoeff = 0.068916889 * element_sqrt( 2*Temp*element_mult(inertia, rbParams->rotDampingCoeff) / timestep );
- }
-}
-
-void RigidBody::addForce(Force f) {
- DebugM(1, "RigidBody "<<key<<" adding f ("<<f<<") to Force " << force << "\n" << endi);
- force += f;
-}
-void RigidBody::addTorque(Force t) {
- DebugM(1, "RigidBody adding t ("<<t<<") to torque " << torque << "\n" << endi);
- torque += t;
-}
-RigidBody::~RigidBody() {}
-
-void RigidBody::addLangevin(Vector w1, Vector w2) {
- // w1 and w2 should be standard normal distributions
-
- /*****************************************************************
- 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
- *****************************************************************/
- // BUT: assume diagonal friction tensor
- // and no Wiener process / stochastic calculus
- // then this is just the same as for translation
- // < T_i(t) T_i(t) > = 2 kT friction inertia
- // friction / kt = Diff
-
- // in RB frame
- Force f = element_mult(transForceCoeff,w1) -
- element_mult(transDampingCoeff, transpose(orientation)*momentum);
-
- Force t = element_mult(rotTorqueCoeff,w2) -
- element_mult(rotDampingCoeff, angularMomentum);
-
- f = orientation * f; // return to lab frame
- t = orientation * t;
-
- addForce(f);
- addTorque(t);
-}
-
-void RigidBody::integrate(Vector *p_trans, Tensor *p_rot, int startFinishAll) {
-// 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
- Vector trans; // = *p_trans;
- Tensor rot = Tensor::identity(); // = *p_rot;
-
-#ifdef DEBUGM
- switch (startFinishAll) {
- case 0: // start
- DebugM(2, "Rigid Body integrating start of cycle" << "\n" << endi);
- case 1: // finish
- DebugM(2, "Rigid Body integrating finish of cycle" << "\n" << endi);
- case 2: // finish and start
- DebugM(2, "Rigid Body integrating finishing last cycle, starting this one" << "\n" << endi);
- }
-#endif
-
- if ( isnan(force.x) || isnan(torque.x) ) // NaN check
- NAMD_die("Rigid Body force was NaN!\n");
-
- // torque = Vector(0,0,10); // debug
- Force tmpTorque = transpose(orientation)*torque; // bring to rigid body frame
-
- DebugM(3, "integrate" <<": force "<<force <<": velocity "<<getVelocity() << "\n" << endi);
- DebugM(3, "integrate" <<": torque "<<tmpTorque <<": orientationalVelocity "<<getAngularVelocity() << "\n" << endi);
-
- if (startFinishAll == 0 || startFinishAll == 1) {
- // propogate momenta by half step
- momentum += 0.5 * timestep * force * impulse_to_momentum;
- angularMomentum += 0.5 * timestep * tmpTorque * impulse_to_momentum;
- } else {
- // propogate momenta by a full timestep
- momentum += timestep * force * impulse_to_momentum;
- angularMomentum += timestep * tmpTorque * impulse_to_momentum;
- }
-
- DebugM(3, " position before: " << position << "\n" << endi);
-
- if (startFinishAll == 0 || startFinishAll == 2) {
- // update positions
- // trans = Vector(0); if (false) {
- trans = timestep * momentum / mass;
- position += trans; // update CoM a full timestep
- // }
-
- // update orientations a full timestep
- Tensor R; // represents a rotation about a principle axis
- R = Rx(0.5*timestep * angularMomentum.x / inertia.x ); // R1
- angularMomentum = R * angularMomentum;
- rot = transpose(R);
- DebugM(1, "R: " << R << "\n" << endi);
- DebugM(1, "Rot 1: " << rot << "\n" << endi);
-
- R = Ry(0.5*timestep * angularMomentum.y / inertia.y ); // R2
- angularMomentum = R * angularMomentum;
- rot = rot * transpose(R);
- DebugM(1, "R: " << R << "\n" << endi);
- DebugM(1, "Rot 2: " << rot << "\n" << endi);
-
- R = Rz( timestep * angularMomentum.z / inertia.z ); // R3
- angularMomentum = R * angularMomentum;
- rot = rot * transpose(R);
- DebugM(1, "R: " << R << "\n" << endi);
- DebugM(1, "Rot 3: " << rot << "\n" << endi);
-
- R = Ry(0.5*timestep * angularMomentum.y / inertia.y ); // R4
- angularMomentum = R * angularMomentum;
- rot = rot * transpose(R);
- DebugM(1, "R: " << R << "\n" << endi);
- DebugM(1, "Rot 4: " << rot << "\n" << endi);
-
- R = Rx(0.5*timestep * angularMomentum.x / inertia.x ); // R5
- angularMomentum = R * angularMomentum;
- rot = rot * transpose(R);
- DebugM(1, "R: " << R << "\n" << endi);
- DebugM(1, "Rot 5: " << rot << "\n" << endi);
-
- // DebugM(3,"TEST: " << Ry(0.01) <<"\n" << endi); // DEBUG
-
- // update actual orientation
- Tensor newOrientation = orientation*rot; // not 100% sure; rot could be in rb frame
- orientation = newOrientation;
- rot = transpose(rot);
-
- DebugM(2, "trans during: " << trans
- << "\n" << endi);
- DebugM(2, "rot during: " << rot
- << "\n" << endi);
-
- clearForce();
- clearTorque();
-
- *p_trans = trans;
- *p_rot = rot;
- }
- DebugM(3, " position after: " << position << "\n" << endi);
-}
-
-// Rotations about axes
-
-// for very small angles 10^-8, cos^2+sin^2 != 1
-// concerned about the accumulation of errors in non-unitary transformations!
-Tensor RigidBody::Rx(BigReal t) {
- BigReal qt = 0.25*t*t; // for approximate calculations of sin(t) and cos(t)
- BigReal cos = (1-qt)/(1+qt);
- BigReal sin = t/(1+qt);
-
- Tensor tmp;
- tmp.xx = 1; tmp.xy = 0; tmp.xz = 0;
- tmp.yx = 0; tmp.yy = cos; tmp.yz = -sin;
- tmp.zx = 0; tmp.zy = sin; tmp.zz = cos;
- return tmp;
-}
-Tensor RigidBody::Ry(BigReal t) {
- BigReal qt = 0.25*t*t; // for approximate calculations of sin(t) and cos(t)
- BigReal cos = (1-qt)/(1+qt);
- BigReal sin = t/(1+qt);
-
- Tensor tmp;
- tmp.xx = cos; tmp.xy = 0; tmp.xz = sin;
- tmp.yx = 0; tmp.yy = 1; tmp.yz = 0;
- tmp.zx = -sin; tmp.zy = 0; tmp.zz = cos;
- return tmp;
-}
-Tensor RigidBody::Rz(BigReal t) {
- BigReal qt = 0.25*t*t; // for approximate calculations of sin(t) and cos(t)
- BigReal cos = (1-qt)/(1+qt);
- BigReal sin = t/(1+qt);
-
- Tensor tmp;
- tmp.xx = cos; tmp.xy = -sin; tmp.xz = 0;
- tmp.yx = sin; tmp.yy = cos; tmp.yz = 0;
- tmp.zx = 0; tmp.zy = 0; tmp.zz = 1;
- return tmp;
-}
-Tensor RigidBody::eulerToMatrix(const Vector e) {
- // convert euler angle input to rotation matrix
- // http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Conversion_formulae_between_formalisms
- return Rz(e.z) * Ry(e.y) * Rx(e.x);
-}
diff --git a/RigidBody.c b/RigidBody.c
new file mode 100644
index 0000000000000000000000000000000000000000..36bd3f85f2c5cea9e96197cc4c3c3ffb3c0626cc
--- /dev/null
+++ b/RigidBody.c
@@ -0,0 +1,265 @@
+/**
+*** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
+*** The Board of Trustees of the University of Illinois.
+*** All rights reserved.
+**/
+
+#ifndef MIN_DEBUG_LEVEL
+#define MIN_DEBUG_LEVEL 5
+#endif
+#define DEBUGM
+#include "Debug.h"
+
+#include <iostream>
+#include <typeinfo>
+
+#include "RigidBody.h"
+#include "Vector.h"
+#include "RigidBodyParams.h"
+
+
+RigidBody::RigidBody(SimParameters *simParams, RigidBodyParams *rbParams)
+: impulse_to_momentum(0.0004184) {
+
+ int len = strlen(rbParams->rigidBodyKey);
+ key = new char[len+1];
+ strncpy(key,rbParams->rigidBodyKey,len+1);
+
+ DebugM(4, "Initializing Rigid Body " << key << "\n" << endi);
+
+ timestep = simParams->dt;
+ Temp = simParams->langevinTemp;
+
+ mass = rbParams->mass;
+ inertia = rbParams->inertia;
+
+ position = rbParams->position;
+ orientation = rbParams->orientation; // orientation a matrix that brings a vector from the RB frame to the lab frame
+
+ momentum = rbParams->velocity * mass; // lab frame
+ DebugM(4, "velocity " << rbParams->velocity << "\n" << endi);
+ DebugM(4, "momentum " << momentum << "\n" << endi);
+
+ angularMomentum = rbParams->orientationalVelocity; // rigid body frame
+ angularMomentum.x *= rbParams->inertia.x;
+ angularMomentum.y *= rbParams->inertia.y;
+ angularMomentum.z *= rbParams->inertia.z;
+
+ isFirstStep = true; // this might not work flawlessly...
+
+ clearForce();
+ clearTorque();
+
+ DebugM(4, "RigidBody initial Force: " << force << "\n" << endi);
+
+ // setup for langevin
+ langevin = rbParams->langevin;
+ if (langevin) {
+ /*=======================================================\
+ | DiffCoeff = kT / dampingCoeff mass |
+ | |
+ | rbParams->DampingCoeff has units of (1/ps) |
+ | |
+ | f = - dampingCoeff * momentum |
+ | |
+ | units "(1/ps) * (amu AA/fs)" "kcal_mol/AA" * 2.3900574 |
+ \=======================================================*/
+ transDampingCoeff = 2.3900574 * rbParams->transDampingCoeff;
+
+ /*================================================================\
+ | < f(t) f(t') > = 2 kT dampingCoeff mass delta(t-t') |
+ | |
+ | units "sqrt( k K (1/ps) amu / fs )" "kcal_mol/AA" * 0.068916889 |
+ \================================================================*/
+ transForceCoeff = 0.068916889 * element_sqrt( 2*Temp*mass*rbParams->transDampingCoeff/timestep );
+
+ // T = - dampingCoeff * angularMomentum
+ rotDampingCoeff = 2.3900574 * rbParams->rotDampingCoeff;
+ // < f(t) f(t') > = 2 kT dampingCoeff inertia delta(t-t')
+ rotTorqueCoeff = 0.068916889 * element_sqrt( 2*Temp*element_mult(inertia, rbParams->rotDampingCoeff) / timestep );
+ }
+}
+
+void RigidBody::addForce(Force f) {
+ DebugM(1, "RigidBody "<<key<<" adding f ("<<f<<") to Force " << force << "\n" << endi);
+ force += f;
+}
+void RigidBody::addTorque(Force t) {
+ DebugM(1, "RigidBody adding t ("<<t<<") to torque " << torque << "\n" << endi);
+ torque += t;
+}
+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 |
+| |
+| |
+| BUT: assume diagonal friction tensor and no Wiener process / stochastic |
+| calculus then this is just the same as for translation |
+| |
+| < T_i(t) T_i(t) > = 2 kT friction inertia |
+| |
+| friction / kt = Diff |
+\===========================================================================*/
+void RigidBody::addLangevin(Vector w1, Vector w2) {
+ // w1 and w2 should be standard normal distributions
+
+ // in RB frame
+ Force f = element_mult(transForceCoeff,w1) -
+ element_mult(transDampingCoeff, transpose(orientation)*momentum);
+
+ Force t = element_mult(rotTorqueCoeff,w2) -
+ element_mult(rotDampingCoeff, angularMomentum);
+
+ f = orientation * f; // return to lab frame
+ t = orientation * t;
+
+ addForce(f);
+ addTorque(t);
+}
+
+/*==========================================================================\
+| 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(Vector *p_trans, Tensor *p_rot, int startFinishAll) {
+ Vector trans; // = *p_trans;
+ Tensor rot = Tensor::identity(); // = *p_rot;
+
+#ifdef DEBUGM
+ switch (startFinishAll) {
+ case 0: // start
+ DebugM(2, "Rigid Body integrating start of cycle" << "\n" << endi);
+ case 1: // finish
+ DebugM(2, "Rigid Body integrating finish of cycle" << "\n" << endi);
+ case 2: // finish and start
+ DebugM(2, "Rigid Body integrating finishing last cycle, starting this one" << "\n" << endi);
+ }
+#endif
+
+ if ( isnan(force.x) || isnan(torque.x) ) // NaN check
+ NAMD_die("Rigid Body force was NaN!\n");
+
+ // torque = Vector(0,0,10); // debug
+ Force tmpTorque = transpose(orientation)*torque; // bring to rigid body frame
+
+ DebugM(3, "integrate" <<": force "<<force <<": velocity "<<getVelocity() << "\n" << endi);
+ DebugM(3, "integrate" <<": torque "<<tmpTorque <<": orientationalVelocity "<<getAngularVelocity() << "\n" << endi);
+
+ if (startFinishAll == 0 || startFinishAll == 1) {
+ // propogate momenta by half step
+ momentum += 0.5 * timestep * force * impulse_to_momentum;
+ angularMomentum += 0.5 * timestep * tmpTorque * impulse_to_momentum;
+ } else {
+ // propogate momenta by a full timestep
+ momentum += timestep * force * impulse_to_momentum;
+ angularMomentum += timestep * tmpTorque * impulse_to_momentum;
+ }
+
+ DebugM(3, " position before: " << position << "\n" << endi);
+
+ if (startFinishAll == 0 || startFinishAll == 2) {
+ // update positions
+ // trans = Vector(0); if (false) {
+ trans = timestep * momentum / mass;
+ position += trans; // update CoM a full timestep
+ // }
+
+ // update orientations a full timestep
+ Tensor R; // represents a rotation about a principle axis
+ R = Rx(0.5*timestep * angularMomentum.x / inertia.x ); // R1
+ angularMomentum = R * angularMomentum;
+ rot = transpose(R);
+ DebugM(1, "R: " << R << "\n" << endi);
+ DebugM(1, "Rot 1: " << rot << "\n" << endi);
+
+ R = Ry(0.5*timestep * angularMomentum.y / inertia.y ); // R2
+ angularMomentum = R * angularMomentum;
+ rot = rot * transpose(R);
+ DebugM(1, "R: " << R << "\n" << endi);
+ DebugM(1, "Rot 2: " << rot << "\n" << endi);
+
+ R = Rz( timestep * angularMomentum.z / inertia.z ); // R3
+ angularMomentum = R * angularMomentum;
+ rot = rot * transpose(R);
+ DebugM(1, "R: " << R << "\n" << endi);
+ DebugM(1, "Rot 3: " << rot << "\n" << endi);
+
+ R = Ry(0.5*timestep * angularMomentum.y / inertia.y ); // R4
+ angularMomentum = R * angularMomentum;
+ rot = rot * transpose(R);
+ DebugM(1, "R: " << R << "\n" << endi);
+ DebugM(1, "Rot 4: " << rot << "\n" << endi);
+
+ R = Rx(0.5*timestep * angularMomentum.x / inertia.x ); // R5
+ angularMomentum = R * angularMomentum;
+ rot = rot * transpose(R);
+ DebugM(1, "R: " << R << "\n" << endi);
+ DebugM(1, "Rot 5: " << rot << "\n" << endi);
+
+ // DebugM(3,"TEST: " << Ry(0.01) <<"\n" << endi); // DEBUG
+
+ // update actual orientation
+ Tensor newOrientation = orientation*rot; // not 100% sure; rot could be in rb frame
+ orientation = newOrientation;
+ rot = transpose(rot);
+
+ DebugM(2, "trans during: " << trans
+ << "\n" << endi);
+ DebugM(2, "rot during: " << rot
+ << "\n" << endi);
+
+ clearForce();
+ clearTorque();
+
+ *p_trans = trans;
+ *p_rot = rot;
+ }
+ DebugM(3, " position after: " << position << "\n" << endi);
+}
+
+// Rotations about axes
+
+// for very small angles 10^-8, cos^2+sin^2 != 1
+// concerned about the accumulation of errors in non-unitary transformations!
+Tensor RigidBody::Rx(BigReal t) {
+ BigReal qt = 0.25*t*t; // for approximate calculations of sin(t) and cos(t)
+ BigReal cos = (1-qt)/(1+qt);
+ BigReal sin = t/(1+qt);
+
+ Tensor tmp;
+ tmp.xx = 1; tmp.xy = 0; tmp.xz = 0;
+ tmp.yx = 0; tmp.yy = cos; tmp.yz = -sin;
+ tmp.zx = 0; tmp.zy = sin; tmp.zz = cos;
+ return tmp;
+}
+Tensor RigidBody::Ry(BigReal t) {
+ BigReal qt = 0.25*t*t; // for approximate calculations of sin(t) and cos(t)
+ BigReal cos = (1-qt)/(1+qt);
+ BigReal sin = t/(1+qt);
+
+ Tensor tmp;
+ tmp.xx = cos; tmp.xy = 0; tmp.xz = sin;
+ tmp.yx = 0; tmp.yy = 1; tmp.yz = 0;
+ tmp.zx = -sin; tmp.zy = 0; tmp.zz = cos;
+ return tmp;
+}
+Tensor RigidBody::Rz(BigReal t) {
+ BigReal qt = 0.25*t*t; // for approximate calculations of sin(t) and cos(t)
+ BigReal cos = (1-qt)/(1+qt);
+ BigReal sin = t/(1+qt);
+
+ Tensor tmp;
+ tmp.xx = cos; tmp.xy = -sin; tmp.xz = 0;
+ tmp.yx = sin; tmp.yy = cos; tmp.yz = 0;
+ tmp.zx = 0; tmp.zy = 0; tmp.zz = 1;
+ return tmp;
+}
+Tensor RigidBody::eulerToMatrix(const Vector e) {
+ // convert euler angle input to rotation matrix
+ // http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Conversion_formulae_between_formalisms
+ return Rz(e.z) * Ry(e.y) * Rx(e.x);
+}
diff --git a/RigidBodyController.C b/RigidBodyController.c
similarity index 68%
rename from RigidBodyController.C
rename to RigidBodyController.c
index bb5030d050faab0f55c21f34cd22ed58c415062d..ad6dcc14bab7c7ab9cffb39620a24b7bd48664c3 100644
--- a/RigidBodyController.C
+++ b/RigidBodyController.c
@@ -1,8 +1,3 @@
-/**
-*** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
-*** The Board of Trustees of the University of Illinois.
-*** All rights reserved.
-**/
#ifndef MIN_DEBUG_LEVEL
#define MIN_DEBUG_LEVEL 5
#endif
@@ -31,6 +26,19 @@
#include <errno.h>
+/*============================\
+| include "RigidBodyParams.h" |
+\============================*/
+
+#include "MStream.h"
+
+#include "DataStream.h"
+#include "InfoStream.h"
+
+#include "Debug.h"
+#include <stdio.h>
+
+
RigidBodyController::RigidBodyController(const NamdState *s, const int reductionTag, SimParameters *sp) : state(s), simParams(sp)
{
DebugM(2, "Rigid Body Controller initializing"
@@ -314,3 +322,196 @@ void RigidBodyController::integrate(int step) {
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;
+}
+
+int RigidBodyParamsList::index_for_key(const char* key)
+{
+ RBElem* cur = head;
+ RBElem* found = NULL;
+ int result = -1;
+
+ int idx = 0;
+ while (found == NULL && cur != NULL) {
+ if (!strcasecmp((cur->elem).rigidBodyKey,key)) {
+ found = cur;
+ } else {
+ cur = cur->nxt;
+ idx++;
+ }
+ }
+ if (found != NULL) {
+ result = idx;
+ }
+ return result;
+}
+
+RigidBodyParams* RigidBodyParamsList::add(const char* key)
+{
+ // If the key is already in the list, we can't add it
+ if (find_key(key)!=NULL) {
+ return NULL;
+ }
+
+ RBElem* new_elem = new RBElem();
+ int len = strlen(key);
+ RigidBodyParams* elem = &(new_elem->elem);
+ elem->rigidBodyKey = new char[len+1];
+ strncpy(elem->rigidBodyKey,key,len+1);
+ elem->mass = NULL;
+ elem->inertia = Vector(NULL);
+ elem->langevin = NULL;
+ elem->temperature = NULL;
+ elem->transDampingCoeff = Vector(NULL);
+ elem->rotDampingCoeff = Vector(NULL);
+ elem->gridList.clear();
+ elem->position = Vector(NULL);
+ elem->velocity = Vector(NULL);
+ elem->orientation = Tensor();
+ elem->orientationalVelocity = Vector(NULL);
+
+ elem->next = NULL;
+ new_elem->nxt = NULL;
+ if (head == NULL) {
+ head = new_elem;
+ }
+ if (tail != NULL) {
+ tail->nxt = new_elem;
+ tail->elem.next = elem;
+ }
+ tail = new_elem;
+ n_elements++;
+
+ return elem;
+}
+
+const void RigidBodyParams::print() {
+ iout << iINFO
+ << "printing RigidBodyParams("<<rigidBodyKey<<"):"
+ <<"\n\t" << "mass: " << mass
+ <<"\n\t" << "inertia: " << inertia
+ <<"\n\t" << "langevin: " << langevin
+ <<"\n\t" << "temperature: " << temperature
+ <<"\n\t" << "transDampingCoeff: " << transDampingCoeff
+ <<"\n\t" << "position: " << position
+ <<"\n\t" << "orientation: " << orientation
+ <<"\n\t" << "orientationalVelocity: " << orientationalVelocity
+ << "\n" << endi;
+
+}
+const void RigidBodyParamsList::print() {
+ iout << iINFO << "Printing " << n_elements << " RigidBodyParams\n" << endi;
+
+ RigidBodyParams *elem = get_first();
+ while (elem != NULL) {
+ elem->print();
+ elem = elem->next;
+ }
+}
+const void RigidBodyParamsList::print(char *s) {
+ iout << iINFO << "("<<s<<") Printing " << n_elements << " RigidBodyParams\n" << endi;
+
+ RigidBodyParams *elem = get_first();
+ while (elem != NULL) {
+ elem->print();
+ elem = elem->next;
+ }
+}
+
+void RigidBodyParamsList::pack_data(MOStream *msg) {
+ DebugM(4, "Packing rigid body parameter list\n");
+ print();
+
+ int i = n_elements;
+ msg->put(n_elements);
+
+ RigidBodyParams *elem = get_first();
+ while (elem != NULL) {
+ DebugM(4, "Packing a new element\n");
+
+ int len;
+ Vector v;
+
+ len = strlen(elem->rigidBodyKey) + 1;
+ msg->put(len);
+ msg->put(len,elem->rigidBodyKey);
+ msg->put(elem->mass);
+
+ // v = elem->
+ msg->put(&(elem->inertia));
+ msg->put( (elem->langevin?1:0) );
+ msg->put(elem->temperature);
+ msg->put(&(elem->transDampingCoeff));
+ msg->put(&(elem->rotDampingCoeff));
+
+ // elem->gridList.clear();
+
+ msg->put(&(elem->position));
+ msg->put(&(elem->velocity));
+ // Tensor data = elem->orientation;
+ msg->put( & elem->orientation );
+ msg->put(&(elem->orientationalVelocity)) ;
+
+ i--;
+ elem = elem->next;
+ }
+ if (i != 0)
+ NAMD_die("MGridforceParams message packing error\n");
+}
+void RigidBodyParamsList::unpack_data(MIStream *msg) {
+ DebugM(4, "Could be unpacking rigid body parameterlist (not used & not implemented)\n");
+
+ int elements;
+ msg->get(elements);
+
+ for(int i=0; i < elements; i++) {
+ DebugM(4, "Unpacking a new element\n");
+
+ int len;
+ msg->get(len);
+ char *key = new char[len];
+ msg->get(len,key);
+ RigidBodyParams *elem = add(key);
+ delete [] key;
+
+ msg->get(&(elem->inertia));
+
+ int j;
+ msg->get(j);
+ elem->langevin = (j != 0);
+
+ msg->get(elem->temperature);
+ msg->get(&(elem->transDampingCoeff));
+ msg->get(&(elem->rotDampingCoeff));
+
+ // elem->gridList.clear();
+
+ msg->get(&(elem->position));
+ msg->get(&(elem->velocity));
+ msg->get( & elem->orientation );
+ msg->get(&(elem->orientationalVelocity)) ;
+
+ elem = elem->next;
+ }
+
+ DebugM(4, "Finished unpacking rigid body parameter list\n");
+ print();
+}
diff --git a/RigidBodyParams.C b/RigidBodyParams.C
deleted file mode 100644
index de6213508b42c57e756ce2700970d97387509743..0000000000000000000000000000000000000000
--- a/RigidBodyParams.C
+++ /dev/null
@@ -1,206 +0,0 @@
-/*
- * RigidBodyParams.C
- */
-
-#define DEBUGM
-#include "RigidBodyParams.h"
-#include "MStream.h"
-
-#include "DataStream.h"
-#include "InfoStream.h"
-
-#include "Debug.h"
-#include <stdio.h>
-
-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;
-}
-
-int RigidBodyParamsList::index_for_key(const char* key)
-{
- RBElem* cur = head;
- RBElem* found = NULL;
- int result = -1;
-
- int idx = 0;
- while (found == NULL && cur != NULL) {
- if (!strcasecmp((cur->elem).rigidBodyKey,key)) {
- found = cur;
- } else {
- cur = cur->nxt;
- idx++;
- }
- }
- if (found != NULL) {
- result = idx;
- }
- return result;
-}
-
-RigidBodyParams* RigidBodyParamsList::add(const char* key)
-{
- // If the key is already in the list, we can't add it
- if (find_key(key)!=NULL) {
- return NULL;
- }
-
- RBElem* new_elem = new RBElem();
- int len = strlen(key);
- RigidBodyParams* elem = &(new_elem->elem);
- elem->rigidBodyKey = new char[len+1];
- strncpy(elem->rigidBodyKey,key,len+1);
- elem->mass = NULL;
- elem->inertia = Vector(NULL);
- elem->langevin = NULL;
- elem->temperature = NULL;
- elem->transDampingCoeff = Vector(NULL);
- elem->rotDampingCoeff = Vector(NULL);
- elem->gridList.clear();
- elem->position = Vector(NULL);
- elem->velocity = Vector(NULL);
- elem->orientation = Tensor();
- elem->orientationalVelocity = Vector(NULL);
-
- elem->next = NULL;
- new_elem->nxt = NULL;
- if (head == NULL) {
- head = new_elem;
- }
- if (tail != NULL) {
- tail->nxt = new_elem;
- tail->elem.next = elem;
- }
- tail = new_elem;
- n_elements++;
-
- return elem;
-}
-
-const void RigidBodyParams::print() {
- iout << iINFO
- << "printing RigidBodyParams("<<rigidBodyKey<<"):"
- <<"\n\t" << "mass: " << mass
- <<"\n\t" << "inertia: " << inertia
- <<"\n\t" << "langevin: " << langevin
- <<"\n\t" << "temperature: " << temperature
- <<"\n\t" << "transDampingCoeff: " << transDampingCoeff
- <<"\n\t" << "position: " << position
- <<"\n\t" << "orientation: " << orientation
- <<"\n\t" << "orientationalVelocity: " << orientationalVelocity
- << "\n" << endi;
-
-}
-const void RigidBodyParamsList::print() {
- iout << iINFO << "Printing " << n_elements << " RigidBodyParams\n" << endi;
-
- RigidBodyParams *elem = get_first();
- while (elem != NULL) {
- elem->print();
- elem = elem->next;
- }
-}
-const void RigidBodyParamsList::print(char *s) {
- iout << iINFO << "("<<s<<") Printing " << n_elements << " RigidBodyParams\n" << endi;
-
- RigidBodyParams *elem = get_first();
- while (elem != NULL) {
- elem->print();
- elem = elem->next;
- }
-}
-
-void RigidBodyParamsList::pack_data(MOStream *msg) {
- DebugM(4, "Packing rigid body parameter list\n");
- print();
-
- int i = n_elements;
- msg->put(n_elements);
-
- RigidBodyParams *elem = get_first();
- while (elem != NULL) {
- DebugM(4, "Packing a new element\n");
-
- int len;
- Vector v;
-
- len = strlen(elem->rigidBodyKey) + 1;
- msg->put(len);
- msg->put(len,elem->rigidBodyKey);
- msg->put(elem->mass);
-
- // v = elem->
- msg->put(&(elem->inertia));
- msg->put( (elem->langevin?1:0) );
- msg->put(elem->temperature);
- msg->put(&(elem->transDampingCoeff));
- msg->put(&(elem->rotDampingCoeff));
-
- // elem->gridList.clear();
-
- msg->put(&(elem->position));
- msg->put(&(elem->velocity));
- // Tensor data = elem->orientation;
- msg->put( & elem->orientation );
- msg->put(&(elem->orientationalVelocity)) ;
-
- i--;
- elem = elem->next;
- }
- if (i != 0)
- NAMD_die("MGridforceParams message packing error\n");
-}
-void RigidBodyParamsList::unpack_data(MIStream *msg) {
- DebugM(4, "Could be unpacking rigid body parameterlist (not used & not implemented)\n");
-
- int elements;
- msg->get(elements);
-
- for(int i=0; i < elements; i++) {
- DebugM(4, "Unpacking a new element\n");
-
- int len;
- msg->get(len);
- char *key = new char[len];
- msg->get(len,key);
- RigidBodyParams *elem = add(key);
- delete [] key;
-
- msg->get(&(elem->inertia));
-
- int j;
- msg->get(j);
- elem->langevin = (j != 0);
-
- msg->get(elem->temperature);
- msg->get(&(elem->transDampingCoeff));
- msg->get(&(elem->rotDampingCoeff));
-
- // elem->gridList.clear();
-
- msg->get(&(elem->position));
- msg->get(&(elem->velocity));
- msg->get( & elem->orientation );
- msg->get(&(elem->orientationalVelocity)) ;
-
- elem = elem->next;
- }
-
- DebugM(4, "Finished unpacking rigid body parameter list\n");
- print();
-
-}
diff --git a/RigidBodyType.cpp b/RigidBodyType.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2dec83a64019a69385681c564f5fe2af620f89f5
--- /dev/null
+++ b/RigidBodyType.cpp
@@ -0,0 +1,59 @@
+#include "RigidBodyType.h"
+
+void RigidBodyType::clear() {
+ num = 0; // TODO: not 100% sure about this
+ if (reservoir != NULL) delete reservoir;
+ reservoir = NULL;
+ // pmf = NULL;
+ mass = 1.0;
+
+ // TODO: make sure that this actually removes grid data
+ potentialGrids.clear();
+ densityGrids.clear();
+ potentialGrids_D.clear();
+ densityGrids_D.clear();
+
+}
+
+// void RigidBodyType::copy(const RigidBodyType& src) {
+// this = new RigidBodyType(src.name);
+// num = src.num;
+// // if (src.pmf != NULL) pmf = new BaseGrid(*src.pmf);
+// if (src.reservoir != NULL) reservoir = new Reservoir(*src.reservoir);
+
+// numPotentialGrid = src.numPotentialGrid;
+// // TODO: fix this: BaseGrid*[]
+// for (int i=0; i < numPotentialGrid; i++) {
+
+// }
+
+// numDensityGrid = src.numDensityGrid;
+// }
+
+// RigidBodyType& RigidBodyType::operator=(const RigidBodyType& src) {
+// clear();
+// copy(src);
+// return *this;
+// }
+
+
+KeyGrid RigidBodyType::createKeyGrid(String s) {
+ // tokenize and return
+ int numTokens = s.tokenCount();
+ if (numTokens != 2) {
+ printf("ERROR: could not add Grid.\n"); // TODO improve this message
+ exit(1);
+ }
+ String* token = new String[numTokens];
+ s.tokenize(token);
+ KeyGrid g;
+ g.key = token[0];
+ g.grid = * new BaseGrid(token[1]);
+ return g;
+}
+void RigidBodyType::addPotentialGrid(String s) {
+ potentialGrids.push_back( createKeyGrid(s) );
+}
+void RigidBodyType::addDensityGrid(String s) {
+ densityGrids.push_back( createKeyGrid(s) );
+}
diff --git a/RigidBodyType.h b/RigidBodyType.h
new file mode 100644
index 0000000000000000000000000000000000000000..571ae8aa5a84268671c77138ab0e57cb7edaa1f2
--- /dev/null
+++ b/RigidBodyType.h
@@ -0,0 +1,68 @@
+// RigidBodyType.h (2015)
+// Author: Chris Maffeo <cmaffeo2@illinois.edu>
+
+#pragma once
+// #include <vector>
+#include <thrust/host_vector.h>
+#include <thrust/device_vector.h>
+#include "Reservoir.h"
+#include "useful.h"
+#include "BaseGrid.h"
+
+// Stores particle type's potential grid and other information
+struct KeyGrid {
+ String key;
+ BaseGrid grid;
+KeyGrid() :
+ key(NULL), grid() { }
+};
+
+class RigidBodyType {
+private:
+ // Deletes all members
+ void clear();
+ // void copy(const RigidBodyType& src);
+
+ KeyGrid createKeyGrid(String s);
+
+public:
+RigidBodyType(const String& name = "") :
+ name(name), num(0),
+ reservoir(NULL), mass(1.0f), inertia(), transDamping(),
+ rotDamping(), potentialGrids(NULL), densityGrids(NULL),
+ potentialGrids_D(NULL), densityGrids_D(NULL) { }
+
+ /* RigidBodyType(const RigidBodyType& src) { copy(src); } */
+ ~RigidBodyType() { clear(); }
+
+ /* RigidBodyType& operator=(const RigidBodyType& src); */
+
+ // crop
+ // Crops all BaseGrid members
+ // @param boundries to crop to (x0, y0, z0) -> (x1, y1, z1);
+ // whether to change the origin
+ // @return success of function (if false nothing was done)
+ /* bool crop(int x0, int y0, int z0, int x1, int y1, int z1, bool keep_origin); */
+
+ void addPotentialGrid(String s);
+ void addDensityGrid(String s);
+
+public:
+ String name;
+ int num; // number of particles of this type
+
+ Reservoir* reservoir;
+ /* BaseGrid* pmf; */
+
+ float mass;
+ Vector3 inertia;
+ Vector3 transDamping;
+ Vector3 rotDamping;
+
+ /* std::vector<KeyGrid> potentialGrids; */
+ /* std::vector<KeyGrid> densityGrids; */
+ thrust::host_vector<KeyGrid> potentialGrids;
+ thrust::device_vector<KeyGrid> potentialGrids_D;
+ thrust::host_vector<KeyGrid> densityGrids;
+ thrust::device_vector<KeyGrid> densityGrids_D;
+};
diff --git a/makefile b/makefile
index ca3ef362444d45d38e0e3d49faf8065d802a8e77..32d4d6d33bae46224c1a89ec9dad2c0b5e0314e9 100644
--- a/makefile
+++ b/makefile
@@ -1,7 +1,7 @@
### Paths, libraries, includes, options
# CUDA_PATH ?= /Developer/NVIDIA/CUDA-6.5
-CUDA_PATH ?= /software/cuda-toolkit-4.1-x86_64/cuda
+CUDA_PATH ?= /software/cuda-toolkit-6.5-x86_64
#CUDA_PATH ?= /usr/local/encap/cuda-5.5
@@ -30,8 +30,9 @@ endif
LD_FLAGS = -L$(LIBRARY) -lcurand -lcudart -Wl,-rpath,$(LIBRARY)
#CODE_10 := -gencode arch=compute_10,code=sm_10
-CODE_12 := -gencode arch=compute_12,code=sm_12
+#CODE_12 := -gencode arch=compute_12,code=sm_12
#CODE_20 := -gencode arch=compute_20,code=sm_20
+CODE_20 := -arch=sm_20
#CODE_30 := -gencode arch=compute_30,code=sm_30
#CODE_35 := -gencode arch=compute_35,code=\"sm_35,compute_35\"