# Script to test all of Aten's forcefield potentials for energy/force consistency

# Variables
int n, npoints = 500;
double delta;
Model m;
Forcefield ff;
double rmse, u[npoints], du[npoints], x[npoints];

#
# van der Waals Potentials
#
npoints = 500;
delta = 0.01;

m = newModel("Buckingham Test");
newAtom("C",0,0,0);
newAtom("C",1.2,0,0);
ff = newFF("Buckingham Test");
ff.addType(1,"C","C",C,"");
ff.addInter("buck",1,0.0,2000000.0,0.2,10000.0);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.vdwEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx +=delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : van der Waals : Buckingham    = %f (~ 0.18)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Lennard-Jones Test");
newAtom("C",0,0,0);
newAtom("C",2.0,0,0);
ff = newFF("Lennard-Jones Test");
ff.addType(1,"C","C",C,"");
ff.addInter("lj",1,0.0,0.65,3.20);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.vdwEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx += delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : van der Waals : LJ 12-6       = %f (~ 0.26)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Lennard-Jones AB Test");
newAtom("C",0,0,0);
newAtom("C",2.0,0,0);
ff = newFF("Lennard-Jones AB Test");
ff.addType(1,"C","C",C,"");
ff.addInter("ljab",1,0.0,2767011.611056,2576.980378);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.vdwEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx += delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : van der Waals : LJ AB         = %f (~ 0.25)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Morse Test");
newAtom("C",0,0,0);
newAtom("C",2.0,0,0);
ff = newFF("Morse Test");
ff.addType(1,"C","C",C,"");
ff.addInter("morse",1,0.0,75.0,1.0,3.0);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.vdwEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx += delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : van der Waals : Morse         = %f (~ 0.01)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Inverse Power Test");
newAtom("C",0,0,0);
newAtom("C",2.0,0,0);
ff = newFF("Inverse Power Test");
ff.addType(1,"C","C",C,"");
ff.addInter("inversepower",1,0.0,3.0,4.0,3.0);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.vdwEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx += delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : van der Waals : Inverse Power = %f (~ 0.00)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

#
# Bond Potentials
#
npoints = 100;
delta = 0.01;

m = newModel("Harmonic Test");
chain("C",0,0,0);
chain("C",1.0,0,0);
ff = newFF("Harmonic Test");
ff.addType(1,"C","C",C,"");
ff.addInter("lj",1,0.0,0.0,0.0);
ff.addBond("harmonic","C","C", 3000.0, 1.5);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.bondEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx += delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : bond          : Harmonic      = %f (~ 0.01)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Morse Test");
chain("C",0,0,0);
chain("C",1.0,0,0);
ff = newFF("Morse Test");
ff.addType(1,"C","C",C,"");
ff.addInter("lj",1,0.0,0.0,0.0);
ff.addBond("morse","C","C", 75.0, 1.0, 1.5);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = m.atoms[2].rx;
	u[n] = m.bondEnergy();
	du[n] = m.atoms[1].fx;
	m.atoms[2].rx += delta; 
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - (u[n+1]-u[n-1])/(x[n+1]-x[n-1]))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], (u[n+1]-u[n-1])/(x[n+1]-x[n-1]));
printf("RMSE for : bond          : Morse         = %f (~ 0.01)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);


#
# Angle Potentials
#
npoints = 500;
delta = 0.05;

m = newModel("Harmonic Test");
chain("C",0,0,0);
chain("C",1.0,0,0);
chain("C",1.0,1.0,0);
ff = newFF("Harmonic Test");
ff.addType(1,"C","C",C,"");
ff.addInter("lj",1,0.0,0.0,0.0);
ff.addBond("harmonic","C","C", 0.0, 0.0);
ff.addAngle("harmonic","C","C","C", 400.0, 109.47);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = geometry(1,2,3);
	u[n] = m.angleEnergy();
	du[n] = m.atoms[1].fy;
	setAngle(1,2,3,x[n]+delta);
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - ((u[n+1]-u[n-1])/(x[n+1]-x[n-1])*DEGRAD))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], ((u[n+1]-u[n-1])/(x[n+1]-x[n-1])*DEGRAD));
printf("RMSE for : bond          : Harmonic      = %f (~ 0.00)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Cosine Test");
chain("C",0,0,0);
chain("C",1.0,0,0);
chain("C",1.0,1.0,0);
ff = newFF("Cosine Test");
ff.addType(1,"C","C",C,"");
ff.addInter("lj",1,0.0,0.0,0.0);
ff.addBond("harmonic","C","C", 0.0, 0.0);
ff.addAngle("cos","C","C","C", 400.0, 3, 109.47, 1.0);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = geometry(1,2,3);
	u[n] = m.angleEnergy();
	du[n] = m.atoms[1].fy;
	setAngle(1,2,3,x[n]+delta);
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - ((u[n+1]-u[n-1])/(x[n+1]-x[n-1])*DEGRAD))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], ((u[n+1]-u[n-1])/(x[n+1]-x[n-1])*DEGRAD));
printf("RMSE for : bond          : Cosine        = %f (~ 0.00)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);

m = newModel("Cos2 Test");
chain("C",0,0,0);
chain("C",1.0,0,0);
chain("C",1.0,1.0,0);
ff = newFF("Cos2 Test");
ff.addType(1,"C","C",C,"");
ff.addInter("lj",1,0.0,0.0,0.0);
ff.addBond("harmonic","C","C", 0.0, 0.0);
ff.addAngle("cos2","C","C","C", 400.0, 1.0, 2.0, 3.0);
for (n=1; n<=npoints; ++n)
{
	modelForces();
	x[n] = geometry(1,2,3);
	u[n] = m.angleEnergy();
	du[n] = m.atoms[1].fy;
	setAngle(1,2,3,x[n]+delta);
}
rmse = 0.0;
for (n=2; n<npoints; ++n) rmse += (du[n] - ((u[n+1]-u[n-1])/(x[n+1]-x[n-1])*DEGRAD))^2;
#for (n=2; n<npoints; ++n) printf(" %f %f %f %f\n", x[n], u[n], du[n], ((u[n+1]-u[n-1])/(x[n+1]-x[n-1])*DEGRAD));
printf("RMSE for : bond          : Cos2          = %f (~ 0.00)\n", sqrt(rmse/(npoints-2)));
deleteModel();
deleteFF(ff);










# Done
quit();
