clear
echo on
format compact
%Solution to Problem 28
%
%In this problem we are asked to calculate
%the internal energy per molecule of three
%methane molecules confined to a box 10
%Angstroms on a side.  We need to do this
%to a precision of 1% (1 sigma error).  We
%can do this by writing a routine which, for
%some chosen number of points in a monte carlo
%integration routine, returns both the estimate
%of the internal energy and the error.  Then,
%by recognizing that the error goes as 1/N^.5,
%we can determine the number of points we have
%to use to get the desired accuracy.  We do all
%this in the attached program mcint.m.
pause
%In order to make our program more general, we
%shall make it dimensionless.  We will scale the
%lengths with the size of the box, and the energy
%with kT where k is the Boltzmann constant and
%T is the absolute temperature.  Thus:
T=273;
k=1.38062e-16; %erg / deg K
s=10; %Angstroms
%And the properties of methane:
eps=147.95; %deg K
sigma=3.73; %Angstroms
%Which yield the normalized values:
estar=eps/T;
sigstar=sigma/s;
pause
%We also pick the number of molecules in the box:
m=3;
%and the number of points:
n=10^3;
%which seems like a nice place to start.  Now we
%are ready to go:
[inte,error]=mcint(estar,sigstar,m,n)
pause
%We want to get the value to 1% accuracy, thus we
%need to increase the number of points:
n=n*(abs(error/inte)/0.01)^2
%So now we do it again:
[inte,error]=mcint(estar,sigstar,m,n)
%which yields the desired error this time:
relerror=abs(error/inte)
pause
%Note that the internal energy is negative, and close
%to zero.  It is negative because the molecules attract
%each other, and it is close to zero because it is a
%very low density gas.  The internal energy would be
%identically zero for an ideal gas.
%
%The actual internal energy per molecule is given by:
dimenergy=inte*k*T
%in ergs.




_________________________



function [inte,error]=mcint(estar,sigstar,m,n)
%This function uses Monte Carlo integration over
%n points to determine the internal energy per
%molecule of m molecules in a box of unit dimension.
%The parameters describing the molecules are given
%by estar and sigstar, and the interaction potential
%is given by the Lennard Jones 6-12 potential.
%We return both the energy and the error in the 
%energy.
sumn=0;
sumd=0;
sumn2=0;
sumd2=0;
sumnd=0; %This is needed for the covariance.
a=4*estar*sigstar^12;
b=4*estar*sigstar^6;
for i=1:n
 r=rand(3,m); %This gets the positions of all m molecules
 u=0;
 for j=1:m-1
  for k=j+1:m
   dist=norm(r(:,j)-r(:,k)); %The relative distance
   u=u+a/dist^12-b/dist^6;
  end
 end
 denom=exp(-u);
 num=u*denom;
 sumn=sumn+num;
 sumd=sumd+denom;
 sumn2=sumn2+num*num;
 sumd2=sumd2+denom*denom;
 sumnd=sumnd+num*denom;
end
inte=sumn/sumd/m;
varsumn=(n*sumn2-sumn*sumn)/(n-1);
varsumd=(n*sumd2-sumd*sumd)/(n-1);
covarsumnd=(n*sumnd-sumn*sumd)/(n-1);
error=inte*(varsumn/sumn^2+varsumd/sumd^2-2*covarsumnd/sumn/sumd)^.5;
%which is the standard formula for a ratio.