matlab

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));

where P(Σ) is a given probability distribution (e.g., the Boltzmann distribution P(Σ) = Z 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.

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

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.