diff --git a/LICENSE b/LICENSE
deleted file mode 100644
index fd81b3557bf18090428847ff438040c8e2248021..0000000000000000000000000000000000000000
--- a/LICENSE
+++ /dev/null
@@ -1,35 +0,0 @@
-University of Illinois Open Source License
-Copyright 2010-2016 Aksimentiev Group,
-All rights reserved.
-
-Developed by: Aksimentiev Group
-University of Illinois at Urbana-Champaign
-http://bionano.phsyics.illinois.edu
-
-Permission is hereby granted, free of charge, to any person obtaining a copy of
-this software and associated documentation files (the Software), to deal with
-the Software without restriction, including without limitation the rights to
-use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
-of the Software, and to permit persons to whom the Software is furnished to
-do so, subject to the following conditions:
-
-- Redistributions of source code must retain the above copyright notice,
-this list of conditions and the following disclaimers.
-
-- Redistributions in binary form must reproduce the above copyright notice,
-this list of conditions and the following disclaimers in the documentation
-and/or other materials provided with the distribution.
-
-- Neither the names of the Aksimentiev Group, University of Illinois at
-Urbana-Champaign, nor the names of its contributors may be used to endorse or
-promote products derived from this Software without specific prior written
-permission.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
-THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
-OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
-ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-OTHER DEALINGS WITH THE SOFTWARE.
-
diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu
index 95b5b23c57bf164ad552462caf666bb9ec915609..4bc38245b85d6e6b49ec8a512fdfec49f5f8bc17 100644
--- a/src/ComputeForce.cu
+++ b/src/ComputeForce.cu
@@ -22,7 +22,7 @@ void runSort(int2 *d1, int *d2, float *key,
int2 *scratch1, int *scratch2, float *scratchKey,
unsigned int count);
-ComputeForce::ComputeForce(const Configuration& c, const int numReplicas) :
+ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) :
num(c.num), numParts(c.numParts), sys(c.sys), switchStart(c.switchStart),
switchLen(c.switchLen), electricConst(c.coulombConst),
cutoff2((c.switchLen + c.switchStart) * (c.switchLen + c.switchStart)),
@@ -30,25 +30,9 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas) :
numBonds(c.numBonds), numTabBondFiles(c.numTabBondFiles),
numExcludes(c.numExcludes), numAngles(c.numAngles),
numTabAngleFiles(c.numTabAngleFiles), numDihedrals(c.numDihedrals),
- numTabDihedralFiles(c.numTabDihedralFiles), numReplicas(numReplicas) {
-
-// int num, const BrownianParticleType part[],
-// int numParts, const BaseGrid* g, float switchStart,
-// float switchLen, float electricConst,
-// int fullLongRange, int numBonds, int numTabBondFiles,
-// int numExcludes, int numAngles, int numTabAngleFiles,
-// int numDihedrals, int numTabDihedralFiles,
-// float pairlistDistance, int numReplicas) :
-// num(num), numParts(numParts), sys(g), switchStart(switchStart),
-// switchLen(switchLen), electricConst(electricConst),
-// cutoff2((switchLen + switchStart) * (switchLen + switchStart)),
-// decomp(g->getBox(), g->getOrigin(), switchStart + switchLen, numReplicas),
-// numBonds(numBonds), numTabBondFiles(numTabBondFiles),
-// numExcludes(numExcludes), numAngles(numAngles),
-// numTabAngleFiles(numTabAngleFiles), numDihedrals(numDihedrals),
-// numTabDihedralFiles(numTabDihedralFiles), numReplicas(numReplicas) {
- // Allocate the parameter tables.
+ numTabDihedralFiles(c.numTabDihedralFiles), numRestraints(c.numRestraints), numReplicas(numReplicas) {
+ // Allocate the parameter tables.
decomp_d = NULL;
pairlistdist2 = (sqrt(cutoff2) + c.pairlistDistance);
@@ -187,6 +171,11 @@ ComputeForce::~ComputeForce() {
gpuErrchk( cudaFree(excludes_d) );
gpuErrchk( cudaFree(excludeMap_d) );
}
+ if (numRestraints > 0) {
+ gpuErrchk( cudaFree(restraintIds_d) );
+ gpuErrchk( cudaFree(restraintLocs_d) );
+ gpuErrchk( cudaFree(restraintSprings_d) );
+ }
}
gpuErrchk(cudaFree(numPairs_d));
@@ -630,6 +619,9 @@ float ComputeForce::computeTabulated(bool get_energy) {
if (dihedralList_d != NULL && tableDihedral_d != NULL)
computeTabulatedDihedrals<<<nb, numThreads>>>(forceInternal_d, pos_d, sys_d, numDihedrals*numReplicas, dihedralList_d, dihedralPotList_d, tableDihedral_d);
+ if (restraintIds_d != NULL )
+ computeHarmonicRestraints<<<1, numThreads>>>(forceInternal_d, pos_d, sys_d, numRestraints*numReplicas, restraintIds_d, restraintLocs_d, restraintSprings_d);
+
// Calculate the energy based on the array created by the kernel
// TODO: return energy
@@ -689,7 +681,7 @@ void ComputeForce::setForceInternalOnDevice(Vector3* f) {
}
-void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals)
+void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints)
{
// type_d
gpuErrchk(cudaMalloc(&type_d, sizeof(int) * num * simNum));
@@ -735,6 +727,28 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
cudaMemcpyHostToDevice));
}
+ if (numRestraints > 0) {
+ int restraintIds[numRestraints];
+ Vector3 restraintLocs[numRestraints];
+ float restraintSprings[numRestraints];
+ for (int i = 0; i < numRestraints; ++i) {
+ restraintIds[i] = restraints[i].id;
+ restraintLocs[i] = restraints[i].r0;
+ restraintSprings[i] = restraints[i].k;
+ }
+
+ gpuErrchk(cudaMalloc(&restraintIds_d, sizeof(int) * numRestraints));
+ gpuErrchk(cudaMalloc(&restraintLocs_d, sizeof(Vector3) * numRestraints));
+ gpuErrchk(cudaMalloc(&restraintSprings_d, sizeof(float) * numRestraints));
+
+ gpuErrchk(cudaMemcpyAsync(restraintIds_d, restraintIds,
+ sizeof(int) * numRestraints, cudaMemcpyHostToDevice));
+ gpuErrchk(cudaMemcpyAsync(restraintLocs_d, restraintLocs,
+ sizeof(Vector3) * numRestraints, cudaMemcpyHostToDevice));
+ gpuErrchk(cudaMemcpyAsync(restraintSprings_d, restraintSprings,
+ sizeof(float) * numRestraints, cudaMemcpyHostToDevice));
+ }
+
gpuErrchk(cudaDeviceSynchronize());
}
diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh
index 3291dab23c8c8166c448c7af16644190ebbec3d2..64ee8b1148dff16463eec3738f4791f152c21e5a 100644
--- a/src/ComputeForce.cuh
+++ b/src/ComputeForce.cuh
@@ -760,3 +760,18 @@ void computeTabulatedDihedrals(Vector3* force, const Vector3* __restrict__ pos,
// }
}
}
+
+__global__
+void computeHarmonicRestraints(Vector3* force, const Vector3* __restrict__ pos,
+ const BaseGrid* __restrict__ sys,
+ int numRestraints, const int* const __restrict__ particleId,
+ const Vector3* __restrict__ r0, const float* __restrict__ k) {
+
+ // Loop over ALL dihedrals in ALL replicas
+ for (int i = threadIdx.x+blockIdx.x*blockDim.x; i < numRestraints; i+=blockDim.x*gridDim.x) {
+ const int& id = particleId[i];
+ const Vector3 dr = sys->wrapDiff(pos[id]-r0[i]);
+ Vector3 f = -k[i]*dr;
+ atomicAdd( &force[ id ], f );
+ }
+}
diff --git a/src/ComputeForce.h b/src/ComputeForce.h
index fb103462e16c3f6c72661968944e56ee4574c37a..623245d4612874fadd5f4c3cef7b854768d873f9 100644
--- a/src/ComputeForce.h
+++ b/src/ComputeForce.h
@@ -15,11 +15,14 @@
#include "BaseGrid.h"
#include "BrownianParticleType.h"
#include "CellDecomposition.h"
-#include "JamesBond.h"
-#include "TabulatedPotential.h"
+
+// Simple classes
+#include "Restraint.h"
#include "useful.h"
#include "Exclude.h"
#include "Angle.h"
+#include "JamesBond.h"
+#include "TabulatedPotential.h"
#include "TabulatedAngle.h"
#include "TabulatedDihedral.h"
@@ -30,6 +33,7 @@
const unsigned int NUM_THREADS = 256;
+// Configuration
class Configuration;
class ComputeForce {
@@ -70,7 +74,7 @@ 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);
+ void copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints);
// void createBondList(int3 *bondList);
void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList);
@@ -185,16 +189,26 @@ private:
Vector3* pos_d;
Vector3* forceInternal_d;
int* type_d;
+
Bond* bonds_d;
int2* bondMap_d;
+
Exclude* excludes_d;
int2* excludeMap_d;
+
Angle* angles_d;
Dihedral* dihedrals_d;
+
int3* bondList_d;
int4* angleList_d;
int4* dihedralList_d;
int* dihedralPotList_d;
+
+ int numRestraints;
+ int* restraintIds_d;
+ Vector3* restraintLocs_d;
+ float* restraintSprings_d;
+
};
#endif
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index 6a090d23b3579ecfbd58f5ccee00fac910848265..750356033869a0b7d261801e4e0d86a8aeb2efaf 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -1,4 +1,8 @@
#include "Configuration.h"
+#include "Angle.h"
+#include "Dihedral.h"
+#include "Restraint.h"
+
#include <cassert>
@@ -114,6 +118,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
if (readExcludesFromFile) readExcludes();
if (readAnglesFromFile) readAngles();
if (readDihedralsFromFile) readDihedrals();
+ if (readRestraintsFromFile) readRestraints();
kT = temperature * 0.0019872065f; // `units "k K" "kcal_mol"`
if (temperatureGridFile.length() != 0) {
@@ -581,6 +586,10 @@ void Configuration::setDefaults() {
dihedrals = NULL;
numTabDihedralFiles = 0;
+ readRestraintsFromFile = false;
+ numRestraints = 0;
+ restraints = NULL;
+
// Hidden parameters
// Might be parameters later
numCapFactor = 5;
@@ -733,7 +742,7 @@ int Configuration::readParameters(const char * config_file) {
numTabBondFiles++;
} else if (param == String("inputParticles")) {
if (readPartsFromFile) {
- printf("WARNING: More than one particle file specified. Discarding new file.\n");
+ printf("WARNING: More than one particle file specified. Ignoring new file.\n");
} else {
partFile = value;
readPartsFromFile = true;
@@ -741,14 +750,14 @@ int Configuration::readParameters(const char * config_file) {
}
} else if (param == String("inputBonds")) {
if (readBondsFromFile) {
- printf("WARNING: More than one bond file specified. Discarding new bond file.\n");
+ printf("WARNING: More than one bond file specified. Ignoring new bond file.\n");
} else {
bondFile = value;
readBondsFromFile = true;
}
} else if (param == String("inputExcludes")) {
if (readExcludesFromFile) {
- printf("WARNING: More than one exclude file specified. Discarding new exclude file.\n");
+ printf("WARNING: More than one exclude file specified. Ignoring new exclude file.\n");
} else {
excludeFile = value;
readExcludesFromFile = true;
@@ -757,7 +766,7 @@ int Configuration::readParameters(const char * config_file) {
excludeRule = value;
} else if (param == String("inputAngles")) {
if (readAnglesFromFile) {
- printf("WARNING: More than one angle file specified. Discarding new angle file.\n");
+ printf("WARNING: More than one angle file specified. Ignoring new angle file.\n");
} else {
angleFile = value;
readAnglesFromFile = true;
@@ -775,7 +784,7 @@ int Configuration::readParameters(const char * config_file) {
numTabAngleFiles++;
} else if (param == String("inputDihedrals")) {
if (readDihedralsFromFile) {
- printf("WARNING: More than one dihedral file specified. Discarding new dihedral file.\n");
+ printf("WARNING: More than one dihedral file specified. Ignoring new dihedral file.\n");
} else {
dihedralFile = value;
readDihedralsFromFile = true;
@@ -791,6 +800,13 @@ int Configuration::readParameters(const char * config_file) {
}
if (readDihedralFile(value, ++currDihedral))
numTabDihedralFiles++;
+ } else if (param == String("inputRestraints")) {
+ if (readRestraintsFromFile) {
+ printf("WARNING: More than one restraint file specified. Ignoring new restraint file.\n");
+ } else {
+ restraintFile = value;
+ readRestraintsFromFile = true;
+ }
} else if (param == String("gridFileScale")) {
partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
} else if (param == String("rigidBodyPotential")) {
@@ -1369,6 +1385,62 @@ void Configuration::readDihedrals() {
// dihedrals[i].print();
}
+void Configuration::readRestraints() {
+ FILE* inp = fopen(restraintFile.val(), "r");
+ char line[256];
+ int capacity = 16;
+ numRestraints = 0;
+ restraints = new Restraint[capacity];
+
+ // If the restraint file cannot be found, exit the program
+ if (inp == NULL) {
+ printf("WARNING: Could not open `%s'.\n", restraintFile.val());
+ printf(" This simulation will not use restraints.\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);
+
+ // inputs have 6 tokens
+ // RESTRAINT | INDEX1 | k | x0 | y0 | z0
+ if (numTokens != 6) {
+ printf("WARNING: Invalid restraint input line: %s\n", line);
+ continue;
+ }
+
+ // Discard any empty line
+ if (tokenList == NULL) continue;
+
+ int id = atoi(tokenList[1].val());
+ float k = (float) strtod(tokenList[2].val(), NULL);
+ float x0 = (float) strtod(tokenList[3].val(), NULL);
+ float y0 = (float) strtod(tokenList[4].val(), NULL);
+ float z0 = (float) strtod(tokenList[5].val(), NULL);
+
+ if (id >= num) continue;
+
+ if (numRestraints >= capacity) {
+ Restraint* temp = restraints;
+ capacity *= 2;
+ restraints = new Restraint[capacity];
+ for (int i = 0; i < numRestraints; ++i)
+ restraints[i] = temp[i];
+ delete[] temp;
+ }
+
+ Restraint tmp(id, Vector3(x0,y0,z0), k);
+ restraints[numRestraints++] = tmp;
+ delete[] tokenList;
+ }
+ // std::sort(restraints, restraints + numRestraints, compare());
+}
+
+
void Configuration::populate() {
int pn = 0;
int p = 0;
diff --git a/src/Configuration.h b/src/Configuration.h
index 24091c993120c2c352a8a9b441b4cf7e27b2744a..64f4acf1407748f3585cd6c0124ca76a308d34e7 100644
--- a/src/Configuration.h
+++ b/src/Configuration.h
@@ -20,13 +20,11 @@
#include "BaseGrid.h"
#include "OverlordGrid.h"
#include "ComputeForce.h"
-#include "Angle.h"
#include "Reader.h"
#include "TrajectoryWriter.h"
#include "TabulatedPotential.h"
#include "TabulatedAngle.h"
#include "GPUManager.h"
-#include "Dihedral.h"
#include "RigidBodyType.h"
#include "RigidBody.h"
@@ -37,9 +35,14 @@
// Energy: kcal/mol (6.947694e-24 kJ)
// Temperature: Kelvin
// Time: nanoseconds
-// Length: nanometers
+// Length: angstroms
+// Forward declerations
+class Angle;
+class Dihedral;
+struct Restraint;
+
class Configuration {
struct compare {
bool operator()(const String& lhs, const String& rhs);
@@ -59,6 +62,7 @@ class Configuration {
void readBonds();
void readExcludes();
void readDihedrals();
+ void readRestraints();
bool readTableFile(const String& value, int currTab);
bool readBondFile(const String& value, int currBond);
@@ -168,17 +172,20 @@ public:
int numExcludes;
int numAngles;
int numDihedrals;
+ int numRestraints;
int* numPartsOfType;
String partFile;
String bondFile;
String excludeFile;
String angleFile;
String dihedralFile;
+ String restraintFile;
bool readPartsFromFile;
bool readBondsFromFile;
bool readExcludesFromFile;
bool readAnglesFromFile;
bool readDihedralsFromFile;
+ bool readRestraintsFromFile;
String* partGridFile;
float* partGridFileScale;
std::vector< std::vector<String> > partRigidBodyGrid;
@@ -208,6 +215,8 @@ public:
String* dihedralTableFile;
int numTabDihedralFiles;
+ Restraint* restraints;
+
// RigidBody parameters.
RigidBodyType* rigidBody;
int numRigidTypes;
diff --git a/src/GrandBrownTown.cu b/src/GrandBrownTown.cu
index 41ff55b1317c4c9f90695d69a0df58dee4a4b6c0..d8321c84a58ced07b5db80190382baa95862cd9c 100644
--- a/src/GrandBrownTown.cu
+++ b/src/GrandBrownTown.cu
@@ -184,7 +184,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);
+ internal -> copyToCUDA(c.simNum, c.type, c.bonds, c.bondMap, c.excludes, c.excludeMap, c.angles, c.dihedrals, c.restraints);
// TODO: check for duplicate potentials
if (c.tabulatedPotential) {
@@ -200,7 +200,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
}
}
}
- printf("Using %d non-bonded exclusions\n",numExcludes/2);
+ printf("Using %d non-bonded exclusions\n",c.numExcludes/2);
if (c.readBondsFromFile) {
printf("Loading %d tabulated bond potentials...\n", numTabBondFiles);