SciComp Project 3 Solved

SciComp Project 3

1        Background
1.1       Importance of Rapid Convergence in Root Finding and Optimization
Finding a minimum energy configuration for a system of mutiple atoms is a widely applicaple for designing molecular systems in solid state physics, chemistry, medicine and other fields. Many applications will require tens to hundreds of atoms in the model. Calculating the energy involves solving the many-body Schrødinger equation, in which all the electrons interact with each other.

With hundreds of atoms and thousands of electrons, solving the many-body Schr¨odinger equation (even approximately) can take weeks of computetime to calculate a single energy point, whereby a full optimization can easily take many years of CPU-time to run. Even on a parallel computer with, say, 1000 CPU cores, optimizing the geometry of a single molecule can easily take months in “human time”.

When every evaluation of the energy function takes, for example, a CPU-week to compute (burning large amounts of CO2), the number of function evaluations matter: a problem that may be extra bad in quantum chemistry, but is common to many fields.

In this assignment, you will try your hand at optimizing a simple system of interacting atoms. However, we do not have months to sit and wait for computations to complete, so we will use an extremely simple approximation to the potential energy, the Lennard-Jones-potential. This approximation is appropriate for noble gas atoms, yielding a crude picture of the formation; yet we will pretend that it is a full ab initio quantum chemical energy calculation, and assume that it takes a week to complete. Thus, you will be asked to count energy function evaluations and report the time spent.

1.2       Minimizing the potential in a cloud of Argon atoms using the LennardJones potential
The Lennard-Jones potential works for systems of neutral atoms or molecules. It assumes that no new electronic bonds nor new molecules are formed, and that the potential can be modeled only with a short-range repulsion force (arising from the Pauli exclusion principle) and a long-range van der Waals attraction. This yields a classical description of the system, in which the neutral atoms or molecules move as classical particles, interacting only through this simple potential.

The Lennard-Jones potential can be written in several ways, but the most common is:

 !

(1)

vLJ(rij) is the potential at distance rij = kxi − xjk2 between two atoms i and j. The repulsive force is stronger than the van der Waals-force, but decreases more rapidly with distance. Thus, for a pair of atoms, there is a distance where the total potential is minimal, the potential well.  is this minimal potential between two atoms, and σ is the inter-particle distance where the potential is zero: vLJ(σ) = 0. The constants  and σ are generally found by fitting to proper ab initio calculations, or using experimental data.

The full potential LJ-energy of N neutral particles is the sum of pair-potentials:

                                                                    N        N

VLJ(X) = X X vLJ(rij)

i=1 j=i+1

where X ∈ RN×3 is the matrix of coordinates for the N particles.

In this assignment, we will be modeling argon, for which
(2)
                                               σ = 3.401˚A and   997kJ/mol                                                                 (3)

1.3       Array Programming
Try to use array programming when programming with Python/Numpy or Matlab, as the speed difference can be up to a factor of hundreds. As an example of how to use array programming, here is how a distance matrix (to be used in the LJ-potential) can be calculated with array operations using NA = newaxis (as I explain further here):

                                                     Dij = kxi − xjk2 = kxij − xijk2                                                                         (4)

In Eq. (4), xi is extended with a newaxis j and xj with a newaxis i to put them on the same footing, so that the pair differences can be calculated as elementwise subtraction xij − xij. In Numpy code, it looks as follows:



# points:                                        (N,3)−array of (x,y,z) coordinates for N points

# distance(points): returns (N,N)−array of inter−point distances:

def distance( points ):

displacement = points[:,NA] − points[NA,:] return sqrt( sum(displacement∗displacement, axis=−1) )



2        Questions for Week 4: Solving Nonlinear Equations
2.1       Solving Nonlinear Equations in 1D
a.   Write a function that computes the Lennard-Jones energy for a collection of particles described by an (N,3)-array points. The resulting function should not take more arguments than this, as it will be called with by your root-finding functions below. In Python, you can do this using a programming construct called closures as follows:



def LJ(sigma,epsilon): def V( points ):

...implementation that depends on sigma and epsilon... return V



Then ) will return a function specialized to σ and  such that V (points) returns the total Lennard-Jones energy for the system.

Demonstrate your solution by making 1) a plot of the potential between two Ar atoms, one placed at x0 = (x,0,0) the other at x1 = (0,0,0) with x ranging from 3 to 11; and 2) a plot of the LJ-potential in the same range with two extra points added: x2 = (14,0,0) and x3 = (7,3.2,0).

b.   Write a bisection root finding function x, ncalls = bisectionroot(f,a,b,tolerance=1e−13) that finds x such that f(x) = 0 given a bracket x ∈ [a;b], and counts the number of calls to the function f. (A reasonable length is 8-12 lines of code).

In this assignment, let the convergence test be on how close we get f(x) to zero. Test it to find the zero of the the LJ-potential between two argon atoms as a function of interatomic distance, and verify that you get x = σ. How many calls to the energy function were needed to get from the start bracket [a,b] = [2,6] to |f(x)| < 10−13?

c.    The derivative of the pair-potential ) is

 

Write a Newton-Rhapson solver x, ncalls = newtonroot(f,df,x0,tolerance,maxiterations), (a reasonable length is 4-8 lines of code), and test it in the same way as above. For simplicity, assume that a call to the derivative df has the same cost as a call to f. How many calls were needed to get from x0 = 2 to |f(x∗)| < 10−12, i.e., 12 decimal digits for x∗ after the comma?

d.   Make a combination of Newton-Rhapson and bisection that is guaranteed to converge, but takes advantage of the quadratic convergence of Newton-Rhapson iteration, and test it on the same example. How many calls to the LJ-energy function was needed to get from x0 = 2, [a,b] = [2,6] to obtain |f(x∗)| < 10−13?

Note: If you have trouble completing this step, simply skip it and move on to the remaining questions, which only requires your bisection root solver to work. You can always return to solve it once you have completed tasks (e) and (f).

2.2       Solving N-dimensional Nonlinear Equations
Using the chain rule for derivatives f(g(x))0 = f0(g(x))g0(x), we can write down the derivative of a pair-potential with respect to the position of one of the particles:1

                                                                                               (5)

This we can then use to find an expression for the gradient of the total LJ-energy:

                                                                                    (6)

The position of all the particles is a point X ∈ R3N, which we organize as a (N,3)-array, so that xi = (xi,yi,zi). The total potential energy is a function VLJ : R3N → R, and the gradient taken at any particular configuration of particle positions is hence a 3N-dimensional vector: VLJ(X) ∈ R3N. Thus, the gradient to VLJ is a function ∇VLJ : R3N → R3N. The negative of the gradient is the force acting on the system, and its direction is that in which the potential decreases most rapidly. But notice that it is a 3N-dimensional direction: it acts on all the particles at once.

We can write it down in Python as follows (in Matlab, it would look very similar):



def LJgradient(sigma, epsilon): def gradV(X): d = X[:,NA] − X[NA,:] # (N,N,3) displacement vectors r = sqrt( sum(d∗d,axis=−1) ) # (N,N) distances fill diagonal(r,1) # Don’t divide by zero

T = 6∗(sigma∗∗6) ∗ (r∗∗−7) − 12∗(sigma∗∗12) ∗ (r∗∗−13) # (N,N)−matrix of r−derivatives

# Using the chain rule, we turn the (N,N)−matrix of r−derivatives into # the (N,3)−array of derivatives to Cartesian coordinate: the gradient.

# (Automatically sets diagonal to (0,0,0) = X[i]−X[i]) u = d/r[:,:,NA]           # u is (N,N,3)−array of unit vectors in direction of X[i] − X[j] return 4∗epsilon∗sum(T[:,:,NA] ∗ u, axis=1)

return gradV



1NB: You don’t need to understand these derivations to solve the problem, as the actual code for the gradient is provided to you. The derivations are here to show you how to do it yourself in the future.

Finding the minima in R3N of the potential involves finding points where its 3N-dimensional gradient is 0. An important component needed in next week’s work is called line search, where we search for a point along a single direction d ∈ R3N at which the gradient along this line is zero. That is, we want to start in a point X0 ∈ R3N, and then find α ∈ [0;b] such that 0 = d · ∇VLJ(X0 + αd), i.e., the gradient has no component in the direction of d. This finds an optimum for the one-dimensional function VLJ(X0 + αd).

e.   Look at the gradient of the 2-particle system in question (a) with x1 = (0,0,0) and x0 = (x,0,0), x ∈ [3;10]. Why are exactly two components nonzero? Why are they equal and opposite? Plot the nonzero component for the derivative of the x-coordinate of x0 (0,0coordinate of gradient) together with the potential, and notice the relationship between the zero of the derivative and the minimum of the potential.

Next look at the gradient for the 4-particle system from (a) at one of the minima of your plot. Why is the gradient not zero?

f.    Write a function x0, ncalls = linesearch(F,X0, d, alphamax, tolerance, maxiterations) that takes a function F : RN×3 → RN×3, a start position X0 ∈ RN×3 and finds the zero along the line-segment X0 + αd of d · F(X0 + αd).[1] Use your bisection solver at this point, as using Newton-Rhapson requires second derivatives when F is the gradient.

Test your function by finding the mimimum along X0 + αd with

X 

and d = −∇VLJ(X0), and with α ∈ [0;1].

[1] You can use either a closure or a lambda expression to define the one-dimensional restriction of F to the line.