Configuration.cpp

		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) {
			BondAngle* temp = bondAngles;
			capacity *= 2;
			bondAngles = new BondAngle[capacity];
			for (int i = 0; i < numBondAngles; i++)
				bondAngles[i] = temp[i];
			delete[] temp;
		}

		BondAngle a(ind1, ind2, ind3, ind4, file_name1, file_name2, file_name3);
		bondAngles[numBondAngles++] = a;
		delete[] tokenList;
	}
	std::sort(bondAngles, bondAngles + numBondAngles, compare());

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

void Configuration::readProductPotentials() {
	FILE* inp = fopen(productPotentialFile.val(), "r");
	char line[256];
	int capacity = 256;
	numProductPotentials = 0;
	productPotentials = new ProductPotentialConf[capacity];

	// If the angle file cannot be found, exit the program
	if (inp == NULL) {
		printf("WARNING: Could not open `%s'.\n", productPotentialFile.val());
		printf("This simulation will not use product potentials.\n");
		return;
	}
	printf("DEBUG: READING PRODUCT POTENTAL FILE\n");
	std::vector<std::vector<int>> indices;
	std::vector<int> tmp;
	std::vector<String> pot_names;

	while(fgets(line, 256, inp)) {
		if (line[0] == '#') continue;
		String s(line);
		int numTokens = s.tokenCount();
		String* tokenList = new String[numTokens];
		s.tokenize(tokenList);

		indices.clear();
		tmp.clear();
		pot_names.clear();		    

		printf("\rDEBUG: reading line %d",numProductPotentials+1);

		// Legitimate ProductPotential inputs have at least 7 tokens
		// BONDANGLE | INDEX1 | INDEX2 | INDEX3 | [TYPE1] | POT_FILENAME1 | INDEX4 | INDEX5 | [TYPE2] POT_FILENAME2 ...
		if (numTokens < 5) {
		    printf("WARNING: Invalid product potential input line (too few tokens %d): %s\n", numTokens, line);
			continue;
		}

		// Discard any empty line
		if (tokenList == NULL)
			continue;

		SimplePotentialType type = BOND; // initialize to suppress warning
		bool type_specified = false;
		for (int i = 1; i < numTokens; ++i) {
		    char *end;
		    // printf("DEBUG: Working on token %d '%s'\n", i, tokenList[i].val());

		    // Try to convert token to integer
		    int index = (int) strtol(tokenList[i].val(), &end, 10);
		    if (tokenList[i].val() == end || *end != '\0' || errno == ERANGE) {
			// Failed to convert token to integer; therefore it must be a potential name or type

			// Try to match a type
			String n = tokenList[i];
			n.lower();
			if (n == "bond") { type = DIHEDRAL; type_specified = true; }
			else if (n == "angle")  { type = DIHEDRAL; type_specified = true; }
			else if (n == "dihedral")  { type = DIHEDRAL; type_specified = true; }
			else if (n == "vecangle") { type = VECANGLE; type_specified = true; }
			else { // Not a type, therefore a path to a potential
			    n = tokenList[i];
			    indices.push_back(tmp);
			    pot_names.push_back( n );
			    // TODO: Key should be tuple of (type,n)
			    std::string n_str = std::string(n.val());
			    if ( simple_potential_ids.find(n_str) == simple_potential_ids.end() ) {
				// Could not find fileName in dictionary, so read and add it
				unsigned int s = tmp.size();
				if (s < 2 || s > 4) {
				    printf("WARNING: Invalid product potential input line (indices of potential %d == %d): %s\n", i, s, line);
				    continue;
				}
				simple_potential_ids[ n_str ] = simple_potentials.size();
				if (not type_specified) type = s==2? BOND: s==3? ANGLE: DIHEDRAL;
				simple_potentials.push_back( SimplePotential(n.val(), type) );
			    }
			    tmp.clear();
			    type_specified = false;

			}
		    } else {
			if (index >= num) {
			    continue;
			}
			tmp.push_back(index);
		    }
		}

		if (numProductPotentials >= capacity) {
			ProductPotentialConf* temp = productPotentials;
			capacity *= 2;
			productPotentials = new ProductPotentialConf[capacity];
			for (int i = 0; i < numProductPotentials; i++)
				productPotentials[i] = temp[i];
			delete[] temp;
		}

		ProductPotentialConf a(indices, pot_names);
		productPotentials[numProductPotentials++] = a;
		delete[] tokenList;
	}
	printf("\nDEBUG: Sorting\n");
	std::sort(productPotentials, productPotentials + numProductPotentials, compare());

	// for(int i = 0; i < numAngles; i++)
	// 	angles[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 + num_rb_attached_particles + numGroupSites) 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) {
	    printf("ERROR: Could not open file for reading: %s\n", file_name);
	    exit(-1);
	    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");

	printf("WARNING: ");
	if (tabulatedPotential) {
		if (fullLongRange) printf("(6) was specified by config file\n");
		else printf("(5) was specified by config file\n");
	} else {
		if (fullLongRange != 0) printf("(%d) was specified by config file\n", fullLongRange);
		else printf("(4) was specified by config file\n");
	}

	char buffer[256];
	int choice;
	while (true) {
		printf("Choose a force computation (1 - 6): ");
		fgets(buffer, 256, stdin);
		bool good = sscanf(buffer, "%d", &choice) && (choice >= 1 && choice <= 6);
		if (good)
			break;
	}
	switch(choice) {
		case 1:
			tabulatedPotential = 0;
			fullLongRange = 1;
			break;
		case 2:
			tabulatedPotential = 0;
			fullLongRange = 2;
			break;
		case 3:
			tabulatedPotential = 0;
			fullLongRange = 3;
			break;
		case 4:
			tabulatedPotential = 0;
			fullLongRange = 0;
			break;
		case 5:
			tabulatedPotential = 1;
			fullLongRange = 0;
			break;
		case 6:
			tabulatedPotential = 1;
			fullLongRange = 1;
			break;
		default:
			tabulatedPotential = 0;
			fullLongRange = 1;
			break;
	}
	printf("\n");
}
//Han-Yi Chou setting boltzman distribution of momentum with a given center of mass velocity
//Before using this code, make sure the array type list and serial list are both already initialized
bool Configuration::Boltzmann(const Vector3& v_com, int N)
{
    int count = 0;
    Vector3 total_momentum = Vector3(0.);

    RandomCPU random = RandomCPU(seed + 2); /* +2 to avoid using same seed elsewhere */

    for(int i = 0; i < N; ++i)
    {
        int typ = type[i];
        double M = part[typ].mass;
        double sigma = sqrt(kT * M) * 2.046167337e4;
   
        Vector3 tmp = random.gaussian_vector() * sigma;
        tmp = tmp * 1e-4;
        total_momentum += tmp;
        momentum[(size_t)count] = tmp;
        ++count;
    }
    if(N > 1)
    {
        total_momentum = total_momentum / (double)N;

        for(int i = 0; i < N; ++i)
        {
            int typ = type[i];
            double M = part[typ].mass;
        
            momentum[i] = momentum[i] - total_momentum + M * v_com;
        }
    }

    
    return true;
} 

//////////////////////////
// Comparison operators //
//////////////////////////
bool Configuration::compare::operator()(const String& lhs, const String& rhs) {
	String* list_lhs = new String[lhs.tokenCount()];
	String* list_rhs = new String[rhs.tokenCount()];
	lhs.tokenize(list_lhs);
	rhs.tokenize(list_rhs);
	int key_lhs = atoi(list_lhs[1].val());
	int key_rhs = atoi(list_rhs[1].val());
	delete[] list_lhs;
	delete[] list_rhs;
	return key_lhs < key_rhs;
}

bool Configuration::compare::operator()(const Bond& lhs, const Bond& rhs) {
	int diff = lhs.ind1 - rhs.ind1;
	if (diff != 0)
		return lhs.ind1 < rhs.ind1;
	return lhs.ind2 < rhs.ind2;
}

bool Configuration::compare::operator()(const Exclude& lhs, const Exclude& rhs) {
	int diff = lhs.ind1 - rhs.ind1;
	if (diff != 0)
		return lhs.ind1 < rhs.ind1;
	return lhs.ind2 < rhs.ind2;
}

bool Configuration::compare::operator()(const Angle& lhs, const Angle& rhs) {
	int diff = lhs.ind1 - rhs.ind1;
	if (diff != 0)
		return lhs.ind1 < rhs.ind1;
	diff = lhs.ind2 - rhs.ind2;
	if (diff != 0)
		return lhs.ind2 < rhs.ind2;
	return lhs.ind3 < rhs.ind3;
}

bool Configuration::compare::operator()(const Dihedral& lhs, const Dihedral& rhs) {
	int diff = lhs.ind1 - rhs.ind1;
	if (diff != 0) 
		return lhs.ind1 < rhs.ind1;
	diff = lhs.ind2 - rhs.ind2;
	if (diff != 0) 
		return lhs.ind2 < rhs.ind2;
	diff = lhs.ind3 - rhs.ind3;
	if (diff != 0) 
		return lhs.ind3 < rhs.ind3;
	return lhs.ind4 < rhs.ind4;
}

bool Configuration::compare::operator()(const BondAngle& lhs, const BondAngle& rhs) {
	int diff = lhs.ind1 - rhs.ind1;
	if (diff != 0)
		return lhs.ind1 < rhs.ind1;
	diff = lhs.ind2 - rhs.ind2;
	if (diff != 0)
		return lhs.ind2 < rhs.ind2;
	diff = lhs.ind3 - rhs.ind3;
	if (diff != 0) 
		return lhs.ind3 < rhs.ind3;
	return lhs.ind4 < rhs.ind4;
}

bool Configuration::compare::operator()(const ProductPotentialConf& lhs, const ProductPotentialConf& rhs) {
    int diff = rhs.indices.size() - lhs.indices.size();
    if (diff != 0) return diff > 0;

    for (unsigned int i = 0; i < lhs.indices.size(); ++i) {
	diff = rhs.indices[i].size() - lhs.indices[i].size();
	if (diff != 0) return diff > 0;
    }

    for (unsigned int i = 0; i < lhs.indices.size(); ++i) {
	for (unsigned int j = 0; j < lhs.indices[i].size(); ++j) {
	    diff = rhs.indices[i][j] - lhs.indices[i][j];
	    if (diff != 0) return diff > 0;
	}
    }
    return true;
}