How do I simulate the interaction between two atoms?

Question

I want to build a rough atom interaction simulator, where pushes two (or more) atoms toward each other and they should behave physically correctly: attract, bond or repel with some force. So that is what I need to calculate.

I have been looking into ionic and covalent bonds theory and that helps.

But still I do don't understand how to calculate the final outcome based on the colliding atoms' types and forces with which they are pushed towards each other.

Input:

• Atom X: and its electronic properties and movement force,
• Atom Y: and its electronic properties and movement force,
• Distance between them

Outcome:

• new force for X and Y + rough changes in electronic properties (bonds etc.)

I got advice to look into the "Morse potential", but reading trough the Wikipedia article is a bit like reading rocket science for me and I can't really 100% associate it with my case. And I think it's mostly about molecules. Is it really also possible to calculate a rough outcome of two colliding molecules (taking some average forces?)?

Any help is appreciated!

Answer 1

EDIT Doing what you want is hard! You will need a full quantum mechanics based simulation. This is unlikely to be something you can build yourself at the current time.

Based on your new additions to the question what you need is Car-Parrinello or Born Oppenheimer MD. These essentially automate the idea of do a quantum electronic structure calculation, take the derivatives of the potential energy landscape with respect to nuclear coordinates to calculate forces, update the nuclear position, rinse and repeat. The difference is BOMD essentially forgets and resolves the entire electronic structure each time step while cpmd tries to calculate how the electronic structure changes (though due to this requires far lower time steps).

Free software does exist out there to do this, but learning how to use it is not always straight forward. A quick google brings up cp2k

More generally the question of computationally modelling exactly how atoms come together is unfortunately a very complex quantum mechanical mess, typically reserved for the fourth or fifth year of a chemistry degree. Forget about two atoms, we are unable to completely solve the problems of one atom with two electrons! While less accurate, non quantum, ways to calculate interactions do exist, these are typically not able to form or break bonds which seems to be what you are interested in. Typically uses are instead protein folding, measuring the energy of crystal structures or modding solvent molecules interacting around a solute (essentially where there are simply too many atoms and the quantum treatment is too slow)

Based on old question

@Camps answer is the most accurate way to go, however from your question I'm assuming complete accuracy is not what you're looking for (if it is something like MP2 or DFT is the way to go, for this you will need a quantum chemistry program, Psi4 is quick to install and free).

If however even the morse potential is too complex can I recommend the Lennard-Jones? Anything simpler (ie using just electrostatics and kinetic energy) and you're just bouncing charged spheres off each other. Essentially this ignores electrostatics (atoms are neutral overall) and only models dispersion and Pauli repulsion. In simple terms dispersion is due to the fact that the electrons moving round each atom aren't always perfectly evenly distributed, at times they will be more on one side or the other which produces an instantaneous dipole. The direction of this dipole is essentially random but as it is lower energy for both atoms to have attracting instantaneous dipoles they will tend to do so on average giving a weak attractive interaction that tends to scale with one over the seperation distance to the sixth power. (a more accurate description is given by inter molecular perturbation theory, I can't find a nice simple description of this though). Pauli repulsion occurs due to the quantum mechanics of electrons. You may have heard of the Pauli exclusion principle which simply states that two electrons can not be in the same state at the same time, this means as we push two atoms together the electrons will tend to shift towards the outsides of the forming bonds, leaving two positively charged hemispheres being pushed towards each other, which of course leads to repulsion. (Again for more accuracy a full quantum description is needed, and falls out of the fact that due to electrons being fermions their wave function must change sign of we swap two of them).

In order to use the LJ potential you will need two constants that determine the strength of the dispersion and repulsion between your atoms/molecules, a list is available in table B2 here. If your two atoms are of different types the geometric mean of the parameters is and often used approximation.

If dealing with ions add the Coulomb force on top of this.

To convert the potential to a force simply take the negative differential. To model its effect you will probably want some sort of discrete time step algorithm such as the Verlet algorithm.

More generally the question of computationally modelling exactly how atoms come together is unfortunately a very complex quantum mechanical mess, typically reserved for the fourth or fifth year of a chemistry degree. Forget about two atoms, we are unable to completely solve the problems of one atom with two electrons! While less accurate, non quantum, ways to calculate interactions do exist, these are typically not able to form or break bonds which seems to be what you are interested in. Typically uses are instead protein folding, measuring the energy of crystal structures or modding solvent molecules interacting around a solute (essentially where there are simply too many atoms and the quantum treatment is too slow) The LJ potential is probably the simplest of these.

Answer 2

How to proceed depends on how accurate you want the outcome. Throughout my answer I will provide blue buttons which demonstrate that there's entire tags in our community to address certain aspects of the simulation.

What you are describing is (essentially) what we call molecular-dynamics even though you're dealing with atoms rather than molecules. In MD (molecular dynamics) simulations, we start with exactly what you have as the input: the initial substances (in your case atoms, but often they are molecules or something else), the distance between them (what we would call "geometry"), and the forces between them (what we would call "forcefields").

A classical MD code like lammps or gromacs or cp2k will take care of this for you, and a software like PACKMOL can help you set up the initial geometry in a format that the MD software packages will understand, and since you're talking about bonding, an appropriate force-field for you might be ReaxFF. With the geometry (perhaps from PACKMOL) and force-field (perhaps ReaxFF) the classical MD software will simply evolve the system using Newton's 2nd law (F = ma) on each individual atom or molecule in your system. You can run the simulation after choosing a temperature, pressure, and amount of time for that you wish to simulate the system (e.g. in picoseconds), and the output will be the new geometry, at which the forces will be different (depending on the force-field you're using). If your output is for example in SMILES format, most software will be able to detect whether a bond has been formed or not based on the geometry, and while this "classical MD" approach involves a lot of empirical data built into the forcefields and the deciding of whether or not bonds have been formed, and a lot of approximations (again in the forcefield and in using Newton's 2nd law rather than the time-dependent Schroedinger equation) it is enough for most cases.

If you want a more accurate simulation, you could do quantum-dynamics or semi-classical dynamics or other things instead of using Newton's 2nd law (which is classical rather than quantum-mechanical), but I don't recommend any of that for a beginner (you said your don't understand the Morse potential so I'm calling you a beginner 😊). You could also do ab-initio calculations or ab-initio-molecular-dynamics instead of empirically determined force-fields if you want to re-calculate the forces at each new geometry during the simulation, but again, I would not recommend this for a beginner. Once you're able to get classical MD done with an empirical force-field, then you can start considering whether or not you want to improve the accuracy of the calculation. You're welcome to ask more follow-up questions if you get stuck along the way, and if you want to learn more about the Morse potential or something else, you can ask separate questions about that!

Answer 3

Generally speaking, most work in molecular dynamics tries to simulate how actual molecules behave (i.e. quantum mechanics) and it doesn't sound like you want to go down that rabbit hole. I completely agree, but I'll begin with a disclaimer that looking up "molecular dynamics" probably won't turn up the kind of results you want.

Since your comments describe wanting a simple simulation based on kinetic energy, forces and collisions (i.e. Newtonian mechanics), it will probably be more helpful to look into more general numerical simulations of Newtonian mechanics.

---

The type of model you are thinking of is very similar to how Newtonian gravity is usually modeled. They are basically identical from a programming perspective, the only difference is that charges are signed (so you can have repulsion as well as attraction).

The simplest way to approach this kind of problem is with numerical integration - set up your initial conditions $x(0)$ and $v(0)$, then update both little bit in each time step $\Delta t$. The most basic form is:

$$ v_x(t + \Delta t) = v_x(t) + a_x(t) \Delta t $$ $$ x(t + \Delta t) = x(t) + v_x(t) \Delta t $$

I used the x-axis here, but you can update the y and z components the same way. This is a really crude numerical integrator (Euler's method) but it's a good start to understand the concept. Once you feel comfortable with it, you can use more advanced, accurate numerical integration schemes.

This acceleration function is really general and you can include as many relevant terms as you see fit. To capture the force of the electric field, you can calculate the acceleration using Coulomb's law, but you can expand it with additional components (e.g. the force generated by the Lennard-Jones potential).

The process to calculate acceleration is typically to calculate the net force on a particle, then divide by its mass to get its acceleration. Again, this is essentially how gravitational acceleration is calculated (so feel free to refer to that material for guidance).

If you want a (really, really) simple model of molecules, you can treat them using rigid body dynamics. You start off by calculating the force on each atom in exactly the same way, then combine those into a net force and torque on each molecule. This can give you a more sophisticated charge distribution than simple spheres (though it will complicate accurate collision detection).

---

Collisions are more complex but can be simplified if you treat them as a collision of two spheres. At that point it's really a 2D ball collision detection problem - I'll refer to other SE questions for how to handle that.

https://stackoverflow.com/questions/345838/ball-to-ball-collision-detection-and-handling?rq=1 https://stackoverflow.com/questions/35211114/2d-elastic-ball-collision-physics

This can be a really tricky problem when you are dealing with discrete jumps in time (what if the particles would have collided in between those steps?). Funnily enough, it helps that Euler's method is so bad - it will force you to use a small $\Delta t$ so you'll catch more collisions, but you'll still miss some unless you implement some clever methods to handle intermediate collisions.

---

Once you've implemented these - congratulations! You've coded a simple model of classical particles. There are plenty of ways to improve it:

1. This numerical integrator is terrible (see: Runge-Kutta 4 or symplectic integration)
2. It won't support too many particles (see: Barnes-Hut or fast multipole methods)
3. The collision detection could be dramatically improved (see: particle-particle collisions)
4. Include simplified approximations of bonding potential (see: Max's answer, Lennard-Jones potential)

To name a few, but you've taken your first step towards a functional model of classical charged particles. Over time you can refine this simulation, better approximate their interactions and perhaps even go down that quantum rabbit hole.

Answer 4

In order to simulate two atom interaction, you have different path to follow.

One is to use Density Functional Theory (DFT) or ab initio, make a script where the distance of the atoms is decreasing, and for each distance, you calculate the system energy. The image below is the result of such calculation but between a metal atom and a carbon nanotube¹. You can use this for atom-atom, atom-system or system-system. This approach is a parameter free method and is independent of other theories.

For your simulator, and to get more accurate and faster results, I recommend to use the semi-empirical method implemented in software like MOPAC (free).

The other way is to use force-field implementation. In this case, you will need o use the "best" potential that represents each type of system. One example is the Morse potential you mention, but there many are others (take a look at the GULP page, for example. The force-field (or potentials) are full of parameters that, in general, are obtained from calculations as described above or when possible, from experimental data.

¹ M. Bastos, I. Camps, First-principles calculations of nickel, cadmium, and lead adsorption on a single-walled (10,0) carbon nanotube. J Mol Model (2014) 20:2094 DOI: 10.1007/s00894-014-2094-y

Answer 5

One thing to be absolutely clear from the beginning is that there is no simple theory of the chemical bonding between two atoms. We do understand a whole lot about the physics that goes into this process, but the physical understanding is not really transformable to formulas that could be reasonably evaluated. Hell, theoretical physicists are still working on fully understanding the energy of a single hydrogen atom.

Since the cutting-edge physics knowledge can not be transformed into real simulations, one always needs to draw the line somewhere and make approximations. The branch of science that deals with these approximations is called quantum chemistry. Quantum chemists know a plethora of different theories and levels of abstractions, and based on the problem in hand they can set up a model. Sadly, there's also no simple recipe on how to choose a model, this all comes down to personal experience (and, in many cases, preference). So this is the second thing to consider: not only the physics can not be accurately described, often the approximations themselves require a good deal of understanding.

Developing the theoretical understanding to apply proper quantum chemical method usually takes years of studying, and includes a lot of pitfalls. To avoid this, I agree with the suggestion of your friend to try using a simplified potential, such as a Morse curve.

There's two main advantages to using a Morse curve as I see:

1. The function that describes the interaction is simple enough that you can focus on the numerical implementation as well as the physics you learn from doing the simulation;
2. The parameters that go into the curve are experimentally determined, so you are poised to get good results even without really understanding the theory.

It is a rather unintuitive, but nevertheless true fact that describing the movement of atoms as classical particles moving in the potential energy field of each other is a good approximation. This means that all you have to do to get a realistic simulation is to take your two atoms, select their initial position and momenta (a lot of science goes into this step too, but I guess you can set them randomly as long as you stay close to the equilibrium energy), then solve the classical equations of motion. This is a relatively simple task to do - I usually assign this to my students as homework. A quick Google search also yielded me this manuscript, which seems to go into a bit more technical details, but I haven't read it myself so I can't vouch for it.

If, for some reason, you don't take my word for the fact that atoms can be modeled classically, you are also free to work with them as quantum mechanical object. Morse potential is also great for this task, since it's eigenvalues and eigenstates are known, so you don't have to bother with implementing anything, you can take a pen and a paper an write down the equation of motion of the system. From this, you can extract all interesting properties - the limit here is your experience with the formalism of quantum mechanics.

EDIT: Reading through the other answers I can see that different molecular dynamics codes had been suggested. While it is always helpful to use existing tools, and they do have wonderful dynamical theories implemented, it is important to keep in mind that out-of-the-box MD usually handles chemical reactivity pretty badly, which seems to be the thing OP is interested in. Because of this, having a better potential (or Hamiltonian, or force field, of however you want to call it) is more important than having good dynamical simulation power. This being said, the popular MD tools sometimes have good reactive force fields implemented, which could work wonderfully here, but they are usually not that straightforward to use.

Answer 6

The simplest possible atomic simulation: a noble gas

As explained by Nike Dattani, what method you use depends on what you want to simulate and why. He gives you a good roadmap for what method to choose based on what you want to do.

I wanted to take a different perspective. If you're familiar with writing code and want to program the simulation yourself, then the easiest atomic system that you can simulate with any reasonable accuracy is a noble gas system. Below is code I wrote in Python for a classical molecular dynamics simulation of liquid neon (using the Lennard-Jones interaction model). The output is a multi-frame PDB file that can be visualized in software such as VMD. The system has periodic boundary conditions. (These are commonly used in simulations and may be confusing to understand at first. You can read about them. They mean that neon atoms exiting the right side of the cubic system will reenter at the left side and so forth for all six sides.)

First, I set the initial positions and velocities of the atoms. And then I run 10000 simulation steps (see the MAIN SIMULATION LOOP). In each step, the net force on each neon atom is computed from the positions of all other atoms using the Lennard-Jones formula. All atoms are then moved forward in time using an approximate time-integration method (explained by the answer of Gumbercules). Each step represents 10 fs of simulated time. The result is a three-dimensional movie of the droplet of neon (in the PDB file).

Note this code is meant to be pedagogical and is not optimized. There are many tricks that can be used to make it run faster for systems of more atoms (and are implemented in codes like LAMMPS and Gromacs). People have spent decades optimizing these codes for molecular simulations, so you are better off using their codes for anything serious.

import math
import random

# Parameters
num_atoms = 20
steps = 10000
out_file = 'sim_neon.pdb'
write_period = 10; # Number of steps between trajectory writing

temper = 26.0; # en K
box = (40.0, 40.0, 40.0); # in Å
timestep = 10.0; # timestep in fs
lj_r = 3.06; # in Å for neon
lj_e = 0.08545; # in kcal/mol for neon
time_factor = 48.88821; # sqrt(angstrom^2*Da/(kcal/mol))
mass = 20.017970; # masa de neon en 
boltz = 0.001987191; # kcal/mol/K


# Pre-computed values
dt = timestep/time_factor
dt2 = dt*dt
v0 = math.sqrt(boltz*temper/mass)
lj_r6 = lj_r**6
lj_r12 = lj_r**12

# -0.5*l <= x <= 0.5*l
def wrap_periodic(x,l):
    imagen = math.floor(x/l)
    x -= imagen*l
    if x >= 0.5*l: x -= l
    return x

def wrap3(position,box_size):
    r=[]
    for c in range(3):
        r.append(wrap_periodic(position[c],box_size[c]))
    return r

def write_frame(output,pos,box):
    name = 'NE'
    resname = 'NE'
    chain = 'A'
    segname = 'A'
    # Write the box size
    output.write('CRYST1 %8.3f %8.3f %8.3f  90.00  90.00  90.00 P 1           1\n' % (box[0], box[1], box[2]))
    for i in range(len(pos)):
        resid = i + 1
        r = wrap3(pos[i],box)
        output.write('ATOM  %5d  %-3s %-4s%s%4d    %8.3f%8.3f%8.3f %5.2f %5.2f      %-4s\n' % (i+1, name, resname, chain, resid, r[0], r[1], r[2], 1.0, 1.0, segname))
        
    output.write('END\n')

    
    
def interact_lj(pi,pj):
    rij = wrap3([ pj[0]-pi[0], pj[1]-pi[1], pj[2]-pi[2] ], box)
    dist2 = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2]

    # force and potential energy for Lennard-Jones interaction (12-6)
    dist6 = dist2**3
    dist12 = dist6**2
    dist8 = dist6*dist2
    dist14 = dist12*dist2
    f_r = -12.0*lj_e*(lj_r12/dist14 - lj_r6/dist8)
    ener = lj_e*(lj_r12/dist12 - 2.0*lj_r6/dist6)

    f=[]
    for c in range(3): f.append(f_r*rij[c])

    # Return energy and force
    return (ener, f)

    
def main():
    # Initial positions
    pos = []
    vel0 = []
    momentum = [0.0, 0.0, 0.0];
    d = 1.0*lj_r
    dx = -0.2*box[0]+0.5*d
    dy = -0.2*box[1]+0.5*d
    dz = -0.2*box[2]+0.5*d
    for i in range(num_atoms):
        dx += d
        if dx > 0.2*box[0]-0.5*d:
            dx = -0.2*box[0]+0.5*d
            dy += d
        if dy > 0.2*box[1]-0.5*d:
            dx = -0.2*box[0]+0.5*d
            dy = -0.2*box[1]+0.5*d
            dz += d

        # Initial position.
        pos.append([dx, dy, dz])
        # Random velocity based on temperature
        vel0.append([random.gauss(0,v0), random.gauss(0,v0), random.gauss(0,v0)])
        # Sum the momentum
        for c in range(3):
            momentum[c] += mass*vel0[i][c]
        
    # Remove any initial net momentum
    residual_mom = [];
    for c in range(3):
        residual_mom.append(momentum[c]/num_atoms)
    for i in range(num_atoms):
        for c in range(3):
            vel0[i][c] -= residual_mom[c]/mass

    # Infer the previous step (use Euler algorithm to start Verlet)
    pos_prev = []
    force = []
    for i in range(num_atoms):
        force.append([0.0, 0.0, 0.0])
        p = []
        for c in range(3):
            p.append(pos[i][c] - vel0[i][c]*dt)
        pos_prev.append(p)
            
    # Write the first frame
    output = open(out_file,'w')
    write_frame(output,pos,box)

    # MAIN SIMULATION LOOP
    for step in range(steps):
        # Zero the forces
        for i in range(num_atoms):
            for c in range(3):
                force[i][c] = 0.0
           
        # Calculate interaction forces
        ener_pot = 0.0
        for i in range(num_atoms):
            for j in range(i+1,num_atoms):
                (e,f) = interact_lj(pos[i],pos[j])
                ener_pot += e
                for c in range(3):
                    force[i][c] += f[c]
                    force[j][c] -= f[c]

        # Do Verlet integration
        pos_next = []
        for i in range(num_atoms):
            p = []
            for c in range(3):
                p.append(2.0*pos[i][c] - pos_prev[i][c] + force[i][c]/mass*dt2)
            pos_next.append(p)
            
        # Calculate the kinetic energy
        ener_kin = 0.0
        for i in range(num_atoms):
            v = []
            for c in range(3):
                v.append((pos_next[i][c] - pos_prev[i][c])/(2.0*dt))
            ener_kin += 0.5*mass*(v[0]**2 + v[1]**2 + v[2]**2)
        ener_total = ener_pot + ener_kin
        mean_temper = 2.0*ener_kin/(3*boltz*(num_atoms-1))

                
        # Cycle the variables
        pos_prev = pos
        pos = pos_next
        
        # Write the current frame
        if step % write_period == 0:
            write_frame(output,pos,box)
            print("step %9d  ener_total %9.4f  mean_temper %8.3f" % (step,ener_total,mean_temper))
    
    output.close()

# main    
if __name__ == '__main__':
    main()