# Create test forcefield (compatible with reference forces)
Forcefield waterff = newFF("Ethanol Test");
waterff.units = "kj";

# Create type definitions
waterff.addType(1,"CM","CM",C, "-C,nh=3", "Methyl carbon");
waterff.addType(2,"CT","CT",C, "-C,-O,nh=2", "Carbon with OH");
waterff.addType(3,"OH_a","OH",O, "", "Alcohol oxygen");
waterff.addType(4,"HC","HC",H, "-C", "Aliphatic hydrogen");
waterff.addType(5,"HO","HO",H, "-O", "Alcoholic hydrogen");
waterff.addInter("ljgeom", 1, -0.18,  0.276144, 3.50);
waterff.addInter("ljgeom", 2,  0.145, 0.276144, 3.50); 
waterff.addInter("ljgeom", 3, -0.683, 0.711280, 3.12);
waterff.addInter("ljgeom", 4,  0.06,  0.125520, 2.50);
waterff.addInter("ljgeom", 5,  0.418, 0.0,      0.0);
waterff.addBond("harmonic", "CM", "CT", 2242.624, 1.529);
waterff.addBond("harmonic", "CM", "HC", 2845.120, 1.090);
waterff.addBond("harmonic", "CT", "HC", 2845.120, 1.090);
waterff.addBond("harmonic", "CT", "OH", 2677.760, 1.410);
waterff.addBond("harmonic", "OH", "HO", 4627.504, 0.945);
waterff.addAngle("harmonic", "CT", "CM", "HC", 313.800, 110.7000);
waterff.addAngle("harmonic", "HC", "CM", "HC", 276.144, 107.8000);
waterff.addAngle("harmonic", "CM", "CT", "OH", 418.400, 109.5000);
waterff.addAngle("harmonic", "CM", "CT", "HC", 313.800, 110.7000);
waterff.addAngle("harmonic", "OH", "CT", "HC", 292.880, 109.5000);
waterff.addAngle("harmonic", "HC", "CT", "HC", 276.144, 107.8000);
waterff.addAngle("harmonic", "CT", "OH", "HO", 460.240, 108.5000);
waterff.addTorsion("cos3", "HC", "C*", "C*", "HC", 0.0, 0.0, 1.2552);
waterff.addTorsion("cos3", "HC", "C*", "C*", "OH", 0.0, 0.0, 1.9581);
waterff.addTorsion("cos3", "CM", "CT", "OH", "HO", -1.4895, -0.7280, 2.0585);
waterff.addTorsion("cos3", "HC", "CT", "OH", "HO", 0.0, 0.0, 1.4744);
waterff.finalise();

# Load reference forces models - Note that all forces in these files are in units of 10J/mol rather than kJ/mol
aten.prefs.zMap = "ff";
loadModel("data/test/ethanol-forces-elec.CONFIG");
Model elecref = aten.model;
loadModel("data/test/ethanol-forces-vdw.CONFIG");
Model vdwref = aten.model;
loadModel("data/test/ethanol-forces-intra.CONFIG");
Model intraref = aten.model;

# Load another copy of one of the reference models so we have the coordinates
loadModel("data/test/ethanol-forces-vdw.CONFIG");

# Check various force components
aten.prefs.elecCutoff = 7.0;
aten.prefs.vdwCutoff = 7.0;
aten.prefs.elecMethod = "ewald";
aten.prefs.ewaldAlpha = 0.46582;
aten.prefs.ewaldKMax = {8,8,8}; 
Vector v;

# Electrostatics (via Ewald sum)
aten.prefs.calculateIntra = FALSE;
aten.prefs.calculateVdw = FALSE;
modelForces();
double rmse_elec = 0.0;
for (int i=1; i<=aten.model.nAtoms; ++i)
{
	v = aten.model.atoms[i].f - elecref.atoms[i].f/100.0;
	rmse_elec += v.x*v.x + v.y*v.y + v.z*v.z;
}
rmse_elec = sqrt(rmse_elec / aten.model.nAtoms);

# Intramolecular terms
aten.prefs.elecMethod = "none";
aten.prefs.calculateIntra = TRUE;
aten.prefs.calculateVdw = FALSE;
modelForces();
double rmse_intra = 0.0;
for (int i=1; i<=aten.model.nAtoms; ++i)
{
	v = aten.model.atoms[i].f - intraref.atoms[i].f/100.0;
	rmse_intra += v.x*v.x + v.y*v.y + v.z*v.z;
}
rmse_intra = sqrt(rmse_intra / aten.model.nAtoms);

# Short-range (vdW)
aten.prefs.elecMethod = "none";
aten.prefs.calculateIntra = FALSE;
aten.prefs.calculateVdw = TRUE;
modelForces();
double rmse_vdw = 0.0;
for (int i=1; i<=aten.model.nAtoms; ++i)
{
	v = aten.model.atoms[i].f - vdwref.atoms[i].f/100.0;
	rmse_vdw += v.x*v.x + v.y*v.y + v.z*v.z;
}
rmse_vdw = sqrt(rmse_vdw / aten.model.nAtoms);

printf("Electrostatic  force RMSE = %f\n", rmse_elec);
printf("Short-range    force RMSE = %f\n", rmse_vdw);
printf("Intramolecular force RMSE = %f\n", rmse_intra);

quit();

