# Script to test UFF parameter generation against various systems

printf("This script tests UFF type assignment and parameter generation against data taken from\nthe UFF implementation in MCCCS Towhee.\n\n");

# Load forcefield and switch to Kelvin for units of energy
loadFF("uff.ff");
aten.prefs.energyUnit = "K";
int n;
FFBound ffb;
int errors = 0;
double tol = 0.2;	# Percentage tolerance

# Define some commonly-used functions
int testBond(double tol, int i, int j, string si, string sj, double d1, double d2, double d3, double d4, double d5, double d6)
{
	int n, result = 0;
	double diff[6], data[6] = { d1, d2, d3, d4, d5, d6 };
	FFBound ffb = generateBond(i,j);
	if (ffb)
	{
		for (n=1; n<=6; ++n) diff[n] =  ffb.data[n] - data[n];
		printf("            Bond %s-%s: generated (expected) parameters are:\n", si, sj);
		for (n=1; n<=ffb.nParams; ++n)
		{
			printf("    %2i %12s %14.5f (%14.5f)   [error = %8.4f%%]", n, ffb.dataName[n], ffb.data[n], data[n], (data[n]/ffb.data[n] * 100.0)-100.0);
			if (abs(data[n]/ffb.data[n]*100.0-100.0) < tol) printf("  Ok\n");
			else { printf("  *** ERROR\n"); result++; }
		}
	}
	else
	{
		printf("   Failed to generate bond!\n");
		result = 1;
	}
	return result;
}
int testAngle(double tol, int i, int j, int k, string si, string sj, string sk, double d1, double d2, double d3, double d4, double d5, double d6)
{
	int n, result = 0;
	double diff[6], data[6] = { d1, d2, d3, d4, d5, d6 };
	FFBound ffb = generateAngle(i,j,k);
	if (ffb)
	{
		for (n=1; n<=6; ++n) diff[n] =  ffb.data[n] - data[n];
		printf("        Angle %s-%s-%s: generated (expected) parameters are:\n", si, sj, sk);
		for (n=1; n<=ffb.nParams; ++n)
		{
			printf("    %2i %12s %14.5f (%14.5f)   [error = %8.4f%%]", n, ffb.dataName[n], ffb.data[n], data[n], (data[n]/ffb.data[n] * 100.0)-100.0);
			if (abs(data[n]/ffb.data[n]*100.0-100.0) < tol) printf("  Ok\n");
			else { printf("  *** ERROR\n"); result++; }
		}
	}
	else
	{
		printf("   Failed to generate angle!\n");
		result = 1;
	}
	return result;
}
int testTorsion(double tol, int i, int j, int k, int l, string si, string sj, string sk, string sl,  double d1, double d2, double d3, double d4, double d5, double d6)
{
	int n, result = 0;
	double diff[6], data[6] = { d1, d2, d3, d4, d5, d6 };
	FFBound ffb = generateTorsion(i,j,k,l);
	if (ffb)
	{
		for (n=1; n<=6; ++n) diff[n] =  ffb.data[n] - data[n];
		printf("   Torsion %s-%s-%s-%s: generated (expected) parameters are:\n", si, sj, sk, sl);
		for (n=1; n<=ffb.nParams; ++n)
		{
			printf("    %2i %12s %14.5f (%14.5f)   [error = %8.4f%%]", n, ffb.dataName[n], ffb.data[n], data[n], (data[n]/ffb.data[n] * 100.0)-100.0);
			if (abs(data[n]/ffb.data[n]*100.0-100.0) < tol) printf("  Ok\n");
			else { printf("  *** ERROR\n"); result++; }
		}
	}
	else
	{
		printf("   Failed to generate torsion!\n");
		result = 1;
	}
	return result;
}
int testType(int id, string name)
{
	int result = 1;
	atom i = aten.model.atoms[id];
	if (i == NULL) printf("Invalid atom ID (%i) passed to testType.", id);
	else if (i.type == NULL) printf("  Atom %2i does not have an assigned type.", id);
	else if (i.type.name != name) printf("  Atom %2i assigned type DOES NOT match expected type (%s)", id, name); 
	else
	{
		result = 0;
		printf("  Atom %2i assigned type matches expected type (%s)", id, name);
	}
	if (result == 0) printf("  Ok\n");
	else printf("  *** ERROR\n");
	return result;
}
int testVdw(double tol, int id, string si, double d1, double d2, double d3, double d4, double d5, double d6)
{
	int n, result = 0;
	double diff[6], data[6] = { d1, d2, d3, d4, d5, d6 };
	ffatom ffa = generateVdw(id);
	if (ffa)
	{
		for (n=1; n<=6; ++n) diff[n] =  ffa.data[n] - data[n];
		printf("          Atomtype %5s: generated (expected) parameters are:\n", si);
		for (n=1; n<=ffa.nParams; ++n)
		{
			printf("    %2i %12s %14.5f (%14.5f)   [error = %8.4f%%]", n, ffa.dataName[n], ffa.data[n], data[n], (data[n]/ffa.data[n]*100.0)-100.0);
			if (abs(data[n]/ffa.data[n]*100.0-100.0) < tol) printf("  Ok\n");
			else { printf("  *** ERROR\n"); result++; }
		}
	}
	else
	{
		printf("  Failed to generate VDW parameters.\n");
		result = 1;
	}
	return result;
}

printf("\n=================================\n");
printf("Test 1) Water\n");
printf("=================================\n");
string water_types[3] = { "O_3", "H_", "H_" };
loadModel("data/test/water.xyz");
typeModel();
# Check assigned atom types
for (n=1; n<=aten.model.nAtoms(); ++n) errors += testType(n, water_types[n]);
# Create and check intramolecular terms
errors += testBond(tol,1,2, water_types[1], water_types[2], 281798.9*2.0, 0.9903, 0, 0, 0, 0 ); 
errors += testAngle(tol,2,1,3, water_types[2], water_types[1], water_types[3], 60637.511, 0.300, 0.267, 0.267, 0.0, 0.0 );
# Create and check VDW types
errors += testVdw(tol,1, water_types[1], 30.193, 3.118, 0.0, 0.0, 0.0, 0.0);
errors += testVdw(tol,2, water_types[2], 22.142, 2.571, 0.0, 0.0, 0.0, 0.0);

printf("\n=================================\n");
printf("Test 2) Methane\n");
printf("=================================\n");
string methane_types[5] = { "C_3", "H_", "H_", "H_", "H_" };
newModel("methane");
newAtom(C,0,0,0);
addHydrogen();
typeModel();
# Check assigned atom types
for (n=1; n<=aten.model.nAtoms(); ++n) errors += testType(n, methane_types[n]);
# Create and check intramolecular terms
errors += testBond(tol,1,2, methane_types[1], methane_types[2], 166599.6*2.0, 1.1094, 0, 0, 0, 0 ); 
errors += testAngle(tol,2,1,3, methane_types[2], methane_types[1], methane_types[3], 37992.232, 0.344, 0.375, 0.281, 0.0, 0.0 );
# Create and check VDW types
errors += testVdw(tol,1, methane_types[1], 52.838, 3.431, 0.0, 0.0, 0.0, 0.0);
errors += testVdw(tol,2, methane_types[2], 22.142, 2.571, 0.0, 0.0, 0.0, 0.0);

printf("\n=================================\n");
printf("Test 3) Benzene\n");
printf("=================================\n");
string benzene_types[12] = { "C_R", "H_", "C_R", "H_", "C_R", "H_", "C_R", "H_", "C_R", "H_", "C_R", "H_" };
loadScript("data/scripts/benzene.txt", "benzene");
runScript("benzene");
typeModel();
# Check assigned atom types
for (n=1; n<=aten.model.nAtoms(); ++n) errors += testType(n, benzene_types[n]);
# Create and check intramolecular terms (N.B. fourth torsion parameter 's' is opposite sign since Towhee uses (1-cos) term in equation)
errors += testBond(tol,1,2, benzene_types[1], benzene_types[2], 179869.9*2.0, 1.0814, 0, 0, 0, 0 ); 
errors += testBond(tol,1,3, benzene_types[1], benzene_types[3], 232815.7*2.0, 1.3793, 0, 0, 0, 0 ); 
errors += testAngle(tol,2,1,3, benzene_types[2], benzene_types[1], benzene_types[3], 57657.568, 0.500, 0.667, 0.333, 0.0, 0.0 );
errors += testAngle(tol,1,3,5, benzene_types[1], benzene_types[3], benzene_types[5], 112013.508, 0.500, 0.667, 0.333, 0.0, 0.0 );
errors += testTorsion(tol,1,3,5,7, benzene_types[1], benzene_types[3], benzene_types[5], benzene_types[7], 6780.452, 2.0, 0.0, -1.0, 0.0, 0.0 );
# Create and check VDW types
errors += testVdw(tol,1, benzene_types[1], 52.838, 3.431, 0.0, 0.0, 0.0, 0.0);
errors += testVdw(tol,2, benzene_types[2], 22.142, 2.571, 0.0, 0.0, 0.0, 0.0);

printf("\n=================================\n");
printf("Results of all tests....\n");
printf("=================================\n");
if (errors == 0) printf("\n *** All parameters generated correctly ***\n\n");
else printf("\n !!! Failed with %i error(s) !!!\n\n", errors);

quit();
