diff --git a/ComputeForce.cu b/ComputeForce.cu
index 81e67264f50ca9dc9469ee735717b985f849bdd9..7f9207ab9452bea77de4c8495d16e2a34bce5e22 100644
--- a/ComputeForce.cu
+++ b/ComputeForce.cu
@@ -556,9 +556,13 @@ float ComputeForce::compute(Vector3 force[], Vector3 pos[], int type[], bool get
return energy;
}
+//MLog: added Bond* bondList to the list of passed in variables.
+/*float ComputeForce::computeTabulated(Vector3* force, Vector3* pos, int* type,
+ Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap,
+ Angle* angles, Dihedral* dihedrals, bool get_energy, Bond* bondList) {*/
float ComputeForce::computeTabulated(Vector3* force, Vector3* pos, int* type,
Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap,
- Angle* angles, Dihedral* dihedrals, bool get_energy) {
+ Angle* angles, Dihedral* dihedrals, bool get_energy, int3* bondList_d) {
float energy = 0.0f;
gridSize = (num * numReplicas) / NUM_THREADS + 1;
@@ -590,13 +594,12 @@ 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)
+ //Mlog: the commented function doesn't use bondList, uncomment for testing.
if(bondMap != NULL && tableBond_d != NULL)
{
- computeTabulatedBonds <<<numBlocks, numThreads>>> ( force, pos, num, numParts, sys_d, bonds, bondMap, numBonds, numReplicas, energies_d, get_energy, tableBond_d);
+ //computeTabulatedBonds <<<numBlocks, numThreads>>> ( force, pos, num, numParts, sys_d, bonds, bondMap, numBonds, numReplicas, energies_d, get_energy, tableBond_d);
+ computeTabulatedBonds <<<numBlocks, numThreads>>> ( force, pos, num, numParts, sys_d, bondList_d, (numBonds/2), numReplicas, energies_d, get_energy, tableBond_d);
}
-/***************************************** 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 debebd7ae376aef363e4de06fcd61760ce61de23..359b879551923b6e1214b902d4a6d4629b48f391 100644
--- a/ComputeForce.cuh
+++ b/ComputeForce.cuh
@@ -644,9 +644,9 @@ void computeAngles(Vector3 force[], Vector3 pos[],
}
}
-/***************************************** computeBonds *****************************************/
-__global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
- int numParts, BaseGrid* sys, Bond bonds[],
+//Mlog: the commented function doesn't use bondList, uncomment for testing.
+/*__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)
{
@@ -660,7 +660,7 @@ __global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
float energy_local = 0.0f;
const int repID = i / num; //TODO: ask why repID is even here, it's unnecessary in this kernel
-
+ */
// BONDS
/* Each particle may have a varying number of bonds
* bondMap is an array with one element for each particle
@@ -671,7 +671,7 @@ __global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
*/
// Get bond_start and bond_end from the bondMap
- const int bond_start = bondMap[i - repID * num].x;
+ /*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)
@@ -699,30 +699,52 @@ __global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
// If the user has specified the ADD option, then add the bond
// force to the tabulated potential value
- /* start testing */
- //cout << "tableBond points to " << tableBond << ".\n"
- //<< "tableBond has a size of " << sizeof(tableBond)<< ".\n";
- //Bond test_bond = bonds[currBond];
- //int test_num = test_bond.tabFileIndex;
- //TabulatedPotential *test_pot = tableBond[test_num];
- //EnergyForce test_energy = test_pot -> compute(dr);
- /* end testing */
-
force_local += tableBond[ bonds[currBond].tabFileIndex ]->computef(dr,dr.length2());
-
- /*EnergyForce ft = tableBond[ bonds[currBond].tabFileIndex ]->compute(dr);
- force_local += ft.f;
-
- if (get_energy && j > i)
- energy_local += ft.e;*/
}
atomicAdd( &force[i], force_local );
if (get_energy)
atomicAdd( &g_energies[i], energy_local);
}
+}*/
+
+__global__ void computeTabulatedBonds( Vector3* force, Vector3* pos, int num,
+ int numParts, BaseGrid* sys, int3* bondList_d,
+ int numBonds, int numReplicas, float* g_energies,
+ bool get_energy, TabulatedPotential** tableBond)
+{
+ // Thread's unique ID.
+ int currBond = blockIdx.x * blockDim.x + threadIdx.x; // first particle ID
+
+ // Loop over ALL bonds in ALL replicas
+ if( currBond < (numBonds * numReplicas) ) //numBonds should be divided by 2 in the kernel call
+ {
+ // Initialize interaction energy (per particle)
+ float energy_local = 0.0f;
+
+ int i = bondList_d[currBond].x;
+ int j = bondList_d[currBond].y;
+
+ // 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] - pos[i]);
+
+ Vector3 force_local = tableBond[ bondList_d[currBond].z ]->computef(dr,dr.length2());
+
+ atomicAdd( &force[i], force_local );
+ atomicAdd( &force[j], -force_local );
+
+ if (get_energy)
+ {
+ //TODO: clarification on energy computation needed, consider changing.
+ atomicAdd( &g_energies[i], energy_local);
+ //atomicAdd( &g_energies[j], energy_local);
+ }
+ }
}
-/***************************************** end *****************************************/
__global__
void computeDihedrals(Vector3 force[], Vector3 pos[],
diff --git a/ComputeForce.h b/ComputeForce.h
index ea1761d9fbec60137a02f36f33a3dc5fb054b041..d1d0b350bc6078d37ea0528a96e1e5027aac39ea 100644
--- a/ComputeForce.h
+++ b/ComputeForce.h
@@ -63,9 +63,13 @@ public:
float compute(Vector3* force, Vector3* pos, int* type, bool get_energy);
float computeFull(Vector3* force, Vector3* pos, int* type, bool get_energy);
+ //MLog: the commented function doesn't use bondList, uncomment for testing.
+ /*float computeTabulated(Vector3* force, Vector3* pos, int* type,
+ Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap,
+ Angle* angles, Dihedral* dihedrals, bool get_energy);*/
float computeTabulated(Vector3* force, Vector3* pos, int* type,
Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap,
- Angle* angles, Dihedral* dihedrals, bool get_energy);
+ Angle* angles, Dihedral* dihedrals, bool get_energy, int3* bondList_d);
float computeTabulatedFull(Vector3* force, Vector3* pos, int* type,
Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap,
Angle* angles, Dihedral* dihedrals, bool get_energy);
diff --git a/Configuration.cpp b/Configuration.cpp
index 492d3968d98a7102369b19c9dc6d9e5ee2f98348..75cb5eb25aad6ff5e1243dac2c7b95e4fe80ad18 100644
--- a/Configuration.cpp
+++ b/Configuration.cpp
@@ -375,6 +375,8 @@ Configuration::~Configuration() {
delete[] partForceZGridFile;
delete[] partDiffusionGridFile;
delete[] partReservoirFile;
+
+ // TODO: plug memory leaks
if (partsFromFile != NULL) delete[] partsFromFile;
if (bonds != NULL) delete[] bonds;
if (bondMap != NULL) delete[] bondMap;
diff --git a/GrandBrownTown.cu b/GrandBrownTown.cu
index e5700e10c9480cd6cc275632637e4c8c39794343..c2db8bf12592c179b8e446f15eb4aeea8a95234a 100644
--- a/GrandBrownTown.cu
+++ b/GrandBrownTown.cu
@@ -218,30 +218,30 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
internal->addDihedralPotential(dihedralTableFile[p].val(), p, dihedrals, dihedrals_d);
}
-/***************************************** init. bondList *****************************************/
- bondList = new int3[numBonds/2];
+ //Mlog: this is where we create the bondList.
+ bondList = new int3[ (numBonds / 2) * numReplicas ];
int j = 0;
-
- for(int i = 0; i < numBonds; ++i)
+
+ for(int k = 0 ; k < numReplicas; k++)
{
- if(bonds[i].ind1 < bonds[i].ind2)
+ for(int i = 0; i < numBonds; ++i)
{
- bondList[j] = make_int3( bonds[i].ind1, bonds[i].ind2, bonds[i].tabFileIndex );
- cout << "Displaying: bondList["<< j <<"].x = " << bondList[j].x << ".\n"
- << "Displaying: bondList["<< j <<"].y = " << bondList[j].y << ".\n"
- << "Displaying: bondList["<< j <<"].z = " << bondList[j].z << ".\n";
- j++;
+ if(bonds[i].ind1 < bonds[i].ind2)
+ {
+ bondList[j] = make_int3( (bonds[i].ind1 + k * num), (bonds[i].ind2 + k * num), bonds[i].tabFileIndex );
+ cout << "Displaying: bondList["<< j <<"].x = " << bondList[j].x << ".\n"
+ << "Displaying: bondList["<< j <<"].y = " << bondList[j].y << ".\n"
+ << "Displaying: bondList["<< j <<"].z = " << bondList[j].z << ".\n";
+ j++;
+ }
}
}
- //TODO: memCopy bondList here to device, pass in list to decompose and change the function signitures. Afterwards edit the function to make use of the list.
createBondList();
// TODO: copy data to device (not here)
// TODO: use bondList in computeTabulatedBonds kernel
-/***************************************** end *****************************************/
-
forceInternal = new Vector3[num * numReplicas];
if (fullLongRange != 0)
@@ -364,11 +364,10 @@ void GrandBrownTown::run() {
internal->decompose(pos_d, type_d);
//internal->updatePairlists(pos_d); // perhaps this way?
}
- energy = internal->computeTabulated(forceInternal_d, pos_d, type_d,
- bonds_d, bondMap_d,
- excludes_d, excludeMap_d,
- angles_d, dihedrals_d,
- get_energy);
+
+ //MLog: added Bond* bondList to the list of passed in variables.
+ /*energy = internal->computeTabulated(forceInternal_d, pos_d, type_d, bonds_d, bondMap_d, excludes_d, excludeMap_d, angles_d, dihedrals_d, get_energy);*/
+ energy = internal -> computeTabulated(forceInternal_d, pos_d, type_d, bonds_d, bondMap_d, excludes_d, excludeMap_d, angles_d, dihedrals_d, get_energy, bondList_d);
break;
default: // [ N^2 ] interactions, no cutoff | decompositions
energy =
@@ -988,11 +987,11 @@ void GrandBrownTown::copyToCUDA() {
void GrandBrownTown::createBondList()
{
- size_t size = (numBonds/2) * sizeof(int3);
- gpuErrchk( cudaMalloc( &bondList_d, sizeof(size) ) );
- gpuErrchk( cudaMemcpyAsync( bondList_d, bondList, sizeof(size), cudaMemcpyHostToDevice) );
+ size_t size = (numBonds / 2) * numReplicas * sizeof(int3);
+ gpuErrchk( cudaMalloc( &bondList_d, size ) );
+ gpuErrchk( cudaMemcpyAsync( bondList_d, bondList, size, cudaMemcpyHostToDevice) );
- for(int i = 0 ; i < numBonds/2 ; i++)
+ for(int i = 0 ; i < (numBonds / 2) * numReplicas ; i++)
{
cout << "Displaying: bondList_d["<< i <<"].x = " << bondList[i].x << ".\n"
<< "Displaying: bondList_d["<< i <<"].y = " << bondList[i].y << ".\n"
diff --git a/makefile b/makefile
index 881f152328b10c8f85d3309e2b181dcd2f5d1bf1..33b13510a78a85c70c141056297bc647edfa09c4 100644
--- a/makefile
+++ b/makefile
@@ -12,7 +12,7 @@ include ./findcudalib.mk
INCLUDE = $(CUDA_PATH)/include
-# DEBUG = -g -G -O0
+DEBUG = -g -G -O0
CC_FLAGS = -Wall -Wno-write-strings -I$(INCLUDE) $(DEBUG) -std=c++0x -pedantic
# NV_FLAGS = --maxrregcount 63 -Xptxas -v # -v,-abi=no
NV_FLAGS = -Xptxas -v # -v,-abi=no