% 10.34 - Fall 2006 % HW Set #10 % Prof. Green and Rob Ashcraft % 11/26/2006 % This problem is a Monte Carlo sampling of molecular geometries within a % given molecular potential. We want to compute a number of quantities, % including the <1/R_HH^6> and for a number of temperatures. % % The parameters in this code are the number of MC steps (N), the % temperature values at which to calculate results (T_vec), and the maximum % step size (del) in the +/- direction in cartesian coordinates. clear; clc; time_in = cputime; N = 100000; % number of MC steps T_vec = [300; 600; 1000; 5000]; % temperature values to run at del = 0.10; %maximum step size, Angstroms % determine the equilibrium structure via minimization of the V(q) q_init = [1;1;1;1;1;1]; options = optimset('display','iter','TolX',1e-5,'TolFun',1e-10); [q_equil] = fminunc(@V,q_init,[]); invR6_equil = RHH_6(q_equil); % calculate the <1/R^6> value for equilibrium using RHH_6 function R_HH_equil = (1/invR6_equil)^(1/6); % since there is only one value, we can use this formula to get R_HH disp(' '); disp(' '); disp('Equilibrium value at 0 K:') disp([' = ',num2str(R_HH_equil),' (Angstroms)']); disp(['<1/R_HH^6> = ',num2str(invR6_equil),' (1/A^6)']); %all x,y,z in Angstroms for the equilibrium geometry xO1 = 0; yO1 = 0; zO1 = 0; xO2 = q_equil(1); yO2 = 0; zO2 = 0; xH1 = q_equil(2); yH1 = q_equil(3); zH1 = 0; xH2 = q_equil(4); yH2 = q_equil(5); zH2 = q_equil(6); % equilibrium structure and plotting X = [xH1,xO1,xO2,xH2]; Y = [yH1,yO1,yO2,yH2]; Z = [zH1,zO1,zO2,zH2]; figure; stem3(X,Y,Z+1,'-','linewidth',2); hold on; plot3(X,Y,Z+1,'.-r','linewidth',4,'markersize',50); grid on; title('Equilibrium Geometry at T = 0 K'); axis([min(X)-.5 max(X)+.5 min(Y)-.5 max(Y)+.5 0 max(Z+1)+.5 ]); text(X(1)+0.2, Y(1)+0.2, Z(1)+1.2,'H_1') text(X(2)+0.2, Y(2)+0.2, Z(2)+1.2,'O_1') text(X(3)+0.2, Y(3)+0.2, Z(3)+1.2,'O_2') text(X(4)+0.2, Y(4)+0.2, Z(4)+1.2,'H_2') disp(' '); disp(' '); disp(['Number of MC steps taken: ',num2str(N)]); % Run the MC simulation for each temperature and generate some plots. for t=1:length(T_vec) T = T_vec(t); % Perform the Monte Carlo steps for each of the temperature values [fAverage(:,t),accepted,q_all,f_val,f_val_sqr,RHH_val,X_all,Y_all,Z_all,V_all] = Metropolis(@RHH_6,T,q_equil,N,del); % estimate the uncertainty from the fluctuations uncertainty(t,1) = sqrt( 1/N*(sum(f_val_sqr)/N - (sum(f_val)/N)^2) ); % estimate the uncertainty by taking the standard deviation of the % <1/R_HH^6> values for the final 50% of the MC points StDev(t,1) = std(fAverage(end/2:end,t)); disp(' '); disp(['<1/R_HH^6> for T = ',num2str(T_vec(t)),' K: ',num2str(fAverage(end,t)),' +/- ',num2str(uncertainty(t))]); disp(['<1/R_HH^6> for T = ',num2str(T_vec(t)),' K: ',num2str(fAverage(end,t)),' with sigma = ',num2str(StDev(t))]); disp([' for T = ',num2str(T_vec(t)),' K: ',num2str(mean(RHH_val))]); disp(['There were ',num2str(accepted),' steps accepted out of ',num2str(N),' total steps']); disp(' '); % Generate some plots, 3-D geometry plots and histograms figure; plot3(X,Y,Z,'.-r','linewidth',4,'markersize',50); grid on; hold on; plot3(X_all,Y_all,Z_all,'.'); hold off; title(['Results for T = ',num2str(T),' K and ',num2str(N),' MC steps']); figure; subplot(1,2,1); hist(RHH_val,linspace(min(RHH_val),max(RHH_val),round((max(RHH_val)-min(RHH_val))/0.01))); axis([0 5 0 1]); axis 'auto y'; ylabel('Frequency'); xlabel('H-H Distance (A)'); title(['H-H Distr. for T = ',num2str(T),' K and ',num2str(N),' MC steps']); subplot(1,2,2); hist(V_all,50); xlabel('Potential Energy (pJ)'); set(get(0,'CurrentFigure'),'position',[288 595 1344 503]); title(['Potential Distr. for T = ',num2str(T),' K and ',num2str(N),' MC steps']); [n_V,V_bin] = hist(V_all,50); % find the most frequent V(q) value from the histogram [val,loc] = max(n_V); V_max = V_bin(loc); disp(['Most probable MC V(q) = ',num2str(V_max),' pJ']); disp(['Value of kB*T = ',num2str(1.3807e-23*T*1e12),' pJ']); disp(['Time up to now = ',num2str(cputime-time_in),' seconds']); disp(' '); pause(2); end % plot the progression of the <1/R_HH^6> value as the # steps increases N_vec = linspace(1,N,N); figure; plot(N_vec,fAverage(:,1),'-',N_vec,fAverage(:,2),':',... N_vec,fAverage(:,3),'-.',N_vec,fAverage(:,4),'--') title(['Results for ',num2str(N),' MC steps']); xlabel('Step Number'); ylabel('Cummulative <1/R_H_H^6> Value (A^-^6)'); legend('T = 300 K','T = 600 K','T = 1000 K','T = 5000 K');