function [fAverage,accepted,q_all,f_val,f_val_sqr,RHH_val,X_all,Y_all,Z_all,V_all] = Metropolis(f,T,q0,N,del)
% computes average f(q) for Boltzmann V(q), temperature T (Kelvin)
% using N points starting from q0
%
% "f" is a function that return the value of interest <1/R_HH^6>
% T is the temperature
% q0 is the initial geometry of the system (6 degrees of freedom)
% N is the number of MC steps
% del is the maximum step size in any one coordinate
% O1 = (0,0,0)
% O2 = (xO2,0,0)
% H1 = (xH1,yH1,0);
% xO2=q(1); %all x,y,z in Angstroms
% xH1=q(2); yH1=q(3);
% xH2=q(4); yH2=q(5); zH2=q(6);
Delta=del*ones(size(q0)); % adjusting Delta changes heuristics
invkT=1/(1.3807e-11*T); %[=] 1/picoJoule
Vcurrent=V(q0); % V must have units pJ/molecule
q=q0;
w_current = q(1)^2*abs(q(3))*exp(-invkT*Vcurrent); %initial weighting factor
sizeq=size(q);
accepted=0; % initialize "accepted" at zero
% allocate vectors to store values
f_val = zeros(N,1);
fAverage = zeros(N,1);
RHH_val = zeros(N,1);
q_all = zeros(N,6);
X_all = zeros(N,4);
Y_all = zeros(N,4);
Z_all = zeros(N,4);
for i=1:N
qtry = q + Delta.*(rand(sizeq)-0.5)*2; % take a step to get a new orientation
[junk,Vtry] = V(qtry); % evaluate the energy in pJ at this new point
w_try = qtry(1)^2*abs(qtry(3))*exp(-invkT*Vtry); % the weight at the new point
if(Vtry rand(1) ) % accept an uphill step based on the ratio of weights
accepted = accepted+1;
f_val(i) = feval(f,qtry);
Vcurrent = Vtry;
q = qtry; % assign the new q value
else
% rejected, don't update the q value (take old one again)
f_val(i) = feval(f,q);
end
w_current = q(1)^2*abs(q(3))*exp(-invkT*Vcurrent); % calculate weight at new point
% store various quantities
q_all(i,:) = q';
V_all(i) = Vcurrent;
fAverage(i) = sum(f_val(1:i))/i; % this average is updated at each step
f_val_sqr(i) = f_val(i)^2; % this is the (1/R_HH^6)^2 value
% store the atom positions at each MC step to use in plots
% atoms: H1 O1 O2 H2
X_all(i,:) = [q(2) ,0 ,q(1) ,q(4)];
Y_all(i,:) = [q(3) ,0 ,0 ,q(5)];
Z_all(i,:) = [0 ,0 ,0 ,q(6)];
RHH_val(i) = (1/f_val(i))^(1/6); % calculate the H-H distance
end