Homework # 7 Numerical Methods in Chemical Engineering

Problem 1. (100 points) Nanoparticle self-assembly Interactions between nanoparticles can lead to self-assembled arrays that resemble the structure of crystals. In this problem, you will use Metropolis Monte Carlo simulations to describe the distribution of a set of non-charged nanoparticles. We will approximate the interaction between two particles using a Lenard-Jones potential: ??(??????) = 4?? [( ?? ?????? ) 12 ? ( ?? ?????? ) 6 ] ??????? ??????= |???? ? ????| where ra and rb are the positions of particle a and b respectively, and the values for ? and ? are provided below. You may assume that the sample at equilibrium follows a Boltzmann distribution with probability of finding the assemble in a state q, ??(??) ? exp(?????????(??)/(??????)), where ????????(??) = ? ? ??(??????) ?? ??=??+1 ?? ??=1 and q is the state vector containing the position of all the nanoparticles in the system. For simplicity, we will consider the self-assembly of a group of N spherical nanoparticles of radius R in two dimensions, and the state vector has the following structure, q = [x1 y1 x2 y2 …. xN yN] T .

i. (25 points) Write a function that computes the value of Utot for a given state vector q.

ii. function z = my_fun(x)

iii. z = zeros(size(x,1),3); % allocate output

iv. z(:,1) = x(:,1).^2 – 2*x(:,1).*x(:,2) + 6*x(:,1) + …

v. 4*x(:,2).^2 – 3*x(:,2);

vi. z(:,2) = x(:,1).*x(:,2) + cos(3*x(:,2)./(2+x(:,1)));

vii. z(:,3) = tanh(x(:,1) + x(:,2));

$displaystyle lim_{N rightarrow infty} frac{N_{xi}}{N} = P(xi),$

where P(?) is a given probability distribution (e.g., the Boltzmann distribution P(?) = Z$ ^{-1}$ exp[-?E(?)]) andN ? is the number of configurations ? (e.g., the number of configurations generated with a particular arrangement of [x1 y1 x2 y2 …. xN yN spins

viii. (40 points) Write a function that performs a Monte Carlo routine with an N_MC Monte Carlo steps, a constant temperature T, and an N number of nanoparticles in the initial state q_0. For each MC step, you should attempt moves of one nanoparticle at a time for all the particles. The function should output a matrix containing the values of the state vector at each MC step in different columns, and a vector containing Utot for each MC step.

$displaystyle P_{xi_n, xi_t} = begin{cases}R=frac{P(xi_t)}{P(xi_n)}, hsp…….2cm} P(xi_t) < P(xi_n), \ 1, hspace{.2cm} text{otherwise}. end{cases}$ The Metropolis Monte Carlo algorithm can be described as follows: function mypi = approxpi(n) % Input: n = number of points to generate % Default is n = 1e6 % Larger values of n should perform better if nargin < 1 n = 1e6; end % Generate uniformly distributed points in % [0, 1] x [0, 1] xy = rand(n, 2); % Compute distance from (0.5, 0.5) r = sqrt((xy(:,1)-0.5).^2+(xy(:,2)-0.5).^2); % Count fraction of points within 1/2 unit of (0.5, 0.5) frac = sum(r <= 0.5) / n; % Since square has side 1, circle has radius (1/2) % and should have area of pi*(1/2)^2 % frac is approximately pi/4 so pi is approximately 4*frac mypi = 4*frac; Step (1): Pick a configuration ?n Step (2): Pick a trial configuration ?t (usually a configuration similar to ?n) and compute the probability ratio R=p(?t)/p(?n) . Pick a random number p with value between 0 and 1. Make ?n+1=?t if p?R . Otherwise, make ?n+1=?t. Step (3): Go to (2) replacing ?n by ?n+1 . Step (3) is repeated N times, where N is a sufficiently large number. Note that, according to step (2), the probability of accepting a trial configuration ?t by making ?n+1=?t from a configuration ?n is $displaystyle P_{xi_n, xi_t} = begin{cases}R=frac{P(xi_t)}{P(xi_n)}, hsp.......2cm} P(xi_t) < P(xi_n), \ 1, hspace{.2cm} text{otherwise}. end{cases}$ ix. (35 points) Use your function to generate histograms containing the distribution of nanoparticle distances, Rab, for a set of 16 nanoparticles. Assume that the particles are initially arranged in a square grid of 4x4 nanoparticles, with center-to-center distance between the particles equal to 2R. Generate graphs for T = 100 K, 500 K, 1000 K. Use a minimum of 100,000 MC steps. Discuss your results. Physical parameters: R = 5 nm, ? = 5 nm, ?/kb = 5000 K.