Configuration.cpp

		else if (param == String("temperature"))
			temperature =  (float) strtod(value.val(),NULL);
		else if (param == String("temperatureGrid"))
			temperatureGridFile = value;
		else if (param == String("numberFluct"))
			numberFluct = atoi(value.val());
		else if (param == String("numberFluctPeriod"))
			numberFluctPeriod = atoi(value.val());
		else if (param == String("interparticleForce"))
			interparticleForce = atoi(value.val());
		else if (param == String("fullLongRange") || param == String("fullElect") )
			fullLongRange = atoi(value.val());
		else if (param == String("coulombConst"))
			coulombConst = (float) strtod(value.val(), NULL);
		else if (param == String("electricField"))
			electricField = (float) strtod(value.val(), NULL);
		else if (param == String("cutoff"))
			cutoff = (float) strtod(value.val(), NULL);
		else if (param == String("switchLen"))
			switchLen = (float) strtod(value.val(), NULL);
		else if (param == String("pairlistDistance"))
			pairlistDistance = (float) strtod(value.val(), NULL);
		else if (param == String("scaleIMDForce"))
			imdForceScale = (float) strtod(value.val(), NULL);		
		else if (param == String("outputPeriod"))
			outputPeriod = atoi(value.val());
		else if (param == String("outputEnergyPeriod"))
			outputEnergyPeriod = atoi(value.val());
		else if (param == String("outputFormat"))
			outputFormat = TrajectoryWriter::getFormatCode(value);
		else if (param == String("currentSegmentZ"))
			currentSegmentZ = (float) strtod(value.val(), NULL);
		else if (param == String("numCap"))
			numCap = atoi(value.val());
		else if (param == String("decompPeriod"))
			decompPeriod = atoi(value.val());

                //Han-Yi Chou
                else if (param == String("ParticleDynamicType"))
                    ParticleDynamicType = value;
                else if (param == String("RigidBodyDynamicType"))
                    RigidBodyDynamicType = value;
                else if (param == String("ParticleLangevinIntegrator"))
                    ParticleLangevinIntegrator = value;
                else if (param == String("ParticleInterpolationType"))
                    ParticleInterpolationType = atoi(value.val());
                else if (param == String("RigidBodyInterpolationType"))
                    RigidBodyInterpolationType = atoi(value.val());
		// PARTICLES
		else if (param == String("particle")) {
		    part[++currPart] = BrownianParticleType(value);
		    currPartClass = partClassPart;
		}
                else if (param == String("mu")) { // for Nose-Hoover Langevin
		    if (currPart < 0) exit(1);
		    part[currPart].mu = (float) strtod(value.val(), NULL);
		} else if (param == String("forceXGridFile")) {
		    if (currPart < 0) exit(1);
		    partForceXGridFile[currPart] = value;
		} else if (param == String("forceYGridFile")) {
		    if (currPart < 0) exit(1);
		    partForceYGridFile[currPart] = value;
		} else if (param == String("forceZGridFile")) {
		    if (currPart < 0) exit(1);
		    partForceZGridFile[currPart] = value;
		} else if (param == String("diffusionGridFile")) {
		    if (currPart < 0) exit(1);
		    partDiffusionGridFile[currPart] = value;
		} else if (param == String("diffusion")) {
		    if (currPart < 0) exit(1);
		    part[currPart].diffusion = (float) strtod(value.val(), NULL);
		} else if (param == String("charge")) {
		    if (currPart < 0) exit(1);
		    part[currPart].charge = (float) strtod(value.val(), NULL);
		} else if (param == String("radius")) {
		    if (currPart < 0) exit(1);
		    part[currPart].radius = (float) strtod(value.val(), NULL);
		} else if (param == String("eps")) {
		    if (currPart < 0) exit(1);
		    part[currPart].eps = (float) strtod(value.val(), NULL);
		} else if (param == String("reservoirFile")) {
		    if (currPart < 0) exit(1);
		    partReservoirFile[currPart] = value;
		}
		else if (param == String("tabulatedPotential"))
			tabulatedPotential = atoi(value.val());
		else if (param == String("tabulatedFile"))
			readTableFile(value, ++currTab);
		else if (param == String("tabulatedBondFile")) {
			if (numTabBondFiles >= btfcap) {
				String* temp = bondTableFile;
				btfcap *= 2;	
				bondTableFile = new String[btfcap];
				for (int j = 0; j < numTabBondFiles; j++)
					bondTableFile[j] = temp[j];
				delete[] temp;
			}
			if (readBondFile(value, ++currBond))
				numTabBondFiles++;
		} else if (param == String("inputParticles")) {
			if (readPartsFromFile) {
				printf("WARNING: More than one particle file specified. Ignoring new file.\n");
			} else {
				partFile = value;
				readPartsFromFile = true;
				loadedCoordinates = true;
			}
		} else if (param == String("inputGroups")) {
			if (readGroupSitesFromFile) {
				printf("WARNING: More than one group file specified. Ignoring new file.\n");
			} else {
				groupSiteFile = value;
				readGroupSitesFromFile = true;
			}
		} else if (param == String("inputBonds")) {
			if (readBondsFromFile) {
				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. Ignoring new exclude file.\n");
			} else {
			    printf("inputExclude %s\n", value.val());
				excludeFile = value;				
				readExcludesFromFile = true;
			}
		} else if (param == String("exclude") or param == String("exclusion")) {
			excludeRule = value; 
		} else if (param == String("inputAngles")) {
			if (readAnglesFromFile) {
				printf("WARNING: More than one angle file specified. Ignoring new angle file.\n");
			} else {
				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("inputCrossPotentials")) {
			if (readBondAnglesFromFile) {
				printf("WARNING: More than one cross potential file specified. Ignoring new file.\n");
			} else {
			        crossPotentialFile = value;
				readCrossPotentialsFromFile = true;
			}
		} else if (param == String("tabulatedAngleFile")) {
			if (numTabAngleFiles >= atfcap) {
				String* temp = angleTableFile;
				atfcap *= 2;	
				angleTableFile = new String[atfcap];
				for (int j = 0; j < numTabAngleFiles; j++)
					angleTableFile[j] = temp[j];
				delete[] temp;
			}
			if (readAngleFile(value, ++currAngle))
				numTabAngleFiles++;
		} else if (param == String("inputDihedrals")) {
			if (readDihedralsFromFile) {
				printf("WARNING: More than one dihedral file specified. Ignoring new dihedral file.\n");
			} else {
				dihedralFile = value;
				readDihedralsFromFile = true;
			}
		} else if (param == String("tabulatedDihedralFile")) {
			if (numTabDihedralFiles >= dtfcap) {
				String * temp = dihedralTableFile;
				dtfcap *= 2;
				dihedralTableFile = new String[dtfcap];
				for (int j = 0; j < numTabDihedralFiles; j++)
					dihedralTableFile[j] = temp[j];
				delete[] temp;
			}
			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")) {
		    if (currPart < 0) exit(1);
			//partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
			  stringToArray<float>(&value, part[currPart].numPartGridFiles, 
                                                      &partGridFileScale[currPart]);
		} else if (param == String("gridFileBoundaryConditions")) {
		    if (currPart < 0) exit(1);
		    register size_t num = value.tokenCount();
		    if (num > 0) {
			String *tokens  = new String[num];
			BoundaryCondition *data = new BoundaryCondition[num];
			value.tokenize(tokens);
			for(size_t i = 0; i < num; ++i) {
			    tokens[i].lower();
			    if (tokens[i] == "dirichlet")
				data[i] = dirichlet;
			    else if (tokens[i] == "neumann")
				data[i] = neumann;
			    else if (tokens[i] == "periodic")
				data[i] = periodic;
			    else {
				fprintf(stderr,"WARNING: Unrecognized gridFile boundary condition \"%s\". Using Dirichlet.\n", tokens[i].val() );
				data[i] = dirichlet;
			    }
			}
			delete[] tokens;
			part[currPart].set_boundary_conditions(num, data);
			delete[] data;
		    }
		} else if (param == String("rigidBodyPotential")) {
		    if (currPart < 0) exit(1);
		    partRigidBodyGrid[currPart].push_back(value);
		}
                //Han-Yi Chou initial COM velocity for total particles
                else if (param == String("COM_Velocity"))
                    COM_Velocity = stringToVector3(value);

		// RIGID BODY
		else if (param == String("rigidBody")) {
			// part[++currPart] = BrownianParticleType(value);
			rigidBody[++currRB] = RigidBodyType(value, this);
			currPartClass = partClassRB;
		}
		else if (param == String("inertia")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].inertia = stringToVector3( value );
		} else if (param == String("rotDamping")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].rotDamping = stringToVector3( value );
		} else if (param == String("attachedParticles")) {
			rigidBody[currRB].append_attached_particle_file(value);
		} else if (param == String("densityGrid")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].addDensityGrid(value);
		} else if (param == String("potentialGrid")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].addPotentialGrid(value);
		} else if (param == String("densityGridScale")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].scaleDensityGrid(value);
		} else if (param == String("potentialGridScale")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].scalePotentialGrid(value);
		} else if (param == String("pmfScale")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].scalePMF(value);
		} else if (param == String("position")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].initPos = stringToVector3( value );
		} else if (param == String("orientation")) {
		    if (currRB < 0) exit(1);
			rigidBody[currRB].initRot = stringToMatrix3( value );
                } else if (param == String("momentum")) {
                        rigidBody[currRB].initMomentum = stringToVector3(value);
                } else if (param == String("angularMomentum")) {
		    if (currRB < 0) exit(1);
                        rigidBody[currRB].initAngularMomentum = stringToVector3(value);
		}
		else if (param == String("inputRBCoordinates"))
			inputRBCoordinates = value;
		else if (param == String("restartRBCoordinates"))
		        restartRBCoordinates = value;
		
		// COMMON
		else if (param == String("num")) {
		    if (currPartClass == partClassPart) {
			if (currPart < 0) exit(1);
			part[currPart].num = atoi(value.val());
		    } else if (currPartClass == partClassRB) {
			if (currRB < 0) exit(1);
			rigidBody[currRB].num = atoi(value.val());
		    }
		}
                //set mass here Han-Yi Chou
                else if (param == String("mass"))
                {
                    if (currPartClass == partClassPart) {
			if (currPart < 0) exit(1);
                        part[currPart].mass    = (float) strtod(value.val(),NULL);
                    } else if (currPartClass == partClassRB) {
			if (currRB < 0) exit(1);
                        rigidBody[currRB].mass = (float) strtod(value.val(),NULL);
		    }
                }
                //set damping here, using anisotropic damping, i.e. data type Vector3 Han-Yi Chou
                else if (param == String("transDamping"))
                {
                    if (currPartClass == partClassPart) {
			if (currPart < 0) exit(1);
                        part[currPart].transDamping    = stringToVector3(value);
		    } else if (currPartClass == partClassRB) {
			if (currRB < 0) exit(1);
                        rigidBody[currRB].transDamping = stringToVector3(value);
		    }
                }
		else if (param == String("gridFile")) {
			if (currPartClass == partClassPart)
                        {
			    if (currPart < 0) exit(1);
                                printf("Applying grid file '%s'\n", value.val());
				stringToArray<String>(&value, part[currPart].numPartGridFiles, 
                                                             &partGridFile[currPart]);
				const int& num = part[currPart].numPartGridFiles;
				partGridFileScale[currPart] = new float[num];
                                for(int i = 0; i < num; ++i) {
                                    // printf("%s ", partGridFile[currPart]->val());
				    partGridFileScale[currPart][i] = 1.0f;
				}

				// Set default boundary conditions for grids
				BoundaryCondition *bc = part[currPart].pmf_boundary_conditions;
				if (bc == NULL) {
				    bc = new BoundaryCondition[num];
				    for(int i = 0; i < num; ++i) {
					bc[i] = dirichlet;
				    }
				    part[currPart].pmf_boundary_conditions = bc;
				}
                        }
			else if (currPartClass == partClassRB) {
			    if (currRB < 0) exit(1);
				rigidBody[currRB].addPMF(value);
			}
		}
		// UNKNOWN
		else {
			printf("ERROR: Unrecognized keyword `%s'.\n", param.val());
			exit(1);
		}
	}

	// extra configuration for RB types
	for (int i = 0; i < numRigidTypes; i++)
		rigidBody[i].setDampingCoeffs(timestep);

        //For debugging purpose Han-Yi Chou
        //Print();
	return numParams;
}
//Han-Yi Chou
void Configuration::Print()
{
    printf("The dynamic type for particle is %s \n", ParticleDynamicType.val());
    for(int i = 0; i < numParts; ++i)
    {
        printf("The type %d has mass %f \n", i,part[i].mass);
        printf("The diffusion coefficient is %f \n", part[i].diffusion);
        printf("The translational damping is %f %f %f \n", part[i].transDamping.x, part[i].transDamping.y, part[i].transDamping.z);
    }
    printf("Done with check for Langevin");
    //assert(1==2);
}

void Configuration::PrintMomentum()
{
    for(int i = 0; i < num; ++i)
    {
        printf("%f %f %f\n", momentum[i].x, momentum[i].y, momentum[i].z);
    }
    //assert(1==2);
}
Vector3 Configuration::stringToVector3(String s) {
	// tokenize and return
	int numTokens = s.tokenCount();
	if (numTokens != 3) {
		printf("ERROR: could not convert input to Vector3.\n"); // TODO improve this message
		exit(1);
	}
	String* token = new String[numTokens];
	s.tokenize(token);
	Vector3 v( (float) strtod(token[0], NULL),
						 (float) strtod(token[1], NULL),
						 (float) strtod(token[2], NULL) );
	return v;
}
Matrix3 Configuration::stringToMatrix3(String s) {
	// tokenize and return
	int numTokens = s.tokenCount();
	if (numTokens != 9) {
		printf("ERROR: could not convert input to Matrix3.\n"); // TODO improve this message
		exit(1);
	}
	String* token = new String[numTokens];
	s.tokenize(token);
	Matrix3 m( (float) strtod(token[0], NULL),
						 (float) strtod(token[1], NULL),
						 (float) strtod(token[2], NULL),
						 (float) strtod(token[3], NULL),
						 (float) strtod(token[4], NULL),
						 (float) strtod(token[5], NULL),
						 (float) strtod(token[6], NULL),
						 (float) strtod(token[7], NULL),
						 (float) strtod(token[8], NULL) );
	return m;
}

void Configuration::readAtoms() {
	// Open the file
	FILE* inp = fopen(partFile.val(), "r");
	char line[256];

	// If the particle file cannot be found, exit the program
	if (inp == NULL) {
		printf("ERROR: Could not open `%s'.\n", partFile.val());
		bool found = true;
		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);
		}
		printf("Using default coordinates file\n");
		return;
	}

	// Our particle array has a starting capacity of 256
	// We will expand this later if we need to.
	int capacity = 256;
	numPartsFromFile = 0;
	partsFromFile = new String[capacity];
	indices = new int[capacity];
	indices[0] = 0;

	// Get and process all lines of input
	while (fgets(line, 256, 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 ATOM input lines have 6 tokens: 
		// ATOM | Index | Name | X-coord | Y-coord | Z-coord
		// A line without exactly six tokens should be discarded.
                if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin")) {
		    if (numTokens != 9) {
			printf("Error: Invalid particle file line: %s\n", line);
			exit(-1);
		    }
		} else {
		    if (numTokens != 6) {
			printf("Error: Invalid particle file line: %s\n", line);
			exit(-1);
		    }
		}

		// Ensure that this particle's type was defined in the config file.
		// If not, discard this line.
		bool found;
		for (int j = 0; j < numParts; j++) {
			// If this particle type exists, add a new one to the list
			if (part[j].name == tokenList[2]) {
				found = true;
				part[j].num++;
			}
		}

		// If the particle's type does not exist according to the config file, discard it.
		if (!found) {
			printf("WARNING Unknown particle type %s found and discarded.\n", tokenList[2].val());
			continue;
		}

		// If we don't have enough room in our particle array, we need to expand it.
		if (numPartsFromFile >= capacity) {
			// Temporary pointers to the old arrays
			String* temp = partsFromFile;	
			int* temp2 = indices;

			// Double the capacity
			capacity *= 2;

			// Create pointers to new arrays which are twice the size of the old ones
			partsFromFile = new String[capacity];
			indices = new int[capacity];
		
			// Copy the old values into the new arrays
			for (int j = 0; j < numPartsFromFile; j++) {
				partsFromFile[j] = temp[j];
				indices[j] = temp2[j];
			}

			// delete the old arrays
			delete[] temp;
			delete[] temp2;
		}
		// Make sure the index of this particle is unique.
		// NOTE: The particle list is sorted by index. 
		bool uniqueID = true;		
		int key = atoi(tokenList[1].val());
		int mid = 0;

		// If the index is greater than the last index in the list, 
		// this particle belongs at the end of the list. Since the 
		// list is kept sorted, we know this is okay.
		if (numPartsFromFile == 0 || key > indices[numPartsFromFile - 1]) {
			indices[numPartsFromFile] = key;
			partsFromFile[numPartsFromFile++] = line;
		}
		// We need to do a binary search to figure out if
		// the index already exists in the list. 
		// The assumption is that input files SHOULD have their indices sorted in 
		// ascending order, so we shouldn't actually use the binary search 
		// or the sort (which is pretty time consuming) very often.
		else {
			int low = 0, high = numPartsFromFile - 1;
			
			while (low <= high) {
				mid = (int)((high - low) / 2 + low);
				int curr = indices[mid];
				if (curr < key) {
					low = mid + 1;
				} else if (curr > key) {
					high = mid - 1;
				} else {
					// For now, particles with non-unique IDs are simply not added to the array
					// Other possible approaches which are not yet implemented:
					// 1: Keep track of these particles and assign them new IDs after you have
					//    already added all of the other particles. 	
					// 2: Get rid of ALL particles with that ID, even the ones that have already 
					//    been added.
					printf("WARNING: Non-unique ID found: %s\n", line);
					uniqueID = false;
					break;
				}
			}
			if (uniqueID) {
				// Add the particle to the end of the array, then sort it. 
				indices[numPartsFromFile] = key;
				partsFromFile[numPartsFromFile++] = line;
				std::sort(indices, indices + numPartsFromFile);
				std::sort(partsFromFile, partsFromFile + numPartsFromFile, compare());		
			}
		}
	}
}
void Configuration::readGroups() {
	// Open the file
    const size_t line_char = 16384;
	FILE* inp = fopen(groupSiteFile.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", partFile.val());
		exit(1);
	}

	// Our particle array has a starting capacity of 256
	// We will expand this later if we need to.
	// int capacity = 256;
	numGroupSites = 0;

	// partsFromFile = new String[capacity];
	// indices = new int[capacity];
	// indices[0] = 0;

	// 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: 
		// GROUP | Atom_1_idx | Atom_2_idx | ...
		// A line without exactly six tokens should be discarded.
		if (numTokens < 3) {
		    printf("Error: Invalid group file line: %s\n", line);
		    exit(-1);
		}

		// Make sure the index of this particle is unique.
		// NOTE: The particle list is sorted by index. 
		std::vector<int> tmp;
		for (int i=1; i < numTokens; ++i) {
		    const int ai = atoi(tokenList[i].val());
		    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 );
		}
		groupSiteData.push_back(tmp);
		numGroupSites++;
	}
}

void Configuration::readBonds() {
	// Open the file
	FILE* inp = fopen(bondFile.val(), "r");
	char line[256];

	// If the particle file cannot be found, exit the program
	if (inp == NULL) {
		printf("WARNING: Could not open `%s'.\n", bondFile.val());
		printf("         This simulation will not use particle bonds.\n");
		return;
	}

	// Our particle array has a starting capacity of 256
	// We will expand this later if we need to.
	int capacity = 256;
	numBonds = 0;
	bonds = new Bond[capacity];

	// Get and process all lines of input
	while (fgets(line, 256, 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 BOND input lines have 4 tokens: 
		// BOND | OPERATION_FLAG | INDEX1 | INDEX2 | FILENAME 
		// A line without exactly five tokens should be discarded.
		if (numTokens != 5) {
			printf("WARNING: Invalid bond file line: %s\n", line);
			continue;
		}

		String op = tokenList[1];
		int ind1 = atoi(tokenList[2].val());
		int ind2 = atoi(tokenList[3].val());
		String file_name = tokenList[4];

		if (ind1 == ind2) {
			printf("WARNING: Invalid bond file line: %s\n", line);
			continue;
		}

		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);
		}

		
		// If we don't have enough room in our bond array, we need to expand it.
		if (numBonds+1 >= capacity) { // "numBonds+1" because we are adding two bonds to array
			// Temporary pointer to the old array
			Bond* temp = bonds;	

			// Double the capacity
			capacity *= 2;

			// Create pointer to new array which is twice the size of the old one
			bonds = new Bond[capacity];
		
			// Copy the old values into the new array
			for (int j = 0; j < numBonds; j++)
				bonds[j] = temp[j];

			// delete the old array
			delete[] temp;
		}
		// Add the bond to the bond array
		// We must add it twice: Once for (ind1, ind2) and once for (ind2, ind1)
		
		// RBTODO: add ind1/2 to exclusion list here iff op == REPLACE

		if (op == "REPLACE")
		    addExclusion(ind1, ind2);

		Bond* b = new Bond(op, ind1, ind2, file_name);
		bonds[numBonds++] = *b;
		b = new Bond(op, ind2, ind1, file_name);
		bonds[numBonds++] = *b;
		delete[] tokenList;
	}	
	// Call compareBondIndex with qsort to sort the bonds by BOTH ind1 AND ind2
	std::sort(bonds, bonds + numBonds, compare());

	/* Each particle may have a varying number of bonds
	 * bondMap is an array with one element for each particle
	 * which keeps track of where a particle's bonds are stored
	 * in the bonds array.
	 * 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+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;
	}
	int currPart = -1;
	int lastPart = -1;
	for (int i = 0; i < numBonds; i++) {
		if (bonds[i].ind1 != currPart) {
			currPart = bonds[i].ind1;
			bondMap[currPart].x = i;
			if (lastPart >= 0) bondMap[lastPart].y = i;
			lastPart = currPart;
		}
	}
	if (bondMap[lastPart].x > 0)
		bondMap[lastPart].y = numBonds;
}

void Configuration::readExcludes()
{
	// Open the file
	FILE* inp = fopen(excludeFile.val(), "r");
	char line[256];

	// If the exclusion file cannot be found, exit the program
	if (inp == NULL) {
		printf("WARNING: Could not open `%s'.\n", excludeFile.val());
		printf("This simulation will not use exclusions.\n");
		return;
	}


	// Get and process all lines of input
	while (fgets(line, 256, 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 EXCLUDE input lines have 3 tokens: 
		// BOND | INDEX1 | INDEX2
		// A line without exactly three tokens should be discarded.
		if (numTokens != 3) {
			printf("WARNING: Invalid exclude file line: %s\n", line);
			continue;
		}
		int ind1 = atoi(tokenList[1].val());
		int ind2 = atoi(tokenList[2].val());
		addExclusion(ind1, ind2);
		delete[] tokenList;
	}
}
void Configuration::addExclusion(int ind1, int ind2) {
    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;
    }
		
    // If we don't have enough room in our bond array, we need to expand it.
    if (numExcludes >= excludeCapacity) {
	// Temporary pointer to the old array
	Exclude* temp = excludes;	

	// Double the capacity
	excludeCapacity *= 2;

	// Create pointer to new array which is twice the size of the old one
	excludes = new Exclude[excludeCapacity];
		
	// Copy the old values into the new array
	for (int j = 0; j < numExcludes; j++)
	    excludes[j] = temp[j];

	// delete the old array
	delete[] temp;
    }

    // Add the bond to the exclude array
    // We must add it twice: Once for (ind1, ind2) and once for (ind2, ind1)
    Exclude ex1(ind1, ind2);
    excludes[numExcludes++] = ex1;
    Exclude ex2(ind2, ind1);
    excludes[numExcludes++] = ex2;
    
}    

void Configuration::buildExcludeMap() {
    // Call compareExcludeIndex with qsort to sort the excludes by BOTH ind1 AND ind2
    std::sort(excludes, excludes + numExcludes, compare());

    /* Each particle may have a varying number of excludes
     * excludeMap is an array with one element for each particle
     * which keeps track of where a particle's excludes are stored
     * in the excludes array.
     * 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+num_rb_attached_particles];
    for (int i = 0; i < num+num_rb_attached_particles; i++) {
	excludeMap[i].x = -1;
	excludeMap[i].y = -1;
    }
    int currPart = -1;
    int lastPart = -1;
    for (int i = 0; i < numExcludes; i++) {
	if (excludes[i].ind1 != currPart) {
	    currPart = excludes[i].ind1;
	    assert(currPart < num+num_rb_attached_particles);
	    excludeMap[currPart].x = i;
	    if (lastPart >= 0)
		excludeMap[lastPart].y = i;
	    lastPart = currPart;
	}
    }
    if (excludeMap[lastPart].x > 0)
	excludeMap[lastPart].y = numExcludes;
}

void Configuration::readAngles() {
	FILE* inp = fopen(angleFile.val(), "r");
	char line[256];
	int capacity = 256;
	numAngles = 0;
	angles = new Angle[capacity];

	// If the angle file cannot be found, exit the program
	if (inp == NULL) {
		printf("WARNING: Could not open `%s'.\n", angleFile.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 ANGLE inputs have 5 tokens
		// ANGLE | INDEX1 | INDEX2 | INDEX3 | FILENAME
		// Any angle input line without exactly 5 tokens should be discarded
		if (numTokens != 5) {
			printf("WARNING: Invalid 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_name = tokenList[4];
		//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)
			continue;

		if (numAngles >= capacity) {
			Angle* temp = angles;
			capacity *= 2;
			angles = new Angle[capacity];
			for (int i = 0; i < numAngles; i++)
				angles[i] = temp[i];
			delete[] temp;
		}

		Angle a(ind1, ind2, ind3, file_name);
		angles[numAngles++] = a;
		delete[] tokenList;
	}
	std::sort(angles, angles + numAngles, compare());	

	// for(int i = 0; i < numAngles; i++)
	// 	angles[i].print();
}

void Configuration::readDihedrals() {
	FILE* inp = fopen(dihedralFile.val(), "r");
	char line[256];
	int capacity = 256;
	numDihedrals = 0;
	dihedrals = new Dihedral[capacity];

	// If the dihedral file cannot be found, exit the program
	if (inp == NULL) {
		printf("WARNING: Could not open `%s'.\n", dihedralFile.val());
		printf("This simulation will not use dihedrals.\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 DIHEDRAL inputs have 6 tokens
		// DIHEDRAL | INDEX1 | INDEX2 | INDEX3 | INDEX4 | FILENAME
		// Any angle input line without exactly 6 tokens should be discarded
		if (numTokens != 6) {
			printf("WARNING: Invalid dihedral 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 (numDihedrals >= capacity) {
			Dihedral* temp = dihedrals;
			capacity *= 2;
			dihedrals = new Dihedral[capacity];
			for (int i = 0; i < numDihedrals; ++i)
				dihedrals[i] = temp[i];
			delete[] temp;
		}

		Dihedral d(ind1, ind2, ind3, ind4, file_name);
		dihedrals[numDihedrals++] = d;
		delete[] tokenList;
	}
	std::sort(dihedrals, dihedrals + numDihedrals, compare());	

	// for(int i = 0; i < numDihedrals; i++)
	// 	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 != 8) {
			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());
		int ind4 = atoi(tokenList[4].val());
		String file_name1 = tokenList[5];
		String file_name2 = tokenList[6];
		String file_name3 = tokenList[7];
		//printf("file_name %s\n", file_name.val());
		if (ind1 >= num or ind2 >= num or ind3 >= num or ind4 >= num)
			continue;

		if (numBondAngles >= capacity) {