diff --git a/src/ComputeGridGrid.cu b/src/ComputeGridGrid.cu
index 1ecfb9fad191cb1e70d787e384198b57b2a87883..f6cb72cc78899b85ca0395e28e9021dedd9c92ed 100644
--- a/src/ComputeGridGrid.cu
+++ b/src/ComputeGridGrid.cu
@@ -171,9 +171,11 @@ void computePartGridForce(const Vector3* __restrict__ pos, Vector3* particleForc
}
}
+
__global__
void createPartlist(const Vector3* __restrict__ pos,
const int numTypeParticles, const int* __restrict__ typeParticles_d,
+ const int attached_particle_start, const int attached_particle_end,
int* numParticles_d, int* particles_d,
const Vector3 gridCenter, const float radius2, BaseGrid* sys_d) {
const int tid = threadIdx.x;
@@ -182,11 +184,13 @@ void createPartlist(const Vector3* __restrict__ pos,
const int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < numTypeParticles) {
int aid = typeParticles_d[i];
- float dist = (sys_d->wrapDiff(pos[aid] - gridCenter)).length2();
-
- if (dist <= radius2) {
+ if (aid < attached_particle_start || aid >= attached_particle_end) {
+ float dist = (sys_d->wrapDiff(pos[aid] - gridCenter)).length2();
+
+ if (dist <= radius2) {
int tmp = atomicAggInc(numParticles_d, warpLane);
particles_d[tmp] = aid;
+ }
}
}
}
diff --git a/src/ComputeGridGrid.cuh b/src/ComputeGridGrid.cuh
index 739aa622e2514840f8586c47d0a277d306427ff4..27fc32ccb6b3b3fd43dd005fd9eda7eb1689c8f8 100644
--- a/src/ComputeGridGrid.cuh
+++ b/src/ComputeGridGrid.cuh
@@ -29,6 +29,7 @@ void computePartGridForce(const Vector3* __restrict__ pos, Vector3* particleForc
extern __global__
void createPartlist(const Vector3* __restrict__ pos,
const int numTypeParticles, const int* __restrict__ typeParticles_d,
+ const int attached_particle_start, const int attached_particle_end,
int* numParticles_d, int* particles_d,
const Vector3 gridCenter, const float radius2, BaseGrid* sys_d);
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index 679daac15e9d595c81ec74848126c6a58aad79bf..e7bdfd71c57e163646c644d6ef3bb85485f4e7a7 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -133,6 +133,20 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
serial = new int[(num+num_rb_attached_particles) * simNum];
posLast = new Vector3[(num+num_rb_attached_particles) * simNum];
+ {
+ int pidx = 0;
+ for (int i = 0; i < numRigidTypes; i++) { // Loop over RB types
+ RigidBodyType &rbt = rigidBody[i];
+ for (int j = 0; j < rbt.num; ++j) { // Loop over RBs
+ for (const int& t: rbt.get_attached_particle_types()) {
+ type[num+pidx] = t;
+ serial[num+pidx] = num+pidx;
+ pidx++;
+ }
+ }
+ }
+ }
+
//Han-Yi Chou
if(ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
momLast = new Vector3[(num+num_rb_attached_particles) * simNum];
@@ -434,10 +448,9 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
// Create exclusions from the exclude rule, if it was specified in the config file
- assert(false); // TODO implement RB exclusions
if (excludeRule != String("")) {
int oldNumExcludes = numExcludes;
- Exclude* newExcludes = makeExcludes(bonds, bondMap, num+num_rb_attached_particles, numBonds, excludeRule, numExcludes);
+ Exclude* newExcludes = makeExcludes(bonds, bondMap, num, numBonds, excludeRule, numExcludes);
if (excludes == NULL) {
excludes = new Exclude[numExcludes];
} else if (numExcludes >= excludeCapacity) {
@@ -445,13 +458,45 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
excludes = new Exclude[numExcludes];
for (int i = 0; i < oldNumExcludes; i++)
excludes[i] = tempExcludes[i];
- delete tempExcludes;
+ delete [] tempExcludes;
}
for (int i = oldNumExcludes; i < numExcludes; i++)
excludes[i] = newExcludes[i - oldNumExcludes];
printf("Built %d exclusions.\n",numExcludes);
}
+ { // Add exclusions for RB attached particles
+ std::vector<Exclude> ex;
+ int start = num;
+ for (int i = 0; i < numRigidTypes; i++) { // Loop over RB types
+ RigidBodyType &rbt = rigidBody[i];
+ const int nap = rbt.num_attached_particles();
+ for (int j = 0; j < rbt.num; ++j) { // Loop over RBs
+ for (int ai = 0; ai < nap-1; ++ai) {
+ for (int aj = ai+1; aj < nap; ++aj) {
+ ex.push_back( Exclude( ai+start, aj+start ) );
+ }
+ }
+ start += nap;
+ }
+ }
+ // copy
+ int oldNumExcludes = numExcludes;
+ numExcludes = numExcludes + ex.size();
+ if (excludes == NULL) {
+ excludes = new Exclude[numExcludes];
+ } else if (numExcludes >= excludeCapacity) {
+ Exclude* tempExcludes = excludes;
+ excludes = new Exclude[numExcludes];
+ for (int i = 0; i < oldNumExcludes; i++)
+ excludes[i] = tempExcludes[i];
+ delete [] tempExcludes;
+ }
+ for (int i = oldNumExcludes; i < numExcludes; i++)
+ excludes[i] = ex[i - oldNumExcludes];
+ }
+
+ printf("Built %d exclusions.\n",numExcludes);
buildExcludeMap();
// Count number of particles of each type
@@ -459,7 +504,6 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
for (int i = 0; i < numParts; ++i) {
numPartsOfType[i] = 0;
}
- assert(false); // First assign type[i] for rb particles
for (int i = 0; i < num+num_rb_attached_particles; ++i) {
++numPartsOfType[type[i]];
}
@@ -1303,7 +1347,7 @@ void Configuration::readAtoms() {
for (int i = 0; i < numParts; i++)
if (part[i].num == 0)
found = false;
- assert(false); // TODO probably relax constraint that particle must be found; could just be in RB
+ // assert(false); // TODO probably relax constraint that particle must be found; could just be in RB
if (!found) {
printf("ERROR: Number of particles not specified in config file.\n");
exit(1);
diff --git a/src/GrandBrownTown.cu b/src/GrandBrownTown.cu
index f6c9d380421fbd63a22ccfe65d9dbcb234ed636d..339a9f37938bdb0481e8d52f022e7b3d6a288f8a 100644
--- a/src/GrandBrownTown.cu
+++ b/src/GrandBrownTown.cu
@@ -503,8 +503,7 @@ void GrandBrownTown::RunNoseHooverLangevin()
// Open the files for recording ionic currents
for (int repID = 0; repID < numReplicas; ++repID)
{
-
- writers[repID]->newFile(pos + (repID * num), name, 0.0f, num); // 'pos + (repID*num)' == array-to-pointer decay
+ writers[repID]->newFile(pos + repID*(num+num_rb_attached_particles), name, 0.0f, num); // 'pos + (repID*num)' == array-to-pointer decay
if(particle_dynamic == String("Langevin") || particle_dynamic == String("NoseHooverLangevin"))
momentum_writers[repID]->newFile((momentum + repID * num), name, 0.0f, num); // 'pos + (repID*num)' == array-to-pointer decay
//random_writers[repID]->newFile(random + (repID * num), name, 0.0f, num);
@@ -606,6 +605,18 @@ void GrandBrownTown::RunNoseHooverLangevin()
internal->clear_force();
internal->clear_energy();
const std::vector<Vector3*>& _pos = internal->getPos_d();
+
+ if (num_rb_attached_particles > 0) {
+ #pragma omp parallel for
+ for(int i = 0; i < numReplicas; ++i) {
+ RBC[i]->update_attached_particle_positions(
+ internal->getPos_d()[0]+i*(num+num_rb_attached_particles),
+ internal->getForceInternal_d()[0]+i*(num+num_rb_attached_particles),
+ internal->getEnergy()+i*(num+num_rb_attached_particles),
+ sys_d);
+ }
+ }
+
if (numGroupSites > 0) updateGroupSites<<<(numGroupSites/32+1),32>>>(_pos[0], groupSiteData_d, num + num_rb_attached_particles, numGroupSites, numReplicas);
#ifdef USE_NCCL
@@ -720,7 +731,7 @@ void GrandBrownTown::RunNoseHooverLangevin()
internal->clear_energy();
gpuman.sync();
- assert(false); // TODO update RB particels
+
if(particle_dynamic == String("Langevin"))
updateKernelBAOAB<<< numBlocks, NUM_THREADS >>>(internal->getPos_d()[0], internal->getMom_d(), internal->getForceInternal_d()[0], internal->getType_d(), part_d, kT, kTGrid_d, electricField, tl, timestep, num, num_rb_attached_particles, sys_d, randoGen_d, numReplicas, ParticleInterpolationType);
else if(particle_dynamic == String("NoseHooverLangevin"))
@@ -760,6 +771,17 @@ void GrandBrownTown::RunNoseHooverLangevin()
}
}
+ if (num_rb_attached_particles > 0) {
+ #pragma omp parallel for
+ for(int i = 0; i < numReplicas; ++i) {
+ RBC[i]->update_attached_particle_positions(
+ internal->getPos_d()[0]+i*(num+num_rb_attached_particles),
+ internal->getForceInternal_d()[0]+i*(num+num_rb_attached_particles),
+ internal->getEnergy()+i*(num+num_rb_attached_particles),
+ sys_d);
+ }
+ }
+
if (s % outputPeriod == 0) {
// Copy particle positions back to CPU
gpuErrchk(cudaDeviceSynchronize());
@@ -842,6 +864,7 @@ void GrandBrownTown::RunNoseHooverLangevin()
for(int i = 0; i < numReplicas; ++i)
RBC[i]->clearForceAndTorque();
+
if (numGroupSites > 0) {
gpuman.sync();
updateGroupSites<<<(numGroupSites/32+1),32>>>(internal->getPos_d()[0], groupSiteData_d, num + num_rb_attached_particles, numGroupSites, numReplicas);
diff --git a/src/Makefile b/src/Makefile
index b0356b86521cc4ca5448f5fdb8b0b35ea106b2dc..4f6bb5779b9248ff6e4ef7e018387e236b18ab58 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -102,11 +102,16 @@ endif
all: $(TARGET)
@echo "Done ->" $(TARGET)
-$(TARGET): $(CU_OBJ) $(CC_OBJ) arbd.cpp vmdsock.c imd.c imd.h
- $(NVCC) $(NVLD_FLAGS) $(CU_OBJ) $(CC_OBJ) -o $(TARGET)_link.o
+$(TARGET)_link.o: $(CU_OBJ)
+# $(NVCC) $(NVLD_FLAGS) $(CU_OBJ) $(CC_OBJ) -o $(TARGET)_link.o
$(NVCC) $(NVLD_FLAGS) $(CU_OBJ) -o $(TARGET)_link.o
+
+$(TARGET): $(TARGET)_link.o $(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) 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)
$(NVCC) $(NV_FLAGS) $(EX_FLAGS) -dc $< -o $@
diff --git a/src/RigidBody.cu b/src/RigidBody.cu
index f65fcc75298816e893a6052045e4815c0ef08527..abd693bdd3895f19af9027655132825501b52c37 100644
--- a/src/RigidBody.cu
+++ b/src/RigidBody.cu
@@ -18,10 +18,10 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t
}
#endif
-RigidBody::RigidBody(String name, const Configuration& cref, const RigidBodyType& tref, RigidBodyController* RBCref)
- : name(name), c(&cref), t(&tref), RBC(RBCref), impulse_to_momentum(4.1867999435271e4) /*impulse_to_momentum(4.184e8f)*/ { init(); }
+RigidBody::RigidBody(String name, const Configuration& cref, const RigidBodyType& tref, RigidBodyController* RBCref, int attached_particle_start, int attached_particle_end)
+ : name(name), c(&cref), t(&tref), RBC(RBCref), attached_particle_start(attached_particle_start), attached_particle_end(attached_particle_end), impulse_to_momentum(4.1867999435271e4) /*impulse_to_momentum(4.184e8f)*/ { init(); }
RigidBody::RigidBody(const RigidBody& rb)
- : name(rb.name), c(rb.c), t(rb.t), RBC(rb.RBC), impulse_to_momentum(4.1867999435271e4)/*impulse_to_momentum(4.184e8f)*/ { init(); }
+ : name(rb.name), c(rb.c), t(rb.t), RBC(rb.RBC), attached_particle_start(rb.attached_particle_start), attached_particle_end(rb.attached_particle_end), impulse_to_momentum(4.1867999435271e4)/*impulse_to_momentum(4.184e8f)*/ { init(); }
void RigidBody::init() {
// units "(kcal_mol/AA) * ns" "dalton AA/ns" * 4.184e+08
timestep = c->timestep;
@@ -103,6 +103,23 @@ int RigidBody::appendNumParticleBlocks( std::vector<int>* blocks ) {
return ret;
}
+__global__
+void update_particle_positions_kernel(Vector3* __restrict__ pos, const int start, const int num,
+ const Vector3* __restrict__ pos_rb,
+ const Vector3 center, const Matrix3 orientation) {
+ const int i = blockIdx.x * blockDim.x + threadIdx.x;
+ if (i < num) {
+ const int aid = i+start;
+ pos[aid] = orientation*pos_rb[i] + center;
+ }
+}
+void RigidBody::update_particle_positions(Vector3* pos_d, Vector3* force_d, float* energy_d) {
+ int num = attached_particle_end - attached_particle_start;
+ int nb = floor(num/NUMTHREADS) + 1;
+ update_particle_positions_kernel<<<nb,NUMTHREADS>>>(pos_d, attached_particle_start, num,
+ t->attached_particle_positions_d, position, orientation);
+}
+
void RigidBody::addForce(Force f) {
force += f;
}
@@ -113,8 +130,142 @@ void RigidBody::addEnergy(float e)
{
energy += e;
}
+
+// TODO move code snippet to more general location
+template<class C, class T>
+auto contains(const C& v, const T& x)
+-> decltype(end(v), true)
+{
+ return end(v) != std::find(begin(v), end(v), x);
+}
+/*
+void RigidBody::initialize_possible_particles()
+{
+ std::vector<int> atomic_ids;
+
+ // Loop over potential grids
+ for (int i = 0; i < t->numPotGrids; ++i) {
+ atomic_ids.clear();
+
+ String& gridName = t->potentialGridKeys[i];
+
+ // Loop over particle types to count the number of particles
+ for (int j = 0; j < conf->numParts; ++j) {
+ if (contains(conf->partRigidBodyGrid[j], gridName)) {
+ // gridNames contained gridName, so add the particles to atomic_ids
+ for (int aid = 0; aid < conf->num + conf->num_rb_attached_particles; ++aid) {
+ if (conf->type[aid] == j && (aid < exclude_start || aid > exclude_end)) {
+ atomic_ids.push_back(aid);
+ }
+ }
+ }
+ }
+
+ // Initialize device data
+ size_t sz = sizeof(int) * atomic_ids.size();
+ gpuErrchk(cudaMalloc( &(possible_particles_d[i]), sz ));
+ gpuErrchk(cudaMemcpyAsync( possible_particles_d[i], &atomic_ids[0], sz, cudaMemcpyHostToDevice))
+
+ // // Add particles attached to OTHER RBs
+ // int rb_particle_offset = conf->num;
+ // for (const auto& rbs: RBC->rigidBodyByType)
+ // {
+ // int rb_type
+ // const RigidBodyType& rb_type = rbs[0].t;
+ // if (rbs.t !=
+ // for (int& ptype: rb_type->attached_particle_types) {
+ // // This could be made much more efficient
+ // if (contains(conf->partRigidBodyGrid[ptype], gridName)) {
+ // // Add particles
+ // rb
+ // rbs[0].sum
+ // }
+ // rb_particle_offset += rb_type.attached_particle_types
+
+ // for (int k = 0; k < t->gridNames.size(); ++k) {
+ // if (t->gridNames[k] == gridName) {
+ // t->numParticles[i] += conf->numPartsOfType[j];
+ // }
+
+
+ // // Loop over rigid body grid names associated with particle type
+ // const std::vector<String>& gridNames = conf->partRigidBodyGrid[j];
+ // for (int k = 0; k < gridNames.size(); ++k) {
+ // if (gridNames[k] == gridName) {
+ // // Copy type j particles to particles[i]
+ // memcpy( &(particles[i][pid]), tmp, sizeof(int)*currId );
+ // assert(currId == conf->numPartsOfType[j]);
+ // pid += conf->numPartsOfType[j];
+ // }
+ // }
+ // }
+
+ // numParticles[i] = 0;
+
+ // Count the particles interacting with potential grid i
+ // Loop over particle types
+ for (int j = 0; j < conf->numParts; ++j) {
+
+ // Loop over rigid body grid names associated with particle type
+ const std::vector<String>& gridNames = conf->partRigidBodyGrid[j];
+ for (int k = 0; k < t->gridNames.size(); ++k) {
+ if (t->gridNames[k] == gridName) {
+ numParticles[i] += conf->numPartsOfType[j];
+ }
+ }
+ }
+
+ // Add RB particles
+ for (const auto& rbv: RBC->rigidBodyByType)
+ {
+ int ptype = rbv[0].t->attached_particle_types;
+ const std::vector<String>& gridNames = conf->partRigidBodyGrid[ptype];
+ for (int k = 0; k < t->gridNames.size(); ++k) {
+ if (t->gridNames[k] == gridName) {
+ t->numParticles[i] += conf->numPartsOfType[j];
+ }
+
+ attached_particle_
+
+ if (numParticles[i] > 0) {
+
+ // allocate array of particle ids for the potential grid
+ particles[i] = new int[numParticles[i]];
+ int pid = 0;
+
+ // Loop over particle types to count the number of particles
+ for (int j = 0; j < conf->numParts; ++j) {
+
+ // Build temporary id array of type j particles
+ int tmp[conf->numPartsOfType[j]];
+ int currId = 0;
+ for (int aid = 0; aid < conf->num; ++aid) {
+ if (conf->type[aid] == j)
+ tmp[currId++] = aid;
+ }
+ if (currId == 0) continue;
+
+ // Loop over rigid body grid names associated with particle type
+ const std::vector<String>& gridNames = conf->partRigidBodyGrid[j];
+ for (int k = 0; k < gridNames.size(); ++k) {
+ if (gridNames[k] == gridName) {
+ // Copy type j particles to particles[i]
+ memcpy( &(particles[i][pid]), tmp, sizeof(int)*currId );
+ assert(currId == conf->numPartsOfType[j]);
+ pid += conf->numPartsOfType[j];
+ }
+ }
+ }
+
+ // Initialize device data
+ size_t sz = sizeof(int) * numParticles[i];
+ gpuErrchk(cudaMalloc( &(particles_d[i]), sz ));
+ gpuErrchk(cudaMemcpyAsync( particles_d[i], particles[i], sz, cudaMemcpyHostToDevice));
+ }
+ }
+}
+*/
void RigidBody::updateParticleList(Vector3* pos_d, BaseGrid* sys_d) {
- assert(false); // TODO update particles, excluding own particles
for (int i = 0; i < t->numPotGrids; ++i) {
numParticles[i] = 0;
int& tnp = t->numParticles[i];
@@ -132,10 +283,12 @@ void RigidBody::updateParticleList(Vector3* pos_d, BaseGrid* sys_d) {
int nb = floor(tnp/NUMTHREADS) + 1;
#if __CUDA_ARCH__ >= 300
createPartlist<<<nb,NUMTHREADS>>>(pos_d, tnp, t->particles_d[i],
+ attached_particle_start, attached_particle_end,
tmp_d, particles_d[i],
gridCenter + position, cutoff*cutoff, sys_d);
#else
createPartlist<<<nb,NUMTHREADS,NUMTHREADS/WARPSIZE>>>(pos_d, tnp, t->particles_d[i],
+ attached_particle_start, attached_particle_end,
tmp_d, particles_d[i],
gridCenter + position, cutoff*cutoff, sys_d);
#endif
diff --git a/src/RigidBody.h b/src/RigidBody.h
index 5854a082a29f975d11281b9ad45ffd5bc8cc0fbf..636c80e94fd0f17475156bbc84e0cd74409d8e94 100644
--- a/src/RigidBody.h
+++ b/src/RigidBody.h
@@ -31,7 +31,9 @@ class RigidBody { // host side representation of rigid bodies
| splitting methods for rigid body molecular dynamics". J Chem Phys. (1997) |
`––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/
public:
- RigidBody(String name, const Configuration& c, const RigidBodyType& t, RigidBodyController* RBC);
+ RigidBody(String name, const Configuration& c, const RigidBodyType& t, RigidBodyController* RBC,
+ int attached_particle_start, int attached_particle_end);
+
RigidBody(const RigidBody& rb);
// RigidBody(const RigidBody& rb) : RigidBody(rb.name, *rb.c, *rb.t) {};
void init();
@@ -40,6 +42,8 @@ class RigidBody { // host side representation of rigid bodies
int appendNumParticleBlocks( std::vector<int>* blocks );
+ void update_particle_positions(Vector3* pos_d, Vector3* force_d, float* energy_d);
+
HOST DEVICE void addForce(Force f);
HOST DEVICE void addTorque(Force t);
HOST DEVICE void addEnergy(float e);
@@ -76,6 +80,7 @@ class RigidBody { // host side representation of rigid bodies
return Vector3( angularMomentum.x, angularMomentum.y, angularMomentum.z);
}
+ void initializeParticleLists();
void updateParticleList(Vector3* pos_d, BaseGrid* sys_d);
void callGridParticleForceKernel(Vector3* pos_d, Vector3* force_d, int s, float* energy, bool get_energy, int scheme, BaseGrid* sys, BaseGrid* sys_d, ForceEnergy* forcestorques_d, const std::vector<int>& forcestorques_offset, int& fto_idx);
void applyGridParticleForces(BaseGrid* sys, ForceEnergy* forcestorques, const std::vector<int>& forcestorques_offset, int& fto_idx);
@@ -122,9 +127,12 @@ private:
bool isFirstStep;
int* numParticles; /* particles affected by potential grids */
+ int** possible_particles_d;
int** particles_d;
const cudaStream_t** particleForceStreams;
-
+
+ int attached_particle_start, attached_particle_end;
+
/*–––––––––––––––––––––––––––––––––––––––––.
| units "kcal_mol/AA * ns" "(AA/ns) * amu" |
`–––––––––––––––––––––––––––––––––––––––––*/
diff --git a/src/RigidBodyController.cu b/src/RigidBodyController.cu
index c5e855ebb34692a168a72198431ff1450d914a53..b16eedbe24610071705a3772919ac3561eef8840 100644
--- a/src/RigidBodyController.cu
+++ b/src/RigidBodyController.cu
@@ -45,15 +45,14 @@ RigidBodyController::RigidBodyController(const Configuration& c, const char* pre
gpuErrchk(cudaDeviceSynchronize()); /* RBTODO: this should be extraneous */
construct_grids();
- for (int i = 0; i < conf.numRigidTypes; i++)
- conf.rigidBody[i].initializeParticleLists();
int numRB = 0;
// grow list of rbs
-
+ int attached_particle_offset = 0;
for (int i = 0; i < conf.numRigidTypes; i++) {
numRB += conf.rigidBody[i].num;
+ int attached_particle_in_type = conf.rigidBody[i].num_attached_particles();
std::vector<RigidBody> tmp;
// RBTODO: change conf.rigidBody to conf.rigidBodyType
const int jmax = conf.rigidBody[i].num;
@@ -64,13 +63,16 @@ RigidBodyController::RigidBodyController(const Configuration& c, const char* pre
snprintf(stmp, 128, "#%d", j);
name.add( stmp );
}
- RigidBody r(name, conf, conf.rigidBody[i], this);
+ RigidBody r(name, conf, conf.rigidBody[i], this,
+ attached_particle_offset, attached_particle_offset+attached_particle_in_type);
+ attached_particle_offset += attached_particle_in_type;
+
int nb = r.appendNumParticleBlocks( &particleForceNumBlocks );
tmp.push_back( r );
}
rigidBodyByType.push_back(tmp);
}
-
+
totalParticleForceNumBlocks = 0;
for (int i=0; i < particleForceNumBlocks.size(); ++i) {
particleForce_offset.push_back(2*totalParticleForceNumBlocks);
@@ -525,6 +527,16 @@ void RigidBodyController::initializeForcePairs() {
}
+void RigidBodyController::update_attached_particle_positions(Vector3* pos_d, Vector3* force_d, float* energy_d, BaseGrid* sys_d) {
+ for (int i = 0; i < rigidBodyByType.size(); i++) {
+ for (int j = 0; j < rigidBodyByType[i].size(); j++) {
+ rigidBodyByType[i][j].update_particle_positions(pos_d, force_d, energy_d);
+ }
+ }
+ assert(false); // clear force and energy
+}
+
+
void RigidBodyController::updateParticleLists(Vector3* pos_d, BaseGrid* sys_d) {
for (int i = 0; i < rigidBodyByType.size(); i++) {
for (int j = 0; j < rigidBodyByType[i].size(); j++) {
diff --git a/src/RigidBodyController.h b/src/RigidBodyController.h
index 0e2dbadbf1e5eb35f9ad4fd007fe1eeb145e8f81..bf8e58637d2b9860a8d1fd132f2c34a108ecfcb2 100644
--- a/src/RigidBodyController.h
+++ b/src/RigidBodyController.h
@@ -110,6 +110,7 @@ public:
void integrateDLM(BaseGrid* sys, int step);
void updateForces(Vector3* pos_d, Vector3* force_d, int s, float* energy, bool get_energy, int scheme, BaseGrid* sys, BaseGrid* sys_d);
void updateParticleLists(Vector3* pos_d, BaseGrid* sys_d);
+ void update_attached_particle_positions(Vector3* pos_d, Vector3* force_d, float* energy_d, BaseGrid* sys_d);
void clearForceAndTorque();
void KineticEnergy();
void print(int step);
diff --git a/src/RigidBodyType.cu b/src/RigidBodyType.cu
index 8f7739ceecc4a15b50dd0258e0f58d82ea6480d8..ec0f8b65c7a25587755326291ce45542af777c44 100644
--- a/src/RigidBodyType.cu
+++ b/src/RigidBodyType.cu
@@ -128,6 +128,9 @@ void RigidBodyType::attach_particles() {
}
fclose(inp);
}
+ size_t sz = sizeof(Vector3)*attached_particle_positions.size();
+ gpuErrchk(cudaMalloc( &(attached_particle_positions_d), sz ));
+ gpuErrchk(cudaMemcpyAsync( attached_particle_positions_d, &attached_particle_positions[0], sz, cudaMemcpyHostToDevice));
}
void RigidBodyType::addGrid(String s, std::vector<String> &keys, std::vector<String> &files) {
diff --git a/src/RigidBodyType.h b/src/RigidBodyType.h
index 7a6edcc1dc6d1294e2a8cb37ce61b4eacd221870..5cc1442052a102da7198e1e6eeda806c4cdaa8b7 100644
--- a/src/RigidBodyType.h
+++ b/src/RigidBodyType.h
@@ -17,9 +17,10 @@ class BaseGrid;
class RigidBodyGrid;
class Configuration;
class RigidBodyController;
+class RigidBody;
class RigidBodyType {
-
+ friend class RigidBody;
public:
RigidBodyType(const String& name = "", const Configuration* conf = NULL ) :
name(name), conf(conf), num(0),
@@ -41,7 +42,8 @@ public:
void append_attached_particle_file(String s) { attached_particle_files.push_back(s); }
void attach_particles();
- size_t num_attached_particles() { return attached_particle_types.size() ;}
+ size_t num_attached_particles() const { return attached_particle_types.size() ;}
+ const std::vector<int>& get_attached_particle_types() const { return attached_particle_types; }
void addPotentialGrid(String s);
void addDensityGrid(String s);
@@ -59,8 +61,10 @@ public:
private:
const Configuration* conf;
std::vector<String> attached_particle_files;
- std::vector<int>attached_particle_types;
+ std::vector<int>attached_particle_types;
+private:
std::vector<Vector3>attached_particle_positions;
+ Vector3* attached_particle_positions_d;
public:
int num; // number of particles of this type