Materials Modeling at the Molecular Level

Lecture 2: Efficiency and Error Estimation

1. Review of Lecture 1 and Lab 1

2. More efficient neighbor finding

3. Estimating errors

Laura Frink

Sandia National Laboratories

Computational Materials Science Dept.

lfrink@cs.sandia.gov

lecture 1: An introduction to molecular modeling

lab 1: Getting to know MD and MC Codes

lecture 2: How to write efficient MC and MD codes

lab 2: Implementation of routines for improved efficiency and error estimation

lecture 3: What to do with an MD or MC Code ????

lab 3: Computing materials properties

Last Time

1. Overview of Materials Modeling in general;

and Molecular simulation in particular.

2. Introduction to Monte Carlo and Molecular Dynamics algorithms.

3. Got a feel for the components of MD and MC codes.

4. Computed pressures along two isotherms

(one sub-critical, one super-critical)

5. Found that the programs run slowly ... especially the MC code.

This Time

1. Improving code performance

2. Estimating errors

Force (MD) /Energy (MC) Computation

90% of the Effort

Loop over all particles in simulation volume

MD -- O(N2) -- have moved all N particles at once.

MC -- O(N) -- have only moved one particle at a time.

MD

loop over i (1,N)

loop over j (1,N)[

if (i¹ j) calculate f_ij

]

MC

loop over j (1,N)[

if (i¹ j) calculate u_ij

]

loop over i (1,N-1)

loop over j (i+1,N)[

]

MC: ???

Reduce the size of the loop over particles by keeping track of particles that are within the cut-off distance of interest.

• Neighbor lists

• Linked Cells

• Hybrid Approach - construct neighbor list from linked cells !!

• Goal: reduce scaling to O(N) for MD and O(1) for MC
• the different methods will yield different prefactors, c, with time = cN

Verlet Neighbor Lists

r_c = cutoff distance

• every M steps, create a list for each particle that contains all the particles within r_n.

• in between steps, use this reduced list to find atoms inside r_c

• neighbor list creation - O(N²) operation

MD: force calculation - O(N) operation

MC: energy calculation - O(1) operation

• 2 data structures
• nlist one long list of neighbors
• nlistpnt pointers (integer offsets) into the list

nlist = [2, 3, 4, 3, 4, 4, 7, 7]

nlistpnt = [1, 4, 6, 7, 7, 8, 9, 9]

particle 1 2 3 4 5 6 7

• sample ‘neighbor’ routine

do I=1, N-1 (N = # of atoms)

do j=i+1, N

compute r² (with periodic boundary conditions)

if (r² < r_n²) then

add "j" to list of atom "i"

endif

enddo

enddo

(1) Main routine :

• call neighbor

(2) Integrate routine :

update v (1/2 a timestep)

update x (full step)

• if (Mth step) call neighbor

call force

(3) Force and energy routines

(old)

do i=1,natoms-1

do j=i+1, natoms

compute distance between atoms i,j

(new)

do i=1,natoms-1

do k=nlistpnt(i), nlistpnt(i+1)-1

j=nlist(k)

compute distance between atoms i,j

Things to think about

What is "correct" reneighboring criterion ?

If you run with and without neighbor lists, will you get the same answer ?

Cautions:

consider the case where the number of neighbors is zero !!

Time = CPU time per step in seconds, (Thompson, 1983)

Neighbor Lists in MC

need to know all information for all the particles because we operate on particles one at a time

• Changes to program:

Linked Cells

• hash atoms into bins every timestep

regular grid of bins ibin_x = Nbins_x*x(i)/Lx

• search only surrounding bins to find atoms within r_c

• data structure : 2D array of bins, each bin containing list of atoms (ibin, atoms_bin)

• Things to think about

• What does the bin size need to be ?

• How many bins to be searched in 2D ? 3D?

• What is the scaling ?

More or less efficient than neighbor lists ?

Linked Cells and neighbor lists

• hash atoms into bins every timestep

regular grid of bins ibin_x = Nbins_x*x(i)/Lx

• use surrounding bins to set up a neighbor list

• use neighbor list to find atoms within r_c

Statistical Errors

• Two problems
• averages over finite time
• successive configurations may be highly correlated

• Want to find the RMS deviation in some average property <A>

For a simulation of infinite length ...

For a simulation of finite length, t_B with correlated samples characterized by a correlation length of

Need to establish ...in the limit t_B gets large

Strategy : Break a simulation up into blocks and take block averages

For a block of length t_B, where we can calculate

Plot vs t_B, and find the long time limit

• uncorrelated configurations occur ever 2650 steps !!

Lab 2 : Computational Efficiency and Error Estimation.

A new MC code will be provided which has neighbor lists already implmemented. This code also has the basic coding for block averages for estimating errors in the potential energy.

• to implement neighbor lists in your molecular dynamics program.
• to study the scaling behavior of the codes with and without neighbor lists
• to implement block averages and error estimation in the MC code for the pressure
• to study the dependence of the errors on the state point of interest and the number of MC moves.

1. First let’s see what we can expect from the neighbor lists. In the new MC code provided, set up the following parameters:

2. Implement neighbor lists in the molecular dynamics program. This will entail 1) writing the subroutine ‘neighbor’ discussed in class, and modifying the loops in the force, energy, and pressure routines. All the other necessary code will be provided.

3. The MC code provided calculates the statistical inefficiency as a function of the length of block averages in the program.

Test the sensitivity of the standard deviation of the energy (i.e. the error in the average), and the statistical inefficiency (the number of steps between uncorrelated moves) to the number of MC moves and the state point.

start with

and vary the density, temperature, number of moves, and equilibration

1. implement linked cells and/or the hybrid linked cell/neighbor list method.

2. implement error calculation in the MD code.

3. implement error calculation for pressure.

4. investigate ‘equilibrium’ with both MC and MD codes.