Calculations on molecular, biological, and physical systems come with a common goal: determining the lowest, most stable energetic state(s), and thus the best example of the system to study computationally. To do this, the values of interest, usually energy, need to converge: i.e., after repeated changes, each result is nearly the same as the one before. Sometimes, this is easy; sometimes, this is one form of a frozen hellscape, a nightmare loop of never-ending calculations always jumping away from a final answer. This is stress-inducing and infuriating. Dealing with this appropriately is a must for assessing the pieces involved, and finally arriving at the “best” answer.
My discussion below skims the surface of the process involved for those unfamiliar with the task, but perhaps familiar with concepts of energy and PESs, discusses the “why” of this traversal difficulty, and then maneuvers into some approaches taken to help a system converge.
For those individuals who work with molecular modeling, or any type of computational chemistry/biology/physics simulation, one key word comes up all the time: convergence. In order to even begin asserting that one’s calculations are “real” or “applicable,” it must converge.
Great… so, what does that mean?
Bluntly, this means that the object under study (a molecule, material, or biological substance) has arrived at a place of stability and minimal change under the desired conditions, and is restrained by the assumptions thrown on top of it during the calculation. If this is a bit befuddling, check out Notes 1 & 2 below for leads to more information and analogies that may be useful for following along.
Convergence can be REALLY HARD to achieve (hence the threat of computer defenestration, of which we would like to avoid). Thus, it may be useful to discuss some of the process, to get a better appreciation for the difficulties involved…
A helpful visualization for anyone unfamiliar with this field of study is a potential energy surface (PES). If you are unfamiliar with this idea or want a refresher on potential vs. kinetic energy, hop on over to my other article that discusses these in greater depth (with analogies and examples), as I will be using PES here to discuss optimization procedures (i.e. finding the “best” structure of a molecule).
When convergence is being discussed, it is a search for the structure (or state) residing in the deepest well, the global minima of the PES based on, for example, the molecular coordinates of the system (bond lengths and angles). However, PESs can be enormous and difficult to search rigorously, so finding the deepest local minima (lowest point in a region of PES) in a searchable PES will suffice.
Question: Why are we searching for the lowest potential energy on a PES?
Answer: During structure optimization, when you have a particular drawing or model of your system that you would like to use, finding the lowest potential energy means finding the structure at which there is minimal stress AND, coincidentally, is the one most often found naturally. In other words, a fully “relaxed” structure is (usually) more physically relevant than any other state of the system.
If you read through my analogy used on the PES blog post, this is equivalent to finding Timmy safely on the ground, perhaps even lain out on the grassy field, versus standing on the roof preparing for an intense cannonball swimming pool bonanza.
Other locations on the PES are just possible structures that the system could have been in (transition states, or spots Timmy could also be, like on the pool deck), or could get to under certain conditions (e.g. by increasing temperature, wherein Timmy may forego the jumping and just plop right on in).
Generally, finding the deepest well uses some variation of an algorithm known as “steepest descent” (aka Gradient Descent). Given a starting geometry (a position on the PES), the algorithm looks around that point and gauges the steepness of the curve, or how quickly any particular direction descends downhill. After identifying the fastest way down from that point, it moves a certain distance towards that direction (known as a “step size”) and checks all directions again, repeating this process until (hopefully) finding the minimum, which is identified by all directions being an uphill climb.
Traveling along the PES is well visualized in the graphic below, and further information on the algorithm can be found here. Warning: a fair bit of mathematics knowledge would be useful reading through that article – however, figures and images are useful guides, and well used within it.
All of this sounds reasonable… until your search gallivants across the PES unbounded, or until you reach a well that is not the minima desired! These situations are known as divergence and “getting stuck,” colloquially and respectively.
How do these situations happen?
Well, divergence occurs when the PES is very flat, or lacks steep curves to easily travel down. Due to the steps that Gradient descent algorithms use, on a flat surface the steps are more likely to either leap over any form of minima, even at smaller step sizes (which then begin to take up much more time to search). Searching this surface is completely analagous to walking along a meadow and looking for a minor depression in the ground that is less-minor then the other fluctuations in the soil… not an easy, nor pleasant, task…
Becoming “stuck” in a well is fairly easy, and depends highly on where you begin searching the PES. Re-examining Fig. 1, starting on an adjacent hillside position results in different paths down, and thus leads the system into a different minimum state.
Hmm, so if those are possible culprits of convergence troubles, how do we spot them?
Divergence is easy: if you were to list the energies being output during an optimization, it would be apparent that the energies were not decreasing to a lower value overall, and may also be jumping up and down. Most likely, there would be wild fluctuations or an overall increase in the energy being tracked. For computationally savvy users (or those on computer clusters, with access to software output files), one can use grep to search for an list energy values as the calculation proceeds in many software: Gaussian, Quantum Espresso, CASINO, MolPro, etc…
When the system has become “stuck,” on the other hand… this involves a little extra effort to identify, as it looks (and technically has) converged. It is common practice to do a sanity check of sorts on calculations, especially those that do not have previously studied experimental or high accuracy energy values to compare to. To do said sanity check, the calculation is run multiple times, with different starting orientations, and the final results (structurally and energetically) are compared, with any “stuck” positions on the PES beginning to stand out if any were accidentally located.
Wait a second – why all this searching? Why travel somewhat blindly along a surface that we have spoken about regularly, and not just jump right into the pit?
A much longer process would involve trying to map out the entire PES of the system under study. Imagine trying to gauge the height of a tree at maturation – what variables are needed? Nitrates, phosphates, mineral levels (soil composition) at the very least, not to mention water availability (rainfall, floods, watering by hand, sunny days evaporating moisture…), sunlight, and air quality. Even applying an assumption that these factors are somehow independent (sounds crazy, right? Scientists due this often to get a starting point to work from though!), try graphing all those factors into one “surface” that described every direction of possible growth patterns the tree could live up to.
Got the picture?
If so, become a scientist and do computer modeling, because that visualization is HARD!
Albeit my example is extreme, and many small molecule systems can be understood by focusing on certain structural parameters (bond lengths and angles) with some foundational physics behind them (bunch of springs and weights stuck together is something we are quite confident working with), the sheer amount of different combinations is nearly insurmountable.
Tricks and clever constraints are used to focus in on the physically relevant combinations (e.g. atoms need to be fairly close to bind, for instance, or won’t sit inside each other), but this attempt is no easy task. Many researchers focus on one particular variable, such as one bond changing or a term known as “reaction coordinate” (trying to look at the process by which something changes, like how far along your mechanic is when changing your car’s tires). However, those more adventurous may attempt to study larger swaths of a surface.
Alternatively, some researchers have even focused on how molecules roam about these surfaces as they react, leading to a fresh understanding of molecular dynamics! Check out Further Reading #1 for more information, written and researched by Drs. Vale Cofer-Shabica and Richard M. Stratt, friend and mentor of mine, respectively.
I digress, however. Now that it is clearer how difficult searching a system’s PES actually is, the expediency at arriving in relaxed states, in the deepest wells, is quite astounding. Algorithmic techniques have improved the ability to find these wells in less time and with more reliability, along with further improvements to theoretical principles enhancing the accuracy of the surfaces (or parts explored) themselves.
But how do we, as researchers, actually help our calculation achieve convergence with these continuously improving techniques?
Simple: freely fidget with your starting conditions, and work with the community.
Now, “simple” is probably going to ignite a few torches, soon to be paired with the pitchforks, but I have found in my studies that constant tweaking is the first step to letting the software in use truly shine. When I still found myself unable to achieve steady convergence, or was finding strange results that didn’t match previous work, I went to the community: my lab mates, my mentor, and online resources such as ResearchGate, CCL.net, or Quantum Espresso’s forums.
Manuals are of course of use in many cases starting out, but as one branches away from the cookbook style calculations, the manuals get a bit… murkier in terms of direct applicability.
I believe the key to keeping your computer (and your sanity) safe is this: share the frustrations with those who understand the same discouragement of blindly searching, and be free to explore and change the way you approach the calculation.
With over 25,000 results on only the American Chemical Society’s journal database when searching for “potential energy surface,” nearly 1,000 of which published in the last year, researchers continue to explore ways of probing the variety of changes systems can undergo.
For the curious, and those interested in other ways of traveling a space of vast options in hopes of finding the ideal state(s), Dr. Brenda Rubenstein’s (my mentor’s) research foci of Quantum Monte Carlo might be of interest (alongside pivotal leaders of the field such as Drs. David Ceperley, Berni Alder, etc… see Further Reading #2). Instead of a path of changes to follow, it uses a random sampling of all possible states a system could be in, and a few clever selection conditions to both expedite the search and even improve accuracy of results.
- Cofer-Shabica, D. V., & Stratt, R. M. (2017). What is special about how roaming chemical reactions traverse their potential surfaces? Differences in geodesic paths between roaming and non-roaming events. The Journal of Chemical Physics, 146(21), 214303. http://doi.org/10.1063/1.4984617
- Ceperley, D., & Alder, B. (1986). Quantum Monte Carlo. Science, 231(4738), 555–560. https://doi.org/10.1126/science.231.4738.555
- Rubenstein B. (2017) Introduction to the Variational Monte Carlo Method in Quantum Chemistry and Physics. In: Wu J. (eds) Variational Methods in Molecular Modeling. Molecular Modeling and Simulation (Applications and Perspectives). Springer, Singapore