Computational Lab 5

If you were able to finish this lab before the start of final projects, you'll get into the hall of fame.

This lab contains 3 parts. First part is high dimensional optimization with Newton's method. The second part is machine learning with gradient descend (non-linear fitting). The third part is multibody dynamics simulation.

Part 1: High dimensional optimization with Newton's method

Consider the system where we have 3 atoms.
Atom2 and atom3 are identical and atom1 has a shorter equilibrium bond length.

V_LJ(r₁₂) = 4 ε₁ [ (σ ₁/ r₁₂)¹² - (σ ₁/ r₁₂)⁶ ]

V_LJ(r₂₃) = 4 ε₂ [ (σ ₂/ r₂₃)¹² - (σ ₂/ r₂₃)⁶ ]

ε₁=5

ε₂=3

σ₁=1

σ₂=1.732

The total energy of the system is V = V_LJ(r₁₂) + V_LJ(r₁₃) + V_LJ(r₂₃)
The starting position for all three atoms are:
x₁ = 10.0 , y₁ = 10.0 , z₁ = 10.0 x₂ = 8.0 , y₂ = 12.0 , z₂ = 10.0 x₃ = 12.0 , y₃ = 12.0 , z₃ = 10.0

Hint: the above set of coordinates will not converge to the LJ potential well. You've encountered this problem in lab 3. If you don't know how to make your NM optimizer converge to the minimum in stead, use the following set of coordinates to start with:
x₁ = 10.0 , y₁ = 10.0 , z₁ = 10.0 x₂ = 9.4 , y₂ = 11.0 , z₂ = 10.0 x₃ = 10.6 , y₃ = 11.0 , z₃ = 10.0

1. Create a python script for this part of the assignment.
2. Write a function that will return the total energy of the system.
3. Write a function that will return the gradient of total energy of the system, given appropriate coordinates. If you treat your variables as being independent of each other when they are not, you will produce wrong optimized structure, energy and force but you can still receive full credit for a functional optimizer.
   A gradient vector in python looks like this:
   
   Extra challenge : Produce the correct gradient using appropriate coordinates/variables. There are two ways to approach:
   First is redefining the variables of the systems using given initial cartesian coordinates. The three atoms defines a triangle. You can use lengths of two of the sides and the angle they form as your variables, so that the length of the third side can be expressed as a function of the length of the other two sides and their angle (Law of cosines). With these variables you can calculate the gradient with respect two these variables correctly. CAUTION: Your gradient with respect to these variables is not directly the force you are familiar with in cartesian coordinates. As you've seen in Lab 3, Newton's method can be unstable or fail to converge to a local minimum. When you use these variables, you are warping the space (you are in a different coordinate system that is neither cartesian nor spherical) such that the curvature of energy surface in this coordinate system can be very far from harmonic (Newton's method can be very unstable if the potential energy surface is far from being harmonic/quadratic). You may need to improvise to make sure that your optimizer converges.
   
   Second, you can express V_LJ using cartesian coordinates of the atoms involved, instead of separations. That is, for example, your r₁₂ is converted to SQRT( (x₁ - x₂)² + ...). In this case, your variables will be all the cartesian coordinates of all atoms. This approach is higher dimensional compared to the previous one, but the energy surface in this coordinate system can be more harmonic, therefore may increase the stability of your optimizer.
   
   
4. Write a function that will return the Hessian matrix given appropriate coordinates.
   The hessian matrix in python looks like this:
   
5. Write a function that will return the magnitude of the given gradient.
6. Write a function that will perform a Newton's method step.
7. Write a function that will keep taking steps until the magnitude of gradient is smaller than 0.01. If it cannot converge, scale your Newton's method stepsize by 1/2. Keep repeating this until it converges. Alternatively, put a cap on how large any single step could be so that the optimizer is stable.
8. EXTRA CREDIT: Use a software (Mathematica or Matlab) to analytically solve for the position that you should converge to.
9. Have your program print out the number of steps it took to converge, your scaling factor ( if you did not have to set this, it's 1), the gradient at the last step, the energy of the minima, the relavant coordinates of the 3 atoms. Creat a lab4.php, link it to your index.php and post your answers there.
10. Examine the positions the 3 atoms are in. What are the ideal bond length and how do they compare to their corresponding sigma? Are they the same? Why should/shoudn't they be the same?

Part 2: Simple machine learning: fitting with Lennard Jones potential.

In this practice we have a collection of images and their corresponding energies.
The system contains 5 atoms, with the first one having stronger interactions with the rest, while all other atoms are identical, similar to the case you studied above.
This means that there are two sets of ε and σ. One for atom 1 and one for the rest.
Our goal is to find out ε₁,ε₂,σ₁,σ₂
To start, copy the following scripts to your lab4 directory.
/home/fri/lab_files_2021/lab4/dataGen.py
/home/fri/lab_files_2021/lab4/myFirstML.py
You will work on myFirstML.py. This script provides you a set of data, generated by dataGen.py.
The data format is :
image1 [ [ [x₁₁,y₁₁,z₁₁],[x₁₂,y₁₂,z₁₂],[x₁₃,y₁₃,z₁₃],[x₁₄,y₁₄,z₁₄],[x₁₅,y₁₅,z₁₅] ],energy₁ ]
image2 [ [ [x₂₁,y₂₁,z₂₁],[x₂₂,y₂₂,z₂₂],[x₂₃,y₂₃,z₂₃],[x₂₄,y₂₄,z₂₄],[x₂₅,y₂₅,z₂₅] ],energy₂ ]
etc...
This part is worth 6 points total. You can do it however you like, or use the following example method with gradient descend:

1. Make a list of 4 variables: epsilon1, sigma1, epsilon2, sigma2. Set their initial values to 1 1 1 1. This list is your initial guess. You can guess any numbe r, but if your initial guess is too different from the true values, the optimization may take a long time.
2. Write a function that takes in two lists of coordinates, and returns a distance.
3. Write a function that reads in the data variable one element at a time. This funciton then calculates the total energy of the system using your specifiec parameters.
4. Write a function that calculates the energy differences squared at every point and sums it up. This is your error.
5. Write a function that varies the 4 variables: epsilon1, sigma1, epsilon2, sigma2. These 4 are your independent variables to your function that produces the error. In other words, treat your error as a function of those 4 variables. Vary each of them by a small initial amount (such as +- 0.01). Calculate the first order derivative of error with respect to each variables independently. This is the function that produces the gradient of your error function.
6. Write a gradient descend optimizer for the error function.
7. Output the parameters your optimizer find.

Bonus challenge : Examine the dataGen.py, and print out the correct answer your optimizer should obtain. You can acquire these points independently from the tasks above.