Configuration.cpp

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

    gsl_rng *gslcpp_rng = gsl_rng_alloc(gsl_rng_default);
    srand(time(NULL));
    gsl_rng_set (gslcpp_rng, rand() % 100000);

    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(gsl_ran_gaussian(gslcpp_rng,sigma),gsl_ran_gaussian(gslcpp_rng,sigma),gsl_ran_gaussian(gslcpp_rng,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;
        }
    }

    gsl_rng_free(gslcpp_rng);
    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;
}