diff --git a/README b/README
index 56081944f6054a7fd057a935da3c78298fcd6278..7a6ec3527e8fd6df97cfb3f32a862ece17c2008b 100644
--- a/README
+++ b/README
@@ -1,6 +1,6 @@
-/============================================\
-| Atomic Resolution Brownian Dynamics (ARBD) |
-\============================================/
+/===========================================================\
+| Atomic Resolution Brownian Dynamics (ARBD) - alpha Nov 16 |
+\===========================================================/
Brownian dynamics (BD) simulation is method for studying biomolecules, ions, and
nanomaterials that balances detail with computational efficiency.
@@ -12,8 +12,10 @@ potentials to be associated with rigid body particles that rotate and
translate to represent larger molecules. Most importantly, the code is
designed to run quickly on modern NVIDIA GPUs.
-Please be aware that ARBD is being actively developed, is offered as alpha
-software without warranty. Expect to run into the occasional bug.
+ARBD is a rewrite of the BrownianMover code, moving almost all computations to
+the GPU and enabling grid-specified particle models. Please be aware that ARBD
+is being actively developed and is offered without warranty.
+
/==============\
| Installation |
@@ -25,21 +27,33 @@ Note that ARBD was developed using CUDA-8.0 and targets NVIDIA GPUs featuring
6.0 compute capability. The code should still run using older NVIDIA hardware,
but no guarantees are being made.
+
/========\
| Citing |
\========/
-If you publish results obtained using ARBD, please cite the following manuscript:
+If you publish results obtained using ARBD, please cite the following manuscripts:
+
+"DNA base-calling from a nanopore using a Viterbi algorithm"
+Winston Timp, Jeffrey Comer, and Aleksei Aksimentiev
+Biophys J 102(10) L37-9 (2012)
+
+"Predicting the DNA sequence dependence of nanopore ion current using atomic-resolution Brownian dynamics"
+Jeffrey Comer and Aleksei Aksimentiev.
+J Phys Chem C Nanomater Interfaces 116:3376-3393 (2012).
"Atoms-to-microns model for small solute transport through sticky nanochannels"
Rogan Carr, Jeffrey Comer, Mark D. Ginsberg, and Aleksei Aksimentiev
Lab Chip 11(22) 3766-73 (2011)
+
/=========\
| Authors |
\=========/
-ARBD is developed by the Aksimentiev group (http://bionano.physics.illinois.edu).
+ARBD is developed by the Aksimentiev group (http://bionano.physics.illinois.edu)
+as a part of the NIH Center for Macromolecular Modeling and Bioinformatics
+(http://www.ks.uiuc.edu/).
Please direct questions or problems to Chris.
diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu
index 92445aa335a0c1307789a7052a58936885a05cd8..66bd8bb1d1817e6cd503a51d35373a916313a815 100644
--- a/src/ComputeForce.cu
+++ b/src/ComputeForce.cu
@@ -95,6 +95,14 @@ ComputeForce::ComputeForce(int num, const BrownianParticleType part[],
}
gpuErrchk(cudaMalloc(&tableDihedral_d, sizeof(TabulatedDihedralPotential*) * numTabDihedralFiles));
+ { // allocate device for pairlists
+ // RBTODO: select maxpairs in better way
+ const int maxPairs = 1<<20;
+ gpuErrchk(cudaMalloc(&numPairs_d, sizeof(int)));
+ gpuErrchk(cudaMalloc(&pairLists_d, sizeof(int2)*maxPairs));
+ gpuErrchk(cudaMalloc(&pairTabPotType_d, sizeof(int)*maxPairs));
+ }
+
//Calculate the number of blocks the grid should contain
gridSize = num / NUM_THREADS + 1;
@@ -109,10 +117,10 @@ ComputeForce::~ComputeForce() {
delete[] tableEps;
delete[] tableRad6;
delete[] tableAlpha;
- cudaFree(tableEps_d);
- cudaFree(tableAlpha_d);
- cudaFree(tableRad6_d);
-
+ gpuErrchk(cudaFree(tableEps_d));
+ gpuErrchk(cudaFree(tableAlpha_d));
+ gpuErrchk(cudaFree(tableRad6_d));
+
for (int j = 0; j < numParts * numParts; ++j)
delete tablePot[j];
delete[] tablePot;
@@ -149,6 +157,11 @@ ComputeForce::~ComputeForce() {
gpuErrchk( cudaFree(dihedrals_d) );
gpuErrchk( cudaFree(bondList_d) );
}
+
+ gpuErrchk(cudaFree(numPairs_d));
+ gpuErrchk(cudaFree(pairLists_d));
+ gpuErrchk(cudaFree(pairTabPotType_d));
+
}
void ComputeForce::updateNumber(int newNum) {
@@ -395,17 +408,6 @@ void ComputeForce::decompose() {
// initializePairlistArrays
int nCells = decomp.nCells.x * decomp.nCells.y * decomp.nCells.z;
int blocksPerCell = 10;
- if (newDecomp) {
- // RBTODO: free memory elsewhere
- // allocate device data
- const int maxPairs = 1<<20;
- gpuErrchk(cudaMalloc(&numPairs_d, sizeof(int)));
-
- gpuErrchk(cudaMalloc(&pairLists_d, sizeof(int2)*maxPairs));
- gpuErrchk(cudaMalloc(&pairTabPotType_d, sizeof(int)*maxPairs));
-
- gpuErrchk(cudaDeviceSynchronize());
- }
/* cuMemGetInfo(&free,&total); */
@@ -451,12 +453,6 @@ CellDecomposition ComputeForce::getDecomp() { return decomp; }
float ComputeForce::decompCutoff() { return decomp.getCutoff(); }
-// TODO: Fix this
-int* ComputeForce::neighborhood(Vector3 r) {
- // return decomp.getCell(r)->getNeighbors();
- return NULL;
-}
-
float ComputeForce::computeFull(bool get_energy) {
float energy = 0.0f;
gridSize = (num * numReplicas) / NUM_THREADS + 1;
diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh
index bd524dc23d8420bcf534fc90931f67a9478a72d8..6a7e4fe0ea4535947f335aa55cc5743a259f51d3 100644
--- a/src/ComputeForce.cuh
+++ b/src/ComputeForce.cuh
@@ -299,7 +299,9 @@ void createPairlists(Vector3* __restrict__ pos, const int num, const int numRepl
/* if (tid == 0) printf("Cell%d: found %d pairs\n",cID,g_numPairs[cID]); */
}
}
- __global__
+
+// TODO: deprecate?
+__global__
void computeKernel(Vector3 force[], Vector3 pos[], int type[],
float tableAlpha[], float tableEps[], float tableRad6[],
int num, int numParts, BaseGrid* sys,
@@ -431,17 +433,14 @@ __global__ void computeTabulatedEnergyKernel(Vector3* force, const Vector3* __re
// RBTODO: implement wrapDiff2, returns dr2 (???)
Vector3 dr = pos[aj] - pos[ai];
dr = sys->wrapDiff(dr);
- // Calculate the force based on the tabulatedPotential file
float d2 = dr.length2();
- /* RBTODO: filter tabpot particles ahead of time */
// RBTODO: order pairs according to distance to reduce divergence // not actually faster
if (tablePot[ind] != NULL && d2 <= cutoff2) {
EnergyForce fe = tablePot[ind]->compute(dr,d2);
- // RBTODO: is this the best approach?
atomicAdd( &force[ai], -fe.f );
atomicAdd( &force[aj], fe.f );
- // RBTODO: why are energies calculated for each atom? Could be reduced
+ // RBTODO: reduce energies
atomicAdd( &(g_energies[ai]), fe.e );
atomicAdd( &(g_energies[aj]), fe.e );
}
@@ -579,70 +578,7 @@ void computeAngles(Vector3 force[], Vector3 pos[],
}
}
-//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)
-{
- // 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; //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
- * 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
-
- force_local += tableBond[ bonds[currBond].tabFileIndex ]->computef(dr,dr.length2());
- }
- atomicAdd( &force[i], force_local );
-
- if (get_energy)
- atomicAdd( &g_energies[i], energy_local);
- }
-}*/
-
+// TODO: add kernels for energy calculations
__global__ void computeTabulatedBonds(Vector3* force,
Vector3* __restrict__ pos,
BaseGrid* __restrict__ sys,
@@ -655,9 +591,6 @@ __global__ void computeTabulatedBonds(Vector3* force,
int i = bondList_d[bid].x;
int j = bondList_d[bid].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]);
@@ -676,7 +609,6 @@ __global__ void computeTabulatedBonds(Vector3* force,
}
}
-// TODO: add kernel for energy calculation
__global__
void computeTabulatedAngles(Vector3* force,
Vector3* __restrict__ pos,
@@ -793,7 +725,6 @@ void computeTabulatedDihedrals(Vector3* force, const Vector3* __restrict__ pos,
// int currDihedral = blockIdx.x * blockDim.x + threadIdx.x; // first particle ID
// Loop over ALL dihedrals in ALL replicas
- // TODO make grid stride loop
for (int i = threadIdx.x+blockIdx.x*blockDim.x; i < numDihedrals; i+=blockDim.x*gridDim.x) {
const int4& ids = dihedralList_d[i];
const int& id = dihedralPotList_d[i];
diff --git a/src/ComputeGridGrid.cu b/src/ComputeGridGrid.cu
index c6d1f27bcf515df601b004ee76b3df20f84897f3..b074fde082cdd327ccb30c29ec11fc98b0042cb2 100644
--- a/src/ComputeGridGrid.cu
+++ b/src/ComputeGridGrid.cu
@@ -30,7 +30,7 @@ void computeGridGridForce(const RigidBodyGrid* rho, const RigidBodyGrid* u,
r_pos = basis_rho.transform( r_pos ) + origin_rho_minus_origin_u; /* real space */
const Vector3 u_ijk_float = basis_u_inv.transform( r_pos );
- // RBTODO What about non-unit delta?
+ // RBTODO: Test for non-unit delta
/* Vector3 tmpf = Vector3(0.0f); */
/* float tmpe = 0.0f; */
/* const ForceEnergy fe = ForceEnergy( tmpf, tmpe); */
@@ -90,7 +90,7 @@ void computePartGridForce(const Vector3* __restrict__ pos, Vector3* particleForc
atomicAdd( &particleForce[id], force[tid] ); // apply force to particle
// Calculate torque about origin_u in the lab frame
- torque[tid] = p.cross(force[tid]); // RBTODO: test if sign is correct!
+ torque[tid] = p.cross(force[tid]); // RBTODO: test sign
}
// Reduce force and torques
diff --git a/src/Makefile b/src/Makefile
index 2450d18bc798877c27baeb0dd6962587aca7646d..220beccaf4d7815c18349ee81f637e6aeee54cda 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -49,7 +49,7 @@ LD_FLAGS += -lcuda
### Sources
CC_SRC := $(wildcard *.cpp)
-CC_SRC := $(filter-out runBrownTown.cpp, $(CC_SRC))
+CC_SRC := $(filter-out arbd.cpp, $(CC_SRC))
CU_SRC := $(wildcard *.cu)
CC_OBJ := $(patsubst %.cpp, %.o, $(CC_SRC))
@@ -57,15 +57,15 @@ CU_OBJ := $(patsubst %.cu, %.o, $(CU_SRC))
### Targets
-TARGET = runBrownCUDA
+TARGET = arbd
all: $(TARGET)
@echo "Done ->" $(TARGET)
-$(TARGET): $(CU_OBJ) $(CC_OBJ) runBrownTown.cpp vmdsock.c imd.c imd.h
+$(TARGET): $(CU_OBJ) $(CC_OBJ) arbd.cpp vmdsock.c imd.c imd.h
$(NVCC) $(NVLD_FLAGS) $(CU_OBJ) $(CC_OBJ) -o $(TARGET)_link.o
$(NVCC) $(NVLD_FLAGS) $(CU_OBJ) -o $(TARGET)_link.o
- $(CC) $(CC_FLAGS) $(EX_FLAGS) runBrownTown.cpp vmdsock.c imd.c $(TARGET)_link.o $(CU_OBJ) $(CC_OBJ) $(LD_FLAGS) -o $(TARGET)
+ $(CC) $(CC_FLAGS) $(EX_FLAGS) arbd.cpp vmdsock.c imd.c $(TARGET)_link.o $(CU_OBJ) $(CC_OBJ) $(LD_FLAGS) -o $(TARGET)
.SECONDEXPANSION:
$(CU_OBJ): %.o: %.cu $$(wildcard %.h) $$(wildcard %.cuh)
diff --git a/src/runBrownTown.cpp b/src/arbd.cpp
similarity index 100%
rename from src/runBrownTown.cpp
rename to src/arbd.cpp