diff --git a/src/RigidBody.cu b/src/RigidBody.cu
index bb498d5b1fcad1351b909053ddf74427064be6eb..7cf80155739ccaf88352e8f03c0228843fcf043d 100644
--- a/src/RigidBody.cu
+++ b/src/RigidBody.cu
@@ -1,6 +1,8 @@
#include <iostream>
#include <typeinfo>
#include "RigidBody.h"
+#include "RigidBodyType.h"
+#include "RigidBodyController.h"
#include "Configuration.h"
#include "ComputeGridGrid.cuh"
@@ -15,10 +17,10 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t
}
-RigidBody::RigidBody(String name, const Configuration& cref, const RigidBodyType& tref)
- : name(name), c(&cref), t(&tref), impulse_to_momentum(4.184e8f) { init(); }
+RigidBody::RigidBody(String name, const Configuration& cref, const RigidBodyType& tref, RigidBodyController* RBCref)
+ : name(name), c(&cref), t(&tref), RBC(RBCref), impulse_to_momentum(4.184e8f) { init(); }
RigidBody::RigidBody(const RigidBody& rb)
- : name(rb.name), c(rb.c), t(rb.t), impulse_to_momentum(4.184e8f) { init(); }
+ : name(rb.name), c(rb.c), t(rb.t), RBC(rb.RBC), impulse_to_momentum(4.184e8f) { init(); }
void RigidBody::init() {
// units "(kcal_mol/AA) * ns" "dalton AA/ns" * 4.184e+08
timestep = c->timestep;
@@ -194,6 +196,7 @@ void RigidBody::callGridParticleForceKernel(Vector3* pos_d, Vector3* force_d, in
| friction / kt = Diff |
\===========================================================================*/
void RigidBody::addLangevin(Vector3 w1, Vector3 w2) {
+ /* commented out for BD
// w1 and w2 should be standard normal distributions
// in RB frame
@@ -211,6 +214,7 @@ void RigidBody::addLangevin(Vector3 w1, Vector3 w2) {
addForce(f);
addTorque(torq);
+ */
}
/*==========================================================================\
@@ -218,8 +222,7 @@ void RigidBody::addLangevin(Vector3 w1, Vector3 w2) {
| 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(int startFinishAll) {
+void RigidBody::integrateDLM(int startFinishAll) {
Vector3 trans; // = *p_trans;
Matrix3 rot = Matrix3(1); // = *p_rot;
@@ -273,7 +276,7 @@ void RigidBody::integrate(int startFinishAll) {
R = Ry(0.5*timestep * angularMomentum.y / t->inertia.y ); // R4
applyRotation(R);
-
+
R = Rx(0.5*timestep * angularMomentum.x / t->inertia.x ); // R5
applyRotation(R);
@@ -286,6 +289,50 @@ void RigidBody::integrate(int startFinishAll) {
}
}
+/* Following:
+Brownian Dynamics Simulation of Rigid Particles of Arbitrary Shape in External Fields
+Miguel X. Fernandes, José GarcÃa de la Torre
+*/
+void RigidBody::integrate(int startFinishAll) {
+ if (startFinishAll == 1) return;
+
+ Matrix3 rot = Matrix3(1); // = *p_rot;
+
+ if ( isnan(force.x) || isnan(torque.x) ) { // NaN check
+ printf("Rigid Body force or torque was NaN!\n");
+ exit(-1);
+ }
+
+ float kT = 1.0f; // temporary
+ Vector3 diffusion = kT / (t->transDamping*t->mass); // TODO: assign diffusion in config file, or elsewhere
+ Vector3 rotDiffusion = kT / (Vector3::element_mult(t->rotDamping,t->inertia));
+
+ // update positions
+ // trans = Vector(0); if (false) {
+ Vector3 rando = getRandomGaussVector();
+ Vector3 offset = Vector3::element_mult( (diffusion / kT), force ) * timestep +
+ // diffGrad term
+ Vector3::element_mult( Vector3::element_sqrt( 2.0f * diffusion*timestep), rando );
+
+ rando = getRandomGaussVector();
+ Vector3 rotationOffset = Vector3::element_mult( (rotDiffusion / kT) , torque * timestep) +
+ // diffGrad term
+ Vector3::element_mult( Vector3::element_sqrt( 2.0f * rotDiffusion*timestep), rando );
+
+ position += offset;
+
+ // Consider whether a DLM-like decomposition of rotations is needed for time-reversibility
+ orientation = orientation *
+ (Rz(rotationOffset.z) * Ry(rotationOffset.y) * Rx(rotationOffset.x));
+
+ // TODO make this periodic
+ // printf("det: %.12f\n", orientation.det());
+ orientation = orientation.normalized();
+ // orientation = orientation/orientation.det();
+ // printf("det2: %.12f\n", orientation.det());
+ // orientation = orientation/orientation.det(); // TODO: see if this can be somehow eliminated (wasn't in original DLM algorithm...)
+}
+
void RigidBody::applyRotation(const Matrix3& R) {
angularMomentum = R * angularMomentum;
diff --git a/src/RigidBody.h b/src/RigidBody.h
index c8ca40e13ced6b121443f8fe9c5d0b462cd01310..0880a405d626822ac73b7b59e723ba8a91bda8ba 100644
--- a/src/RigidBody.h
+++ b/src/RigidBody.h
@@ -4,8 +4,7 @@
#pragma once
#include "useful.h"
-#include "RigidBodyType.h"
-
+#include "RandomCPU.h" /* for BD integration; RBTODO: fix this */
#ifdef __CUDACC__
#define HOST __host__
@@ -15,8 +14,12 @@
#define DEVICE
#endif
+#include "RigidBodyType.h"
+#include "RigidBodyController.h"
+
class Configuration;
+
typedef float BigReal; /* strip this out later */
typedef Vector3 Force;
@@ -28,7 +31,7 @@ class RigidBody { // host side representation of rigid bodies
| splitting methods for rigid body molecular dynamics". J Chem Phys. (1997) |
`––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/
public:
- RigidBody(String name, const Configuration& c, const RigidBodyType& t);
+ RigidBody(String name, const Configuration& c, const RigidBodyType& t, RigidBodyController* RBC);
RigidBody(const RigidBody& rb);
// RigidBody(const RigidBody& rb) : RigidBody(rb.name, *rb.c, *rb.t) {};
void init();
@@ -44,8 +47,9 @@ class RigidBody { // host side representation of rigid bodies
// HOST DEVICE void integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll);
// HOST DEVICE void integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll);
- void integrate(int startFinishAll);
-
+ void integrateDLM(int startFinishAll);
+ void integrate(int startFinishAll);
+
// HOST DEVICE inline String getKey() const { return key; }
// HOST DEVICE inline String getKey() const { return t->name; }
HOST DEVICE inline String getKey() const { return name; }
@@ -69,6 +73,12 @@ class RigidBody { // host side representation of rigid bodies
Vector3 torque; // lab frame (except in integrate())
private:
+
+ RigidBodyController* RBC;
+ inline Vector3 getRandomGaussVector() {
+ return RBC->getRandomGaussVector();
+ }
+
// String key;
String name;
/* static const SimParameters * simParams; */
diff --git a/src/RigidBodyController.cu b/src/RigidBodyController.cu
index dd2a5e6600bc66bed2b232f98976769665d2a6ad..d2caf094de173d7164a1b55f210cc1af2fae5a99 100644
--- a/src/RigidBodyController.cu
+++ b/src/RigidBodyController.cu
@@ -45,7 +45,7 @@ RigidBodyController::RigidBodyController(const Configuration& c, const char* out
snprintf(tmp, 128, "#%d", j);
name.add( tmp );
}
- RigidBody r(name, conf, conf.rigidBody[i]);
+ RigidBody r(name, conf, conf.rigidBody[i], this);
tmp.push_back( r );
}
rigidBodyByType.push_back(tmp);
diff --git a/src/RigidBodyController.h b/src/RigidBodyController.h
index 35ffe73c1d24bb5d737059d0b706e6b99b642121..50e4373a905f0dd55ba566d8fefffc074bb4178d 100644
--- a/src/RigidBodyController.h
+++ b/src/RigidBodyController.h
@@ -14,8 +14,8 @@
class RigidBodyType;
class RigidBody;
class Configuration;
-class RandomCPU;
-
+// class RandomCPU;
+#include "RandomCPU.h"
// TODO: performance: create RigidBodyGridPair so pairlistdist check is done per grid pair, not per RB pair
class RigidBodyForcePair {
@@ -106,7 +106,12 @@ private:
void printLegend(std::ofstream &file);
void printData(int step, std::ofstream &file);
public:
- RigidBodyType** rbType_d;
+ RigidBodyType** rbType_d;
+
+ inline Vector3 getRandomGaussVector() {
+ return random->gaussian_vector();
+ }
+ /* RequireReduction *gridReduction; */
private:
std::ofstream trajFile;
diff --git a/src/useful.h b/src/useful.h
index 25d2951e7e296b1e33a12b91f11e5a5cd9bdca97..addd00f750ed3cdeba3cfcd24e03bfffcd71a331 100644
--- a/src/useful.h
+++ b/src/useful.h
@@ -257,6 +257,14 @@ HOST DEVICE inline Vector3 operator/(Vector3 v, float s) {
return v*sinv;
}
+HOST DEVICE inline Vector3 operator/(float s, Vector3 v) {
+ v.x = s / v.x;
+ v.y = s / v.y;
+ v.z = s / v.z;
+ return v;
+}
+
+
// class Matrix3
// Operations on 3D float matrices
class Matrix3 {