diff --git a/src/Angle.cu b/src/Angle.cu
index acc2a09dc389f96861b2ca3b1ebc5d7acc0f4d2b..28facb728927e782b5cb50697fdb2a80618557aa 100644
--- a/src/Angle.cu
+++ b/src/Angle.cu
@@ -15,3 +15,15 @@ String Angle::toString()
return String("ANGLE ") + ind1 + " " + ind2 + " " + ind3 + " " + fileName;
}
+void BondAngle::print()
+{
+// printf("about to print fileName %p\n", fileName.val());
+// fileName.print();
+ printf("BONDANGLE (%d %d %d) %s; %s; %s\n", ind1, ind2, ind3, angleFileName.val(), bondFileName1.val(), bondFileName2.val());
+}
+
+String BondAngle::toString()
+{
+ return String("BONDANGLE ") + ind1 + " " + ind2 + " " + ind3 + " " + angleFileName + " " + bondFileName1 + " " + bondFileName2;
+}
+
diff --git a/src/Angle.h b/src/Angle.h
index a4e9542c0878b6b195092bfef5b3ea8a4633f1df..a8794e89c099d9b7d18d294c7065246586685ebf 100644
--- a/src/Angle.h
+++ b/src/Angle.h
@@ -59,4 +59,59 @@ public:
String toString();
void print();
};
+
+// TODO consolidate with Angle using inheritence
+class BondAngle
+{
+public:
+ BondAngle() {}
+ BondAngle(int ind1, int ind2, int ind3, String angleFileName, String bondFileName1, String bondFileName2) :
+ ind1(ind1), ind2(ind2), ind3(ind3), angleFileName(angleFileName), bondFileName1(bondFileName1), bondFileName2(bondFileName2), tabFileIndex1(-1), tabFileIndex2(-1), tabFileIndex3(-1) { }
+
+ int ind1, ind2, ind3;
+
+ String angleFileName;
+ String bondFileName1;
+ String bondFileName2;
+ // tabFileIndex will be assigned after ComputeForce loads the
+ // TabulatedAnglePotentials. The tabefileIndex is used by ComputeForce to
+ // discern which TabulatedAnglePotential this Angle uses.
+ int tabFileIndex1;
+ int tabFileIndex2;
+ int tabFileIndex3;
+
+ inline BondAngle(const BondAngle& a) : ind1(a.ind1), ind2(a.ind2), ind3(a.ind3),
+ angleFileName(a.angleFileName),
+ bondFileName1(a.bondFileName1),
+ bondFileName2(a.bondFileName2),
+ tabFileIndex1(a.tabFileIndex1),
+ tabFileIndex2(a.tabFileIndex2),
+ tabFileIndex3(a.tabFileIndex3) { }
+
+ HOST DEVICE inline float calcAngle(Vector3* pos, BaseGrid* sys) {
+ const Vector3& posa = pos[ind1];
+ const Vector3& posb = pos[ind2];
+ const Vector3& posc = pos[ind3];
+ const float distab = sys->wrapDiff(posa - posb).length();
+ const float distbc = sys->wrapDiff(posb - posc).length();
+ const float distac = sys->wrapDiff(posc - posa).length();
+ float cos = (distbc * distbc + distab * distab - distac * distac)
+ / (2.0f * distbc * distab);
+ if (cos < -1.0f) cos = -1.0f;
+ else if (cos > 1.0f) cos = 1.0f;
+ float angle = acos(cos);
+ return angle;
+ }
+
+ HOST DEVICE inline int getIndex(int index) {
+ if (index == ind1) return 1;
+ if (index == ind2) return 2;
+ if (index == ind3) return 3;
+ return -1;
+ }
+
+ String toString();
+ void print();
+};
+
#endif
diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu
index 7ccf22acdab6420843b3266085cef3dc5a51f005..c66fb14b3bd4f76ddc2c5651485f3bd4d0064bf0 100644
--- a/src/ComputeForce.cu
+++ b/src/ComputeForce.cu
@@ -49,6 +49,7 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) :
numExcludes(c.numExcludes), numAngles(c.numAngles),
numTabAngleFiles(c.numTabAngleFiles), numDihedrals(c.numDihedrals),
numTabDihedralFiles(c.numTabDihedralFiles), numRestraints(c.numRestraints),
+ numBondAngles(c.numBondAngles),
numGroupSites(c.numGroupSites),
numReplicas(numReplicas) {
@@ -233,6 +234,7 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) :
}
restraintIds_d = NULL;
+ bondAngleList_d = NULL;
//Calculate the number of blocks the grid should contain
gridSize = (num+num_rb_attached_particles) / NUM_THREADS + 1;
@@ -401,7 +403,7 @@ bool ComputeForce::addTabulatedPotential(String fileName, int type0, int type1)
return true;
}
-bool ComputeForce::addBondPotential(String fileName, int ind, Bond bonds[])
+bool ComputeForce::addBondPotential(String fileName, int ind, Bond bonds[], BondAngle bondAngles[])
{
// TODO: see if tableBond_addr can be removed
if (tableBond[ind] != NULL) {
@@ -415,6 +417,14 @@ bool ComputeForce::addBondPotential(String fileName, int ind, Bond bonds[])
if (bonds[i].fileName == fileName)
bonds[i].tabFileIndex = ind;
+ for (int i = 0; i < numBondAngles; i++)
+ {
+ if (bondAngles[i].bondFileName1 == fileName)
+ bondAngles[i].tabFileIndex2 = ind;
+ if (bondAngles[i].bondFileName2 == fileName)
+ bondAngles[i].tabFileIndex3 = ind;
+ }
+
gpuErrchk(cudaMemcpyAsync(bonds_d, bonds, sizeof(Bond) * numBonds, cudaMemcpyHostToDevice));
tableBond_addr[ind] = tableBond[ind]->copy_to_cuda();
@@ -423,7 +433,7 @@ bool ComputeForce::addBondPotential(String fileName, int ind, Bond bonds[])
return true;
}
-bool ComputeForce::addAnglePotential(String fileName, int ind, Angle* angles) {
+bool ComputeForce::addAnglePotential(String fileName, int ind, Angle* angles, BondAngle* bondAngles) {
if (tableAngle[ind] != NULL) {
delete tableAngle[ind];
gpuErrchk(cudaFree(tableAngle_addr[ind]));
@@ -452,6 +462,10 @@ bool ComputeForce::addAnglePotential(String fileName, int ind, Angle* angles) {
if (angles[i].fileName == fileName)
angles[i].tabFileIndex = ind;
+ for (int i = 0; i < numBondAngles; i++)
+ if (bondAngles[i].angleFileName == fileName)
+ bondAngles[i].tabFileIndex1 = ind;
+
gpuErrchk(cudaMemcpy(angles_d, angles, sizeof(Angle) * numAngles,
cudaMemcpyHostToDevice));
return true;
@@ -714,6 +728,11 @@ float ComputeForce::computeTabulated(bool get_energy) {
//Mlog: the commented function doesn't use bondList, uncomment for testing.
//if(bondMap_d != NULL && tableBond_d != NULL)
+ if(bondAngleList_d != NULL && tableBond_d != NULL && tableAngle_d != NULL)
+ {
+ computeTabulatedBondAngles <<<nb, numThreads, 0, gpuman.get_next_stream()>>> ( forceInternal_d, pos_d, sys_d, numReplicas*numBondAngles, bondAngleList_d, tableAngle_d, tableBond_d, energies_d, get_energy);
+ }
+
if(bondList_d != NULL && tableBond_d != NULL)
{
@@ -865,7 +884,7 @@ void ComputeForce::setForceInternalOnDevice(Vector3* f) {
gpuErrchk(cudaMemcpy(forceInternal_d[0], f, sizeof(Vector3) * tot_num, cudaMemcpyHostToDevice));
}
-void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints)
+void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints, const BondAngle* const bondAngles)
{
assert(simNum == numReplicas); // Not sure why we have both of these things
int tot_num_with_rb = (num+num_rb_attached_particles) * simNum;
@@ -936,6 +955,12 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
sizeof(float) * numRestraints, cudaMemcpyHostToDevice));
}
+ if (numBondAngles > 0) {
+ gpuErrchk(cudaMalloc(&bondAngles_d, sizeof(BondAngle) * numBondAngles));
+ gpuErrchk(cudaMemcpyAsync(bondAngles_d, bondAngles, sizeof(BondAngle) * numBondAngles,
+ cudaMemcpyHostToDevice));
+ }
+
gpuErrchk(cudaDeviceSynchronize());
}
@@ -954,7 +979,7 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
// }
// }
-void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList) {
+void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int2* bondAngleList) {
size_t size;
@@ -980,4 +1005,11 @@ void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *d
gpuErrchk( cudaMalloc( &dihedralPotList_d, size ) );
gpuErrchk( cudaMemcpyAsync( dihedralPotList_d, dihedralPotList, size, cudaMemcpyHostToDevice) );
}
+
+ if (numBondAngles > 0) {
+ size = 3*numBondAngles * numReplicas * sizeof(int2);
+ gpuErrchk( cudaMalloc( &bondAngleList_d, size ) );
+ gpuErrchk( cudaMemcpyAsync( bondAngleList_d, bondAngleList, size, cudaMemcpyHostToDevice) );
+ }
+
}
diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh
index eb619c91cf054511913d6e0bb9756b12765f2bc1..16de7427b2678600fe554002217cfe4503598317 100644
--- a/src/ComputeForce.cuh
+++ b/src/ComputeForce.cuh
@@ -865,6 +865,28 @@ void computeTabulatedAngles(Vector3* force,
}
}
+__global__
+void computeTabulatedBondAngles(Vector3* force,
+ Vector3* __restrict__ pos,
+ BaseGrid* __restrict__ sys,
+ int numBondAngles, int2* __restrict__ bondAngleList_d, TabulatedAnglePotential** tableAngle,
+ TabulatedPotential** tableBond,
+ float* energy, bool get_energy) {
+ // Loop over ALL angles in ALL replicas
+ for (int i = threadIdx.x+blockIdx.x*blockDim.x; i<numBondAngles; i+=blockDim.x*gridDim.x) {
+ int atom1 = bondAngleList_d[3*i ].x;
+ int atom2 = bondAngleList_d[3*i ].y;
+ int atom3 = bondAngleList_d[3*i+1].x;
+
+ int angleInd = bondAngleList_d[3*i+1].y;
+ int bondInd1 = bondAngleList_d[3*i+2].x;
+ int bondInd2 = bondAngleList_d[3*i+2].y;
+
+ computeBondAngle(tableAngle[ angleInd ], tableBond[ bondInd1 ], tableBond[ bondInd2 ], sys, force, pos, atom1, atom2, atom3, energy, get_energy);
+ }
+}
+
+
__global__
void computeDihedrals(Vector3 force[], Vector3 pos[],
diff --git a/src/ComputeForce.h b/src/ComputeForce.h
index b3bc258541f9603f387221363f3c39fe7fbae881..08269d511d7aec24af27700c092dc484fbab7845 100644
--- a/src/ComputeForce.h
+++ b/src/ComputeForce.h
@@ -46,8 +46,8 @@ public:
void makeTables(const BrownianParticleType* part);
bool addTabulatedPotential(String fileName, int type0, int type1);
- bool addBondPotential(String fileName, int ind, Bond* bonds);
- bool addAnglePotential(String fileName, int ind, Angle* angles);
+ bool addBondPotential(String fileName, int ind, Bond* bonds, BondAngle* bondAngles);
+ bool addAnglePotential(String fileName, int ind, Angle* angles, BondAngle* bondAngles);
bool addDihedralPotential(String fileName, int ind, Dihedral* dihedrals);
void decompose();
@@ -75,12 +75,12 @@ public:
//MLog: new copy function to allocate memory required by ComputeForce class.
void copyToCUDA(Vector3* forceInternal, Vector3* pos);
- void copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints);
+ void copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints, const BondAngle* const bondAngles);
void copyToCUDA(Vector3* forceInternal, Vector3* pos, Vector3* mom);
void copyToCUDA(Vector3* forceInternal, Vector3* pos, Vector3* mom, float* random);
// void createBondList(int3 *bondList);
- void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList);
+ void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int2 *bondAngleList);
//MLog: because of the move of a lot of private variables, some functions get starved necessary memory access to these variables, below is a list of functions that return the specified private variable.
std::vector<Vector3*> getPos_d()
@@ -249,6 +249,10 @@ private:
Angle* angles_d;
Dihedral* dihedrals_d;
+ int numBondAngles;
+ BondAngle* bondAngles_d;
+ int2* bondAngleList_d;
+
int3* bondList_d;
int4* angleList_d;
int4* dihedralList_d;
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index 8ee8754485ade2edd9a8124c684e7a68cb3d73b4..a6f06b8e9cc1078a9b3492287420903391d94783 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -271,6 +271,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
if (readAnglesFromFile) readAngles();
if (readDihedralsFromFile) readDihedrals();
if (readRestraintsFromFile) readRestraints();
+ if (readBondAnglesFromFile) readBondAngles();
if (temperatureGridFile.length() != 0) {
printf("\nFound temperature grid file: %s\n", temperatureGridFile.val());
@@ -630,6 +631,7 @@ Configuration::~Configuration() {
if (excludeMap != NULL) delete[] excludeMap;
if (angles != NULL) delete[] angles;
if (dihedrals != NULL) delete[] dihedrals;
+ if (bondAngles != NULL) delete[] bondAngles;
delete[] numPartsOfType;
@@ -861,6 +863,10 @@ void Configuration::setDefaults() {
dihedrals = NULL;
numTabDihedralFiles = 0;
+ readBondAnglesFromFile = false;
+ numBondAngles = 0;
+ bondAngles = NULL;
+
readRestraintsFromFile = false;
numRestraints = 0;
restraints = NULL;
@@ -1122,6 +1128,13 @@ int Configuration::readParameters(const char * config_file) {
angleFile = value;
readAnglesFromFile = true;
}
+ } else if (param == String("inputBondAngles")) {
+ if (readBondAnglesFromFile) {
+ printf("WARNING: More than one angle file specified. Ignoring new angle file.\n");
+ } else {
+ bondAngleFile = value;
+ readBondAnglesFromFile = true;
+ }
} else if (param == String("tabulatedAngleFile")) {
if (numTabAngleFiles >= atfcap) {
String* temp = angleTableFile;
@@ -1605,7 +1618,7 @@ void Configuration::readBonds() {
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);
@@ -1925,6 +1938,68 @@ void Configuration::readDihedrals() {
// dihedrals[i].print();
}
+void Configuration::readBondAngles() {
+ FILE* inp = fopen(bondAngleFile.val(), "r");
+ char line[256];
+ int capacity = 256;
+ numBondAngles = 0;
+ bondAngles = new BondAngle[capacity];
+
+ // If the angle file cannot be found, exit the program
+ if (inp == NULL) {
+ printf("WARNING: Could not open `%s'.\n", bondAngleFile.val());
+ printf("This simulation will not use angles.\n");
+ return;
+ }
+
+ while(fgets(line, 256, inp)) {
+ if (line[0] == '#') continue;
+ String s(line);
+ int numTokens = s.tokenCount();
+ String* tokenList = new String[numTokens];
+ s.tokenize(tokenList);
+
+ // Legitimate BONDANGLE inputs have 7 tokens
+ // BONDANGLE | INDEX1 | INDEX2 | INDEX3 | ANGLE_FILENAME | BOND_FILENAME1 | BOND_FILENAME2
+ // Any angle input line without exactly 5 tokens should be discarded
+ if (numTokens != 7) {
+ printf("WARNING: Invalid bond_angle input line: %s\n", line);
+ continue;
+ }
+
+ // Discard any empty line
+ if (tokenList == NULL)
+ continue;
+
+ int ind1 = atoi(tokenList[1].val());
+ int ind2 = atoi(tokenList[2].val());
+ int ind3 = atoi(tokenList[3].val());
+ String file_name1 = tokenList[4];
+ String file_name2 = tokenList[5];
+ String file_name3 = tokenList[6];
+ //printf("file_name %s\n", file_name.val());
+ if (ind1 >= num or ind2 >= num or ind3 >= num)
+ continue;
+
+ if (numBondAngles >= capacity) {
+ BondAngle* temp = bondAngles;
+ capacity *= 2;
+ bondAngles = new BondAngle[capacity];
+ for (int i = 0; i < numBondAngles; i++)
+ bondAngles[i] = temp[i];
+ delete[] temp;
+ }
+
+ BondAngle a(ind1, ind2, ind3, file_name1, file_name2, file_name3);
+ bondAngles[numBondAngles++] = a;
+ delete[] tokenList;
+ }
+ std::sort(bondAngles, bondAngles + numBondAngles, compare());
+
+ // for(int i = 0; i < numAngles; i++)
+ // angles[i].print();
+}
+
void Configuration::readRestraints() {
FILE* inp = fopen(restraintFile.val(), "r");
char line[256];
@@ -2524,3 +2599,13 @@ bool Configuration::compare::operator()(const Dihedral& lhs, const Dihedral& rhs
return lhs.ind3 < rhs.ind3;
return lhs.ind4 < rhs.ind4;
}
+
+bool Configuration::compare::operator()(const BondAngle& lhs, const BondAngle& rhs) {
+ int diff = lhs.ind1 - rhs.ind1;
+ if (diff != 0)
+ return lhs.ind1 < rhs.ind1;
+ diff = lhs.ind2 - rhs.ind2;
+ if (diff != 0)
+ return lhs.ind2 < rhs.ind2;
+ return lhs.ind3 < rhs.ind3;
+}
diff --git a/src/Configuration.h b/src/Configuration.h
index 6199dd8b24bc019dd7044b2a3a6031312142f13d..66d653984943200bb06b1575a8000ea3292004c7 100644
--- a/src/Configuration.h
+++ b/src/Configuration.h
@@ -48,6 +48,7 @@ class Configuration {
bool operator()(const Exclude& lhs, const Exclude& rhs);
bool operator()(const Angle& lhs, const Angle& rhs);
bool operator()(const Dihedral& lhs, const Dihedral& rhs);
+ bool operator()(const BondAngle& lhs, const BondAngle& rhs);
};
void setDefaults();
@@ -64,6 +65,7 @@ class Configuration {
void buildExcludeMap();
void readDihedrals();
void readRestraints();
+ void readBondAngles();
bool readTableFile(const String& value, int currTab);
@@ -71,6 +73,8 @@ class Configuration {
bool readAngleFile(const String& value, int currAngle);
bool readDihedralFile(const String& value, int currDihedral);
+ bool readBondAngleFile(const String& value, const String& bondfile1, const String& bondfile2, int currBondAngle);
+
// Given the numbers of each particle, populate the type list.
void populate();
@@ -201,6 +205,7 @@ public:
int numExcludes;
int numAngles;
int numDihedrals;
+ int numBondAngles;
int numRestraints;
int* numPartsOfType;
String partFile;
@@ -209,12 +214,14 @@ public:
String angleFile;
String dihedralFile;
String restraintFile;
+ String bondAngleFile;
bool readPartsFromFile;
bool readGroupSitesFromFile;
bool readBondsFromFile;
bool readExcludesFromFile;
bool readAnglesFromFile;
bool readDihedralsFromFile;
+ bool readBondAnglesFromFile;
bool readRestraintsFromFile;
//String* partGridFile;
String **partGridFile;
@@ -254,6 +261,8 @@ public:
String* dihedralTableFile;
int numTabDihedralFiles;
+ BondAngle* bondAngles;
+
Restraint* restraints;
//Han-Yi Chou
diff --git a/src/GrandBrownTown.cu b/src/GrandBrownTown.cu
index aa5295f4271f2ac169033e139e2fa15418654562..7ea4cd27f498458cbd1e2be93bd76cd0a961ad6b 100644
--- a/src/GrandBrownTown.cu
+++ b/src/GrandBrownTown.cu
@@ -204,6 +204,8 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
partTableIndex0 = c.partTableIndex0;
partTableIndex1 = c.partTableIndex1;
+ numBondAngles = c.numBondAngles;
+
numTabBondFiles = c.numTabBondFiles;
bondMap = c.bondMap;
// TODO: bondList = c.bondList;
@@ -220,6 +222,8 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
dihedrals = c.dihedrals;
numTabDihedralFiles = c.numTabDihedralFiles;
+ bondAngles = c.bondAngles;
+
// Device parameters
//type_d = c.type_d;
part_d = c.part_d;
@@ -258,7 +262,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
internal = new ComputeForce(c, numReplicas);
//MLog: I did the other halve of the copyToCUDA function from the Configuration class here, keep an eye on any mistakes that may occur due to the location.
- internal -> copyToCUDA(c.simNum, c.type, c.bonds, c.bondMap, c.excludes, c.excludeMap, c.angles, c.dihedrals, c.restraints);
+ internal -> copyToCUDA(c.simNum, c.type, c.bonds, c.bondMap, c.excludes, c.excludeMap, c.angles, c.dihedrals, c.restraints, c.bondAngles );
if (numGroupSites > 0) init_cuda_group_sites();
// TODO: check for duplicate potentials
@@ -283,11 +287,11 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
if (bondTableFile[p].length() > 0) {
//MLog: make sure to add to all GPUs
// printf("...loading %s\n",bondTableFile[p].val());
- internal->addBondPotential(bondTableFile[p].val(), p, bonds);
+ internal->addBondPotential(bondTableFile[p].val(), p, bonds, bondAngles);
// printf("%s\n",bondTableFile[p].val());
} else {
printf("...skipping %s (\n",bondTableFile[p].val());
- internal->addBondPotential(bondTableFile[p].val(), p, bonds);
+ internal->addBondPotential(bondTableFile[p].val(), p, bonds, bondAngles);
}
}
@@ -298,7 +302,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
if (angleTableFile[p].length() > 0)
{
//MLog: make sure to do this for every GPU
- internal->addAnglePotential(angleTableFile[p].val(), p, angles);
+ internal->addAnglePotential(angleTableFile[p].val(), p, angles, bondAngles);
}
}
@@ -373,7 +377,32 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
}
}
}
- internal->copyBondedListsToGPU(bondList,angleList,dihedralList,dihedralPotList);
+
+ if (numBondAngles > 0) {
+ bondAngleList = new int2[ (numBondAngles*3) * numReplicas ];
+ for(int k = 0 ; k < numReplicas; k++) {
+ for(int i = 0; i < numBondAngles; ++i) {
+ if (bondAngles[i].tabFileIndex1 == -1) {
+ fprintf(stderr,"Error: bondanglefile '%s' was not read with tabulatedAngleFile command.\n", bondAngles[i].angleFileName.val());
+ exit(1);
+ }
+ if (bondAngles[i].tabFileIndex2 == -1) {
+ fprintf(stderr,"Error: bondanglefile1 '%s' was not read with tabulatedBondFile command.\n", bondAngles[i].bondFileName1.val());
+ exit(1);
+ }
+ if (bondAngles[i].tabFileIndex3 == -1) {
+ fprintf(stderr,"Error: bondanglefile2 '%s' was not read with tabulatedBondFile command.\n", bondAngles[i].bondFileName2.val());
+ exit(1);
+ }
+ int idx = i+k*numBondAngles;
+ bondAngleList[idx*3] = make_int2( bondAngles[i].ind1+k*num, bondAngles[i].ind2+k*num );
+ bondAngleList[idx*3+1] = make_int2( bondAngles[i].ind3+k*num, bondAngles[i].tabFileIndex1 );
+ bondAngleList[idx*3+2] = make_int2( bondAngles[i].tabFileIndex2, bondAngles[i].tabFileIndex3 );
+ }
+ }
+ }
+
+ internal->copyBondedListsToGPU(bondList,angleList,dihedralList,dihedralPotList,bondAngleList);
forceInternal = new Vector3[(num+num_rb_attached_particles+numGroupSites)*numReplicas];
if (fullLongRange != 0)
diff --git a/src/GrandBrownTown.h b/src/GrandBrownTown.h
index 56d4c5a3b28c7b765542ba87eb5efa2445768c82..c74ca2b59f243a8ff3867e3581de46f338af2b64 100644
--- a/src/GrandBrownTown.h
+++ b/src/GrandBrownTown.h
@@ -275,6 +275,10 @@ private:
int4 *dihedralList;
int *dihedralPotList;
+ int numBondAngles;
+ BondAngle* bondAngles;
+ int2 *bondAngleList;
+
//Han-Yi Chou
String particle_dynamic;
String rigidbody_dynamic;
diff --git a/src/TabulatedMethods.cuh b/src/TabulatedMethods.cuh
index d92809c99ce0206d8b4b380a6cc08abf3b692112..7965c7736c793bbd7b0cfcf7d6229d8734247f4b 100644
--- a/src/TabulatedMethods.cuh
+++ b/src/TabulatedMethods.cuh
@@ -77,6 +77,88 @@ __device__ inline void computeAngle(const TabulatedAnglePotential* __restrict__
atomicAdd( &force[k], force3 );
}
+__device__ inline void computeBondAngle(const TabulatedAnglePotential* __restrict__ a,
+ const TabulatedPotential* __restrict__ b1, const TabulatedPotential* __restrict__ b2,
+ const BaseGrid* __restrict__ sys, Vector3* force, const Vector3* __restrict__ pos,
+ const int& i, const int& j, const int& k, float* energy, bool get_energy) {
+
+ // Particle's type and position
+ Vector3 posa = pos[i];
+ Vector3 posb = pos[j];
+ Vector3 posc = pos[k];
+
+ // The vectors between each pair of particles
+ const Vector3 ab = sys->wrapDiff(posa - posb);
+ const Vector3 bc = sys->wrapDiff(posb - posc);
+ const Vector3 ac = sys->wrapDiff(posc - posa);
+
+ // Find the distance between each pair of particles
+ float distab = ab.length2();
+ float distbc = bc.length2();
+ EnergyForce fe_b1 = b1->compute(ab,distab);
+ EnergyForce fe_b2 = b2->compute(bc,distbc);
+
+ const float distac2 = ac.length2();
+
+ // Find the cosine of the angle we want - <ABC
+ float cos = (distab + distbc - distac2);
+
+ distab = 1.0f/sqrt(distab); //TODO: test other functiosn
+ distbc = 1.0f/sqrt(distbc);
+ cos *= 0.5f * distbc * distab;
+
+ // If the cosine is illegitimate, set it to 1 or -1 so that acos won't fail
+ if (cos < -1.0f) cos = -1.0f;
+ if (cos > 1.0f) cos = 1.0f;
+
+ // Find the sine while we're at it.
+
+ // Now we can use the cosine to find the actual angle (in radians)
+ float angle = acos(cos);
+
+ // transform angle to units of tabulated array index
+ angle *= a->angle_step_inv;
+
+ // tableAngle[0] stores the potential at angle_step
+ // tableAngle[1] stores the potential at angle_step * 2, etc.
+ // 'home' is the index after which 'convertedAngle' would appear if it were stored in the table
+ int home = int(floorf(angle));
+ home = (home >= a->size) ? (a->size)-1 : home;
+ //assert(home >= 0);
+ //assert(home+1 < a->size);
+
+ // // Make angle the distance from [0,1) from the first index in the potential array index
+ // angle -= home;
+
+ // Linearly interpolate the potential
+ float U0 = a->pot[home];
+ float dUdx = (a->pot[(((home+1)==(a->size)) ? (a->size)-1 : home+1)] - U0) * a->angle_step_inv;
+ float a_e = (dUdx * (angle-home)) + U0;
+ if(get_energy)
+ {
+ float e = a_e * fe_b1.e * fe_b2.e * 0.3333333333f;
+ atomicAdd( &energy[i], e);
+ atomicAdd( &energy[j], e);
+ atomicAdd( &energy[k], e);
+ }
+ float sin = sqrtf(1.0f - cos*cos);
+ dUdx /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
+
+ // Calculate the forces
+ Vector3 force1 = -(dUdx*distab) * (ab * (cos*distab) + bc * distbc); // force on particle 1
+ Vector3 force3 = (dUdx*distbc) * (bc * (cos*distbc) + ab * distab); // force on particle 3
+
+ force1 = (force1 * fe_b1.e - a_e * fe_b1.f) * fe_b2.e;
+ force3 = (force3 * fe_b2.e + a_e * fe_b2.f) * fe_b1.e;
+
+ // assert( force1.length() < 10000.0f );
+ // assert( force3.length() < 10000.0f );
+
+ atomicAdd( &force[i], force1 );
+ atomicAdd( &force[j], -(force1 + force3) );
+ atomicAdd( &force[k], force3 );
+}
+
__device__ inline void computeDihedral(const TabulatedDihedralPotential* __restrict__ d,
const BaseGrid* __restrict__ sys, Vector3* forces, const Vector3* __restrict__ pos,
diff --git a/src/TabulatedPotential.h b/src/TabulatedPotential.h
index e32068a6fa7dc0b8cbada59cb7f6ea8d40933b22..2bf31620f25c73bba9bc20f528b9ff4e3f5aadc3 100644
--- a/src/TabulatedPotential.h
+++ b/src/TabulatedPotential.h
@@ -153,7 +153,7 @@ public:
Vector3 force = (du*drInv/d)*r;
return EnergyForce(energy,force);
}
- HOST DEVICE inline Vector3 computef(Vector3 r, float d) {
+ HOST DEVICE inline Vector3 computef(const Vector3 r, float d) const {
d = sqrt(d);
// float d = r.length();
// RBTODO: precompute so that initial blocks are zero; reduce computation here