From a07ceae20245b64c46e95998b64461684b292900 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Wed, 19 Oct 2016 11:08:16 -0500
Subject: [PATCH] Fixed issues with coordinate wrapping during integration;
improved implementation of select matrix operations
---
BaseGrid.h | 63 ++++++++++++++++++++-------------
ComputeForce.cuh | 4 +--
Configuration.cpp | 4 +++
Configuration.h | 1 +
GrandBrownTown.cuh | 20 +++++------
RigidBody.cu | 31 ++++++++---------
notes.org | 12 +------
useful.cu | 78 -----------------------------------------
useful.h | 86 +++++++++++++++++++++++++++++++++++++---------
9 files changed, 139 insertions(+), 160 deletions(-)
diff --git a/BaseGrid.h b/BaseGrid.h
index 4eca85d..3df015c 100644
--- a/BaseGrid.h
+++ b/BaseGrid.h
@@ -644,6 +644,14 @@ public:
}
// Wrap coordinate: 0 <= x < l
+ HOST DEVICE inline int quotient(float x, float l) const {
+#ifdef __CUDA_ARCH__
+ return int(floorf( __fdividef(x,l) ));
+#else
+ return int(floor(x/l));
+#endif
+ }
+
HOST DEVICE inline float wrapFloat(float x, float l) const {
int image = int(floor(x/l));
x -= image*l;
@@ -651,14 +659,7 @@ public:
}
// Wrap distance: -0.5*l <= x < 0.5*l
- DEVICE static inline float d_wrapDiff(float x, float l) {
- // int image = int(floor(x/l));
- int image = int(floorf( __fdividef(x,l) ));
- x -= image*l;
- x = (x >= 0.5f * l) ? x-l : x;
- return x;
- }
- HOST DEVICE static inline float wrapDiff(float x, float l) {
+ HOST DEVICE static inline float wrapDiff(float x, float l) {
int image = int(floor(x/l));
x -= image*l;
if (x >= 0.5f * l)
@@ -666,15 +667,39 @@ public:
return x;
}
+ // TODO: implement simpler approach when basis is diagonal
+ // TODO: implement device version using __fdividef and floorf
+ // TODO: make BaseGrid an abstract class that diagGrid and nonDiagGrid inherit from
+ // TODO: test wrap and wrapDiff for non-diagonal basis
// Wrap vector, 0 <= x < lx && 0 <= y < ly && 0 <= z < lz
HOST DEVICE inline Vector3 wrap(Vector3 r) const {
- Vector3 l = basisInv.transform(r - origin);
- l.x = wrapFloat(l.x, nx);
- l.y = wrapFloat(l.y, ny);
- l.z = wrapFloat(l.z, nz);
- return basis.transform(l) + origin;
+ Vector3 l = basisInv.transform(r - origin);
+ if ( basis.isDiagonal() ) {
+ r = r - Vector3(quotient(l.x,nx) * nx*basis.exx,
+ quotient(l.y,ny) * ny*basis.eyy,
+ quotient(l.z,nz) * nz*basis.ezz);
+ } else {
+ r = r - quotient(l.x,nx) * nx*basis.ex();
+ r = r - quotient(l.y,ny) * ny*basis.ey();
+ r = r - quotient(l.z,nz) * nz*basis.ez();
+ }
+ return r;
}
-
+
+ HOST DEVICE inline Vector3 wrapDiff(Vector3 r) const {
+ Vector3 l = basisInv.transform(r);
+ if ( basis.isDiagonal() ) {
+ r = r - Vector3(quotient(l.x+0.5f*nx,nx) * nx*basis.exx,
+ quotient(l.y+0.5f*ny,ny) * ny*basis.eyy,
+ quotient(l.z+0.5f*nz,nz) * nz*basis.ezz);
+ } else {
+ r = r - quotient(l.x+0.5f*nx,nx) * nx*basis.ex();
+ r = r - quotient(l.y+0.5f*ny,ny) * ny*basis.ey();
+ r = r - quotient(l.z+0.5f*nz,nz) * nz*basis.ez();
+ }
+ return r;
+ }
+
// Wrap vector distance, -0.5*l <= x < 0.5*l && ...
/* HOST DEVICE inline Vector3 wrapDiff(Vector3 r) const { */
/* Vector3 l = basisInv.transform(r); */
@@ -683,21 +708,13 @@ public:
/* l.z = wrapDiff(l.z, nz); */
/* return basis.transform(l); */
/* } */
- HOST DEVICE inline Vector3 wrapDiff(Vector3 r) const {
+ HOST DEVICE inline Vector3 wrapDiffOrig(Vector3 r) const {
Vector3 l = basisInv.transform(r);
l.x = wrapDiff(l.x, nx);
l.y = wrapDiff(l.y, ny);
l.z = wrapDiff(l.z, nz);
return basis.transform(l);
}
- DEVICE inline Vector3 d_wrapDiff(Vector3 r) const {
- Vector3 l = basisInv.transform(r);
- l.x = d_wrapDiff(l.x, nx);
- l.y = d_wrapDiff(l.y, ny);
- l.z = d_wrapDiff(l.z, nz);
- return basis.transform(l);
- }
-
Vector3 wrapDiffNearest(Vector3 r) const;
// Includes the home node.
diff --git a/ComputeForce.cuh b/ComputeForce.cuh
index d9fe31b..906b169 100644
--- a/ComputeForce.cuh
+++ b/ComputeForce.cuh
@@ -471,7 +471,7 @@ __global__ void computeTabulatedKernel(
/* printf("tid,bid,pos[ai(%d)]: %d %d %f %f %f\n", ai, threadIdx.x, blockIdx.x, pos[ai].x, pos[ai].y, pos[ai].z); //*/
Vector3 dr = pos[aj] - pos[ai];
- dr = sys->d_wrapDiff(dr);
+ dr = sys->wrapDiff(dr);
// Calculate the force based on the tabulatedPotential file
float d2 = dr.length2();
@@ -502,7 +502,7 @@ __global__ void computeTabulatedEnergyKernel(Vector3* force, const Vector3* __re
// RBTODO: implement wrapDiff2, returns dr2 (???)
Vector3 dr = pos[aj] - pos[ai];
- dr = sys->d_wrapDiff(dr);
+ dr = sys->wrapDiff(dr);
// Calculate the force based on the tabulatedPotential file
float d2 = dr.length2();
/* RBTODO: filter tabpot particles ahead of time */
diff --git a/Configuration.cpp b/Configuration.cpp
index 42b9f27..a52b56b 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -649,6 +649,10 @@ int Configuration::readParameters(const char * config_file) {
else if (param == String("kT"))
kT = (float) strtod(value.val(), NULL);
// else if (param == String("kTGridFile")) kTGridFile = value;
+ 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;
else if (param == String("numberFluct"))
diff --git a/Configuration.h b/Configuration.h
index 9d23bb8..7fc75f1 100644
--- a/Configuration.h
+++ b/Configuration.h
@@ -128,6 +128,7 @@ public:
int tabulatedPotential;
int fullLongRange;
float kT;
+ float temperature;
float coulombConst;
float electricField;
float cutoff;
diff --git a/GrandBrownTown.cuh b/GrandBrownTown.cuh
index 736970e..157978d 100644
--- a/GrandBrownTown.cuh
+++ b/GrandBrownTown.cuh
@@ -82,6 +82,7 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
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); */
@@ -90,20 +91,20 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
diffusion = diff.e;
diffGrad = diff.f;
}
+
+ // if (idx == 0) {
+ // printf("force: "); force.print();
+ // }
- /* printf("atom %d: pos: %f %f %f\n", idx, p.x, p.y, p.z); */
- /* printf("atom %d: force: %f %f %f\n", idx, force.x, force.y, force.z); */
- /* printf("atom %d: kTlocal, diffusion, timestep: %f, %f, %f\n", */
- /* idx, kTlocal, diffusion, timestep); */
Vector3 tmp = step(p, kTlocal, force, diffusion, -diffGrad, timestep, sys, randoGen, num);
- assert( tmp.length() < 10000.0f );
+ // assert( tmp.length() < 10000.0f );
pos[idx] = tmp;
}
}
__device__
-Vector3 step(Vector3 r0, float kTlocal, Vector3 force, float diffusion,
+inline Vector3 step(Vector3 r0, float kTlocal, Vector3 force, float diffusion,
Vector3 diffGrad, float timestep, BaseGrid *sys,
Random *randoGen, int num) {
const int idx = blockIdx.x * blockDim.x + threadIdx.x;
@@ -111,16 +112,11 @@ Vector3 step(Vector3 r0, float kTlocal, Vector3 force, float diffusion,
Vector3 rando = randoGen->gaussian_vector(idx, num);
diffusion *= timestep;
-
Vector3 r = r0 + (diffusion / kTlocal) * force
+ timestep * diffGrad
+ sqrtf(2.0f * diffusion) * rando;
// Wrap about periodic boundaries
- Vector3 l = sys->getInverseBasis().transform(r - sys->getOrigin());
- l.x = sys->wrapFloat(l.x, sys->getNx());
- l.y = sys->wrapFloat(l.y, sys->getNy());
- l.z = sys->wrapFloat(l.z, sys->getNz());
- return sys->getBasis().transform(l) + sys->getOrigin();
+ return sys->wrap(r);
}
__global__ void devicePrint(RigidBodyType* rb[]) {
diff --git a/RigidBody.cu b/RigidBody.cu
index b18be2a..bd614d1 100644
--- a/RigidBody.cu
+++ b/RigidBody.cu
@@ -98,8 +98,8 @@ void RigidBody::integrate(int startFinishAll) {
if (startFinishAll == 0 || startFinishAll == 1) {
// propogate momenta by half step
- momentum += 0.5 * timestep * force * impulse_to_momentum;
- angularMomentum += 0.5 * timestep * orientation.transpose()*torque * impulse_to_momentum;
+ momentum += 0.5f * timestep * force * impulse_to_momentum;
+ angularMomentum += 0.5f * timestep * (orientation.transpose()*torque) * impulse_to_momentum;
} else {
// propogate momenta by a full timestep
momentum += timestep * force * impulse_to_momentum;
@@ -156,33 +156,30 @@ Matrix3 RigidBody::Rx(BigReal t) {
BigReal cos = (1-qt)/(1+qt);
BigReal sin = t/(1+qt);
- Matrix3 tmp;
- tmp.exx = 1; tmp.exy = 0; tmp.exz = 0;
- tmp.eyx = 0; tmp.eyy = cos; tmp.eyz = -sin;
- tmp.ezx = 0; tmp.ezy = sin; tmp.ezz = cos;
- return tmp;
+ return Matrix3(
+ 1.0f, 0.0f, 0.0f,
+ 0.0f, cos, -sin,
+ 0.0f, sin, cos);
}
Matrix3 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);
- Matrix3 tmp;
- tmp.exx = cos; tmp.exy = 0; tmp.exz = sin;
- tmp.eyx = 0; tmp.eyy = 1; tmp.eyz = 0;
- tmp.ezx = -sin; tmp.ezy = 0; tmp.ezz = cos;
- return tmp;
+ return Matrix3(
+ cos, 0.0f, sin,
+ 0.0f, 1.0f, 0.0f,
+ -sin, 0.0f, cos);
}
Matrix3 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);
- Matrix3 tmp;
- tmp.exx = cos; tmp.exy = -sin; tmp.exz = 0;
- tmp.eyx = sin; tmp.eyy = cos; tmp.eyz = 0;
- tmp.ezx = 0; tmp.ezy = 0; tmp.ezz = 1;
- return tmp;
+ return Matrix3(
+ cos, -sin, 0.0f,
+ sin, cos, 0.0f,
+ 0.0f, 0.0f, 1.0f);
}
Matrix3 RigidBody::eulerToMatrix(const Vector3 e) {
// convert euler angle input to rotation matrix
diff --git a/notes.org b/notes.org
index 4ba69a8..c748a8d 100644
--- a/notes.org
+++ b/notes.org
@@ -1,16 +1,6 @@
* TODO active
-
-** re-implement exclusions
-how should this be done?
-
-*** in-thread approach, with particles i,j?
- 1) Look at whether i*j is in an array? --> potentially very memory intensive; could run out of integers!
- 2) Look at find i,j in array? --> potentially slow
- 3) Look at
-*** post-pairlist filtering step
-
-
+** Fix wrapping routines (trouble with machine precision
** make units and kT more uniform (references to 295 K are hardcoded into RB code and Temperature code)
*** I would prefer to specify temperature in K, rather than specifying a scaling factor to be applied to all energies to make them kT
diff --git a/useful.cu b/useful.cu
index 9e1a118..b729cac 100644
--- a/useful.cu
+++ b/useful.cu
@@ -410,54 +410,6 @@ Matrix3::Matrix3(const float* d) {
setIsDiag();
}
-const Matrix3 Matrix3::operator*(float s) const {
- Matrix3 m;
- m.exx = s*exx;
- m.exy = s*exy;
- m.exz = s*exz;
- m.eyx = s*eyx;
- m.eyy = s*eyy;
- m.eyz = s*eyz;
- m.ezx = s*ezx;
- m.ezy = s*ezy;
- m.ezz = s*ezz;
- m.isDiag = isDiag;
- return m;
-}
-
-
-const Matrix3 Matrix3::operator*(const Matrix3& m) const {
- Matrix3 ret;
- ret.exx = exx*m.exx + exy*m.eyx + exz*m.ezx;
- ret.eyx = eyx*m.exx + eyy*m.eyx + eyz*m.ezx;
- ret.ezx = ezx*m.exx + ezy*m.eyx + ezz*m.ezx;
-
- ret.exy = exx*m.exy + exy*m.eyy + exz*m.ezy;
- ret.eyy = eyx*m.exy + eyy*m.eyy + eyz*m.ezy;
- ret.ezy = ezx*m.exy + ezy*m.eyy + ezz*m.ezy;
-
- ret.exz = exx*m.exz + exy*m.eyz + exz*m.ezz;
- ret.eyz = eyx*m.exz + eyy*m.eyz + eyz*m.ezz;
- ret.ezz = ezx*m.exz + ezy*m.eyz + ezz*m.ezz;
- ret.setIsDiag();
- return ret;
-}
-
-const Matrix3 Matrix3::operator-() const {
- Matrix3 m;
- m.exx = -exx;
- m.exy = -exy;
- m.exz = -exz;
- m.eyx = -eyx;
- m.eyy = -eyy;
- m.eyz = -eyz;
- m.ezx = -ezx;
- m.ezy = -ezy;
- m.ezz = -ezz;
- m.isDiag = isDiag;
- return m;
-}
-
Matrix3 Matrix3::inverse() const {
Matrix3 m;
if (isDiag) {
@@ -483,11 +435,6 @@ float Matrix3::det() const {
exx*(eyy*ezz-eyz*ezy) - exy*(eyx*ezz-eyz*ezx) + exz*(eyx*ezy-eyy*ezx);
}
-
-Vector3 Matrix3::ex() const { return Vector3(exx,eyx,ezx); }
-Vector3 Matrix3::ey() const { return Vector3(exy,eyy,ezy); }
-Vector3 Matrix3::ez() const { return Vector3(exz,eyz,ezz); }
-
String Matrix3::toString() const {
char s[128];
sprintf(s, "%2.8f %2.8f %2.8f\n%2.8f %2.8f %2.8f\n%2.8f %2.8f %2.8f",
@@ -502,31 +449,6 @@ String Matrix3::toString1() const {
return String(s);
}
-Matrix3 operator*(float s, Matrix3 m) {
- m.exx *= s;
- m.exy *= s;
- m.exz *= s;
- m.eyx *= s;
- m.eyy *= s;
- m.eyz *= s;
- m.ezx *= s;
- m.ezy *= s;
- m.ezz *= s;
- return m;
-}
-
-Matrix3 operator/(Matrix3 m, float s) {
- m.exx /= s;
- m.exy /= s;
- m.exz /= s;
- m.eyx /= s;
- m.eyy /= s;
- m.eyz /= s;
- m.ezx /= s;
- m.ezy /= s;
- m.ezz /= s;
- return m;
-}
// class IndexList
// A growable list of integers.
diff --git a/useful.h b/useful.h
index 75a5d78..02f0e52 100644
--- a/useful.h
+++ b/useful.h
@@ -26,6 +26,7 @@ bool isInt(char c);
int firstSpace(const char* s, int max);
+
/*class int2 {
public:
int2(int x, int y) : x(x), y(y) {}
@@ -235,6 +236,10 @@ public:
return ret;
}
+ HOST DEVICE inline void print() {
+ // printf("%0.3f %0.3f %0.3f\n", x,y,z);
+ printf("%0.12f %0.12f %0.12f\n", x,y,z);
+ }
String toString() const;
float x, y, z;
@@ -257,6 +262,9 @@ HOST DEVICE inline Vector3 operator/(Vector3 v, float s) {
// class Matrix3
// Operations on 3D float matrices
class Matrix3 {
+ friend class TrajectoryWriter;
+ friend class BaseGrid;
+ friend class RigidBodyController; /* for trajectory writing */
public:
HOST DEVICE inline Matrix3() {}
HOST DEVICE Matrix3(float s);
@@ -265,13 +273,57 @@ public:
Matrix3(const Vector3& ex, const Vector3& ey, const Vector3& ez);
Matrix3(const float* d);
- const Matrix3 operator*(float s) const;
- HOST DEVICE const Vector3 operator*(const Vector3& v) const { return this->transform(v); }
+ // Operators
+ HOST DEVICE inline const Matrix3 operator*(float s) const {
+ Matrix3 m;
+ m.exx = s*exx; m.exy = s*exy; m.exz = s*exz;
+ m.eyx = s*eyx; m.eyy = s*eyy; m.eyz = s*eyz;
+ m.ezx = s*ezx; m.ezy = s*ezy; m.ezz = s*ezz;
+ m.isDiag = isDiag;
+ return m;
+ }
+ HOST DEVICE friend inline Matrix3 operator*(float s, Matrix3 m) { return m*s; }
+ HOST DEVICE friend inline Matrix3 operator/(Matrix3 m, float s) {
+ m.exx /= s; m.exy /= s; m.exz /= s;
+ m.eyx /= s; m.eyy /= s; m.eyz /= s;
+ m.ezx /= s; m.ezy /= s; m.ezz /= s;
+ return m;
+ }
- const Matrix3 operator*(const Matrix3& m) const;
+ HOST DEVICE inline const Vector3 operator*(const Vector3& v) const { return this->transform(v); }
- const Matrix3 operator-() const;
+ HOST DEVICE inline Matrix3 operator*(const Matrix3& m) const {
+ Matrix3 ret;
+ ret.exx = exx*m.exx + exy*m.eyx + exz*m.ezx;
+ ret.eyx = eyx*m.exx + eyy*m.eyx + eyz*m.ezx;
+ ret.ezx = ezx*m.exx + ezy*m.eyx + ezz*m.ezx;
+
+ ret.exy = exx*m.exy + exy*m.eyy + exz*m.ezy;
+ ret.eyy = eyx*m.exy + eyy*m.eyy + eyz*m.ezy;
+ ret.ezy = ezx*m.exy + ezy*m.eyy + ezz*m.ezy;
+
+ ret.exz = exx*m.exz + exy*m.eyz + exz*m.ezz;
+ ret.eyz = eyx*m.exz + eyy*m.eyz + eyz*m.ezz;
+ ret.ezz = ezx*m.exz + ezy*m.eyz + ezz*m.ezz;
+ ret.setIsDiag();
+ return ret;
+ }
+
+ HOST DEVICE inline Matrix3 operator-() const {
+ Matrix3 m;
+ m.exx = -exx;
+ m.exy = -exy;
+ m.exz = -exz;
+ m.eyx = -eyx;
+ m.eyy = -eyy;
+ m.eyz = -eyz;
+ m.ezx = -ezx;
+ m.ezy = -ezy;
+ m.ezz = -ezz;
+ m.isDiag = isDiag;
+ return m;
+ }
HOST DEVICE inline Matrix3 transpose() const {
Matrix3 m;
@@ -288,11 +340,6 @@ public:
return m;
}
- HOST DEVICE
- Matrix3 inverse() const;
-
- float det() const;
-
HOST DEVICE inline Vector3 transform(const Vector3& v) const {
Vector3 w;
if (isDiag) {
@@ -307,7 +354,6 @@ public:
return w;
}
-
HOST DEVICE inline Matrix3 transform(const Matrix3& m) const {
Matrix3 ret;
ret.exx = exx*m.exx + exy*m.eyx + exz*m.ezx;
@@ -325,6 +371,12 @@ public:
return ret;
}
+
+ HOST DEVICE
+ Matrix3 inverse() const;
+
+ float det() const;
+
HOST DEVICE void setIsDiag() {
isDiag = (exy == 0 && exz == 0 &&
eyx == 0 && eyz == 0 &&
@@ -333,14 +385,17 @@ public:
- Vector3 ex() const;
- Vector3 ey() const;
- Vector3 ez() const;
-
+ HOST DEVICE inline Vector3 ex() const { return Vector3(exx,eyx,ezx); }
+ HOST DEVICE inline Vector3 ey() const { return Vector3(exy,eyy,ezy); }
+ HOST DEVICE inline Vector3 ez() const { return Vector3(exz,eyz,ezz); }
String toString() const;
String toString1() const;
+ HOST DEVICE inline bool isDiagonal() const { return isDiag; }
+
+
+private:
float exx, exy, exz;
float eyx, eyy, eyz;
float ezx, ezy, ezz;
@@ -348,9 +403,6 @@ public:
};
-Matrix3 operator*(float s, Matrix3 m);
-Matrix3 operator/(Matrix3 m, float s);
-
// class IndexList
// A growable list of integers.
class IndexList {
--
GitLab