matlab

Homework # 7 Numerical Methods in Chemical Engineering

Submission due: Monday, December 6th, 2019 at 12:25 pm through NYU Classes

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. (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.

iii. (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.