#Read data and integrate them
interface(quiet=true);
with(stats):with(linalg):
data:=readdata("data1.tmp",4):
data1:=map(u->[u[1],u[3]],data):
nn:=nops(data):
y0:=0.1+data[1][3]:
 for i from 1 to nn do
     data[i][3]:=data[i][3]-y0:
 od:
data2:=map(u->[u[1],log(-u[3])],data):
data3:=convert(transpose(convert(data2,array)),listlist):
fit2:=fit[leastsquare[[x,y], y=c0+c1*x+c2*x^2+c3*x^3+c4*x^4+c5*c^5, {c0,c1,c2,c3,c4,c5}]](data3):
f2:=unapply(rhs(fit2),x):
result:=evalf(int(-exp(f2(x))+y0,x=0..1)):
writeto("FreeE.res"):
printf("Integrated Free Energy: %15.10f\n",result):
writeto(terminal):
interface(quiet=false):


####data1.tmp
 0.0   -2.373772442322081     0.019110573465948  0.708009021872868
 0.1   -2.374423260444426     0.014947792481552  0.708009021872868
 0.4   -2.374655047081285    -0.002250677002858  0.708009021872868
 0.7   -2.375215606030175    -0.019147495800667  0.708009021872868
 0.8   -2.373858276444851    -0.033666054870812  0.708009021872868
 0.9   -2.374074142484935    -0.045489139855569  0.708009021872868
 1.00   -2.373570137950830    -0.064315275727469  0.708009021872868