Configuration.cpp

		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;
		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 (numTokens != 6) {
			printf("Warning: Invalid particle file line: %s\n", line);
			return;
		}

		// 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::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 || ind2 < 0 || ind2 >=num) {
			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];
	for (int i = 0; i < num; 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 || ind2 >= num) {
	printf("WARNING: Attempted to add an exclusion for an out-of-range particle index (%d or %d >= %d).\n", ind1, ind2, num);
	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];
    for (int i = 0; i < num; 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);
	    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 or ind2 >= num or ind3 >= num)
			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 or ind2 >= num
				or ind3 >= num or ind4 >= num)
			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::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());
}

//populate the type list and serial list
void Configuration::populate() {
    for (int repID = 0; repID < simNum; ++repID) {
                const int offset = repID * num;
                int pn = 0;
                int p = 0;
                for (int i = 0; i < num; ++i) {
                        type[i + offset] = p;
                        serial[i + offset] = currSerial++;

                        if (++pn >= part[p].num) {
                                p++;
                                pn = 0;
                        }
                }
        }
}

bool Configuration::readBondFile(const String& value, int currBond) {
	int numTokens = value.tokenCount();
	if (numTokens != 1) {
		printf("ERROR: Invalid tabulatedBondFile: %s, numTokens = %d\n", value.val(), numTokens);
		return false;
	}

	String* tokenList = new String[numTokens];
	value.tokenize(tokenList);
	if (tokenList == NULL) {
		printf("ERROR: Invalid tabulatedBondFile: %s; tokenList is NULL\n", value.val());
		return false;
	}

	bondTableFile[currBond] = tokenList[0];

	// printf("Tabulated Bond Potential: %s\n", bondTableFile[currBond].val() );

	return true;
}

bool Configuration::readAngleFile(const String& value, int currAngle) {
	int numTokens = value.tokenCount();
	if (numTokens != 1) {
		printf("ERROR: Invalid tabulatedAngleFile: %s, numTokens = %d\n", value.val(), numTokens);
		return false;
	}

	String* tokenList = new String[numTokens];
	value.tokenize(tokenList);
	if (tokenList == NULL) {
		printf("ERROR: Invalid tabulatedAngleFile: %s; tokenList is NULL\n", value.val());
		return false;
	}

	angleTableFile[currAngle] = tokenList[0];

	// printf("Tabulated Angle Potential: %s\n", angleTableFile[currAngle].val() );

	return true;
}

bool Configuration::readDihedralFile(const String& value, int currDihedral) {
	int numTokens = value.tokenCount();
	if (numTokens != 1) {
		printf("ERROR: Invalid tabulatedDihedralFile: %s, numTokens = %d\n", value.val(), numTokens);
		return false;
	}

	String* tokenList = new String[numTokens];
	value.tokenize(tokenList);
	if (tokenList == NULL) {
		printf("ERROR: Invalid tabulatedDihedralFile: %s; tokenList is NULL\n", value.val());
		return false;
	}

	dihedralTableFile[currDihedral] = tokenList[0];

	// printf("Tabulated Dihedral Potential: %s\n", dihedralTableFile[currDihedral].val() );

	return true;
}
//Load the restart coordiantes only
void Configuration::loadRestart(const char* file_name) {
	char line[STRLEN];
	FILE* inp = fopen(file_name, "r");

	if (inp == NULL) {
		printf("GrandBrownTown:loadRestart File `%s' does not exist\n", file_name);
		exit(-1);
	}

	int count = 0;
	while (fgets(line, STRLEN, inp) != NULL) {
		// Ignore comments.
		int len = strlen(line);
		if (line[0] == '#') continue;
		if (len < 2) continue;

		String s(line);
		int numTokens = s.tokenCount();
		if (numTokens != 4) {
			printf("GrandBrownTown:loadRestart Invalid coordinate file line: %s\n", line);
			fclose(inp);	
			exit(-1);
		}

		String* tokenList = new String[numTokens];
		s.tokenize(tokenList);
		if (tokenList == NULL) {
			printf("GrandBrownTown:loadRestart Invalid coordinate file line: %s\n", line);
			fclose(inp);
			exit(-1);
		}

		int typ = atoi(tokenList[0]);
		float x = (float) strtod(tokenList[1],NULL);
		float y = (float) strtod(tokenList[2],NULL);
		float z = (float) strtod(tokenList[3],NULL);

		pos[count] = Vector3(x,y,z);
		type[count] = typ;
		serial[count] = currSerial;
		currSerial++;
		if (typ < 0 || typ >= numParts) {
			printf("GrandBrownTown:countRestart Invalid particle type: %d\n", typ);
			fclose(inp);
			exit(-1);
		}

		count++;
		delete[] tokenList;
	}

	fclose(inp);    
}
//Han-Yi Chou
//First the resart coordinates should be loaded
void Configuration::loadRestartMomentum(const char* file_name) 
{
    char line[STRLEN];
    FILE* inp = fopen(file_name, "r");

    if (inp == NULL) 
    {
        printf("GrandBrownTown:loadRestart File `%s' does not exist\n", file_name);
        exit(-1);
    }
    if(!loadedCoordinates)
    {
        printf("First load the restart coordinates\n");
        assert(1==2);
    }
    int count = 0;
    while (fgets(line, STRLEN, inp) != NULL) 
    {
        // Ignore comments.
        int len = strlen(line);
        if (line[0] == '#') continue;
        if (len < 2) continue;

        String s(line);
        int numTokens = s.tokenCount();
        if (numTokens != 4) 
        {
            printf("GrandBrownTown:loadRestart Invalid momentum file line: %s\n", line);
            fclose(inp);
            exit(-1);
        }

        String* tokenList = new String[numTokens];
        s.tokenize(tokenList);
        if (tokenList == NULL) 
        {
            printf("GrandBrownTown:loadRestart Invalid momentum file line: %s\n", line);
            fclose(inp);
            exit(-1);
        }

        int typ = atoi(tokenList[0]);
        float x = (float) strtod(tokenList[1],NULL);
        float y = (float) strtod(tokenList[2],NULL);
        float z = (float) strtod(tokenList[3],NULL);

        if (typ < 0 || typ >= numParts) 
        {
            printf("GrandBrownTown:countRestart Invalid particle type : %d\n", typ);
            fclose(inp);
            exit(-1);
        }

        if(typ != type[count])
        {
            printf("Inconsistent in momentum file with the position file\n");
            fclose(inp);
            exit(-1);
        }
        momentum[count] = Vector3(x,y,z);
        ++count;
        delete[] tokenList;
    }
    fclose(inp);
}

bool Configuration::loadCoordinates(const char* file_name) {
	char line[STRLEN];
	FILE* inp = fopen(file_name, "r");

	if (inp == NULL) return false;

	int count = 0;
	while (fgets(line, STRLEN, inp) != NULL) {
		// Ignore comments.
		int len = strlen(line);
		if (line[0] == '#') continue;
		if (len < 2) continue;

		String s(line);
		int numTokens = s.tokenCount();
		if (numTokens != 3) {
			printf("ERROR: Invalid coordinate file line: %s\n", line);
			fclose(inp);	
			return false;
		}

		String* tokenList = new String[numTokens];
		s.tokenize(tokenList);
		if (tokenList == NULL) {
			printf("ERROR: Invalid coordinate file line: %s\n", line);
			fclose(inp);
			return false;
		}

		if (count >= num*simNum) {
			printf("WARNING: Too many coordinates in coordinate file %s.\n", file_name);
			fclose(inp);
			return true;
		}

		float x = (float) strtod(tokenList[0],NULL);
		float y = (float) strtod(tokenList[1],NULL);
		float z = (float) strtod(tokenList[2],NULL);
		pos[count] = Vector3(x,y,z);
		count++;

		delete[] tokenList;
	}
	fclose(inp);

	if (count < num) {
		printf("ERROR: Too few coordinates in coordinate file.\n");
		return false;
	}
	return true;
}
//Han-Yi Chou The function populate should be called before entering this function
bool Configuration::loadMomentum(const char* file_name) 
{
    char line[STRLEN];
    FILE* inp = fopen(file_name, "r");

    if (inp == NULL) 
        return false;

    int count = 0;
    while (fgets(line, STRLEN, inp) != NULL) 
    {
        // Ignore comments.
        int len = strlen(line);
        if (line[0] == '#') 
            continue;
        if (len < 2) 
            continue;

        String s(line);
        int numTokens = s.tokenCount();
        if (numTokens != 3) 
        {
            printf("ERROR: Invalid momentum file line: %s\n", line);
            fclose(inp);
            return false;
        }

        String* tokenList = new String[numTokens];
        s.tokenize(tokenList);
        if (tokenList == NULL) 
        {
            printf("ERROR: Invalid momentum file line: %s\n", line);
            fclose(inp);
            return false;
        }

        if (count >= num) 
        {
            printf("WARNING: Too many momentum in momentum file %s.\n", file_name);
            fclose(inp);
            return false;
        }

        float x = (float) strtod(tokenList[0],NULL);
        float y = (float) strtod(tokenList[1],NULL);
        float z = (float) strtod(tokenList[2],NULL);
        momentum[count] = Vector3(x,y,z);
        ++count;
        delete[] tokenList;
    }
    fclose(inp);

    if (count < num) 
    {
        printf("ERROR: Too few momentum in momentum file.\n");
        return false;
    }
    return true;
}

// Count the number of atoms in the restart file.
int Configuration::countRestart(const char* file_name) {
	char line[STRLEN];
	FILE* inp = fopen(file_name, "r");

	if (inp == NULL) {
		printf("ERROR: countRestart File `%s' does not exist\n", file_name);
		exit(-1);
	}

	int count = 0;
	while (fgets(line, STRLEN, inp) != NULL) {
		int len = strlen(line);
		// Ignore comments.
		if (line[0] == '#') continue;
		if (len < 2) continue;

		String s(line);
		int numTokens = s.tokenCount();
		if (numTokens != 4) {
			printf("ERROR: countRestart Invalid coordinate file line: %s\n", line);
			fclose(inp);	
			exit(-1);
		}

		String* tokenList = new String[numTokens];
		s.tokenize(tokenList);
		if (tokenList == NULL) {
			printf("ERROR: countRestart Invalid coordinate file line: %s\n", line);
			fclose(inp);
			exit(-1);
		}

		int typ = atoi(tokenList[0]);
		// float x = strtod(tokenList[1],NULL);
		// float y = strtod(tokenList[2],NULL);
		// float z = strtod(tokenList[3],NULL);
		if (typ < 0 || typ >= numParts) {
			printf("ERROR: countRestart Invalid particle type: %d\n", typ);
			fclose(inp);
			exit(-1);
		}

		count++;
		delete[] tokenList;
	}

	fclose(inp);    
	return count;
}

bool Configuration::readTableFile(const String& value, int currTab) {
	int numTokens = value.tokenCount('@');
	if (numTokens != 3) {
		printf("ERROR: Invalid tabulatedFile: %s\n", value.val());
		return false;
	}

	String* tokenList = new String[numTokens];
	value.tokenize(tokenList, '@');
	if (tokenList == NULL) {
		printf("ERROR: Invalid tabulatedFile: %s\n", value.val());
		return false;
	}

	if (currTab >= numParts*numParts) {
	    printf("ERROR: Number of tabulatedFile entries exceeded %d*%d particle types.\n", numParts,numParts);
	    exit(1);
	}

	partTableIndex0[currTab] = atoi(tokenList[0]);
	partTableIndex1[currTab] = atoi(tokenList[1]);
	partTableFile[currTab] = tokenList[2];

	// printf("Tabulated Potential: %d %d %s\n", partTableIndex0[currTab],
	// 		partTableIndex1[currTab], partTableFile[currTab].val() );
	delete[] tokenList;
	return true;
}

void Configuration::getDebugForce() {
	// Allow the user to choose which force computation to use
	printf("\n");
	printf("(1) ComputeFull [Default]          (2) ComputeSoftcoreFull\n");
	printf("(3) ComputeElecFull                (4) Compute (Decomposed)\n");
	printf("(5) ComputeTabulated (Decomposed)  (6) ComputeTabulatedFull\n");