From f891f45dd2c8948d3d667bd375b1d2091db24cc3 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Mon, 4 Jan 2016 11:37:23 -0600
Subject: [PATCH] code appears to work; much to be done; force evaluation,
integration need to be tested
---
Configuration.cpp | 5 ++
Random.h | 2 +-
RandomCPU.h | 200 +++++++++++++++++++++++++++++++++++++++++
RigidBody.cu | 2 +
RigidBody.h | 5 +-
RigidBodyController.cu | 9 +-
RigidBodyController.h | 4 +-
RigidBodyType.cu | 31 +++----
RigidBodyType.h | 3 +-
notes.org | 7 ++
10 files changed, 242 insertions(+), 26 deletions(-)
create mode 100644 RandomCPU.h
diff --git a/Configuration.cpp b/Configuration.cpp
index 01278da..3a75992 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -792,6 +792,11 @@ int Configuration::readParameters(const char * config_file) {
printf("WARNING: Unrecognized keyword `%s'.\n", param.val());
}
}
+
+ // extra configuration for RB types
+ for (int i = 0; i < numRigidTypes; i++)
+ rigidBody[i].setDampingCoeffs(timestep);
+
return numParams;
}
Vector3 Configuration::stringToVector3(String s) {
diff --git a/Random.h b/Random.h
index 378a09d..4dcd231 100644
--- a/Random.h
+++ b/Random.h
@@ -154,7 +154,7 @@ public:
v1 = 2.0 * uniform() - 1.0;
v2 = 2.0 * uniform() - 1.0;
r = v1*v1 + v2*v2;
- printf("r %f", r);
+ // printf("r %f", r);
}
fac = sqrt(-2.0 * log(r)/r);
// now make the Box-Muller transformation to get two normally
diff --git a/RandomCPU.h b/RandomCPU.h
new file mode 100644
index 0000000..c8d3b9b
--- /dev/null
+++ b/RandomCPU.h
@@ -0,0 +1,200 @@
+/**
+*** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
+*** The Board of Trustees of the University of Illinois.
+*** All rights reserved.
+**/
+
+/*
+ * Copyright (c) 1993 Martin Birgmeier
+ * All rights reserved.
+ *
+ * You may redistribute unmodified or modified versions of this source
+ * code provided that the above copyright notice and this and the
+ * following conditions are retained.
+ *
+ * This software is provided ``as is'', and comes with no warranties
+ * of any kind. I shall in no event be liable for anything that happens
+ * to anyone/anything when using this software.
+ */
+
+/*––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––.
+| Redundant with Random.h but with RandomCPU class since Random class is also |
+| declared in RandomCUDA.h, and I want to use both for rigidBody code |
+`––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/
+// RBTODO: make Random.h include RandomCPU.h and "typedef" this class
+
+#ifndef RANDOM_CPU_H
+#define RANDOM_CPU_H
+
+#ifdef __CUDACC__
+#define HOST __host__
+#define DEVICE __device__
+#else
+#define HOST
+#define DEVICE
+#endif
+
+#include "common.h"
+#include "useful.h"
+
+#ifdef _MSC_VER
+#define INT64_LITERAL(X) X ## i64
+#else
+#define INT64_LITERAL(X) X ## LL
+#endif
+
+#define RAND48_SEED INT64_LITERAL(0x00001234abcd330e)
+#define RAND48_MULT INT64_LITERAL(0x00000005deece66d)
+#define RAND48_ADD INT64_LITERAL(0x000000000000000b)
+#define RAND48_MASK INT64_LITERAL(0x0000ffffffffffff)
+
+class RandomCPU {
+
+private:
+
+ float second_gaussian;
+ int64 second_gaussian_waiting;
+ int64 rand48_seed;
+ int64 rand48_mult;
+ int64 rand48_add;
+
+public:
+
+ // default constructor
+ RandomCPU(void) {
+ init(0);
+ rand48_seed = RAND48_SEED;
+ }
+
+ // constructor with seed
+ RandomCPU(unsigned long seed) {
+ init(seed);
+ }
+
+ // reinitialize with seed
+ HOST DEVICE inline void init(unsigned long seed) {
+ second_gaussian = 0;
+ second_gaussian_waiting = 0;
+ rand48_seed = seed & INT64_LITERAL(0x00000000ffffffff);
+ rand48_seed = rand48_seed << 16;
+ rand48_seed |= RAND48_SEED & INT64_LITERAL(0x0000ffff);
+ rand48_mult = RAND48_MULT;
+ rand48_add = RAND48_ADD;
+ print("INIT\n");
+ }
+
+ HOST DEVICE inline void print(char* string)
+ {
+ printf(string);
+ printf("RAND48_SEED = %d, RAND48_MULT = %d, RAND48_ADD = %d, RAND48_MASK = %d\n", RAND48_SEED, RAND48_MULT, RAND48_ADD, RAND48_MASK);
+ printf("rand48_seed = %d, rand48_mult = %d, rand48_add = %d\n", rand48_seed, rand48_mult, rand48_add);
+ }
+
+ // advance generator by one (seed = seed * mult + add, to 48 bits)
+ HOST DEVICE inline void skip(void) {
+ rand48_seed = ( rand48_seed * rand48_mult + rand48_add ) & RAND48_MASK;
+ }
+
+ // split into numStreams different steams and take stream iStream
+ void split(int iStream, int numStreams) {
+ int i;
+
+ // make sure that numStreams is odd to ensure maximum period
+ numStreams |= 1;
+
+ // iterate to get to the correct stream
+ for ( i = 0; i < iStream; ++i ) skip();
+
+ // save seed and add so we can use skip() for our calculations
+ int64 save_seed = rand48_seed;
+
+ // calculate c *= ( 1 + a + ... + a^(numStreams-1) )
+ rand48_seed = rand48_add;
+ for ( i = 1; i < numStreams; ++i ) skip();
+ int64 new_add = rand48_seed;
+
+ // calculate a = a^numStreams
+ rand48_seed = rand48_mult;
+ rand48_add = 0;
+ for ( i = 1; i < numStreams; ++i ) skip();
+ rand48_mult = rand48_seed;
+
+ rand48_add = new_add;
+ rand48_seed = save_seed;
+
+ second_gaussian = 0;
+ second_gaussian_waiting = 0;
+ print("END SPLIT\n");
+ }
+
+ // return a number uniformly distributed between 0 and 1
+ HOST DEVICE inline BigReal uniform(void) {
+ skip();
+ const float exp48 = ( 1.0 / (float)(INT64_LITERAL(1) << 48) );
+ return ( (float) rand48_seed * exp48 );
+ }
+
+ long poisson(BigReal lambda) {
+ const BigReal l = exp(-lambda);
+ long k = 0;
+ BigReal p = uniform();
+
+ while (p >= l) {
+ p *= uniform();
+ k = k + 1;
+ }
+
+ return k;
+ }
+
+ // return a number from a standard gaussian distribution
+ HOST DEVICE inline BigReal gaussian(void) {
+ BigReal fac, r, v1, v2;
+
+ if (second_gaussian_waiting) {
+ second_gaussian_waiting = 0;
+ return second_gaussian;
+ } else {
+ r = 2.; // r >= 1.523e-8 ensures abs result < 6
+ while (r >=1. || r < 1.523e-8) { // make sure we are within unit circle
+ v1 = 2.0 * uniform() - 1.0;
+ v2 = 2.0 * uniform() - 1.0;
+ r = v1*v1 + v2*v2;
+ // printf("r %f\n", r);
+ }
+ fac = sqrt(-2.0 * log(r)/r);
+ // now make the Box-Muller transformation to get two normally
+ // distributed random numbers. Save one and return the other.
+ second_gaussian_waiting = 1;
+ second_gaussian = v1 * fac;
+ return v2 * fac;
+ }
+
+ }
+
+ // return a vector of gaussian random numbers
+ HOST DEVICE inline Vector3 gaussian_vector(void) {
+ return Vector3( gaussian(), gaussian(), gaussian() );
+ }
+
+ // return a random long
+ HOST DEVICE inline long integer(void) {
+ skip();
+ return ( ( rand48_seed >> 17 ) & INT64_LITERAL(0x000000007fffffff) );
+ }
+
+ // randomly order an array of whatever
+ template <class Elem> void reorder(Elem *a, int n) {
+ for ( int i = 0; i < (n-1); ++i ) {
+ int ie = i + ( integer() % (n-i) );
+ if ( ie == i ) continue;
+ const Elem e = a[ie];
+ a[ie] = a[i];
+ a[i] = e;
+ }
+ }
+
+};
+
+#endif // RANDOM_CPU_H
+
diff --git a/RigidBody.cu b/RigidBody.cu
index 01739b8..9a4b1a4 100644
--- a/RigidBody.cu
+++ b/RigidBody.cu
@@ -92,6 +92,8 @@ void RigidBody::integrate(int startFinishAll) {
Vector3 trans; // = *p_trans;
Matrix3 rot = Matrix3(1); // = *p_rot;
+ printf("Rigid Body force\n");
+
#ifdef DEBUGM
switch (startFinishAll) {
case 0: // start
diff --git a/RigidBody.h b/RigidBody.h
index cebbf0f..46c827a 100644
--- a/RigidBody.h
+++ b/RigidBody.h
@@ -42,7 +42,8 @@ class RigidBody { // host side representation of rigid bodies
// HOST DEVICE void integrate(Vector3& old_trans, Matrix3& old_rot, int startFinishAll);
HOST DEVICE void integrate(int startFinishAll);
- HOST DEVICE inline String getKey() const { return key; }
+ // HOST DEVICE inline String getKey() const { return key; }
+ HOST DEVICE inline String getKey() const { return t->name; }
HOST DEVICE inline Vector3 getPosition() const { return position; }
HOST DEVICE inline Matrix3 getOrientation() const { return orientation; }
@@ -78,7 +79,7 @@ private:
// integration
const Configuration* c;
RigidBodyType* t; /* RBTODO: const? */
- int timestep;
+ float timestep;
Vector3 force; // lab frame
Vector3 torque; // lab frame (except in integrate())
diff --git a/RigidBodyController.cu b/RigidBodyController.cu
index 1fe941d..253ec54 100644
--- a/RigidBodyController.cu
+++ b/RigidBodyController.cu
@@ -13,7 +13,7 @@
// #include <vector>
#include "Debug.h"
-#include "Random.h" /* RBTODO: fix this? */
+#include "RandomCPU.h" /* RBTODO: fix this? */
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, String file, int line, bool abort=true) {
@@ -46,7 +46,7 @@ RigidBodyController::RigidBodyController(const Configuration& c, const char* out
rigidBodyByType.push_back(tmp);
}
- random = new Random(conf.seed + 1); /* +1 to avoid using same seed as RandomCUDA */
+ random = new RandomCPU(conf.seed + 1); /* +1 to avoid using same seed as RandomCUDA */
initializeForcePairs();
}
@@ -151,6 +151,7 @@ void RigidBodyController::integrate(int step) {
RigidBody& rb = rigidBodyByType[i][j];
// thermostat
+ // rb.addLangevin( random->gaussian_vector(), random->gaussian_vector() );
rb.addLangevin( random->gaussian_vector(), random->gaussian_vector() );
}
}
@@ -387,7 +388,7 @@ void RigidBodyController::print(int step) {
printf("Error opening RigidBody trajectory file %s",fname);
exit(1);
}
- trajFile << "# NAMD RigidBody trajectory file" << std::endl;
+ trajFile << "# RigidBody trajectory file" << std::endl;
printLegend(trajFile);
}
printData(step,trajFile);
@@ -417,7 +418,7 @@ void RigidBodyController::print(int step) {
printf("Error opening rigid body restart file %s",fname);
exit(1); // NAMD_err(err_msg);
}
- restartFile << "# NAMD rigid body restart file" << std::endl;
+ restartFile << "# RigidBody restart file" << std::endl;
printLegend(restartFile);
printData(step,restartFile);
if (!restartFile) {
diff --git a/RigidBodyController.h b/RigidBodyController.h
index 79ed18b..5bfde4e 100644
--- a/RigidBodyController.h
+++ b/RigidBodyController.h
@@ -8,7 +8,7 @@
#define NUMTHREADS 256
class Configuration;
-class Random;
+class RandomCPU;
class RigidBodyForcePair {
friend class RigidBodyController;
@@ -87,7 +87,7 @@ private:
const Configuration& conf;
const char* outArg;
- Random* random;
+ RandomCPU* random;
/* RequireReduction *gridReduction; */
Vector3* trans; // would have made these static, but
diff --git a/RigidBodyType.cu b/RigidBodyType.cu
index eca2b46..2d9ae68 100644
--- a/RigidBodyType.cu
+++ b/RigidBodyType.cu
@@ -1,7 +1,7 @@
#include "RigidBodyType.h"
void RigidBodyType::clear() {
- num = 0; // TODO: not 100% sure about this
+ num = 0; // RBTODO: not 100% sure about this
if (reservoir != NULL) delete reservoir;
reservoir = NULL;
// pmf = NULL;
@@ -59,19 +59,17 @@ void RigidBodyType::clear() {
// }
-void RigidBodyType::setDampingCoeffs(float timestep, float tmp_mass, Vector3 tmp_inertia, float tmp_transDamping, float tmp_rotDamping) {
- mass = tmp_mass;
- inertia = tmp_inertia;
- /*–––––––––––––––––––––––––––––––––––––––––––––––––––––––––.
+// void RigidBodyType::setDampingCoeffs(float timestep, float tmp_mass, Vector3 tmp_inertia, float tmp_transDamping, float tmp_rotDamping) {
+void RigidBodyType::setDampingCoeffs(float timestep) { /* MUST ONLY BE CALLED ONCE!!! */
+ /*–––––––––––––––––––––––––––––––––––––––––––––––––––––––.
| DiffCoeff = kT / dampingCoeff mass |
- | |
+ | |
| type->DampingCoeff has units of (1/ps) |
- | |
+ | |
| f[kcal/mol AA] = - dampingCoeff * momentum[amu AA/fs] |
- | |
+ | |
| units "(1/ps) * (amu AA/fs)" "kcal_mol/AA" * 2.3900574 |
- `–––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/
- transDamping = 2.3900574 * tmp_transDamping;
+ `–––––––––––––––––––––––––––––––––––––––––––––––––––––––*/
/*––––––––––––––––––––––––––––––––––––––––––––––––––––.
| < f(t) f(t') > = 2 kT dampingCoeff mass delta(t-t') |
@@ -79,19 +77,22 @@ void RigidBodyType::setDampingCoeffs(float timestep, float tmp_mass, Vector3 tmp
| units "sqrt( k K (1/ps) amu / fs )" "kcal_mol/AA" |
| * 0.068916889 |
`––––––––––––––––––––––––––––––––––––––––––––––––––––*/
- float Temp = 295; /* RBTODO: Fix!!!! */
- transForceCoeff = 0.068916889 * Vector3::element_sqrt( 2*Temp*mass*tmp_transDamping/timestep );
+ float Temp = 295; /* RBTODO: temperature should be read from grid? Or set in uniformly in config file */
+ transForceCoeff = 0.068916889 * Vector3::element_sqrt( 2*Temp*mass*transDamping/timestep );
// setup for langevin
// langevin = rbParams->langevin;
// if (langevin) {
// T = - dampingCoeff * angularMomentum
- rotDamping = 2.3900574 * tmp_rotDamping;
// < f(t) f(t') > = 2 kT dampingCoeff inertia delta(t-t')
rotTorqueCoeff = 0.068916889 *
- Vector3::element_sqrt( 2*Temp* Vector3::element_mult(inertia,tmp_rotDamping) / timestep );
- // }
+ Vector3::element_sqrt( 2*Temp* Vector3::element_mult(inertia,rotDamping) / timestep );
+
+
+ transDamping = 2.3900574 * transDamping;
+ rotDamping = 2.3900574 * rotDamping;
+
}
void RigidBodyType::addPotentialGrid(String s) {
diff --git a/RigidBodyType.h b/RigidBodyType.h
index 04eb55a..707b3d1 100644
--- a/RigidBodyType.h
+++ b/RigidBodyType.h
@@ -38,8 +38,7 @@ RigidBodyType(const String& name = "") :
void addPotentialGrid(String s);
void addDensityGrid(String s);
void updateRaw();
- void setDampingCoeffs(float timestep, float tmp_mass, Vector3 tmp_inertia,
- float tmp_transDamping, float tmp_rotDamping);
+ void setDampingCoeffs(float timestep);
public:
diff --git a/notes.org b/notes.org
index 28b02e9..b64defb 100644
--- a/notes.org
+++ b/notes.org
@@ -2,6 +2,13 @@
** make rigid body object device pointer
** branch rb-device: have RBType objects hold device pointers for raw RigidBodyGrids
+* TODO eventually
+** read restart file
+** multiple "PMF" grids associated with each RB type
+** each RB gets temperature from grid
+** optionally don't eval forces every ts
+
+
* organization
** RgidBodyController (RBC) holds device pointers, manages force evaluation and integration
--
GitLab