src/ComputeForce.cu

@@ -48,7 +48,8 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) :
     numBonds(c.numBonds), numTabBondFiles(c.numTabBondFiles),
     numExcludes(c.numExcludes), numAngles(c.numAngles),
     numTabAngleFiles(c.numTabAngleFiles), numDihedrals(c.numDihedrals),
-    numTabDihedralFiles(c.numTabDihedralFiles), numRestraints(c.numRestraints),
+    numTabDihedralFiles(c.numTabDihedralFiles), numVecangles(c.numVecangles),
+    numTabVecangleFiles(c.numTabVecangleFiles), numRestraints(c.numRestraints),
     numBondAngles(c.numBondAngles), numProductPotentials(c.numProductPotentials),
     numGroupSites(c.numGroupSites),
     numReplicas(numReplicas) {
@@ -143,6 +144,17 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) :
 	}
 	gpuErrchk(cudaMalloc(&tableDihedral_d, sizeof(TabulatedDihedralPotential*) * numTabDihedralFiles));
 
+	// Create the vecangle table
+	tableVecangle = new TabulatedVecanglePotential*[numTabVecangleFiles];
+	tableVecangle_addr = new TabulatedVecanglePotential*[numTabVecangleFiles];
+	vecangleList_d = NULL;
+	tableVecangle_d = NULL;
+	for (int i = 0; i < numTabVecangleFiles; i++) {
+		tableVecangle_addr[i] = NULL;
+		tableVecangle[i] = NULL;
+	}
+	gpuErrchk(cudaMalloc(&tableVecangle_d, sizeof(TabulatedVecanglePotential*) * numTabVecangleFiles));
+
 	{	// allocate device for pairlists
 		// RBTODO: select maxpairs in better way; add assertion in kernel to avoid going past this
 		const int maxPairs = MAX_NLIST_PAIRS;
@@ -303,6 +315,11 @@ ComputeForce::~ComputeForce() {
 			gpuErrchk( cudaFree(dihedralList_d) );
 			gpuErrchk( cudaFree(dihedralPotList_d) );
 		}
+		if (numVecangles > 0) {
+			gpuErrchk( cudaFree(vecangles_d) );
+			gpuErrchk( cudaFree(vecangleList_d) );
+			gpuErrchk( cudaFree(vecanglePotList_d) );
+		}
 		if (numExcludes > 0) {
 			gpuErrchk( cudaFree(excludes_d) );
 			gpuErrchk( cudaFree(excludeMap_d) );
@@ -507,6 +524,41 @@ bool ComputeForce::addDihedralPotential(String fileName, int ind, Dihedral dihed
 	t.pot = NULL;
 	return true;
 }
+bool ComputeForce::addVecanglePotential(String fileName, int ind, Vecangle vecangles[])
+{
+	for (int i = 0; i < numVecangles; i++)
+		if (vecangles[i].fileName == fileName)
+			vecangles[i].tabFileIndex = ind;
+
+	gpuErrchk(cudaMemcpyAsync(vecangles_d, vecangles, sizeof(Vecangle) * numVecangles,
+			cudaMemcpyHostToDevice));
+
+	if (tableVecangle[ind] != NULL) {
+		delete tableVecangle[ind];
+		gpuErrchk(cudaFree(tableVecangle_addr[ind]));
+		tableVecangle[ind] = NULL;
+		tableVecangle_addr[ind] = NULL;
+	}
+
+	tableVecangle[ind] = new TabulatedVecanglePotential(fileName,VECANGLE);
+	TabulatedVecanglePotential t = TabulatedVecanglePotential(*tableVecangle[ind]);
+
+	// Copy tableAngle[ind] to the device
+	float *pot;
+	int size = tableVecangle[ind]->size;
+	gpuErrchk(cudaMalloc(&pot, sizeof(float) * size));
+	gpuErrchk(cudaMemcpyAsync(pot, tableVecangle[ind]->pot,
+			sizeof(float) * size, cudaMemcpyHostToDevice));
+	t.pot = pot;
+
+	gpuErrchk(cudaMalloc(&tableVecangle_addr[ind], sizeof(TabulatedVecanglePotential)));
+	gpuErrchk(cudaMemcpyAsync(tableVecangle_addr[ind], &t,
+			sizeof(TabulatedVecanglePotential), cudaMemcpyHostToDevice));
+	gpuErrchk(cudaMemcpy(tableVecangle_d, tableVecangle_addr,
+			sizeof(TabulatedVecanglePotential*) * numTabVecangleFiles, cudaMemcpyHostToDevice));
+	t.pot = NULL;
+	return true;
+}
 
 void ComputeForce::decompose() {
 	//gpuErrchk( cudaProfilerStart() );
@@ -763,6 +815,18 @@ float ComputeForce::computeTabulated(bool get_energy) {
 	    computeTabulatedDihedrals<<<nb, numThreads, 0, gpuman.get_next_stream()>>>(forceInternal_d[0], pos_d[0], sys_d[0], numDihedrals*numReplicas, 
                 dihedralList_d, dihedralPotList_d, tableDihedral_d, energies_d, get_energy);
         }
+	if (vecangleList_d != NULL && tableVecangle_d != NULL)
+        {
+            //if(get_energy)
+		//computeTabulatedVecangles<<<nb, numThreads>>>(forceInternal_d[0], pos_d[0], sys_d[0], numVecangles*numReplicas, vecangleList_d, vecanglePotList_d, tableVecangle_d);
+	    if (get_energy) {
+		computeTabulatedVecangles<true><<<nb, numThreads, 0, gpuman.get_next_stream()>>>(forceInternal_d[0], pos_d[0], sys_d[0], numVecangles*numReplicas,
+                vecangleList_d, vecanglePotList_d, tableVecangle_d, energies_d);
+	    } else {
+		computeTabulatedVecangles<false><<<nb, numThreads, 0, gpuman.get_next_stream()>>>(forceInternal_d[0], pos_d[0], sys_d[0], numVecangles*numReplicas,
+                vecangleList_d, vecanglePotList_d, tableVecangle_d);
+	}
+        }
 
 	// TODO: Sum energy
 	if (restraintIds_d != NULL )
@@ -812,6 +876,16 @@ float ComputeForce::computeTabulatedFull(bool get_energy) {
 																							  tableDihedral_d, numDihedrals,
 																								num+num_rb_attached_particles, sys_d[0], energies_d,
 																								get_energy);
+	gpuErrchk(cudaDeviceSynchronize());
+	if (get_energy) {
+		computeTabulatedVecangles<true><<<numBlocks, numThreads, 0, gpuman.get_next_stream()>>>(forceInternal_d[0], pos_d[0], sys_d[0], numVecangles*numReplicas,
+                vecangleList_d, vecanglePotList_d, tableVecangle_d, energies_d);
+	    } else {
+		computeTabulatedVecangles<false><<<numBlocks, numThreads, 0, gpuman.get_next_stream()>>>(forceInternal_d[0], pos_d[0], sys_d[0], numVecangles*numReplicas,
+                vecangleList_d, vecanglePotList_d, tableVecangle_d);
+	}
+
+
 	// Calculate the energy based on the array created by the kernel
 	if (get_energy) {
 		gpuErrchk(cudaDeviceSynchronize());
@@ -891,7 +965,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, const BondAngle* const bondAngles, const XpotMap simple_potential_map, const std::vector<SimplePotential> simple_potentials, const ProductPotentialConf* const product_potential_confs)
+void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, Vecangle* vecangles, const Restraint* const restraints, const BondAngle* const bondAngles, const XpotMap simple_potential_map, const std::vector<SimplePotential> simple_potentials, const ProductPotentialConf* const product_potential_confs)
 {
     assert(simNum == numReplicas); // Not sure why we have both of these things
     int tot_num_with_rb = (num+num_rb_attached_particles) * simNum;
@@ -939,6 +1013,13 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
 												 		  sizeof(Dihedral) * numDihedrals,
 														 	cudaMemcpyHostToDevice));
 	}
+	if (numVecangles > 0) {
+		// Vecangles_d
+		gpuErrchk(cudaMalloc(&vecangles_d, sizeof(Vecangle) * numVecangles));
+		gpuErrchk(cudaMemcpyAsync(vecangles_d, vecangles,
+					  sizeof(Vecangle) * numVecangles,
+					  cudaMemcpyHostToDevice));
+	}
 
 	if (numRestraints > 0) {
 	    int restraintIds[numRestraints];
@@ -1079,7 +1160,7 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
 // 	}
 // }
 
-void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int4* bondAngleList, int2* restraintList) {
+void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int4 *vecangleList, int *vecanglePotList, int4* bondAngleList, int2* restraintList) {
 	size_t size;
 
 	if (numBonds > 0) {
@@ -1104,6 +1185,16 @@ void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *d
     gpuErrchk( cudaMemcpyAsync( dihedralPotList_d, dihedralPotList, size, cudaMemcpyHostToDevice) );
 	}
 
+	if (numVecangles > 0) {
+    size = numVecangles * numReplicas * sizeof(int4);
+    gpuErrchk( cudaMalloc( &vecangleList_d, size ) );
+    gpuErrchk( cudaMemcpyAsync( vecangleList_d, vecangleList, size, cudaMemcpyHostToDevice) );
+
+    size = numVecangles * numReplicas * sizeof(int);
+    gpuErrchk( cudaMalloc( &vecanglePotList_d, size ) );
+    gpuErrchk( cudaMemcpyAsync( vecanglePotList_d, vecanglePotList, size, cudaMemcpyHostToDevice) );
+	}
+
 	if (numBondAngles > 0) {
 	    size = 2*numBondAngles * numReplicas * sizeof(int4);
 	    gpuErrchk( cudaMalloc( &bondAngleList_d, size ) );

src/ComputeForce.cuh

@@ -1064,6 +1064,31 @@ void computeTabulatedDihedrals(Vector3* force, const Vector3* __restrict__ pos,
 	// }
     }
 }
+template<bool get_energy>
+__global__
+void computeTabulatedVecangles(Vector3* force, const Vector3* __restrict__ pos,
+			       const BaseGrid* __restrict__ sys,
+			       int numVecangles, const int4* const __restrict__ vecangleList_d,
+			       const int* __restrict__ vecanglePotList_d, TabulatedVecanglePotential** tableVecangle, float* energy=nullptr) {
+
+	// int currVecangle = blockIdx.x * blockDim.x + threadIdx.x; // first particle ID
+
+    // Loop over ALL vecangles in ALL replicas
+	for (int i = threadIdx.x+blockIdx.x*blockDim.x; i < numVecangles; i+=blockDim.x*gridDim.x) {
+		const int4& ids = vecangleList_d[i];
+		const int& id = vecanglePotList_d[i];
+		TabulatedVecanglePotential* __restrict__ p = tableVecangle[ id ];
+		float tmp = p->compute_value( pos, sys, ids);
+		float2 ef = p->compute_energy_and_deriv(tmp);
+		p->apply_force(pos,sys,force, ids, ef.y);
+		if (get_energy) {
+		    atomicAdd( &energy[ids.x], ef.x*0.25 );
+		    atomicAdd( &energy[ids.y], ef.x*0.25 );
+		    atomicAdd( &energy[ids.z], ef.x*0.25 );
+		    atomicAdd( &energy[ids.w], ef.x*0.25 );
+		}
+	}
+}
 
 __global__
 void computeHarmonicRestraints(Vector3* force, const Vector3* __restrict__ pos,

src/ComputeForce.h

@@ -69,6 +69,7 @@ public:
 	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);
+	bool addVecanglePotential(String fileName, int ind, Vecangle* vecangles);
 
 	void decompose();
 	
@@ -95,7 +96,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, const Restraint* const restraints, const BondAngle* const bondAngles,
+	void copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, Vecangle* vecangles, const Restraint* const restraints, const BondAngle* const bondAngles,
 			const XpotMap simple_potential_map,
 			const std::vector<SimplePotential> simple_potentials,
 			const ProductPotentialConf* const product_potential_confs);
@@ -103,7 +104,7 @@ public:
         void copyToCUDA(Vector3* forceInternal, Vector3* pos, Vector3* mom, float* random);
 	
 	// void createBondList(int3 *bondList);
-    void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int4 *bondAngleList, int2 *restraintList);
+    void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int4 *vecangleList, int *vecanglePotList, int4 *bondAngleList, int2 *restraintList);
 	    
 	//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()
@@ -160,6 +161,11 @@ public:
 		return dihedrals_d;
 	}
 
+	Vecangle* getVecangles_d()
+	{
+		return vecangles_d;
+	}
+
 	int3* getBondList_d()
 	{
 		return bondList_d;
@@ -207,6 +213,8 @@ private:
 	int numTabAngleFiles;
 	int numDihedrals;
 	int numTabDihedralFiles;
+	int numVecangles;
+	int numTabVecangleFiles;
 
 	int numGroupSites;
 	int* comSiteParticles;
@@ -217,6 +225,7 @@ private:
 	TabulatedPotential **tableBond;
 	TabulatedAnglePotential **tableAngle;
 	TabulatedDihedralPotential **tableDihedral;
+	TabulatedVecanglePotential **tableVecangle;
 	const BaseGrid* sys;
 	float switchStart, switchLen, electricConst, cutoff2;
 	CellDecomposition decomp;
@@ -237,6 +246,7 @@ private:
 	TabulatedPotential **tableBond_d, **tableBond_addr;
 	TabulatedAnglePotential **tableAngle_d, **tableAngle_addr;
 	TabulatedDihedralPotential **tableDihedral_d, **tableDihedral_addr;
+	TabulatedVecanglePotential **tableVecangle_d, **tableVecangle_addr;
 
 	// Pairlists
 	float pairlistdist2;
@@ -273,6 +283,7 @@ private:
 
 	Angle* angles_d;
 	Dihedral* dihedrals_d;
+	Vecangle* vecangles_d;
 
 	int numBondAngles;
 	BondAngle* bondAngles_d;
@@ -290,6 +301,8 @@ private:
 	int4* angleList_d;
 	int4* dihedralList_d;
 	int* dihedralPotList_d;
+	int4* vecangleList_d;
+	int* vecanglePotList_d;
         int2* restraintList_d;
     
 	int numRestraints;

src/Configuration.cpp

@@ -285,6 +285,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
 	if (readBondsFromFile) readBonds();
 	if (readAnglesFromFile) readAngles();
 	if (readDihedralsFromFile) readDihedrals();
+	if (readVecanglesFromFile) readVecangles();
 	if (readRestraintsFromFile) readRestraints();
 	if (readBondAnglesFromFile) readBondAngles();
 	if (readProductPotentialsFromFile) readProductPotentials();
@@ -382,12 +383,14 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
 	    }
 		if (partForceXGridFile[i].length() != 0) {
 			part[i].forceXGrid = new BaseGrid(partForceXGridFile[i].val());
+			if (partForceGridScale[i] != nullptr) part[i].forceXGrid->scale( partForceGridScale[i][0] );
 			printf("Loaded `%s'.\n", partForceXGridFile[i].val());
 			printf("Grid size %s.\n", part[i].forceXGrid->getExtent().toString().val());
 		}
 
 		if (partForceYGridFile[i].length() != 0) {
 			part[i].forceYGrid = new BaseGrid(partForceYGridFile[i].val());
+			if (partForceGridScale[i] != nullptr) part[i].forceYGrid->scale( partForceGridScale[i][1] );
 			printf("Loaded `%s'.\n", partForceYGridFile[i].val());
 			printf("Grid size %s.\n", part[i].forceYGrid->getExtent().toString().val());
 		}
@@ -396,6 +399,10 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) :
 			part[i].forceZGrid = new BaseGrid(partForceZGridFile[i].val());
 			printf("Loaded `%s'.\n", partForceZGridFile[i].val());
 			printf("Grid size %s.\n", part[i].forceZGrid->getExtent().toString().val());
+			if (partForceGridScale[i] != nullptr) {
+			    printf("Scaling forceGridZ `%s' by %f.\n", partForceZGridFile[i].val(), partForceGridScale[i][2] );
+			    part[i].forceZGrid->scale( partForceGridScale[i][2] );
+			}
 		}
 
 		if (partDiffusionGridFile[i].length() != 0) {
@@ -654,6 +661,7 @@ Configuration::~Configuration() {
 	if (excludeMap != NULL) delete[] excludeMap;
 	if (angles != NULL) delete[] angles;
 	if (dihedrals != NULL) delete[] dihedrals;
+	if (vecangles != NULL) delete[] vecangles;
 	if (bondAngles != NULL) delete[] bondAngles;
 	if (productPotentials != NULL) delete[] productPotentials;
 
@@ -669,6 +677,7 @@ Configuration::~Configuration() {
 	delete[] angleTableFile;
 
 	delete[] dihedralTableFile;
+	delete[] vecangleTableFile;
 
 	//if (type_d != NULL) {
 		//gpuErrchk(cudaFree(type_d));
@@ -755,6 +764,17 @@ void Configuration::copyToCUDA() {
 		    b->diffusionGrid = NULL;
 		}
 		
+		// Copy the diffusion grid
+		if (part[i].forceXGrid != nullptr) {
+		    b->forceXGrid = part[i].forceXGrid->copy_to_cuda();
+		}
+		if (part[i].forceYGrid != nullptr) {
+		    b->forceYGrid = part[i].forceYGrid->copy_to_cuda();
+		}
+		if (part[i].forceZGrid != nullptr) {
+		    b->forceZGrid = part[i].forceZGrid->copy_to_cuda();
+		}
+
 		//b->pmf = pmf;
 		gpuErrchk(cudaMalloc(&part_addr[i], sizeof(BrownianParticleType)));
 		gpuErrchk(cudaMemcpyAsync(part_addr[i], b, sizeof(BrownianParticleType),
@@ -887,6 +907,11 @@ void Configuration::setDefaults() {
 	dihedrals = NULL;
 	numTabDihedralFiles = 0;
 
+	readVecanglesFromFile = false;
+	numVecangles = 0;
+	vecangles = NULL;
+	numTabVecangleFiles = 0;
+
 	readBondAnglesFromFile = false;
 	numBondAngles = 0;
 	bondAngles = NULL;
@@ -936,6 +961,8 @@ int Configuration::readParameters(const char * config_file) {
 	partForceXGridFile = new String[numParts];
 	partForceYGridFile = new String[numParts];
 	partForceZGridFile = new String[numParts];
+        partForceGridScale  = new float*[numParts];
+
 	partDiffusionGridFile = new String[numParts];
 	partReservoirFile = new String[numParts];
 	partRigidBodyGrid.resize(numParts);
@@ -958,6 +985,7 @@ int Configuration::readParameters(const char * config_file) {
         {
             partGridFile[i] = NULL;
             partGridFileScale[i] = NULL;
+            partForceGridScale[i] = nullptr;
             //part[i].numPartGridFiles = -1;
         }
         //for(int i = 0; i < numParts; ++i)
@@ -972,11 +1000,15 @@ int Configuration::readParameters(const char * config_file) {
 	int dtfcap = 10;
 	dihedralTableFile = new String[dtfcap];
 
+	int vatfcap = 10;
+	vecangleTableFile = new String[vatfcap];
+
 	int currPart = -1;
 	int currTab = -1;
 	int currBond = -1;
 	int currAngle = -1;
 	int currDihedral = -1;
+	int currVecangle = -1;
 	int currRB = -1;
 
 	int partClassPart =  0;
@@ -1085,6 +1117,14 @@ int Configuration::readParameters(const char * config_file) {
 		} else if (param == String("forceZGridFile")) {
 		    if (currPart < 0) exit(1);
 		    partForceZGridFile[currPart] = value;
+		} else if (param == String("forceGridScale")) {
+		    if (currPart < 0) exit(1);
+		    int tmp;
+		    stringToArray<float>(&value, tmp, &partForceGridScale[currPart]);
+		    if (tmp != 3) {
+			printf("ERROR: Expected three floating point scale values for x,y,z, but got `%s'.\n", param.val());
+			exit(1);
+		    }
 		} else if (param == String("diffusionGridFile")) {
 		    if (currPart < 0) exit(1);
 		    partDiffusionGridFile[currPart] = value;
@@ -1201,6 +1241,24 @@ int Configuration::readParameters(const char * config_file) {
 			}
 			if (readDihedralFile(value, ++currDihedral))
 				numTabDihedralFiles++;
+		} else if (param == String("inputVecangles")) {
+			if (readVecanglesFromFile) {
+				printf("WARNING: More than one vecangle file specified. Ignoring new vecangle file.\n");
+			} else {
+				vecangleFile = value;
+				readVecanglesFromFile = true;
+			}
+		} else if (param == String("tabulatedVecangleFile")) {
+			if (numTabVecangleFiles >= dtfcap) {
+				String * temp = vecangleTableFile;
+				dtfcap *= 2;
+				vecangleTableFile = new String[dtfcap];
+				for (int j = 0; j < numTabVecangleFiles; j++)
+					vecangleTableFile[j] = temp[j];
+				delete[] temp;
+			}
+			if (readVecangleFile(value, ++currVecangle))
+				numTabVecangleFiles++;
 		} else if (param == String("inputRestraints")) {
 			if (readRestraintsFromFile) {
 				printf("WARNING: More than one restraint file specified. Ignoring new restraint file.\n");
@@ -1975,6 +2033,70 @@ void Configuration::readDihedrals() {
 	// 	dihedrals[i].print();
 }
 
+void Configuration::readVecangles() {
+	FILE* inp = fopen(vecangleFile.val(), "r");
+	char line[256];
+	int capacity = 256;
+	numVecangles = 0;
+	vecangles = new Vecangle[capacity];
+
+	// If the vecangle file cannot be found, exit the program
+	if (inp == NULL) {
+		printf("WARNING: Could not open `%s'.\n", vecangleFile.val());
+		printf("This simulation will not use vecangles.\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 VECANGLE inputs have 6 tokens
+		// VECANGLE | INDEX1 | INDEX2 | INDEX3 | INDEX4 | FILENAME
+		// Any angle input line without exactly 6 tokens should be discarded
+		if (numTokens != 6) {
+			printf("WARNING: Invalid vecangle 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());
+		int ind4 = atoi(tokenList[4].val());
+		String file_name = tokenList[5];
+		//printf("file_name %s\n", file_name.val());
+		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 (numVecangles >= capacity) {
+			Vecangle* temp = vecangles;
+			capacity *= 2;
+			vecangles = new Vecangle[capacity];
+			for (int i = 0; i < numVecangles; ++i)
+				vecangles[i] = temp[i];
+			delete[] temp;
+		}
+
+		Vecangle d(ind1, ind2, ind3, ind4, file_name);
+		vecangles[numVecangles++] = d;
+		delete[] tokenList;
+	}
+	std::sort(vecangles, vecangles + numVecangles, compare());
+
+	// for(int i = 0; i < numVecangles; i++)
+	// 	vecangles[i].print();
+}
+
 void Configuration::readBondAngles() {
 	FILE* inp = fopen(bondAngleFile.val(), "r");
 	char line[256];
@@ -2093,7 +2215,7 @@ void Configuration::readProductPotentials() {
 			// Try to match a type
 			String n = tokenList[i];
 			n.lower();
-			if (n == "bond") { type = ΒΟΝD; type_specified = true; }
+			if (n == "bond") { type = BOND; type_specified = true; }
 			else if (n == "angle")  { type = ANGLE; type_specified = true; }
 			else if (n == "dihedral")  { type = DIHEDRAL; type_specified = true; }
 			else if (n == "vecangle") { type = VECANGLE; type_specified = true; }
@@ -2282,6 +2404,26 @@ bool Configuration::readDihedralFile(const String& value, int currDihedral) {
 
 	return true;
 }
+bool Configuration::readVecangleFile(const String& value, int currVecangle) {
+	int numTokens = value.tokenCount();
+	if (numTokens != 1) {
+		printf("ERROR: Invalid tabulatedVecangleFile: %s, numTokens = %d\n", value.val(), numTokens);
+		return false;
+	}
+
+	String* tokenList = new String[numTokens];
+	value.tokenize(tokenList);
+	if (tokenList == NULL) {
+		printf("ERROR: Invalid tabulatedVecangleFile: %s; tokenList is NULL\n", value.val());
+		return false;
+	}
+
+	vecangleTableFile[currVecangle] = tokenList[0];
+
+	// printf("Tabulated Vecangle Potential: %s\n", vecangleTableFile[currVecangle].val() );
+
+	return true;
+}
 //Load the restart coordiantes only
 void Configuration::loadRestart(const char* file_name) {
 	char line[STRLEN];

src/Configuration.h

@@ -40,6 +40,7 @@
 // Forward declerations
 class Angle;
 class Dihedral;
+using Vecangle = Dihedral;
 struct Restraint;
 
 class Configuration {
@@ -49,6 +50,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 Vecangle& lhs, const Vecangle& rhs);
 		bool operator()(const BondAngle& lhs, const BondAngle& rhs);
 		bool operator()(const ProductPotentialConf& lhs, const ProductPotentialConf& rhs);
 	};
@@ -66,6 +68,7 @@ class Configuration {
 	void addExclusion(int ind1, int ind2);
 	void buildExcludeMap();
 	void readDihedrals();
+	void readVecangles();
 	void readRestraints();
 	void readBondAngles();
 
@@ -74,6 +77,7 @@ class Configuration {
 	bool readBondFile(const String& value, int currBond);
 	bool readAngleFile(const String& value, int currAngle);
 	bool readDihedralFile(const String& value, int currDihedral);
+	bool readVecangleFile(const String& value, int currVecangle);
 
 	bool readBondAngleFile(const String& value, const String& bondfile1, const String& bondfile2, int currBondAngle);
 
@@ -207,6 +211,7 @@ public:
 	int numExcludes;
 	int numAngles;
 	int numDihedrals;
+	int numVecangles;
 	int numBondAngles;
 	int numRestraints;
 	int* numPartsOfType;
@@ -215,6 +220,7 @@ public:
 	String excludeFile;
 	String angleFile;
 	String dihedralFile;
+	String vecangleFile;
 	String restraintFile;
 	String bondAngleFile;
 	bool readPartsFromFile;
@@ -223,6 +229,7 @@ public:
 	bool readExcludesFromFile;
 	bool readAnglesFromFile;
 	bool readDihedralsFromFile;
+	bool readVecanglesFromFile;
 	bool readBondAnglesFromFile;
 	bool readRestraintsFromFile;
 	//String* partGridFile;
@@ -237,6 +244,7 @@ public:
 	String* partForceXGridFile;
 	String* partForceYGridFile;
 	String* partForceZGridFile;
+	float **partForceGridScale;
 	String* partTableFile;
 	String* partReservoirFile;
 	int* partTableIndex0;
@@ -263,6 +271,10 @@ public:
 	String* dihedralTableFile;
 	int numTabDihedralFiles;
 
+        Vecangle* vecangles;
+	String* vecangleTableFile;
+	int numTabVecangleFiles;
+
 	BondAngle* bondAngles;
 
 	Restraint* restraints;

src/Dihedral.h

@@ -28,4 +28,6 @@ public:
 	void print();
 };
 
+using Vecangle = Dihedral;
+
 #endif

src/GrandBrownTown.cu

@@ -188,6 +188,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
 	bondTableFile = c.bondTableFile;
 	angleTableFile = c.angleTableFile;
 	dihedralTableFile = c.dihedralTableFile;
+	vecangleTableFile = c.vecangleTableFile;
 
 	// Other parameters.
 	switchStart = c.switchStart;
@@ -201,6 +202,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
 	numExcludes = c.numExcludes;
 	numAngles = c.numAngles;
 	numDihedrals = c.numDihedrals;
+	numVecangles = c.numVecangles;
 	partTableIndex0 = c.partTableIndex0;
 	partTableIndex1 = c.partTableIndex1;
 
@@ -220,7 +222,9 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
 	numTabAngleFiles = c.numTabAngleFiles;
 
 	dihedrals = c.dihedrals;
+	vecangles = c.vecangles;
 	numTabDihedralFiles = c.numTabDihedralFiles;
+	numTabVecangleFiles = c.numTabVecangleFiles;
 
 	bondAngles = c.bondAngles;
 
@@ -262,7 +266,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, c.bondAngles, c.simple_potential_ids, c.simple_potentials, c.productPotentials );
+	internal -> copyToCUDA(c.simNum, c.type, c.bonds, c.bondMap, c.excludes, c.excludeMap, c.angles, c.dihedrals, c.vecangles, c.restraints, c.bondAngles, c.simple_potential_ids, c.simple_potentials, c.productPotentials );
 	if (numGroupSites > 0) init_cuda_group_sites();
 
 	// TODO: check for duplicate potentials 
@@ -313,6 +317,13 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
 				internal->addDihedralPotential(dihedralTableFile[p].val(), p, dihedrals);
 	}
 
+	if (c.readVecanglesFromFile) {
+		printf("Loading %d tabulated vecangle potentials...\n", numTabVecangleFiles);
+		for (int p = 0; p < numTabVecangleFiles; p++)
+			if (vecangleTableFile[p].length() > 0)
+				internal->addVecanglePotential(vecangleTableFile[p].val(), p, vecangles);
+	}
+
 	auto _get_index = [this](int idx, int replica) {
 	    // Convenient lambda function to deal with increasingly complicated indexing
 	    auto num = this->num;
@@ -377,6 +388,20 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
 	    }
 	}
 	}
+	if (numVecangles > 0) {
+	vecangleList = new int4[ (numVecangles) * numReplicas ];
+	vecanglePotList = new  int[ (numVecangles) * numReplicas ];
+	for(int k = 0 ; k < numReplicas; k++) {
+	    for(int i = 0; i < numVecangles; ++i) {
+			if (vecangles[i].tabFileIndex == -1) {
+				fprintf(stderr,"Error: vecanglefile '%s' was not read with tabulatedVecangleFile command.\n", vecangles[i].fileName.val());
+				exit(1);
+			}
+			vecangleList[i+k*numVecangles] = make_int4( _get_index(vecangles[i].ind1,k), _get_index(vecangles[i].ind2,k), _get_index(vecangles[i].ind3,k), _get_index(vecangles[i].ind4,k) );
+		vecanglePotList[i+k*numVecangles] = vecangles[i].tabFileIndex;
+	    }
+	}
+	}
 
 	if (numBondAngles > 0) {
 	bondAngleList = new int4[ (numBondAngles*2) * numReplicas ];
@@ -422,7 +447,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
 	    }
 	}
 	
-	internal->copyBondedListsToGPU(bondList,angleList,dihedralList,dihedralPotList,bondAngleList,restraintList);
+	internal->copyBondedListsToGPU(bondList,angleList,dihedralList,dihedralPotList,vecangleList,vecanglePotList,bondAngleList,restraintList);
 	if (c.numRestraints > 0) delete[] restraintList;
 	
 	forceInternal = new Vector3[(num+num_rb_attached_particles+numGroupSites)*numReplicas];
@@ -482,6 +507,10 @@ GrandBrownTown::~GrandBrownTown() {
 		delete[] dihedralList;
 		delete[] dihedralPotList;
 	}
+	if (numVecangles > 0) {
+		delete[] vecangleList;
+		delete[] vecanglePotList;
+	}
         if(randoGen->states != NULL)
             {
                 gpuErrchk(cudaFree(randoGen->states));

src/GrandBrownTown.h

@@ -37,6 +37,7 @@
 #include "Angle.h"
 #include "Configuration.h"
 #include "Dihedral.h"
+#include "TabulatedPotential.h"
 /* #include "RigidBody.h" */
 /* #include "RigidBodyType.h" */
 /* #include "RigidBodyGrid.h" */
@@ -228,6 +229,7 @@ private:
 	int numExcludes;
 	int numAngles;
 	int numDihedrals;
+	int numVecangles;
 
     int num_rb_attached_particles;
 
@@ -239,11 +241,13 @@ private:
 	String excludeFile;
 	String angleFile;
 	String dihedralFile;
+	String vecangleFile;
 	bool readPartsFromFile;
 	bool readBondsFromFile;
 	bool readExcludesFromFile;
 	bool readAnglesFromFile;
 	bool readDihedralsFromFile;
+	bool readVecanglesFromFile;
 	String* partGridFile;
 	String* partDiffusionGridFile;
 	String* partForceXGridFile;
@@ -275,6 +279,12 @@ private:
 	int4 *dihedralList;
 	int  *dihedralPotList;
 
+	Vecangle* vecangles;
+	String* vecangleTableFile;
+	int numTabVecangleFiles;
+	int4 *vecangleList;
+	int  *vecanglePotList;
+
 	int numBondAngles;
 	BondAngle* bondAngles;
 	int4 *bondAngleList;

src/TabulatedPotential.h

@@ -254,6 +254,20 @@ public:
 	    val = compute_vecangle(pos, sys, particles[0], particles[1], particles[2], particles[3]);
 	return val;
     }
+    HOST DEVICE inline float compute_value(const Vector3* __restrict__ pos,
+					   const BaseGrid* __restrict__ sys,
+					   const int4& particles) const {
+	float val;
+	if (type == BOND)
+	    val = compute_bond(pos, sys, particles.x, particles.y);
+	else if (type == ANGLE)
+	    val = compute_angle(pos, sys, particles.x, particles.y, particles.z);
+	else if (type == DIHEDRAL)
+	    val = compute_dihedral(pos, sys, particles.x, particles.y, particles.z, particles.w);
+	else if (type == VECANGLE)
+	    val = compute_vecangle(pos, sys, particles.x, particles.y, particles.z, particles.w);
+	return val;
+    }
 
     HOST DEVICE inline float2 compute_energy_and_deriv(float value) {
 	float2 ret;
@@ -285,8 +299,15 @@ public:
 					      int i, int j, int k, int l) const {
 	const Vector3 ab = sys->wrapDiff(pos[j] - pos[i]);
 	const Vector3 bc = sys->wrapDiff(pos[l] - pos[k]);
-	const Vector3 ac = bc+ab;
-	return compute_angle( ab.length2(), bc.length2(), ac.length2() );
+	const Vector3 ac = bc-ab;
+	float tmp = compute_angle( ab.length2(), bc.length2(), ac.length2() );
+	// printf("i,j,r_ij: %d %d %f %f %f\n",
+	//        i,j, ab.x,ab.y,ab.z);
+	// printf("k,l,r_kl: %d %d %f %f %f\n",
+	//        k,l, bc.x,bc.y,bc.z);
+	// printf("i,j,k,l, theta: %d %d %d %d %f\n",
+	//        i,j,k,l, tmp);
+	return tmp;
     }
 
     HOST DEVICE inline float compute_angle(float distab2, float distbc2, float distac2) const {
@@ -358,6 +379,22 @@ public:
 	    apply_vecangle_force(pos, sys, forces, particles[0], particles[1],
 				 particles[2], particles[3], energy_deriv);
     }
+    DEVICE inline void apply_force(const Vector3* __restrict__ pos,
+				   const BaseGrid* __restrict__ sys,
+				   Vector3* __restrict__ forces,
+				   const int4& particles, float energy_deriv) const {
+	if (type == BOND)
+	    apply_bond_force(pos, sys, forces, particles.x, particles.y, energy_deriv);
+	else if (type == ANGLE)
+	    apply_angle_force(pos, sys, forces, particles.x, particles.y,
+			     particles.z, energy_deriv);
+	else if (type == DIHEDRAL)
+	    apply_dihedral_force(pos, sys, forces, particles.x, particles.y,
+				 particles.z, particles.w, energy_deriv);
+	else if (type == VECANGLE)
+	    apply_vecangle_force(pos, sys, forces, particles.x, particles.y,
+				 particles.z, particles.w, energy_deriv);
+    }
 
     __device__ inline void apply_bond_force(const Vector3* __restrict__ pos,
 					const BaseGrid* __restrict__ sys,
@@ -397,9 +434,10 @@ public:
 
 	float sin = sqrtf(1.0f - cos*cos);
 	energy_deriv /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
-	if (abs(sin) < 1e-3) {
-	    printf("BAD ANGLE: sin, cos, energy_deriv, distab, distbc, distac2: (%f %f %f %f %f)\n",
-		   sin,cos,energy_deriv,distab,distbc);	}
+	if (abs(sin) <= 1e-3) {
+	    printf("BAD ANGLE: sin, cos, energy_deriv, distab, distbc, distac2: (%f %f %f %f %f %f)\n",
+		   sin,cos,energy_deriv,ab.length(),bc.length(),sqrt(distac2));
+	}
 
 	// Calculate the forces
 	TwoVector3 force;
@@ -504,21 +542,25 @@ public:
 
 #ifdef __CUDA_ARCH__
 
-	const Vector3 ab = -sys->wrapDiff(pos[j] - pos[i]);
-	const Vector3 bc = -sys->wrapDiff(pos[k] - pos[j]);
+	const Vector3 ab = sys->wrapDiff(pos[j] - pos[i]);
+	const Vector3 cd = -sys->wrapDiff(pos[l] - pos[k]);
 	// const Vector3 ac = sys->wrapDiff(pos[k] - pos[i]);
 
-	TwoVector3 f = get_angle_force(ab,bc, energy_deriv);
+	TwoVector3 f = get_angle_force(ab,cd, energy_deriv);
 
+	// printf("i,j,k,l, df, f1, f2: %d %d %d %d %f (%f %f %f) (%f %f %f)\n",
+	//        i,j,k,l, energy_deriv, f.v1.x, f.v1.y, f.v1.z, f.v2.x, f.v2.y, f.v2.z);
 	atomicAdd( &force[i], f.v1 );
 	atomicAdd( &force[j], -f.v1 );
-	atomicAdd( &force[k], -f.v2 );
-	atomicAdd( &force[l], f.v2 );
+	atomicAdd( &force[k], f.v2 );
+	atomicAdd( &force[l], -f.v2 );
 #endif
     }
 
 };
 
+using TabulatedVecanglePotential = SimplePotential;
+
 #ifndef delgpuErrchk
 #undef  delgpuErrchk
 #undef  gpuErrchk(code)