1

​

CL 413 Fundamentals of molecular Simulation

​

​

​

​

​

​

​

​

​

Chemical Engineering Department

Birla Institute of Technology Mesra, Ranchi

Monte Carlo Simulations

Molecular Simulation

2

• In the broadest sense, molecular modeling and simulation can be defined as the use of computational methods to describe the behavior of matter at the atomistic or molecular level.
• A molecular simulation generally consists of a computer realization of a system in which actual molecular configurations are used to extract structural, thermodynamic and dynamic information.
• The term “configuration” denotes a set of Cartesian coordinates (and momenta in the case of dynamic simulations) for all the atoms or molecules that constitute a system.
• Molecular simulation is based on classical physics and treats interactions between atoms with empirical potential functions (force fields).

Molecular Simulation

3

• Molecular dynamics and Monte Carlo Simulation.

​

• In molecular dynamics, Newton’s equations of motion are integrated to generate a trajectory for the system of interest.
• In Monte Carlo methods, probability distribution functions dictated by statistical mechanics are sampled to generate distinct configurations.
• MC algorithms vary widely depending on the system of interest, and they are not well suited for study of time-dependent processes.
• Macroscopic properties are found as the average of microscopic quantities.

Monte Carlo Simulation

4

• Monte Carlo (MC) simulation is a stochastic simulation method, which relies on the probability distribution of certain system parameters. There is some probability associated with that particular event and the selection of a particular probable event is decided by the random numbers.
• This is different from the Molecular Dynamics (MD) simulation method, which follows a deterministic approach (viz., using Newton's law of motion).

Metropolis Monte Carlo

5

• Consider a system which only contains O₂ molecules at 300 K.
• What would be the energy of each O₂ molecule present in the system?

​

Metropolis Monte Carlo

6

E

E₁

E₂

E₃

E₄

Metropolis Monte Carlo

7

E

E₁

E₂

E₃

E₄

Distribution of Molecular Speeds (Maxwell-Boltzmann Distribution)

8

Metropolis Monte Carlo

9

Assuming that a molecule is at an initial energy state, E_i and is supposed to be moved to final energy state E_f. What is the probability of moving from E_i to E_f?

​

Case I: E_f < E_i 🡪 Since the final energy is lower than or equal to the initial energy, the transition would take place.

​

Case II: E_f > E_i 🡪 if may appear that, since the final energy is higher than the initial energy, the transition would not take place. However, considering the energy distribution described above, there is a possibility to move from E_i to E_f.

Metropolis Monte Carlo

10

The probability is given by : p (i 🡪 f) = exp(-ΔE/k_BT), where ΔE = E_f – E_i

k_B= Boltzmann constant, T= temperature, K

Now, generate a random number, r and

check if r ≤ p (i 🡪 f) 🡪 move will be accepted,

Otherwise, the move will be rejected.

​

This technique is known as the Metropolis Monte Carlo.

​

Simulation Strategy

11

To simulate a system, we need to follow certain steps systematically to perform the simulation smoothly.

The steps are:

0) Understand the process/phenomena (at the macroscopic level) that are to

be simulated.

1. Sample set-up: Size of simulation box, number of particles in the box.
2. Framework: on-lattice or off-lattice system.
3. PBC or no PBC.
4. Modeling the potential energy function.
5. Simulation: algorithm to move molecules in the simulation box.
6. Calculation of properties and data analysis.

MC Simulation in various ensembles

12

• Canonical ensemble (NVT)
• Micro canonical ensemble (NVE)
• Isothermal-isobaric ensemble (NPT)
• Grand canonical ensemble (μVT)
• Gibbs Ensemble (phase equilibria)

​

Monte Carlo Moves

13

• Displacement/Translation move.
• Rotational move about the center-of-mass
• Volume change move.
• Interbox molecule transfer move.
• Insertion or deletion move.

Simulation �in�canonical (NVT) ensemble

14

NVT Ensemble:

15

• Canonical ensemble is the most widely used ensemble to carry out MC simulation.
• Simulation will begin with the pre-defined value of the number of molecules (N), box size (V = L³ ) and temperature (T).
• Once N is decided, generate initial coordinate file, which will contain coordinates ( x, y, z ) of all the molecules. It can be generated manually or by using some algorithm.

NVT Ensemble:

16

• For a generic model system, one can place the molecules sequentially in a cubic box maintaining the distance between any two molecules greater than the molecular diameter, σ.

​

• We can also randomly pick a coordinate (x, y, z) and assign to a particle.  Continue till all the particles are placed and inter-particle distance between any two particles is greater than or equal to the particle diameter. 

Steps:

17

1. Randomly select one molecule from N.
• Generate a random number, r(0,1)
• Multiply by N.
2. Get the initial co-ordinate (x_i, y_i, z_i) of M^th molecule.
3. Use the predefined value of Δr_max (maximum displacement along x, y and z direction) to get the final position. Δr to be chosen from -Δr_max to + Δr_max.
• Generate a random number, r_x (0,1)
• Get Δx from Δx = 2(r_x - 0.5) Δx_max 🡪 for x-coordinate.
• Similarly, for Δy and Δz. (Usually, Δx_max Δy_max and Δz_max are same).

Steps:

18

4. Get the final coordinates: x_f = x_i + Δx; y_f = y_i + Δy; z_f = z_i + Δz

5. Transform (x_f, y_f, z_f) into periodic cell by applying PBC.

6A. Calculate the distance r_m-j (where j = 1 to N, ≠ M between M^th and all other molecules.

6B. If all r_m-j is greater than or equal to σ, then calculate E_i and E_f.

6C. Calculate p ( i -> f) = exp(-βΔE).

6D. Generate a random number, m.

6E. If m ≤ p ( I -> f) 🡺 accept the move and update the co-ordinate for the M^th molecule.

Steps:

19

6F. Else reject the move and retain the initial co-ordinate of M^th molecule.

6G. Go to step (1) to select another M and repeat the whole process.

​

6H. If any of r_m-j is less that σ, reject the final co-ordinate and retain the initial co-ordinate of the M^th molecule and go to step (1).

​

Simulation in NPT ensemble

20

NPT Ensemble:

21

• This is one of the widely used ensembles.
• The simulation is carried out at a constant temperature and pressure, where volume is subject to vary.
• It is usually done by one volume change move coupled with the particle displacement moves in an MC step, as we did in NVT ensemble simulation.

Steps:

22

Steps:

23

Steps:

24

g) If the move accepted, update the coordinates and go to step 4 for the next MC step.

h) If the move rejected, keep the present the coordinates, and go to step 4 for the next MC step.

Grand canonical Monte Carlo (GCMC)

25

GCMC:

26

• In GCMC, the chemical potential (μ), volume (V), and temperature (T) are kept constant. Since N can vary during the simulation, molecules can leave the simulation box and new molecules can enter the simulation box.
• The possible MC moves are:
• Displacement of molecules.
• Insertion of molecules into the simulation box (Total number of molecules in the box will be incremented by one for each successful move).
• Deletion of a molecule from the simulation box (Total number of molecules in the box will be decreased by one for each successful move).

Steps:

27

1. Initialization: Assign the values for N, V(L), T, μ.
2. Get the initial coordinates of all the molecules and Calculate E_i.
3. Select the move from displacement, Insertion and Deletion. Each MC move will consist one of the above move.
1. Generate a random number, r.
2. If (r ≤ 1/3) 🡪 Displacement move.
3. Else if (r > 1/3 && r ≤ 2/3) 🡪 Insertion move.
4. Else 🡪 Deletion move.

Displacement move:

28

This follows the similar formalism as we do for NVT- ensemble simulation.

Insertion move:

29

Deletion move:

30

Estimation of the chemical potential (µ)

31

Steps:

32

• We can estimate the value of µ from an independent simulation (may be with NVT or NPT ensemble) for a given system. The method is known Widom particle insertion method.
• In this method, a “test” particle is inserted into the simulation box with N molecules and the change in total energy due to this inserted particle is calculated. The change in energy is equal to this inserted particle. The change in energy is equal to the interaction energy of the new particle with the existing N particle. Estimate the Boltzmann factor and a weighted average of the Boltzmann factor would give the value of µ.

Steps (NVT Ensemble):

33

1. Fix the value of N, V and T.
2. Read the initial coordinates of N particles.
3. Equilibrate using standard MC move (e.g., Displacement), used in simulating NVT ensemble.
4. Randomly choose a coordinate (x_t, y_t, z_t) for the test particle.
5. Calculate the distance between the test particle and all the N particles in the simulation 🡪 r_tj (j = 1 to N).
6. If r_tj < σ reject and try for new trial, go to step 4.
7. If r_tj ≥ σ, calculate the interaction energy between the test particle and all the N particles (E_t).

Steps (NVT Ensemble):

34

8. Choose another coordinate and estimate E_t and repeat for say, M-times i.e., M number of trial insertion.

9. Estimate the Boltzmann factor, exp(-βE_t) and sum them

🡪 sum = sum + exp(-βE_t)

10. Estimate, µ_ex = -K_BT ln(sum/M)

(µ = µ_id + µ_ex)

​

Note:

1. During M trial insertion, equilibration of the system at regular interval is required for a better result.
2. Test particle never being accepted. They are used to calculate the interaction energy only and not stored.
3. Limitation: For dense system, rate of successful insertion is very low.

Gibbs Ensemble

35

Gibbs Ensemble Monte Carlo Simulation

36

Gibbs Ensemble Monte Carlo Simulation

37

• Assume, we have N₁ and N₂ number of particles in box 1 and box 2, with volume V₁ and V₂ ( V₁ + V₂ = V = Constant).
• In one MC step, on average every particle should be attempted to move (via displacement); few volume changes and a few particle swapping moves.
• Let, the number of displacement move = (N₁ + N₂) including both the boxes.
• The number of volume change move = n_vol.
• The number of particle transfer move = n_trans.
• The total number of moves in a single MC step, n_tot = N₁ + N₂ + n_vol + n_trans.

Gibbs Ensemble Monte Carlo Simulation

38

Selection rule is as follows:

Generate a random number r ( 0 to 1).

If ( r ≤ N₁/n_tot)

🡪 Displace the particle in box 1.

else if ( r > N₁/n_tot && r ≤ (N₁ + N₂)/n_tot)

🡪 Displace the particle in box 2.

else if ( r > (N₁+N₂)/n_tot && r ≤ (N₁ + N₂ + n_vol)/n_tot)

🡪 Do the volume change.

else

🡪 Do the particle transfer.

Usually, n_vol = 1 is used per MC step while N_trans ~ 3% of (N₁ + N₂).

Implementation of moves: Initialization

39

1. Assign the value of N₁, N₂, V₁ and V₂.
2. Model the interaction energy, E(r).
3. Decide the cut-off distance for the calculation of energy.
4. PBC need to be implemented on both boxes individually 🡪 only for the displacement moves.
5. Decide the value of n_vol and n_trans.

Implementation details: Particle displacement in box 1

40

1. The procedure is exactly identical to the NVT simulation.
2. Randomly select one particle out of N₁.
3. Calculate E _i.
4. Give a random displacement Δr.
5. Get the final coordinate (use PBC with L1 as box size).
6. Calculate final energy, E_f.
7. Calculate the acceptance probability,
• p (i 🡪 f) = exp(-βΔE)
• ΔE = E_f(N₁) – E_i(N₁)
8. Generate a random number, r
• If r ≤ p (i → f)  →   Accept the move and update the coordinate.

9. Else reject the move.

Volume change move

41

Steps:

42

1. Randomly select which box to go for volume change, with equal probability.
2. If r ≤ 0.5 🡪 volume change in box 1

else volume change in box 2.

3. If volume change is positive in the box 1, then the volume change will be

negative in the box 2 and vice versa.

4. Randomly choose ΔV from - ΔV_max to +ΔV_max._.

5. Assume volume change in box 1:

a) V₁^f = V₁ⁱ + ΔV and,

b) for box 2: V₂^f = V₂ⁱ + ΔV; to keep the total volume constant.

Steps:

43

Particle Transfer move

44

Particle transfer move

45

1. Randomly select the box (1 or 2) from which the particle will be removed,

and will be inserted into the other box.

• Generate a random number (r).
• If r ≤ 0.5 🡪 select box 1.
• Else 🡪 select box 2.

2. Assume box 1 is selected.
3. Calculate the initial energy E_i (N).

Particle transfer move

46

Particle transfer move

47

9. Generate a random number, m.
10. If m ≤ p ( i -> f) 🡺 accept the move and update N₁ (N₁ - 1) and N₂ (N₂ + 1).
11. Else reject the move and retain the previous values.

​

Note:

1. This move is equivalent to the combined move of particle deletion and insertion in µVT ensemble.
2. Limitation: For strongly interacting system and dense system, particle transfer has very low acceptance probability.

​

​

​

Thank You

48

List of software for Monte Carlo Molecular Modeling

49

1. MCCCS Towhee (Monte Carlo for Complex Chemical Systems) – Prof. J. Ijja Siepmann, University of Minnesota.
2. GPU Optimized Monte Carlo (GOMC) – Prof. Jeffrey potoff, Wayne State University.
3. Cassandra – Prof. Edward J. Maginn, University of Notre Dame.
4. RASPA - Adsorption and Diffusion In Flexible Nanoporous Materials, David Dubbeldam (University of Amsterdam, The Netherlands), Sofia Calero (University Pablo de Olavide, Sevilla, Spain), Thijs Vlugt (Delft University of Technology, The Netherlands), and Randall Q. Snurr (Northwestern University, Evanston, USA).

Key Concepts

50

Force Field

51

The extent to which a classical molecular simulation accurately predicts thermophysical properties depends on the quality of the force field used to model the interactions in the fluid. 

Force Field

52

The forcefield contains the following building blocks for the calculations of energy:

​

1. A list of atoms.
2. A list of atom types.
3. functional forms for the components of the energy expression.
4. Parameters for the functional terms. 

Force Field

53

The forcefield contains the following building blocks for the calculations of energy:

​

1. OPLS (Optimized Potential for Liquid Simulations) – Yale University.
2. AMBER (Assisted Model Building and Energy Refinement) – for proteins.
3. CHARMM (chemistry at HARvard Molecular Mechanics) – Harvard University.
4. GROMOS (GROningen MOlecular Simulation).
5. TraPPE (Transferable Potentials for Phase Equilibria) – University of Minnesota.

​

Potential Energy Function

54