Computational Chemistry 3.2 - Energy Minimization

Short lecture on energy minimization of molecules.
A set of atomic coordinates where the potential energy gradient is equal to zero is called a stationary point. All minimum energy structures are stationary points and local minima. The lowest energy local minimum is the global minimum, i.e. the lowest energy configuration.
Notes Slide: i.imgur.com/Q06ZW88.png
-- About TMP Chem --
All TMP Chem content is free for everyone, everywhere, and created independently by Trent Parker.
Email: [e-mail]
-- Video Links --
Course Playlist: • Computational Chemistry
Other Courses: • PChem Course Intros
Channel Info: • About TMP Chem
-- Social Links --
Facebook: / tmpchem
Twitter: / tmpchem
LinkedIn: / tmpchem
Imgur: tmpchem.imgur.com
GitHub: www.github.com/tmpchem
-- Equipment --
Microphone: Blue Yeti USB Microphone
Drawing Tablet: Wacom Intuos Pen and Touch Small
Drawing Program: Autodesk Sketchbook Express
Screen Capture: Corel Visual Studio Pro X8

[@alichem8]

Thanks for your clear explanation, I really enjoy watching your videos 🙌

[@Learning432]

Dear Prof., hi. You are working nicely and enthusiastically. Make a video explaining how to choose the appropriate hybrid functional & basis set for specific transition metals for inorganic clusters like [XM12O40]n-, Polyoxometalates (POMs) are clusters of units of oxoanions of transition metals, such as Mo, W, V and Nb.

[@sagarkamble2947]

Nice explanation.. I would like to know why we need to do an energy minimization of the system?

[@tusharkush7]

It is widely used in molecular modelling. Most biological molecules like proteins function in their native state which is a local or mostly global energy minima. Or it could be applied to a non scientific statistical function.

[@yhk6187]

Even if I downloaded all files and run it as you showed, still error messages coming out..C:\Users\scripts\molecular_mechanics\mmlib\optimize.py:311: RuntimeWarning: invalid value encountered in double_scalars self.mol.atoms[i].coords[j] += disp_mag * disp_vector[i][j]why does this occur?

[@lyamlikeji]

I also meet the same problem. Did you solve this problem?

[@sagarhm2237]

How to energy minimize pdb in python

[@ThiagoBenazziMaia]

Hello, I can't stop looking at the code and playing with it. One thing that I realize for 3 atoms on the opt.py. If we set the h2o.xyzq to xyz position O(0,0,0) H1(0,1,0) H2(0,-1,0) we got the minimization locked on the same position. Should we add a very small random number on the function to avoid this kind of lock?

[@TMPChem]

What you have inadvertently discovered is the concept of a "saddle point" or "transition state" in a potential energy surface. This is a point in multidimensional space which is a minimum with respect to all but N dimensions, and a maximum in N dimensions. Such a point is an Nth order saddle point, with 1st order saddle points such as this one being by far the most common. In this case the saddle point arises due to symmetry constraints. Under conditions of perfect D_infH symmetry, any deterministic energy minimization algorithm is incapable of breaking that symmetry (except accidentally due to numerical precision limits / errors). However, as you noted, as soon as the symmetry is broken even gradient descent will quickly push away from this structure and down towards the bent C_2V minimum energy structure. This example highlights the importance of the selection of the initial structure for geometry optimization. Often, care should be taken such that any symmetry elements in the structure should be broken through small, random perturbations in the coordinates. If the high symmetry structure is a true minimum, the structure will find its way back too it. If not, it will quickly displace away from it.

[@ThiagoBenazziMaia]

Thank you so much for the complete explanation.

[@johnleed876]

To get more accurate energy minimization values, would apporoximations like Hartree-Fock or DFT be used?

[@TMPChem]

Hi John. Perhaps. Just because a method is more computationally intensive doesn't necessarily mean it will lead to superior outcomes for a particular modeling problem, but if judiciously chosen, more often than not, that will likely be the general direction things tend to move. Sometimes mathematically simple models like molecular mechanics can be highly optimized to particular subdomains, and within that niche may produce surprisingly good results, perhaps even superior to that achieve with HF or DFT. Some methods, like Full CI with a large basis set, will be extremely accurate for almost all chemical scenarios, but will be intractible for all but extremely small systems (such as a single water molecule, or two helium atoms). For most systems in a closed shell ground state (and all atomic numbers being below about Iron), some reasonable post-HF method like MP2 with a triple-zeta basis set will generally do very well.

[@gauravjoshi8194]

I think second derivative zero is not a saddle point, it's inflection point. The saddle point is when the second derivative is negative (depending upon no. of coordinate with the negative second derivative, it is termed as first, second-order saddle point). Correct me, if I am wrong.

[@TMPChem]

Hi Guarav. You're correct that a saddle point is referred to by the number of second derivatives which are negative at a particular location. As for when a second derivative is zero, that situation is still ambiguous. We would need to examine even higher-order derivatives to make a conclusive claim about the type of point that exists whenever a particular second derivative is zero (in addition to all of the first derivatives). In the function f(x) = x^3, x = 0 (where d^2 / dx^2 = 0) is an inflection point, whereas in f(x) = x^4 it is a local minimum.

[@gauravjoshi8194]

Thanks... So, It can be minima as well as maxima.

[@TMPChem]

Yes, a stationary point with a zero second derivative can be a local minimum, a local maximum, or neither.

[@Vincent8234]

You are so good. what the hell!

[@TMPChem]

I know right? It's not fair. ;)

[@bowang475]

@@TMPChem Yea, much better than "distinguished" professors!

[@mrknowmyself]

God level of chemistry lecturer, huge respect!

[@TMPChem]

Thanks Jeju

[@ThiagoBenazziMaia]

Hello, it's me again -). Thank you for making a very clear code, full of comments. It make it so easy to follow. I was trying to understand how the gradient works and I was going through the "gradient.py" on the function "get_g_bonds(mol)". I couldn't find if the mass was used to calculate the vector. I made a test using the h2o.xyzq. I left there just two atoms O (0,0,0) and H (0,1,0). When I run python opt.py ../../input/opt/h2o.opt I got the oxygen moving 0.00216943 and Hydrogen moving 0.00216942. Almost the same distance. Shouldn't the final vector be proportional to the mass between the pair of atoms? Thanks a lot

[@TMPChem]

Welcome back. Minor nuance here. In molecular dynamics, we obey Newton's second law, F=ma. The force comes from the negative gradient of the potential energy surface, and dividing by mass gives you acceleration, which is the derivative of velocity. In this case, from a zero initial velocity with equal and opposite forces (Newton's third law), the displacement of the oxygen should be ~1/16 times that of hydrogen. This leads to all sorts of invariants being obeyed, such as conservation of total energy, constant velocity of the center of mass, etc. Optimization is different. It does not obey laws of motion. It only tries to find the nearest local minimum on the potential energy surface. When displacing atoms in between steps, you won't find a place where displacements are inverse-mass-weighted (as in MD) because I didn't put one in there. Thus the center of mass is allowed to displace during the optimization. You *could* do optimizations in a way that maintains the center of mass position, but it is not done here. Doing so could be beneficial, but it is not necessary for the purpose of finding local minima. As noted earlier, the set of 3N cartesian coordinates are not unique, as there are 3 translations and 3 rotations (for non-linear systems) which are energy invariant. The resulting structure after optimization is one of these redundant structures, and could be translated and/or rotated to the center-of-mass invariant structure if desired after the fact. The method chosen here does not affect our ability to correctly optimize to one of these equivalent structures.

[@ThiagoBenazziMaia]

Thank you for replying. I see now, I should have been using the md.py. So from what I understand so far from your explanation, the videos and the code. If I want to simulate a protein folding I should. 1) Add the aminoacid chain 2) Add water 3) Run md.py, and once is folded 4) Run opt.py

[@TMPChem]

Sort of. I would also highly recommend against attempting to use this program for large-scale or research-quality production runs. It is meant primarily as a pedagogical exercise to assist in understanding what's going on under the hood of more sophisticated and thoroughly developed / tested / optimized packages. A typical workflow in such a program for protein folding simulations would look like this: 1) Generate an initial structure for the amino acid chain (PDB, CCDC, internal script, hand-drawn, etc.) 2) Surround the protein with solvent water molecules in a box (also automatic tools for this) 3) Optimize the structure a little bit, but not all the way 4) Run MD for as long as practical or desired 5) Analyze results for properties of interest Step 3 is a bit tricky to explain, but basically you don't want too high energy or too low energy of a structure to begin with. Your system has thermal energy, and you want your starting coordinates and momenta to be as reflective of a typical timestep as possible (otherwise you waste a lot of timesteps until it has reached such a state, and your earliest steps are all junk). To do this, researchers have developed elaborate sequences of optimization, annealing, dynamics, etc. to begin simulations such that the first "official" timestep has as much statistical significance in the results as any other. But that gets way into the details and quirks of statistical mechanics and molecular biology far more than I'd like to in this playlist.

[@ThiagoBenazziMaia]

Hello, thank you again for the replying. I'm learning in two weeks on your youtube channel, more than I learned for more than one year using another ways to study and I think I'm not the only one who thinks that. Please, don't stop this project =) . Thanks a lot.

[@TMPChem]

I don't expect to stop, we're just hitting a two-fold issue on the timeline: 1) I have an orthogonal, demanding full-time job, and 2) The material is getting much more demanding to produce. The comp chem playlist now starts going into electronic structure theory, which to do as well as the MM / MD sections requires an enormous time investment not only learning the material well enough to present, but having analogous helpful coding demos and exercises. The next thing to be released will be all of Hartree-Fock, which I estimate will be ~40 videos split over two chapters. Hoping for end of year 2017 release. That may be optimistic, but why have goals if we're not going to push ourselves?