Molecular Mechanics: Stretch, Bend & Torsion Terms

[00:00:00]this is my cat onyx sometimes she likes to come up in my home office and give me a hand in the afternoon as i'm doing my work i think right now she wants to take a nap so maybe we'll talk about something exciting it might wake her up like molecular mechanics okay so today let's talk about molecular mechanics these are also sometimes called force field methods and also sometimes this word is used synonymously with molecular dynamics although technically that's not exactly the same thing molecular mechanics has to do with what's the energy or what are the forces on the atoms and molecular dynamics is then based on those forces where would the atoms want to move and how do they move around in time so we will set aside discussion of dynamics for now and just talk about for a molecule in a certain geometry how would i know what the energy of that is and how would i know what the forces on the atoms are so that's what we'll talk about today these approaches use classical models like explain what i mean by that in a second to predict an energy of a molecule as a function of its conformation so based on the bond lengths and the bond angles and all of these kind of things if you are able to predict that then you can predict equilibrium geometries those are just the geometries that minimize the energy or you can predict the geometries of transition states with a few caveats sometimes these molecular mechanics approaches are not as good at deformed geometries like you see sometimes in transition states it could predict relative energies between conformers so cis versus trans conformers of a molecule you could predict the energy difference of that okay and so this is generally useful you can get the potential energy you would need from like the dynamics computations as i just mentioned um but i'll let you know at the outset most force field methods don't really work correctly if you literally break a covalent bond there are a few methods that can do this sometimes they're called reactive force fields and those exist but they're not as well known or well used as the standard for sealed methods so let's start off with the easiest part of a force field method the easiest thing to talk about first is the stretching interaction so if i've got a bond and i stretch the bond or i compress the bond how does that change the energy of a molecule what you could do is just assume some kind of taylor series and in fact that'll be the theme for a lot of molecular mechanics we'll suspect that the energy as a function of a bond length r is going to be some taylor series where i have some force constant k2 times the difference between the bond length and its equilibrium bond length or its preferred bond length r sub e and then i could tack on higher order terms like cubic terms or quartic terms etc we would expect if the taylor series works that the first term is dominant the second term is a small correction and additional terms would become smaller and smaller corrections so that in principle i shouldn't need those now how would i get these parameters like k1 or k2 and k3 and so on um well i i could do some experiments and infer that from experiment or i could do some high level quantum chemistry calculation and extract it out of that now the central idea of force field methods is that i shouldn't have to go do this for every single molecule otherwise that becomes an impossible game because there are too many molecules the idea instead is that for a ch bond that bond length is going to be pretty darn similar no matter what molecule we're in it's going to be something like 1.08 angstroms plus or minus a little bit and that's going to be pretty typical for any molecule so if i see a molecule i haven't seen before i could guess 1.08 angstroms and be in the right ballpark similarly this the stretching frequencies which is related to the force constant k2 for a ch bond are between 2 900 and 3300 wave numbers and so uh if i had to guess i could guess 3 000 wave numbers 3100 wave numbers have something and know that i'm in the ballpark so the fact that these values are similar across many different molecules means of those parameters what we call transferable and that's a good thing and it saves us from having to catalog every single situation in fact what we do is we look at the different classes of situations and then just fix parameters for those so you might imagine that a a ch bond you know maybe is a little different if you're in an aromatic molecule than in a non-aromatic molecule so you might have different parameters for those situations and maybe you would even have different situations if you're in an sp3 versus an sp2 or an sp hybridized carbon but if you distinguish into those categories then hopefully everything within that category can be lumped together with similar parameters that's the idea and it turns out this idea works pretty well as i hope i'm just motivated with this example of ch and the typical bond lengths and vibrational frequencies so in one of the earlier molecular mechanics force fields called mm2 molecular mechanics version number two there were a number of so-called atom types which basically just means what are the different chemical environments that an atom can be in that i need to worry about when i'm trying to catalog all these parameters like k2 so here's an example see there's a bunch of different carbons there's sp3 carbon there's an sp2 carbon when it's an alkene there's an sp2 carbon when it's a carbonyl or an amine there's an sp carbon etc etc there's quite a few examples and similarly nitrogens it depends on the hybridization um sometimes whether or not there's an oxygen nearby or whether you're in a nitro group or that kind of thing there's several different oxygens there's a handful of hydrogens you notice there's not a whole lot of different distinct hydrogen atom types so it matters if you're in a carboxyl group or alcohol or an amine but otherwise there's not a whole lot of different ones there atom type 5 hydrogen is except on in an o is kind of a catch-all hydrogen so they act pretty much the same in most molecules and then you see there's a few sulfurs notice atom type 11 there's only one uh fluoride uh i have type 12 there's only one chloride so some of the halogens don't particularly care what they're bonded to in a lot of molecules and that's kind of good and you notice there's a few noble gases they're listed in a couple different nickel types depending on the oxidation state whether it's oxidation state two or three a couple different cobalts etc some of your favorite atoms may not be listed on this list it might mean that the designers of this force field never got around to developing parameters for that atom and that actually can happen that's a drawback of forcefield situations so the force field in general as i've said maps a geometry of a molecule to an energy and here's how we do that we write the energy as a sum of stretching terms bending terms torsional terms i'll explain what that means in a minute van der waals terms i have to define that too electrostatic terms so that has to do with interactions between charges and then finally depending on which force field we're using sometimes we talk about cross terms and if we separate out the van der waals terms and the electrostatics terms then that tends to make the other terms more transferable so this was a key advance in the early days of force field methods was to pluck out those van der waals and electrostatic terms and make them separate and then that made the other terms more well-behaved and systematic so that was a good thing that allowed folks to move on from what used to be called spectroscopic force fields where you had a force field but it was specific to every single different molecule and that was not as easy and general as these more transferable force fields with equations like equation two for the energy just a few words about where this all came from way back in 1930 andrews proposed extending the spectroscopic force fields i just mentioned to doing a more general molecular mechanics in 1940 westheimer did the only micro mechanics calculation actually done by hand to find a transition state 1961 hendrickson did conformational analysis on rings with more than six atoms and this was a big deal in the 60s and 70s there was a lot of strong interest in conformational analysis and force field methods were handy to use for those applications weiburg and 65 published the first general molecular mechanics program before that everything had been some one-off application to this molecule or that molecule but weiburg had the first general ability to find the minimum energy of any molecule if you had the parameters and then allenger in 76 came along and published uh the first in his mm series of four seals which included mm1 mm2 and then ultimately mm3 and mm4 and then since that time this idea really caught on and there's been a lot of interest in development and force fields over the time since that all right so back to the stretch energy let me see if this is clear how we're going to do this for any pair of atoms that are bonded a and b the energy as a function of the distance between atoms a and b is some quadratic term k2 times the difference between the bond length and the natural bond length i called it re before here i'm calling it r0 a b and notice that term is that distance is squared so if i stretch the bond uh then that energy goes up quadratically or if i compress the bond then our a b is less than r a b zero i get a negative number but i still square it so i get a positive energy whether i stretch or compress and then i could have if i want cubic terms coordinate terms and so on a lot of times honestly in force field methods it's folks decide this is too complicated and they just take that first term that quadratic term and say it's good enough if you do that though you will get wrong behavior at long bond lengths when in reality you would start to break the bond and the energy should start to tail off to some asymptotic value that won't happen if i take a quadratic form the energy would just keep going up and up and up as i stretch which is correct if you want to do better one way to do that would be a function like this which still only has essentially one term but it's more accurate term this is the so-called morse potential so i have like a well depth or dissociation energy d and then i have one minus exponential and some factors there and i square the term in brackets this is more accurate in terms of what the energy would really be like however if i'm way out at long distances then if i computed the force associated with this energy i'd get a really small restoring force pushing me back towards equilibrium because the energy has kind of flattened out at long distances in a morse potential where i have a nice high repulsion at short distances and then a well and then i come up and i go towards dissociation limit the restoring forces that are near the association limit are very small and so if i were to try to optimize a molecule with these forces and i guessed a bond length that was too big that optimization might be slow because the force is so small so for that reason people say yeah just use the quadratic term because then that'll push it down really fast because i'll have a big big force at long distances and so this form although more accurate is is practically speaking not very popular in the community and i can illustrate this with the figure i'm borrowing from frank jensen's excellent textbooks introduction to computational chemistry which has a really nice chapter on monkey mechanics methods that i'm kind of following along here so what uh is done in this figure is this is the stretch energy for methane where i'm uh stretching a bond and i'm comparing the exact energy which is kind of the solid curve with a few approximations so p2 is the dashed curve the heavy dashed curve with a polynomial to order two or basically just that quadratic term so it's a parabola and you can see it looks like a parabola and if you're at very close to the equilibrium and you don't go very far in the delta r direction delta r negative meaning i compress the bond delta r positive meaning i stretch the bond if i compress or stretch by a tiny amount p2 is a fine approximation and it mimics the solid curve very well but if i go very far on the compressed side or very far on the stretch side then it becomes a terrible approximation p4 is a fourth order polynomial and that's the closer dashed curve and you can see for that one it looks pretty darn good for compressed geometries it mimics the solid curve very well but as i go to uh stretched bonds it uh it does better than p2 but eventually it just fails and it keeps going up whereas the real curve starts to level off the morse curve is this other dashed line and it follows the exact potential very very closely you'll see a little difference at long bond lengths but the difference is tiny compared to the errors of the other ones so if you wanted to look at really stretched bonds then the morse potential is your ticket but a lot of times people only really want to know what's equilibrium geometry or they have vibrations close to the equilibrium geometry and in that case as a cost savings measure p2 is actually fine and a lot more convenient those morse potentials involve exponentials and it doesn't take long to to do an exponential calculation on a computer computer is very fast but if you were doing that for many many many many many pairs of atoms maybe millions of pairs of atoms it's a lot easier to take that bond stretch or compression and square it and get that quickly compared to doing an exponential calculation if you had to do that millions of times then that actually can slow down the calculation so that's another practical reason people prefer the quadratic form over the morse potential okay so let's talk about the bend energy i could have a bond angle like in water say and if i bend the angle in water i would expect that to cause the energy to change and we can reflect that with an equation like this one which is just analogous to the stretching energy except now it involves bond angles theta instead of bond links so the simplest possible form for how does the energy change if i compress a bond angle or expand it would be some force constant which depends on the identities of three atoms so i'll call it a kabc for a bond angle a b c and the difference between what the bond angle is in its desired equilibrium bond angle theta abc 0. so i just get that difference and square it multiply it by the force constant and there you go so it's analogous to this stretch but it's a bend i do need to be a little careful uh with this if i approach 180 degrees because at 180 degrees i then suddenly become linear i can't stretch my hands out totally linear but you get a line and then if you keep going so i approach 180 and then i keep going on the other side then i need to have it symmetric so 175 degrees ought to be the same thing as 185 degrees and that won't necessarily happen unless i somehow tweak my functional form to try to enforce that so again a figure from jensen where you can see the exact bending potential for water is the solid line again and you can see at 180 degrees it's leveled off and if i were to go to 185 190 and keep going i get a mirror image of what's to the left of me so the derivative of the energy at 180 goes to zero because i'm at a local maximum and that won't necessarily happen if i do a quadratic term like i had on the previous slide which is called p2 here p2 is that parabola that just goes up up up up so at large angles approaching 180 that has just totally the wrong behavior if i jazz it up by adding a cubic term i get p3 which is the find dashed line and that goes has an energy too high near 180 degrees and i can fix it uh to get to a derivative of zero at 180 degrees by this p3 prime where i've gone in and tweaked it and enforced that derivative so there i would get a mirror image as i go to angles bigger than 180 but you can see the actual energy in this case is not quite right because i still don't have enough parameters to match that okay now let me talk about a different kind of coordinates so far i've talked about stretches i've talked about bends there's another situation called out of plane bending so if you think about the ammonia molecule nh3 i've got a nitrogen on the top of a pyramid and then i've got three hydrogens down below it and i could talk about how does the energy change as i flatten out that pyramid sometimes this is called an umbrella motion you could define that coordinate various ways in this slide we're defining it as this angle chi which you can see there and i could associate an energy associated with that if i wanted and maybe i just make it a force constant times however big chi is squared that's one example of how you could do that some force fields worry about this term this way some handle it other ways it depends on the force field but certainly the energy will change if i flatten out something like ammonia so i may need to account for that okay the torsion energy this is an interesting one this is something like if i had hydrogen peroxide h-o-o-h i could form here's a hydrogen there's an oxygen here's another oxygen here's a hydrogen this thing the energy changes if i rotate one of the hydrogens around with one of the oh bonds with respect to the other oh bond and i need to capture that in my force field and that's done with this torsion angle which we've depicted here this kind of torsion angle if i look down the bond axis of the two atoms that are bonded say the two oxygens and hydrogen peroxide if i look down that axis if the two end atoms line up exactly that's called zero degrees and then if i rotate one around i get to things like 90 and then if i keep rotating and they're pointed in opposite directions then i get to 180 and then i keep going around 360 is the same as 0.

[00:23:49]all right now you might expect me to say the torsion potential is another taylor series no in this case we don't use a taylor series why is that well in this case we're much less certain that we're going to do small perturbations around the preferred angle if i have a bond stretch it really normally unless i heat up the molecule a lot or something probably won't stretch or compress very much bond angle in water probably doesn't compress too much or stretch too much unless i heat up the system a lot but a torsion can freely rotate in most molecules unless i have like a double bond here or something rotation around these is pretty free and can just spin around so i need a functional form that doesn't assume that i stay near the equilibrium geometry which the stretch and the bend taylor series do so in this case i'll use instead of a taylor series a fourier series and fourier series is perfectly good no matter what angle i'm at whereas the taylor series is only good kind of near equilibrium so what does this look like well equation 7 says i'm going to have a bunch of cosine terms and i take my torsional angle omega sometimes they're called tau depending on what you're reading and there's an n in there n is the order of the term in the taylor series uh 1 2 3 4 or something like that so i have a sum of these terms cosine of omega cosine of 2 omega cosine of 3 omega cosine of 4 omega and each of those terms gets multiplied by a coefficient v a b c d sub n and i add all those up and i get a torsional potential now sometimes this is kind of rewritten in a way to make sure that the energy never goes negative uh and frequently you chop it off at about three terms n equals 3 is the maximum so that could get you something like equation 8 which is just kind of a rewritten version of equation 7. so that's kind of typical now if i have a molecule like ethylene ethylene is a c2h4 so i've got a double bond of two carbons and then i've got four hydrogens around that that needs to be periodic because if i start rotating ethylene because i'd break a double bond so that's a little bit high energy but if i start rotating this around by the time i get to 180 i'm back where i started so i need to be periodic there and that means it turns out only even terms in the fourier series can contribute like n equals 2 n equals 4 etc molecules like ethane where i've got 3 hydrogens on either carbon and just a single bond have odd terms like n equals three n equals six et cetera so uh uh out of jensen's bookie illustrates you know what these fourier terms look like so you can see the top panel i've got cosine 3 omega and what that 3 omega means is if i go from 0 degrees to 360 degrees i get 3 maxima or three minima you can see that it's a maximum at zero degrees and i hit another maximum at 120 and another maximum at 240 and then 360 is the same as zero so there's a total of three maxima or three minima there if i look at cosine two omega then i get two maxima or two minima whichever one you're looking at and again zero is the same as 360. so those don't count as different ones okay so i think you get the idea of the periodicity of these functional forms if i just had n equals one then that's simple there's one maximum and one minimum you