# How may I calculate the bond length between two atoms?

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.

• I am slightly confused as to why you need atom energies for a sandbox game? Can you explain what exactly you want about atoms? If it needs to be very accurate, then you have to use QM, otherwise you can just use Morse potential type equations with fitted parameters. – Shoubhik R Maiti Apr 27 at 11:41
• I explained in topic what I exactly want. I do not have to give away my entire project idea. I do not understand Morse very well, but I assume just like Lennard-Jones Potential it needs those "known" parameters. – John T Apr 27 at 16:31
• All combinations of all atoms. And I think I need both ϵ and σ to use Lennard-Jones Potential. Look up table would be awesome. – John T Apr 27 at 18:32
• @J... , JohnT if you need to talk a lot more about this question, I would recommend this chat room about potential energy surfaces: chat.stackexchange.com/rooms/117652/spectroscopy-potentiology – Nike Dattani Apr 27 at 19:11
• @J.. I totally agree! I just want to make sure the system doesn't try to ask us to move this conversation to chat, because I'd like to avoid this: meta.stackexchange.com/q/353643/391772 – Nike Dattani Apr 27 at 19:14

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).

• Using 3rd party software for my program is not an option. – John T Apr 26 at 19:49
• Please share Ian Miller work – John T Apr 26 at 19:50
• Ian's work is not yet published and I'll have to ask him permission first. However, what do you mean by "using 3rd party software is not an option"? You're using a 3rd party operating system and 3rd party hardware and 3rd party software for programming right? Also, can you please edit your first comment to contain the second one, then delete the second comment? If we get too many comments on a post, it can get us into trouble where we have to move the conversation to a different place, and I'd like to try to avoid that. – Nike Dattani Apr 26 at 19:52
• If it were possible to do it without a 3rd party software, 95% of us would be out of our job immediately – user430 Apr 27 at 5:37
• @Nike I have edited the input file in your answer, hope that's okay with you. – Shoubhik R Maiti Apr 27 at 9:01

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.

• This not for Nobel prize, neither it needs to be very precise. I mentioned that. Please pay attention. I just need aprox calculations and outcomes. Lookup table is not a bad idea. Not sure how to use CSD to get A and B, X,Y. – John T Apr 27 at 16:37
• The problem is to get anything even remotely realistic you need quantum. More simplistic models are not predictive in the slightest as far as bond formation is concerned. That is the key problem you have. If you want bonds to form you need quantum. – Max Apr 27 at 16:41
• Fine. Can you elaborate on lookup table idea (to get A and B, X and Y)? – John T Apr 27 at 16:43
• I mean I like the idea of a lookup table of all possible pairs, but where and what do I use to build such table? – John T Apr 27 at 17:08
• Yes, which is why we keep trying to tell you this will be very hard to do. – Max Apr 27 at 18:38

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.

• I do not really want to share everything about my project. let's keep focused on the topic. – John T Apr 27 at 16:27
• @JohnT It's fine if you don't want to share the details, but the question you have asked is very vague, so without knowing what it is you want to do, I cannot give you a focused answer. – Shoubhik R Maiti Apr 27 at 17:03
• My question is very concrete and focused. Check @Nike Dattani answer. Its on POINT. But sadly doesn't include a solution. Like Max idea of lookup table. – John T Apr 27 at 17:07
• @JohnT Are you amenable to using GNU public licenced open source codes in your software? Then I can attempt to write about using XTB for md. – Shoubhik R Maiti Apr 28 at 8:53
• Yes please....... – John T Apr 28 at 10:36

## 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.

• I seem to have no access to UFF. However even I get it, how do I use these parameters to calculate stuff for 2 atom interaction? All I know for now is Lennard-Jones Potential. – John T Apr 27 at 18:22
• The RDKit code implements UFF. I now mention this in my answer. You can find the original UFF paper here: lammpstube.com/wp-content/uploads/2019/10/UFF.pdf – WaterMolecule Apr 27 at 19:21
• If you want to get structures for organic molecules, you can get them from PubChem in SDF format: pubchem.ncbi.nlm.nih.gov Then you can use openbabel to convert them into formats used by VMD or LAMMPS or whatever. – WaterMolecule Apr 27 at 19:28

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.

• Some of the elements in table have no Radii – John T Apr 27 at 18:17
• Can you be more precise? Like where I get a full table? Paste some example of 2 atoms? Bond length is not all I need - they edited my topic. I must be able to use Morse or Lennard-Jones potential to calculate Energy at distance, etc. – John T Apr 27 at 18:41
• What a simulation person would call van der Waals radii are usually much bigger than the covalent bond distance. For instance, the bond distance parameter between benzene carbons is 1.375 Å in the CHARMM force field, while the van der Waals distance Rmin/2 = 3.9848 Å. – WaterMolecule Apr 27 at 19:12

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).

• How do I prepare the table. Where do I get the data? – John T Apr 27 at 18:21