From 64468c9171e7e0acb947560eb157f41e877e8936 Mon Sep 17 00:00:00 2001
From: Han-Yi Chou <hchou10@illinois.edu>
Date: Sun, 7 Jan 2018 06:08:36 -0600
Subject: [PATCH] add some testing directories and add feature supporting
multiple grid files
Conflicts:
src/Configuration.cpp
src/GrandBrownTown.cu
---
src/BaseGrid.cu | 2 +-
src/BaseGrid.h | 2 +
src/BrownianParticleType.cpp | 33 +++--
src/BrownianParticleType.h | 14 ++-
src/Configuration.cpp | 225 +++++++++++++++++++++++++++--------
src/Configuration.h | 7 +-
src/GrandBrownTown.cu | 4 +-
src/GrandBrownTown.cuh | 69 +++++++----
src/OverlordGrid.h | 2 +-
src/TabulatedMethods.cuh | 2 +-
10 files changed, 270 insertions(+), 90 deletions(-)
diff --git a/src/BaseGrid.cu b/src/BaseGrid.cu
index e70cbf1..75ccfbe 100644
--- a/src/BaseGrid.cu
+++ b/src/BaseGrid.cu
@@ -146,7 +146,7 @@ BaseGrid BaseGrid::mult(const BaseGrid& g) {
}
BaseGrid& BaseGrid::operator=(const BaseGrid& g) {
- delete[] val;
+ if(val != NULL) delete[] val;
val = NULL;
nx = g.nx;
ny = g.ny;
diff --git a/src/BaseGrid.h b/src/BaseGrid.h
index 01419a0..5523284 100644
--- a/src/BaseGrid.h
+++ b/src/BaseGrid.h
@@ -37,6 +37,8 @@ class ForceEnergy {
public:
DEVICE ForceEnergy(Vector3 &f, float &e) :
f(f), e(e) {};
+ DEVICE ForceEnergy(float f, float e) :
+ f(f), e(e) {};
Vector3 f;
float e;
};
diff --git a/src/BrownianParticleType.cpp b/src/BrownianParticleType.cpp
index acec0c0..0b56413 100644
--- a/src/BrownianParticleType.cpp
+++ b/src/BrownianParticleType.cpp
@@ -4,15 +4,16 @@
// BrownianParticleType Implementations //
//////////////////////////////////////////
void BrownianParticleType::clear() {
- if (pmf != NULL) delete pmf;
+ if (pmf != NULL) delete [] pmf;
if (diffusionGrid != NULL) delete diffusionGrid;
if (forceXGrid != NULL) delete forceXGrid;
if (forceYGrid != NULL) delete forceYGrid;
if (forceZGrid != NULL) delete forceZGrid;
if (reservoir != NULL) delete reservoir;
+ if (meanPmf != NULL) delete [] meanPmf;
pmf = NULL, diffusionGrid = NULL;
forceXGrid = NULL, forceYGrid = NULL, forceZGrid = NULL;
- reservoir = NULL;
+ reservoir = NULL; meanPmf = NULL;
}
void BrownianParticleType::copy(const BrownianParticleType& src) {
@@ -23,10 +24,23 @@ void BrownianParticleType::copy(const BrownianParticleType& src) {
charge = src.charge;
radius = src.radius;
eps = src.eps;
+ pmf = src.pmf;
meanPmf = src.meanPmf;
+ numPartGridFiles = src.numPartGridFiles;
//Han-Yi Chou
transDamping = src.transDamping;
mu = src.mu;
+ diffusionGrid = NULL;
+ forceXGrid = NULL, forceYGrid = NULL, forceZGrid = NULL;
+ reservoir = NULL;
+ //if (src.pmf != NULL) pmf = new BaseGrid(*src.pmf);
+ if (src.diffusionGrid != NULL) diffusionGrid = new BaseGrid(*src.diffusionGrid);
+ if (src.forceXGrid != NULL) forceXGrid = new BaseGrid(*src.forceXGrid);
+ if (src.forceYGrid != NULL) forceYGrid = new BaseGrid(*src.forceYGrid);
+ if (src.forceZGrid != NULL) forceZGrid = new BaseGrid(*src.forceZGrid);
+ if (src.reservoir != NULL) reservoir = new Reservoir(*src.reservoir);
+
+ /*
pmf = NULL, diffusionGrid = NULL;
forceXGrid = NULL, forceYGrid = NULL, forceZGrid = NULL;
reservoir = NULL;
@@ -35,15 +49,18 @@ void BrownianParticleType::copy(const BrownianParticleType& src) {
if (src.forceXGrid != NULL) forceXGrid = new BaseGrid(*src.forceXGrid);
if (src.forceYGrid != NULL) forceYGrid = new BaseGrid(*src.forceYGrid);
if (src.forceZGrid != NULL) forceZGrid = new BaseGrid(*src.forceZGrid);
- if (src.reservoir != NULL) reservoir = new Reservoir(*src.reservoir);
+ if (src.reservoir != NULL) reservoir = new Reservoir(*src.reservoir);*/
}
BrownianParticleType& BrownianParticleType::operator=(const BrownianParticleType& src) {
- clear();
- copy(src);
+ if(&src != this)
+ {
+ clear();
+ copy(src);
+ }
return *this;
}
-
+/*
bool BrownianParticleType::crop(int x0, int y0, int z0,
int x1, int y1, int z1, bool keep_origin) {
bool success = true;
@@ -104,7 +121,8 @@ bool BrownianParticleType::crop(int x0, int y0, int z0,
return success;
}
-
+*/
+/*
///////////////////////////////////////
// TypeDecomposition Implementations //
///////////////////////////////////////
@@ -146,3 +164,4 @@ const BrownianParticleType* TypeDecomposition::at(size_t i) const {
}
return parts_[i];
}
+*/
diff --git a/src/BrownianParticleType.h b/src/BrownianParticleType.h
index da959ee..67e6346 100644
--- a/src/BrownianParticleType.h
+++ b/src/BrownianParticleType.h
@@ -29,9 +29,9 @@ class BrownianParticleType {
public:
BrownianParticleType(const String& name = "") :
name(name), num(0),
- diffusion(0.0f), radius(1.0f), charge(0.0f), eps(0.0f), meanPmf(0.0f),
+ diffusion(0.0f), radius(1.0f), charge(0.0f), eps(0.0f), meanPmf(NULL),
reservoir(NULL), pmf(NULL), diffusionGrid(NULL),
- forceXGrid(NULL), forceYGrid(NULL), forceZGrid(NULL) { }
+ forceXGrid(NULL), forceYGrid(NULL), forceZGrid(NULL), numPartGridFiles(-1) { }
BrownianParticleType(const BrownianParticleType& src) { copy(src); }
@@ -44,7 +44,7 @@ class BrownianParticleType {
// @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);
+ //bool crop(int x0, int y0, int z0, int x1, int y1, int z1, bool keep_origin);
public:
String name;
@@ -55,7 +55,9 @@ public:
float radius;
float charge;
float eps;
- float meanPmf;
+ //float meanPmf;
+ float *meanPmf;
+ int numPartGridFiles;
float mu; //for Nose-Hoover Langevin dynamics
Reservoir* reservoir;
@@ -66,7 +68,7 @@ public:
BaseGrid* forceZGrid;
};
-
+/*
// Spatially decomposes BrownianParticleTypes
class TypeDecomposition {
private:
@@ -91,5 +93,5 @@ class TypeDecomposition {
int num_cells() const { return num_cells_; }
int num_parts() const { return num_parts_; }
};
-
+*/
#endif
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index c33274b..8dfbb22 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -6,8 +6,6 @@
#include <cassert>
#include <stdlib.h> /* srand, rand */
#include <time.h> /* time */
-
-
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) {
if (code != cudaSuccess) {
@@ -15,7 +13,44 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t
if (abort) exit(code);
}
}
+namespace
+{
+ template<class T>
+ void convertString(const String& token, void* data)
+ {
+ exit(1);
+ }
+
+ template<>
+ void convertString<float>(const String& token, void* data)
+ {
+ float* tmp = (float*)data;
+ *tmp = atof(token);
+ }
+
+ template<>
+ void convertString<String>(const String& token, void* data)
+ {
+ String* tmp = (String*)data;
+ *tmp = token;
+ }
+ template<class T>
+ void stringToArray(String* str, int& size, T** array)
+ {
+ register int num;
+ String *token;
+ num = str->tokenCount();
+ size = num;
+ *array = new T[num];
+ token = new String[num];
+ str->tokenize(token);
+
+ for(int i = 0; i < num; ++i)
+ convertString<T>(token[i], (*array)+i);
+ delete [] token;
+ }
+}
Configuration::Configuration(const char* config_file, int simNum, bool debug) :
simNum(simNum) {
// Read the parameters.
@@ -213,33 +248,69 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
printf("\nFound %d particle types.\n", numParts);
// Load the potential grids.
printf("Loading the potential grids...\n");
- for (int i = 0; i < numParts; i++) {
- // Decide which type of grid is given.
- String map = partGridFile[i];
- int len = map.length();
- if (len >= 3 && map[len-3]=='.' && map[len-2]=='d' && map[len-1]=='x') {
- // A dx file. Load the old-fashioned way.
- part[i].pmf = new BaseGrid(map.val());
- if (partGridFileScale[i] != 1.0f) part[i].pmf->scale(partGridFileScale[i]);
-
- part[i].meanPmf = part[i].pmf->mean();
- printf("Loaded dx potential grid `%s'.\n", map.val());
- printf("Grid size %s.\n", part[i].pmf->getExtent().toString().val());
- } else if (len >= 4 && map[len-4]=='.' && map[len-3]=='d' && map[len-2]=='e' && map[len-1]=='f') {
- // A system definition file.
- String rootGrid = OverlordGrid::readDefFirst(map);
- OverlordGrid* over = new OverlordGrid(rootGrid.val());
- int count = over->readDef(map);
- printf("Loaded system def file `%s'.\n", map.val());
- printf("Found %d unique grids.\n", over->getUniqueGridNum());
- printf("Linked %d subgrids.\n", count);
-
- part[i].pmf = static_cast<BaseGrid*>(over);
- part[i].meanPmf = part[i].pmf->mean();
- } else {
- printf("WARNING: Unrecognized gridFile extension. Must be *.def or *.dx.\n");
- exit(-1);
- }
+ for (int i = 0; i < numParts; i++)
+ {
+ String map = partGridFile[i][0];
+ int len = map.length();
+
+ if (len >= 3 && map[len-3]=='.' && map[len-2]=='d' && map[len-1]=='x')
+ {
+ part[i].pmf = new BaseGrid[part[i].numPartGridFiles];
+ part[i].meanPmf = new float[part[i].numPartGridFiles];
+ for(int j = 0; j < part[i].numPartGridFiles; ++j)
+ {
+ map = partGridFile[i][j];
+ len = map.length();
+ if (!(len >= 3 && map[len-3]=='.' && map[len-2]=='d' && map[len-1]=='x'))
+ {
+ cout << "currently do not support different format " << endl;
+ exit(1);
+ }
+ // Decide which type of grid is given.
+ //String map = partGridFile[i];
+
+ // A dx file. Load the old-fashioned way.
+ //part[i].pmf[j] = new BaseGrid(map.val());
+ BaseGrid tmp(map.val());
+ part[i].pmf[j] = tmp;
+ if (partGridFileScale[i][j] != 1.0f)
+ part[i].pmf[j].scale(partGridFileScale[i][j]);
+
+ //part[i].meanPmf = part[i].pmf->mean();
+ part[i].meanPmf[j] = part[i].pmf[j].mean();
+ printf("Loaded dx potential grid `%s'.\n", map.val());
+ printf("Grid size %s.\n", part[i].pmf[j].getExtent().toString().val());
+ }
+ }
+ else if (len >= 4 && map[len-4]=='.' && map[len-3]=='d' && map[len-2]=='e' && map[len-1]=='f')
+ {
+ OverlordGrid* over = new OverlordGrid[part[i].numPartGridFiles];
+ part[i].meanPmf = new float[part[i].numPartGridFiles];
+ for(int j = 0; j < part[i].numPartGridFiles; ++j)
+ {
+ map = partGridFile[i][j];
+ len = map.length();
+ if (!(len >= 4 && map[len-4]=='.' && map[len-3]=='d' && map[len-2]=='e' && map[len-1]=='f'))
+ {
+ cout << "currently do not support different format " << endl;
+ exit(1);
+ }
+
+ String rootGrid = OverlordGrid::readDefFirst(map);
+ over[j] = OverlordGrid(rootGrid.val());
+ int count = over->readDef(map);
+ printf("Loaded system def file `%s'.\n", map.val());
+ printf("Found %d unique grids.\n", over->getUniqueGridNum());
+ printf("Linked %d subgrids.\n", count);
+ part[i].meanPmf[j] = part[i].pmf[j].mean();
+ }
+ part[i].pmf = static_cast<BaseGrid*>(over);
+ }
+ else
+ {
+ printf("WARNING: Unrecognized gridFile extension. Must be *.def or *.dx.\n");
+ exit(-1);
+ }
if (partForceXGridFile[i].length() != 0) {
part[i].forceXGrid = new BaseGrid(partForceXGridFile[i].val());
@@ -277,6 +348,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
// part[i].diffusionGrid->write(outFile, comment);
}
}
+
}
// Load reservoir files if any
@@ -449,8 +521,24 @@ Configuration::~Configuration() {
// Particle parameters
delete[] part;
- delete[] partGridFile;
- delete[] partGridFileScale;
+ //delete[] partGridFile;
+ //delete[] partGridFileScale;
+ for(int i = 0; i < numParts; ++i)
+ {
+ if(partGridFile[i] != NULL)
+ {
+ delete[] partGridFile[i];
+ partGridFile[i] = NULL;
+ }
+ if(partGridFileScale[i] != NULL)
+ {
+ delete[] partGridFileScale[i];
+ partGridFileScale[i] = NULL;
+ }
+ }
+ delete partGridFile;
+ delete partGridFileScale;
+ //delete numPartGridFiles;
delete[] partForceXGridFile;
delete[] partForceYGridFile;
delete[] partForceZGridFile;
@@ -503,22 +591,35 @@ void Configuration::copyToCUDA() {
gpuErrchk(cudaMalloc(&part_d, sizeof(BrownianParticleType*) * numParts));
// TODO: The above line fails when there is not enough memory. If it fails, stop.
- for (int i = 0; i < numParts; i++) {
+ for (int i = 0; i < numParts; i++)
+ {
BaseGrid *pmf = NULL, *diffusionGrid = NULL;
BrownianParticleType *b = new BrownianParticleType(part[i]);
// Copy PMF
- if (part[i].pmf != NULL) {
- float *val = NULL;
- size_t sz = sizeof(float) * part[i].pmf->getSize();
- gpuErrchk(cudaMalloc(&pmf, sizeof(BaseGrid)));
- gpuErrchk(cudaMalloc(&val, sz));
- gpuErrchk(cudaMemcpyAsync(val, part[i].pmf->val, sz, cudaMemcpyHostToDevice));
- BaseGrid *pmf_h = new BaseGrid(*part[i].pmf);
- pmf_h->val = val;
- gpuErrchk(cudaMemcpy(pmf, pmf_h, sizeof(BaseGrid), cudaMemcpyHostToDevice));
- pmf_h->val = NULL;
+ if (part[i].pmf != NULL)
+ {
+ float *tmp;
+ gpuErrchk(cudaMalloc(&pmf, sizeof(BaseGrid)*part[i].numPartGridFiles));
+ gpuErrchk(cudaMalloc(&tmp, sizeof(float)*part[i].numPartGridFiles));
+ gpuErrchk(cudaMemcpy(tmp, part[i].meanPmf, sizeof(float)*part[i].numPartGridFiles,
+ cudaMemcpyHostToDevice));
+ b->meanPmf = tmp;
+
+ for(int j = 0; j < part[i].numPartGridFiles; ++j)
+ {
+ float *val = NULL;
+ size_t sz = sizeof(float) * part[i].pmf[j].getSize();
+ //gpuErrchk(cudaMalloc(pmf, sizeof(BaseGrid)));
+ gpuErrchk(cudaMalloc(&val, sz));
+ gpuErrchk(cudaMemcpyAsync(val, part[i].pmf[j].val, sz, cudaMemcpyHostToDevice));
+ BaseGrid *pmf_h = new BaseGrid(part[i].pmf[j]);
+ pmf_h->val = val;
+ gpuErrchk(cudaMemcpy(pmf+j, pmf_h, sizeof(BaseGrid), cudaMemcpyHostToDevice));
+ pmf_h->val = NULL;
+ }
+
}
-
+ b->pmf = pmf;
// Copy the diffusion grid
if (part[i].diffusionGrid != NULL) {
float *val = NULL;
@@ -533,7 +634,7 @@ void Configuration::copyToCUDA() {
diffusionGrid_h->val = NULL;
}
- b->pmf = pmf;
+ //b->pmf = pmf;
b->diffusionGrid = diffusionGrid;
gpuErrchk(cudaMalloc(&part_addr[i], sizeof(BrownianParticleType)));
gpuErrchk(cudaMemcpyAsync(part_addr[i], b, sizeof(BrownianParticleType),
@@ -688,8 +789,13 @@ int Configuration::readParameters(const char * config_file) {
// Allocate the particle variables.
part = new BrownianParticleType[numParts];
- partGridFile = new String[numParts];
- partGridFileScale = new float[numParts];
+ //partGridFile = new String[numParts];
+ //partGridFileScale = new float[numParts];
+ partGridFile = new String*[numParts];
+ //partGridFileScale = new float[numParts];
+ partGridFileScale = new float*[numParts];
+ //int numPartGridFiles = new int[numParts];
+
partForceXGridFile = new String[numParts];
partForceYGridFile = new String[numParts];
partForceZGridFile = new String[numParts];
@@ -706,9 +812,19 @@ int Configuration::readParameters(const char * config_file) {
rigidBody = new RigidBodyType[numRigidTypes];
// Set a default
+ /*
for (int i = 0; i < numParts; ++i) {
partGridFileScale[i] = 1.0f;
- }
+ }*/
+
+ for(int i = 0; i < numParts; ++i)
+ {
+ partGridFile[i] = NULL;
+ partGridFileScale[i] = NULL;
+ //part[i].numPartGridFiles = -1;
+ }
+ //for(int i = 0; i < numParts; ++i)
+ // cout << part[i].numPartGridFiles << endl;
int btfcap = 10;
bondTableFile = new String[btfcap];
@@ -758,7 +874,6 @@ int Configuration::readParameters(const char * config_file) {
inputCoordinates = value;
else if (param == String("restartCoordinates"))
restartCoordinates = value;
-
//Han-Yi Chou
else if (param == String("inputMomentum"))
inputMomentum = value;
@@ -919,7 +1034,9 @@ int Configuration::readParameters(const char * config_file) {
readRestraintsFromFile = true;
}
} else if (param == String("gridFileScale")) {
- partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
+ //partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
+ stringToArray<float>(&value, part[currPart].numPartGridFiles,
+ &partGridFileScale[currPart]);
} else if (param == String("rigidBodyPotential")) {
partRigidBodyGrid[currPart].push_back(value);
}
@@ -985,7 +1102,15 @@ int Configuration::readParameters(const char * config_file) {
}
else if (param == String("gridFile")) {
if (currPartClass == partClassPart)
- partGridFile[currPart] = value;
+ {
+ printf("The grid file name %s\n", value.val());
+ //partGridFile[currPart] = value;
+ cout << part[currPart].numPartGridFiles << endl;
+ stringToArray<String>(&value, part[currPart].numPartGridFiles,
+ &partGridFile[currPart]);
+ for(int i = 0; i < part[currPart].numPartGridFiles; ++i)
+ printf("%s ", partGridFile[currPart]->val());
+ }
else if (currPartClass == partClassRB)
rigidBody[currRB].addPMF(value);
}
diff --git a/src/Configuration.h b/src/Configuration.h
index 19404a4..6bd6996 100644
--- a/src/Configuration.h
+++ b/src/Configuration.h
@@ -201,8 +201,11 @@ public:
bool readAnglesFromFile;
bool readDihedralsFromFile;
bool readRestraintsFromFile;
- String* partGridFile;
- float* partGridFileScale;
+ //String* partGridFile;
+ String **partGridFile;
+ //float* partGridFileScale;
+ float **partGridFileScale;
+ //int *numPartGridFiles;
std::vector< std::vector<String> > partRigidBodyGrid;
String* partDiffusionGridFile;
String* partForceXGridFile;
diff --git a/src/GrandBrownTown.cu b/src/GrandBrownTown.cu
index c237838..26df336 100644
--- a/src/GrandBrownTown.cu
+++ b/src/GrandBrownTown.cu
@@ -46,7 +46,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
if(particle_dynamic == String("Langevin") || particle_dynamic == String("NoseHooverLangevin"))
{
std::stringstream restart_file_p, out_momentum_prefix, out_force_prefix;
- restart_file_p << outArg << '.' << i << ".momentum.resart";
+ restart_file_p << outArg << '.' << i << ".momentum.restart";
out_momentum_prefix << outArg << ".momentum." << i;
//out_force_prefix << outArg << ".force." << i;
@@ -1437,7 +1437,7 @@ void GrandBrownTown::writeRestart(int repID) const
}
-
+/*the center is defined by the first pmf*/
void GrandBrownTown::initialCondCen() {
for (int i = 0; i < num; i++)
pos[i] = part[ type[i] ].pmf->getCenter();
diff --git a/src/GrandBrownTown.cuh b/src/GrandBrownTown.cuh
index 7c5741b..2ff3085 100644
--- a/src/GrandBrownTown.cuh
+++ b/src/GrandBrownTown.cuh
@@ -97,9 +97,11 @@ __global__ void updateKernelABOBA(Vector3* pos, Vector3* momentum, Vector3* __re
////The kernel is for Nose-Hoover Langevin dynamics
-__global__ void updateKernelNoseHooverLangevin(Vector3* pos, Vector3* momentum, float* random, Vector3* __restrict__ forceInternal,
- int type[], BrownianParticleType* part[], float kT, BaseGrid* kTGrid, float electricField,
- int tGridLength, float timestep, int num, BaseGrid* sys, Random* randoGen, int numReplicas)
+__global__ void
+updateKernelNoseHooverLangevin(Vector3* __restrict__ pos, Vector3* __restrict__ momentum, float* random,
+ Vector3* __restrict__ forceInternal, int type[], BrownianParticleType* part[],
+ float kT, BaseGrid* kTGrid, float electricField, int tGridLength, float timestep,
+ int num, BaseGrid* sys, Random* randoGen, int numReplicas)
{
const int idx = static_cast<int>(blockIdx.x * blockDim.x + threadIdx.x);
if (idx < num * numReplicas)
@@ -112,7 +114,13 @@ __global__ void updateKernelNoseHooverLangevin(Vector3* pos, Vector3* momentum,
const BrownianParticleType& pt = *part[t];
Vector3 forceExternal = Vector3(0.0f, 0.0f, pt.charge * electricField);
- ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(r0);
+ ForceEnergy fe(0.f, 0.f);
+ for(int i = 0; i < pt.numPartGridFiles; ++i)
+ {
+ ForceEnergy tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(r0);
+ fe.f += tmp.f;
+ }
+ //ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(r0);
#ifndef FORCEGRIDOFF
// Add a force defined via 3D FORCE maps (not 3D potential maps)
@@ -183,8 +191,9 @@ __global__ void updateKernelNoseHooverLangevin(Vector3* pos, Vector3* momentum,
//which is not possible in GPU code.
//Han-Yi Chou
__global__ void updateKernelBAOAB(Vector3* pos, Vector3* momentum, Vector3* __restrict__ forceInternal,
- int type[], BrownianParticleType* part[], float kT, BaseGrid* kTGrid, float electricField,
- int tGridLength, float timestep, int num, BaseGrid* sys, Random* randoGen, int numReplicas)
+ int type[], BrownianParticleType* part[], float kT, BaseGrid* kTGrid,
+ float electricField,int tGridLength, float timestep, int num, BaseGrid* sys,
+ Random* randoGen, int numReplicas)
{
const int idx = static_cast<int>(blockIdx.x * blockDim.x + threadIdx.x);
@@ -198,7 +207,13 @@ __global__ void updateKernelBAOAB(Vector3* pos, Vector3* momentum, Vector3* __re
const BrownianParticleType& pt = *part[t];
Vector3 forceExternal = Vector3(0.0f, 0.0f, pt.charge * electricField);
- ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(r0);
+ //ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(r0);
+ ForceEnergy fe(0.f, 0.f);
+ for(int i = 0; i < pt.numPartGridFiles; ++i)
+ {
+ ForceEnergy tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(r0);
+ fe.f += tmp.f;
+ }
#ifndef FORCEGRIDOFF
// Add a force defined via 3D FORCE maps (not 3D potential maps)
@@ -255,8 +270,9 @@ __global__ void updateKernelBAOAB(Vector3* pos, Vector3* momentum, Vector3* __re
//update momentum in the last step of BAOAB integrator for the Langevin dynamics. Han-Yi Chou
__global__ void LastUpdateKernelBAOAB(Vector3* pos,Vector3* momentum, Vector3* __restrict__ forceInternal,
- int type[], BrownianParticleType* part[], float kT, BaseGrid* kTGrid, float electricField,
- int tGridLength, float timestep, int num, BaseGrid* sys, Random* randoGen, int numReplicas)
+ int type[], BrownianParticleType* part[], float kT, BaseGrid* kTGrid,
+ float electricField, int tGridLength, float timestep, int num,
+ BaseGrid* sys, Random* randoGen, int numReplicas)
{
const int idx = static_cast<int>(blockIdx.x * blockDim.x + threadIdx.x);
@@ -269,7 +285,14 @@ __global__ void LastUpdateKernelBAOAB(Vector3* pos,Vector3* momentum, Vector3* _
const BrownianParticleType& pt = *part[t];
Vector3 forceExternal = Vector3(0.0f, 0.0f, pt.charge * electricField);
- ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(r0);
+ //ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(r0);
+ ForceEnergy fe(0.f, 0.f);
+ for(int i = 0; i < pt.numPartGridFiles; ++i)
+ {
+ ForceEnergy tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(r0);
+ fe.f += tmp.f;
+ }
+
#ifndef FORCEGRIDOFF
// Add a force defined via 3D FORCE maps (not 3D potential maps)
if (pt.forceXGrid != NULL) fe.f.x += pt.forceXGrid->interpolatePotentialLinearly(r0);
@@ -353,12 +376,11 @@ __global__ void LastUpdateKernelABOBA(Vector3* pos, Vector3* momentum, Vector3*
//Update kernel for Brownian dynamics
__global__
-void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
- int type[], BrownianParticleType* part[],
- float kT, BaseGrid* kTGrid,
- float electricField, int tGridLength,
- float timestep, int num, BaseGrid* sys,
- Random* randoGen, int numReplicas) {
+void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal, int type[],
+ BrownianParticleType* part[],float kT, BaseGrid* kTGrid, float electricField,
+ int tGridLength, float timestep, int num, BaseGrid* sys,
+ Random* randoGen, int numReplicas)
+{
// Calculate this thread's ID
const int idx = static_cast<int>(blockIdx.x * blockDim.x + threadIdx.x);
@@ -377,7 +399,13 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
// Compute PMF
// TODO: maybe make periodic and nonPeriodic versions of this kernel
- ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(p);
+ //ForceEnergy fe = pt.pmf->interpolateForceDLinearlyPeriodic(p);
+ ForceEnergy fe(0.f, 0.f);
+ for(int i = 0; i < pt.numPartGridFiles; ++i)
+ {
+ ForceEnergy tmp = pt.pmf[i].interpolateForceDLinearlyPeriodic(p);
+ fe.f += tmp.f;
+ }
#ifndef FORCEGRIDOFF
// Add a force defined via 3D FORCE maps (not 3D potential maps)
@@ -414,8 +442,8 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
Vector3 gridCenter = pt.diffusionGrid->origin +
pt.diffusionGrid->basis.transform( Vector3(0.5*pt.diffusionGrid->nx,
- 0.5*pt.diffusionGrid->ny,
- 0.5*pt.diffusionGrid->nz));
+ 0.5*pt.diffusionGrid->ny,
+ 0.5*pt.diffusionGrid->nz));
Vector3 p2 = p - gridCenter;
p2 = sys->wrapDiff( p2 ) + gridCenter;
/* p2 = sys->wrap( p2 ); */
@@ -430,7 +458,8 @@ void updateKernel(Vector3* pos, Vector3* __restrict__ forceInternal,
// printf("force: "); force.print();
// }
- Vector3 tmp = step(p, kTlocal, force, diffusion, -diffGrad, timestep, sys, randoGen, num * numReplicas);
+ Vector3 tmp = step(p, kTlocal, force, diffusion, -diffGrad, timestep, sys, randoGen,
+ num * numReplicas);
// assert( tmp.length() < 10000.0f );
pos[idx] = tmp;
diff --git a/src/OverlordGrid.h b/src/OverlordGrid.h
index 060f909..af2a415 100644
--- a/src/OverlordGrid.h
+++ b/src/OverlordGrid.h
@@ -22,7 +22,7 @@ public:
initSubgrids();
initUniqueGrids();
}
-
+ OverlordGrid() {};
/*OverlordGrid(const char* systemDefFile) : BaseGrid(readDefFirst(systemDefFile)) {
printf("size: %d\n", size);
diff --git a/src/TabulatedMethods.cuh b/src/TabulatedMethods.cuh
index bbd9b82..90a7839 100644
--- a/src/TabulatedMethods.cuh
+++ b/src/TabulatedMethods.cuh
@@ -43,7 +43,7 @@ __device__ inline void computeAngle(const TabulatedAnglePotential* __restrict__
// tableAngle[0] stores the potential at angle_step
// tableAngle[1] stores the potential at angle_step * 2, etc.
// 'home' is the index after which 'convertedAngle' would appear if it were stored in the table
- int home = int(floor(angle));
+ int home = int(floorf(angle));
home = (home >= a->size) ? (a->size)-1 : home;
//assert(home >= 0);
//assert(home+1 < a->size);
--
GitLab