From 4d7b1a995b03acfdbc5759f59e09991a97a54e85 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Sun, 25 Apr 2021 11:19:23 -0500
Subject: [PATCH] Initial implementation of .bd file keywords to attach point
particles to rigid bodies
---
src/Configuration.cpp | 88 ++++++++++++++++++++++++++++++-------------
src/Configuration.h | 11 +++++-
src/RigidBodyType.cu | 49 ++++++++++++++++++++++++
src/RigidBodyType.h | 8 ++++
4 files changed, 129 insertions(+), 27 deletions(-)
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index b1a89c4..04f2a94 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -97,8 +97,24 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
}
} // end result: variable "num" is set
+ // Count particles associated with rigid bodies
+ num_rb_attached_particles = 0;
+ if (numRigidTypes > 0) {
+ // grow list of rbs
+ for (int i = 0; i < numRigidTypes; i++) {
+ RigidBodyType &rbt = rigidBody[i];
+ rbt.attach_particles();
+ num_rb_attached_particles += rbt.num * rbt.num_attached_particles();
+ }
+ }
+ assert( num_rb_attached_particles == 0 || simNum == 1 ); // replicas not yet implemented
+ // num = num+num_rb_attached_particles;
+
+
// Set the number capacity
printf("\n%d particles\n", num);
+ printf("%d particles attached to RBs\n", num_rb_attached_particles);
+
if (numCap <= 0) numCap = numCapFactor*num; // max number of particles
if (numCap <= 0) numCap = 20;
@@ -106,18 +122,23 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
printf("%d groups\n", numGroupSites);
// Allocate particle variables.
- pos = new Vector3[num * simNum];
+ // First num*simNum entries are for point particles, next num_rb_attached_particles*simNum are for RB particles
+
+
+ pos = new Vector3[ (num+num_rb_attached_particles) * simNum];
+
//Han-Yi Chou
if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
- momentum = new Vector3[num * simNum];
+ momentum = new Vector3[(num+num_rb_attached_particles) * simNum];
+
+ type = new int[(num+num_rb_attached_particles) * simNum];
+ serial = new int[(num+num_rb_attached_particles) * simNum];
+ posLast = new Vector3[(num+num_rb_attached_particles) * simNum];
- type = new int[num * simNum];
- serial = new int[num * simNum];
- posLast = new Vector3[num * simNum];
//Han-Yi Chou
if(ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
- momLast = new Vector3[num * simNum];
- name = new String[num * simNum];
+ momLast = new Vector3[(num+num_rb_attached_particles) * simNum];
+ name = new String[(num+num_rb_attached_particles) * simNum];
currSerial = 0;
@@ -164,7 +185,12 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
String* tokenList = new String[numTokens];
partsFromFile[i].tokenize(tokenList);
- int currType = 0;
+ int currType = find_particle_type(tokenList[2]);
+ if (currType == -1) {
+ printf("Error: Unable to find particle type %s\n", tokenList[2].val());
+ exit(1);
+
+ }
for (int j = 0; j < numParts; j++)
if (tokenList[2] == part[j].name)
currType = j;
@@ -216,7 +242,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
printf("done!\n");
}
else
- loadedMomentum = Boltzmann(COM_Velocity, num * simNum);
+ loadedMomentum = Boltzmann(COM_Velocity, (num * simNum));
}
}
}
@@ -410,9 +436,10 @@ 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, numBonds, excludeRule, numExcludes);
+ Exclude* newExcludes = makeExcludes(bonds, bondMap, num+num_rb_attached_particles, numBonds, excludeRule, numExcludes);
if (excludes == NULL) {
excludes = new Exclude[numExcludes];
} else if (numExcludes >= excludeCapacity) {
@@ -434,7 +461,8 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
for (int i = 0; i < numParts; ++i) {
numPartsOfType[i] = 0;
}
- for (int i = 0; i < num; ++i) {
+ assert(false); // First assign type[i] for rb particles
+ for (int i = 0; i < num+num_rb_attached_particles; ++i) {
++numPartsOfType[type[i]];
}
@@ -679,7 +707,7 @@ void Configuration::copyToCUDA() {
gpuErrchk(cudaMalloc(&sys_d, sizeof(BaseGrid)));
gpuErrchk(cudaMemcpyAsync(sys_d, sys, sizeof(BaseGrid), cudaMemcpyHostToDevice));
/*gpuErrchk(cudaMalloc(&type_d, sizeof(int) * num * simNum));
- gpuErrchk(cudaMemcpyAsync(type_d, type, sizeof(int) * num * simNum, cudaMemcpyHostToDevice));
+ gpuErrchk(cudaMemcpyAsync(type_d, type, sizeof(int+num_rb_attached_particles) * num * simNum, cudaMemcpyHostToDevice));
if (numBonds > 0) {
// bonds_d
@@ -698,7 +726,7 @@ void Configuration::copyToCUDA() {
cudaMemcpyHostToDevice));
// excludeMap_d
- gpuErrchk(cudaMalloc(&excludeMap_d, sizeof(int2) * num));
+ gpuErrchk(cudaMalloc(&excludeMap_d, sizeof(int2) * (num));
gpuErrchk(cudaMemcpyAsync(excludeMap_d, excludeMap, sizeof(int2) * num,
cudaMemcpyHostToDevice));
}
@@ -1117,6 +1145,8 @@ int Configuration::readParameters(const char * config_file) {
else if (param == String("rotDamping"))
rigidBody[currRB].rotDamping = stringToVector3( value );
+ else if (param == String("attachedParticles"))
+ rigidBody[currRB].append_attached_particle_file(value);
else if (param == String("densityGrid"))
rigidBody[currRB].addDensityGrid(value);
else if (param == String("potentialGrid"))
@@ -1275,6 +1305,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
if (!found) {
printf("ERROR: Number of particles not specified in config file.\n");
exit(1);
@@ -1454,9 +1485,11 @@ void Configuration::readGroups() {
std::vector<int> tmp;
for (int i=1; i < numTokens; ++i) {
const int ai = atoi(tokenList[i].val());
- if (ai >= num) {
+ if (ai >= num+num_rb_attached_particles) {
printf("Error: Attempted to include invalid particle in group: %s\n", line);
exit(-1);
+ } else if (ai >= num) {
+ printf("WARNING: including RB particles in group with line: %s\n", line);
}
tmp.push_back( ai );
}
@@ -1514,7 +1547,8 @@ void Configuration::readBonds() {
continue;
}
- if (ind1 < 0 || ind1 >= num+numGroupSites || ind2 < 0 || ind2 >= num+numGroupSites) {
+ if (ind1 < 0 || ind1 >= num+num_rb_attached_particles+numGroupSites ||
+ ind2 < 0 || ind2 >= num+num_rb_attached_particles+numGroupSites) {
printf("ERROR: Bond file line '%s' includes invalid index\n", line);
exit(1);
}
@@ -1562,8 +1596,8 @@ void Configuration::readBonds() {
* 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
*/
- bondMap = new int2[num+numGroupSites];
- for (int i = 0; i < num+numGroupSites; i++) {
+ bondMap = new int2[num+num_rb_attached_particles+numGroupSites];
+ for (int i = 0; i < num+num_rb_attached_particles+numGroupSites; i++) {
bondMap[i].x = -1;
bondMap[i].y = -1;
}
@@ -1621,8 +1655,8 @@ void Configuration::readExcludes()
}
}
void Configuration::addExclusion(int ind1, int ind2) {
- if (ind1 >= num || ind2 >= num) {
- printf("WARNING: Attempted to add an exclusion for an out-of-range particle index (%d or %d >= %d).\n", ind1, ind2, num);
+ if (ind1 >= num+num_rb_attached_particles || ind2 >= num+num_rb_attached_particles) {
+ printf("WARNING: Attempted to add an exclusion for an out-of-range particle index (%d or %d >= %d).\n", ind1, ind2, num+num_rb_attached_particles);
return;
}
@@ -1665,8 +1699,8 @@ void Configuration::buildExcludeMap() {
* excludeMap[i].x is the index in the excludes array where the ith particle's excludes begin
* excludeMap[i].y is the index in the excludes array where the ith particle's excludes end
*/
- excludeMap = new int2[num];
- for (int i = 0; i < num; i++) {
+ excludeMap = new int2[num+num_rb_attached_particles];
+ for (int i = 0; i < num+num_rb_attached_particles; i++) {
excludeMap[i].x = -1;
excludeMap[i].y = -1;
}
@@ -1675,7 +1709,7 @@ void Configuration::buildExcludeMap() {
for (int i = 0; i < numExcludes; i++) {
if (excludes[i].ind1 != currPart) {
currPart = excludes[i].ind1;
- assert(currPart < num);
+ assert(currPart < num+num_rb_attached_particles);
excludeMap[currPart].x = i;
if (lastPart >= 0)
excludeMap[lastPart].y = i;
@@ -1724,7 +1758,7 @@ void Configuration::readAngles() {
int ind3 = atoi(tokenList[3].val());
String file_name = tokenList[4];
//printf("file_name %s\n", file_name.val());
- if (ind1 >= num+numGroupSites or ind2 >= num+numGroupSites or ind3 >= num+numGroupSites)
+ if (ind1 >= num+num_rb_attached_particles+numGroupSites or ind2 >= num+num_rb_attached_particles+numGroupSites or ind3 >= num+num_rb_attached_particles+numGroupSites)
continue;
if (numAngles >= capacity) {
@@ -1785,8 +1819,10 @@ void Configuration::readDihedrals() {
int ind4 = atoi(tokenList[4].val());
String file_name = tokenList[5];
//printf("file_name %s\n", file_name.val());
- if (ind1 >= num+numGroupSites or ind2 >= num+numGroupSites
- or ind3 >= num+numGroupSites or ind4 >= num+numGroupSites)
+ if (ind1 >= num+num_rb_attached_particles+numGroupSites or
+ ind2 >= num+num_rb_attached_particles+numGroupSites or
+ ind3 >= num+num_rb_attached_particles+numGroupSites or
+ ind4 >= num+num_rb_attached_particles+numGroupSites)
continue;
if (numDihedrals >= capacity) {
@@ -1845,7 +1881,7 @@ void Configuration::readRestraints() {
float y0 = (float) strtod(tokenList[4].val(), NULL);
float z0 = (float) strtod(tokenList[5].val(), NULL);
- if (id >= num+numGroupSites) continue;
+ if (id >= num + num_rb_attached_particles + numGroupSites) continue;
if (numRestraints >= capacity) {
Restraint* temp = restraints;
diff --git a/src/Configuration.h b/src/Configuration.h
index fc70e2a..8bc023d 100644
--- a/src/Configuration.h
+++ b/src/Configuration.h
@@ -65,6 +65,7 @@ class Configuration {
void readDihedrals();
void readRestraints();
+
bool readTableFile(const String& value, int currTab);
bool readBondFile(const String& value, int currBond);
bool readAngleFile(const String& value, int currAngle);
@@ -89,6 +90,13 @@ public:
Configuration(const char * config_file, int simNum = 0, bool debug=false);
~Configuration();
+ int find_particle_type(const char* s) const {
+ for (int j = 0; j < numParts; j++)
+ if (s == part[j].name)
+ return j;
+ return -1;
+ }
+
void copyToCUDA();
// Output variables
@@ -122,7 +130,8 @@ public:
Bond* bonds;
int numCap; // max number of particles
int num; // current number of particles
- Vector3* pos; // position of each particle
+ int num_rb_attached_particles;
+ Vector3* pos; // position of each particle
Vector3* momentum; //momentum of each brownian particles Han-Yi Chou
Vector3 COM_Velocity; //center of mass velocity Han-Yi Chou
int* type; // type of each particle
diff --git a/src/RigidBodyType.cu b/src/RigidBodyType.cu
index 969a16a..8f7739c 100644
--- a/src/RigidBodyType.cu
+++ b/src/RigidBodyType.cu
@@ -81,6 +81,55 @@ void RigidBodyType::setDampingCoeffs(float timestep) { /* MUST ONLY BE CALLED ON
}
+void RigidBodyType::attach_particles() {
+ for (const auto& filename: attached_particle_files) {
+ const size_t line_char = 256;
+ FILE* inp = fopen(filename.val(), "r");
+ char line[line_char];
+
+ // If the particle file cannot be found, exit the program
+ if (inp == NULL) {
+ printf("ERROR: Could not open `%s'.\n", filename.val());
+ fclose(inp);
+ exit(1);
+ }
+
+ // Get and process all lines of input
+ while (fgets(line, line_char, inp) != NULL) {
+ // Lines in the particle file that begin with # are comments
+ if (line[0] == '#') continue;
+
+ String s(line);
+ int numTokens = s.tokenCount();
+
+ // Break the line down into pieces (tokens) so we can process them individually
+ String* tokenList = new String[numTokens];
+ s.tokenize(tokenList);
+
+ // Legitimate GROUP input lines have at least 3 tokens:
+ // Particle_type | x | y | z
+ // A line without exactly six tokens should be discarded.
+ if (numTokens != 4) {
+ printf("Error: Invalid attached particle file line: %s\n", line);
+ fclose(inp);
+ exit(-1);
+ }
+
+ // Make sure the index of this particle is unique.
+ // NOTE: The particle list is sorted by index.
+ int type_idx = conf->find_particle_type( tokenList[0].val() );
+ if (type_idx < 0) {
+ printf("Error: Unrecognized particle type: %s\n", line);
+ fclose(inp);
+ exit(-1);
+ }
+ attached_particle_types.push_back( type_idx );
+ attached_particle_positions.push_back( Vector3(atof(tokenList[1].val()), atof(tokenList[2].val()), atof(tokenList[3].val())) );
+ }
+ fclose(inp);
+ }
+}
+
void RigidBodyType::addGrid(String s, std::vector<String> &keys, std::vector<String> &files) {
// tokenize and return
int numTokens = s.tokenCount();
diff --git a/src/RigidBodyType.h b/src/RigidBodyType.h
index 575a313..7a6edcc 100644
--- a/src/RigidBodyType.h
+++ b/src/RigidBodyType.h
@@ -39,6 +39,10 @@ public:
/* RigidBodyType& operator=(const RigidBodyType& src); */
void copyGridsToDevice();
+ 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() ;}
+
void addPotentialGrid(String s);
void addDensityGrid(String s);
void addPMF(String s);
@@ -54,6 +58,10 @@ public:
String name;
private:
const Configuration* conf;
+ std::vector<String> attached_particle_files;
+ std::vector<int>attached_particle_types;
+ std::vector<Vector3>attached_particle_positions;
+
public:
int num; // number of particles of this type
--
GitLab