diff --git a/ComputeForce.cu b/ComputeForce.cu
index d1cfc0bea80ac19cd8c6a06d54341ab068e01c9b..7dce0591d060686c9d6cd8bc5f68da3901f6ad3a 100644
--- a/ComputeForce.cu
+++ b/ComputeForce.cu
@@ -590,6 +590,13 @@ float ComputeForce::computeTabulated(Vector3* force, Vector3* pos, int* type,
computeAngles<<<numBlocks, numThreads>>>(force, pos, angles,
tableAngle_d, numAngles, num, sys_d, energies_d, get_energy);
+/***************************************** call computeTabulatedBonds *****************************************/
+ if(bondMap != NULL && tableBond_d != NULL)
+ {
+ computeTabulatedBonds <<<numBlocks, numThreads>>> ( force, pos, num, numParts, sys_d, bonds, bondMap, numBonds, numReplicas, energies_d, get_energy, tableBond);
+
+ }
+/***************************************** end *****************************************/
computeDihedrals<<<numBlocks, numThreads>>>(force, pos, dihedrals,
tableDihedral_d, numDihedrals, num, sys_d, energies_d, get_energy);
diff --git a/ComputeForce.cuh b/ComputeForce.cuh
index 715df18a342d6b82566e3ff11c94628671f191f8..43628781aa52ae4805f0407c76a5aff85d4bd286 100644
--- a/ComputeForce.cuh
+++ b/ComputeForce.cuh
@@ -299,21 +299,17 @@ void createPairlists(Vector3* __restrict__ pos, int num, int numReplicas,
const int split = 32; /* numblocks should be divisible by split */
/* const int blocksPerCell = gridDim.x/split; */
- const CellDecomposition::cell_t* __restrict__ pairs = decomp->getCells();
+ const CellDecomposition::cell_t* __restrict__ cellInfo = decomp->getCells();
for (int cID = 0 + (blockIdx.x % split); cID < nCells; cID += split) {
// for (int cID = blockIdx.x/blocksPerCell; cID < nCells; cID += split ) {
for (int repID = 0; repID < numReplicas; repID++) {
const CellDecomposition::range_t rangeI = decomp->getRange(cID, repID);
for (int ci = rangeI.first + blockIdx.x/split; ci < rangeI.last; ci += gridDim.x/split) {
- /* for (int ci = rangeI.first + (blockIdx.x % blocksPerCell); */
- /* ci < rangeI.last; */
- /* ci += blocksPerCell) { */
-
// ai - index of the particle in the original, unsorted array
- const int ai = pairs[ci].particle;
+ const int ai = cellInfo[ci].particle;
// const CellDecomposition::cell_t celli = decomp->getCellForParticle(ai);
- const CellDecomposition::cell_t celli = pairs[ci];
+ const CellDecomposition::cell_t celli = cellInfo[ci];
// Vector3 posi = pos[ai];
for (int x = -1; x <= 1; ++x) {
@@ -325,7 +321,7 @@ void createPairlists(Vector3* __restrict__ pos, int num, int numReplicas,
const CellDecomposition::range_t range = decomp->getRange(nID, repID);
for (int n = range.first + tid; n < range.last; n+=blockDim.x) {
- const int aj = pairs[n].particle;
+ const int aj = cellInfo[n].particle;
if (aj <= ai) continue;
// skip ones that are too far away
@@ -441,7 +437,7 @@ __global__ void computeTabulatedKernel(
const Bond* __restrict__ bonds, const int2* __restrict__ bondMap, int numBonds,
float cutoff2, const int* __restrict__ g_numPairs,
const int2* __restrict__ g_pair, const int* __restrict__ g_pairTabPotType) {
-
+//remove int* type. remove bond references
const int numPairs = *g_numPairs;
for (int i = threadIdx.x+blockIdx.x*blockDim.x; i < numPairs; i+=blockDim.x*gridDim.x) {
@@ -648,6 +644,76 @@ void computeAngles(Vector3 force[], Vector3 pos[],
}
}
+/***************************************** computeBonds *****************************************/
+__global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
+ int numParts, BaseGrid* sys, Bond* bonds,
+ int2* bondMap, int numBonds, int numReplicas,
+ float* g_energies, bool get_energy, TabulatedPotential** tableBond)
+{
+ // Thread's unique ID.
+ int i = blockIdx.x * blockDim.x + threadIdx.x; // first particle ID
+
+ // Loop over ALL bonds in ALL replicas
+ if(i < num * numReplicas)
+ {
+ // Initialize interaction energy (per particle)
+ float energy_local = 0.0f;
+
+ const int repID = i / num;
+
+ // BONDS
+ /* Each particle may have a varying number of bonds
+ * bondMap is an array with one element for each particle
+ * which keeps track of where a particle's bonds are stored
+ * in the bonds array.
+ * bondMap[i].x is the index in the bonds array where the ith particle's bonds begin
+ * bondMap[i].y is the index in the bonds array where the ith particle's bonds end
+ */
+
+ // Get bond_start and bond_end from the bondMap
+ const int bond_start = bondMap[i - repID * num].x;
+ const int bond_end = bondMap[i - repID * num].y;
+
+ // Initialize force_local - force on a particle (i)
+ Vector3 force_local(0.0f);
+
+ for (int currBond = bond_start; currBond < bond_end; ++currBond)
+ {
+ // currBond is the index in the bonds array that we should look at next
+ // currBond is initialized to bond_start because that is the first index of the
+ // bonds array where this particle's bonds are stored
+
+ int j = bonds[currBond].ind2; // other particle ID
+
+ // Particle's type and position
+ Vector3 posi = pos[i];
+
+ // Find the distance between particles i and j,
+ // wrapping this value if necessary
+ const Vector3 dr = sys->wrapDiff(pos[j] - posi);
+
+ // Calculate the force on the particle from the bond
+ // If the user has specified the REPLACE option for this bond,
+ // then overwrite the force we calculated from the regular
+ // tabulated potential
+ // If the user has specified the ADD option, then add the bond
+ // force to the tabulated potential value
+ EnergyForce ft = tableBond[ bonds[currBond].tabFileIndex ]->compute(dr);
+
+ force_local += ft.f;
+
+ if (get_energy && j > i)
+ energy_local += ft.e;
+ }
+
+ force[i] = force_local;
+
+ if (get_energy)
+ g_energies[i] = energy_local;
+ }
+}
+/***************************************** end *****************************************/
+
__global__
void computeDihedrals(Vector3 force[], Vector3 pos[],
Dihedral dihedrals[],
diff --git a/Configuration.cpp b/Configuration.cpp
index c74fac5c527d1fb3c397606def72c684be41e2a4..492d3968d98a7102369b19c9dc6d9e5ee2f98348 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -1030,6 +1030,11 @@ void Configuration::readBonds() {
int ind1 = atoi(tokenList[2].val());
int ind2 = atoi(tokenList[3].val());
String file_name = tokenList[4];
+
+ if (ind1 == ind2) {
+ printf("WARNING: Invalid bond file line: %s\n", line);
+ continue;
+ }
// If we don't have enough room in our bond array, we need to expand it.
if (numBonds >= capacity) {
@@ -1051,6 +1056,15 @@ void Configuration::readBonds() {
}
// Add the bond to the bond array
// We must add it twice: Once for (ind1, ind2) and once for (ind2, ind1)
+
+ // RBTODO: add ind1/2 to exclusion list here iff op == REPLACE
+ /*if (op == Bond::REPLACE)
+ if( (int)(op) == 1)
+ {
+ printf("WARNING: Bond exclusions not implemented\n");
+ continue;
+ }*/
+
Bond* b = new Bond(op, ind1, ind2, file_name);
bonds[numBonds++] = *b;
b = new Bond(op, ind2, ind1, file_name);
diff --git a/RigidBodyController.cu b/RigidBodyController.cu
index 85e75561630f394a58ff6d4dd53bb15fb9a1d9fd..9db8d2d7261d19d5b9a45a868068e662c65208bf 100644
--- a/RigidBodyController.cu
+++ b/RigidBodyController.cu
@@ -432,7 +432,6 @@ void RigidBodyController::copyGridsToDevice() {
// RigidBodyGrid* g_d = rb.rawDensityGrids_d[gid]; // convenience
int len = g->getSize();
float* tmpData;
- // tmpData = new float[len];
size_t sz = sizeof(RigidBodyGrid);
gpuErrchk(cudaMalloc((void **) &(rb.rawDensityGrids_d[gid]), sz));
@@ -450,9 +449,6 @@ void RigidBodyController::copyGridsToDevice() {
gpuErrchk(cudaMemcpy( tmpData, rb.rawDensityGrids[gid].val, sz, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy( &(rb.rawDensityGrids_d[gid]->val), &tmpData,
sizeof(float*), cudaMemcpyHostToDevice));
-
- // RBTODO: why can't this be deleted?
- // delete[] tmpData;
}
}
diff --git a/makefile b/makefile
index 5c7ecc0dd31b86145a66f325737fa7f28456ade1..381a292dc6affd336b14fdd651389e9a25c6a9ed 100644
--- a/makefile
+++ b/makefile
@@ -12,7 +12,7 @@ include ./findcudalib.mk
INCLUDE = $(CUDA_PATH)/include
-# DEBUG = -g -O0
+DEBUG = -g -O0
CC_FLAGS = -Wall -Wno-write-strings -I$(INCLUDE) $(DEBUG) -std=c++0x -pedantic -fopenmp
# NV_FLAGS = --maxrregcount 63 -Xptxas -v # -v,-abi=no
NV_FLAGS = -Xptxas -v # -v,-abi=no
@@ -51,7 +51,8 @@ CODE_20 := -arch=sm_20
# NV_FLAGS += $(CODE_10) $(CODE_12) $(CODE_20) $(CODE_30) $(CODE_35)
# NV_FLAGS += -arch=sm_35
-SM=52
+# SM=52
+SM ?= 35
NV_FLAGS += -gencode arch=compute_$(SM),code=compute_$(SM)