How may I calculate the bond length between two atoms?

Question

I have been looking into Lennard-Jones Potential and Morse potential. And what the formulas do is to calculate a potential and the 'well depth' and 'distance at which the intermolecular potential between the two particles is zero' and these are known input parameters. So I guess I am looking at the wrong place. Because I want to take two atoms (with they properties) and calculate bonding length and as well the energy at a given distance. How do I do that? It doesn't have to be very precise.

Edit:

P.S. I am building a sandbox-game-like application with c++. All interaction calculations happen on fly (at least ones that can't be cached). Using 3rd party executables is therefore not an option.

Answer 1

The only proper way to do it ab initio is to calculate the energy (E) for various internuclear distances (R) and then to approximate the equilibrium bond length by choosing the distance which has the lowest energy, or by fitting the E(R) points to something like a Morse potential (which is decent for short-distances), Lennard-Jones potential (which is decent for large distances), or the one that I helped invent: Morse/Long-range potential (which works for both short and long-distances because mathematically it can be a Morse potential for short distances and Lennard-Jones potential for long ones!).

People have attempted to come up with "heuristic" ways to get the bond length in a more "direct" way (for example Ian Miller in some unpublished work that I can ask him to share with you), but there is not yet any successful theory that will work even qualitatively for most of the elements in the periodic table. You are best off with the approach described in paragraph 1, and thanks to a lot of easy-to-use free/open-source software, it's not too hard now, for example in ORCA:

! UHF OPT def2-QZVPP 
%geom Scan 
B 0 1 = 0.3, 1.3, 30
end
end

* xyz 0 1
   H       -4.61685        1.79381        0.00000
   H       -4.14986        1.26166        0.00000
*

will give you the $\ce{H2}$ potential energy surface, from which you can find the equilibrium bond length as described in my firt paragraph (see Shoubhik's answer here for more detail).

Answer 2

Unfortunately the answer to this question is much the same as your last one. What you are trying to do is hard. Very hard. Creating this would be PhD if not postdoc work and doing it well could well win you a nobel prize. From the comments you've made in this and your other question you do not appear to be a specialist in either physics or chemistry, and thus as much as I hat to say this this is not something you can do at the moment. If this is really something you want to do getting a physical chemistry and/or a physics degree is the best way to go.

So writing something yourself aside let's have a quick look at the practicalities. To predict bonds you need to a accurately model electron density of wavefavefunctions. You can not just use simple Lewis Structure type models. So we need an initio or DFT calculations. How practical would it be to use these in the type of program you are describing?

Unfortunately not at all. In order to predict bond lengths these systems work by optimising an initial system to find its lowest energy, performing a full quantum electronic structure calculation at each step. They are far to slow. For just one or two atoms they are reasonably fast, on my reasonably high end consumer grade PC they will run in a few seconds. The problem is that the scaling of these types of calculations is quite horrendous. The best scaling is DFT of Hartree Fock. These scale like the size of the system to the fourth power. So let's go from two molecules to say 500. Now instead of a few seconds were looking at just under 2000 years. Of course we can get round this using supercomputers, but that's hardly practical for a game (and there is still a limit, for example we can't really run full QM calculations on proteins). In order to model a visible amount of matter you would need on the order of ten to the twenty fourth power atoms. This is simply impossible to model.

Using classical force fields are also an option. These are much faster but will not model bond breaking and forming. Also by the time you get to visible quantities of atoms it will still be too slow.

Your best approach is probably to just forget about trying to predict bonds and instead just use a lookup table of something like A bonded to B has length X, energy Y, forms at distance Z with probability P and has probability P' of breaking. Obviously you can make A and B multi atom clusters to improve accuracy, but this will also make your lookup table much larger. Trawling the CSD would probably be the easiest way to get started.

Download a large number of CIF files, these contain information about what atoms are present in a crystal, where they are and how that are bonded. It's all just text so reasonably understandable. Iterate through them and calculate the average bond lengths for atom pairs. Bond energies tend to be on the order of a few hundred kJ/mol. You can probably approximate them as <a few hundred kJ/mol> plus some function of bond length. A 1/r type dependence seems reasonable but you'll have to play with it and see what works well (Incase it's not obvious we've left accurate modeling well behind now). As for capture distance, the sum of the Van Der walls radii is probably reasonable. Bond breaking probability can probably be approximated with boltzman style terms (that is $P_{break} = e^\frac{-E_{bond}}{RT}$, again the constants will require tuning to produce something that looks reasonable, -(E$_{bond}$ - E$_k)$ might work well on the numerator instead). Bond formation could be something like 1-P$_{break}$. You could also run quantum calculations for all possible pairs (or higher order clusters) of atoms or hit the literature to look up experimental bond strengths for greater accuracy.

Note the heavy use of probably, this will not produce anything remotely accurate or predictive, but with a bit of tweaking could look reasonable. You'll obviously need some kind of long range interactions as well (here an LJ potential plus some numeric integration to update positions and speeds each time step as in the other answers is probably what you'll want) and the check all pairs of atoms against the table to see if they bond. Tuning this will be all about making something that looks good in your game and is something you will have to do yourself.

Answer 3

I don't know why you need atom-atom interactions in a sandbox game. My only guess is that you want to simulate or visualize the movement of molecules of some kind in real time. If your system is large i.e. more than ~20 atoms, you have to resort to molecular dynamics to get reasonable speed. This has already been covered in detail in Max's answer, so I will not go into that.

Instead, I will try to write about how you can implement that in your code. There are already some open-source codes that have C++ API that you can use for molecular dynamics; for example, simple-molecular-dynamics.

If you need a robust library for molecular dynamics (MD) simulation, you can use LAMMPS which can be compiled as a C library and called from another software (see LAMMPS library interfaces). However, LAMMPS is a dedicated MD code with a large number of features, so it might be overkill for your purposes.

Finally, an option that I think might be the best for your purpose is to the software XTB. It is an open-source software that exposes a C API for running calculations. The software has the option to use polarizable GFN-FF forcefield which is quite fast. The documentation of the API and an example use is given here. If the API is too complicated to work, then just take the easy route and run the compiled executable with system().

These are all open source codes under GNU public license, so there should be no issue with you including them in your program, as long as you acknowledge the original author.

Answer 4

Get your parameters from an existing force field (UFF? ReaxFF?)

People in the molecular dynamics world use established force fields to describe forces between different types of atoms. These force fields are mathematical functions that have parameters like bond lengths, bond angles, dihedrals, van der Waals radii, etc. Usually these parameters are obtained by optimization to match quantum mechanical calculations or experimental data. Making a good force field requires decades of work!

UFF

The simplest force field to implement that might do what you want is the UFF. The mathematical form is fairly simple and each type of atom has only a few parameters. However, it's generally regarded as pretty inaccurate, but it covers almost all elements and probably would be good enough for a sandbox game.

UFF is implemented in the RDKit cheminformatics Python/C++ code. You can find the code implementing it here.

ReaxFF

I think the best thing for your sandbox game would be to use ReaxFF, which would allow you to represent chemical reactions with a lot of different molecules and substances. The mathematical formulas describing the interatomic potential energies in this force field are very complex (and then you have to differentiate them to get the forces). See the supporting information of Chenoweth, van Duin and Goddard. Luckily there is a free and open source (GPL'd) implementation of this force field. See the LAMMPS code. The implementation of the force field terms are included in the reaxc_*.cpp files. Various parameters for this force field exist and some are included with LAMMPS in the examples/reaxc directory. Parameters for compounds of H, C, O, and N can be found in the supporting information of this paper. Here's the thing though, this force field consists of hundreds of parameters. There isn't just one "bond length", but many different values that describe bonding and nonbonded interactions in a rather complex way. There are quite a few different ReaxFF parameter sets floating around and they encode a lot of information about how chemistry works.

Answer 5

The most trivial model to estimate bond lengths between two atoms is probably a simple sum of Van-der-Waals Radii. Just load a table with Van-der-Waals Radii and take the equilibrium distance for two atoms as the sum of the Van-der-Waals Radii of both atoms. This is by no means accurate but for the intention of the poster it might suffice.

If this is still too complicated, you can take as a poor thumb of rule (with many exceptions in reality) that the Van der Waals radius decreases for atoms in a period and while it increases in a group and build a function that maps the atomic number to a Van der Waals radius following this rule of thumb.

I endorse this of course in no way for actual quantitative scientific work. And other answers have already explained the proper scientific approach to this question which are preferable by any means.

https://en.wikipedia.org/wiki/Van_der_Waals_radius

Answer 6

I don't think you are looking for anything too physically correct so I think for your purposes, you can get away with a simple morse model.

Consider which elements you would like to use, hopefully you can restrict it to basic organic elements such as C,O,N,H,Cl,B and consider which elements that need to be able to bond together. Prepare a table such that you can look up bonds in the format of "bond_distance["C"]["O"]" with some bond distance calculated or looked up. You can probably look up common bond distances or just find a structure with that type of bond and call it good enough.

The problem comes when you want to model a correcting force such that you create something that looks like a morse potential. The bottom of the well will be found at that bond_distance discussed before, but the shape of the potential is harder to determine without proper calculations. You will need to fit a well depth and an alpha parameter which determines the well width.

To be honest, if you do not need a physically correct system you might just play around with these parameters and find something that works well. You might be able to assume all elements have similar values and it works for your game. If you want to get more accurate numbers you may need to find a computational researcher to run a dataset for you (but this might cost you).