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