I want to build a rough atom interaction simulator, where is push 2 (or more) atom towards each other and they should behave some have physically correct - 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 colliding atoms type and force they are pushed towards each other.

Input (Atom X - it's electronic properties and movement force, Atom Y and it's electronic properties and movement force, distance between) Outcome => (new force for X and Y + rough changes in electronic properties (bonds etc))

Any help appreciated.

  • 2
    $\begingroup$ You need to decide on a model. How exactly do you want them to interact ? Purely classical like charged point particles ? Or do you also want to allow bonding which would require a quantum mechanical description. $\endgroup$
    – Hans Wurst
    Apr 21, 2021 at 18:02
  • $\begingroup$ This is a non trivial task, depending on the sophistication/accuracy that you want to achieve. One approach to describe dynamic bonding of atoms/molecules are molecular dynamics with reactive force fields. $\endgroup$
    – Hans Wurst
    Apr 21, 2021 at 18:28
  • $\begingroup$ As I wrote a "rough" simulator. I was hoping to make some function that takes 2 atom flow direction, speed, atom type and could calculate the outcome (bond, repeal, etc.). Can't we just take into account some average charge diff between atoms + the bounding and done? $\endgroup$
    – John T
    Apr 21, 2021 at 18:35
  • $\begingroup$ You could simply take a Morse potential if you want to model two atoms interacting. The potential provides you forces via its negative gradient that you can plug into the integration of motion for your atoms. $\endgroup$
    – Hans Wurst
    Apr 21, 2021 at 18:47
  • 2
    $\begingroup$ Crossposted to mattermodeling.stackexchange.com/q/4811 $\endgroup$
    – Qmechanic
    Apr 27, 2021 at 20:40

1 Answer 1


We can model the interaction of two atoms or more by setting up a potential. The negative gradient of the potential defines the force that acts on the atoms. The forces can then be used in your simulation of atoms "moving around" and possibly "bonding". This approach treats nuclei as classical particles and is typically called a MM/MD(molecular mechanics/ molecular dynamics) simulation. The potential needs to be modeled and you are unlikely to find a simple potential that describes all physical effects correctly. You can also define a forces/force field instead of the potential, since you need only the forces to integrate the equation of motion. Finding good potentials and force fields is active research and a comprehensive discussion goes beyond a simple Q&A.

A way to obtain the potentials/forces would be a potential energy surface scan with electronic structure calculations. Such a scan would allow you to obtain a potential numerically, based on ab-initio calculations that take only atomic parameters as input. This approach is computational expensive and requires access and likewise importunately, the ability to use electronic structure programs properly. Another possibility are on-the-fly dynamics where you couple your integration of motion to function calls to electronic structure programs to provide forces for the current system at each step. This avoids the need to calculate the full potential.

Another possible way is to define the potentials/forces by yourself. A simple potential to model bonding between two atoms is the Morse potential. You could do a very simplistic approach by defining a Morse potential for each pair of atom types that you have in your system and the total potential as a simple sum of all pair potentials. The forces are provided by the negative gradient of this potential and can be calculated analytically. You could model the Morse potential defining parameters based on the atom type pairs. I.e. the potentials parameters could then be functions of two atomic numbers.

This would be a very simplistic approach since the real physical potential is not simply a sum of pair potentials but such a model can certainly serve as first step to get an idea how to proceed.

You could also try to ask for more details on https://mattermodeling.stackexchange.com/ for help. That site is pretty much all about modelling.


I have added some horrible spaghetti code in Python how a very simplistic model with 3 particles and the Morse potential could look like. The code is pretty bad since i was to lazy to generalize the construction of all possible pairs. I have also made up the parameters manually instead of building them from atomic numbers. The particles have all mass 1 and i did not base the parameters on any actual atoms. To generalize this one would have to build a logic that handles the pairwise interactions nicely in the gradient calculation and make a function to create the potential parameters based on atomic properties like the atomic number.

The model includes dampening of the momentum. The behavior of the particles
can be changed by fiddling with the parameters that define the pairwise potentials.

It should work with any standard python installation of Version 3 that has the common packages numpy and matplotlib. To start/pause the animation left click into the matplotlib window.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

import time

def humanize_time_hours(secs):
     mins, secs  = divmod(secs, 60)
     hours, mins = divmod(mins, 60)
     return "%02d:%02d:%02d" % (hours, mins, secs)

def DoPri45Step(f,t,x,h):
    Ref : https://stackoverflow.com/questions/54494770/how-to-set-fixed-step-size-with-scipy-integrate
    k1 = f(t,x)
    k2 = f(t + 1./5*h, x + h*(1./5*k1) )
    k3 = f(t + 3./10*h, x + h*(3./40*k1 + 9./40*k2) )
    k4 = f(t + 4./5*h, x + h*(44./45*k1 - 56./15*k2 + 32./9*k3) )
    k5 = f(t + 8./9*h, x + h*(19372./6561*k1 - 25360./2187*k2 + 64448./6561*k3 - 212./729*k4) )
    k6 = f(t + h, x + h*(9017./3168*k1 - 355./33*k2 + 46732./5247*k3 + 49./176*k4 - 5103./18656*k5) )

    v5 = 35./384*k1 + 500./1113*k3 + 125./192*k4 - 2187./6784*k5 + 11./84*k6
    k7 = f(t + h, x + h*v5)
    v4 = 5179./57600*k1 + 7571./16695*k3 + 393./640*k4 - 92097./339200*k5 + 187./2100*k6 + 1./40*k7;

    return v4,v5

def DoPri45integrate(f, t, x0, show_progress=True):
    Ref : https://stackoverflow.com/questions/54494770/how-to-set-fixed-step-size-with-scipy-integrate
    N = len(t)
    x = [x0]
    for k in range(N-1):
        if show_progress==True:
            if k == 0:
                prog_n = 1
                start_time = time.time()
                block_time0 = time.time()
            if k > prog_n*N/100*10:
                block_time1 = time.time()
                delta_time = block_time1 - block_time0
                print("{:12.0f}/{:} = ".format(prog_n*N/100*10,N),
                      "{:3.0f} %".format(k/N*100),
                      " wall_time[sec] = {:4.4e}".format(delta_time),
                      "HH:MM:SS = ", humanize_time_hours(delta_time) )
                prog_n +=1
                block_time0 = block_time1
        v4, v5 = DoPri45Step(f,t[k],x[k],t[k+1]-t[k])
        x.append(x[k] + (t[k+1]-t[k])*v5)

    if show_progress == True:
        block_time1 = time.time()
        delta_time = block_time1 - block_time0
        print("{:12.0f}/{:} = ".format(N,N),
                      "{:3.0f} %".format(N/N*100),
                      " wall_time[sec] = {:4.4e}".format(delta_time),
                      "HH:MM:SS = ",time.strftime('%H:%M:%S', time.gmtime(delta_time)))
        total_wall_time = time.time() - start_time
        print("Total time[sec] = {:4.4e}".format(total_wall_time), "HH:MM:SS = ",
              " for a total of N={:} steps".format(N))
        single_step_time = total_wall_time/N
        print("Average time for a single step = {:4.4e}\n".format(single_step_time),
              "Average steps per second = {:4.4e}\n".format(single_step_time**-1))
    return np.array(x)

def morse_potential(r1,r2, d0, omega, de, v0=0.0):
    d = np.linalg.norm(r1-r2)
    a = omega*np.sqrt(2*de)**-1
    v = de*( 1.0 - np.exp(-a*(d-d0)) )**2 + v0
    return v

def morse_potential_derivative(ra, rb, d0, omega, de, v0=0.0):
    a = omega*np.sqrt(2*de)**-1
    d = np.linalg.norm(ra-rb)
    ef = np.exp(-a*(d-d0))
    x1a, x2a = ra
    x1b, x2b = rb
    dvdx1a = 2*de*a*(1-ef)*ef * (x1a-x1b)/d
    dvdx2a = 2*de*a*(1-ef)*ef * (x2a-x2b)/d
    dvdx1b = 2*de*a*(1-ef)*ef * (x1b-x1a)/d
    dvdx2b = 2*de*a*(1-ef)*ef * (x2b-x2a)/d
    return np.array([[dvdx1a, dvdx2a], [dvdx1b, dvdx2b]])

# Initial conditions
r1 = np.array([0.0, 0.])
r2 = np.array([10.0, 0.])
r3 = np.array([5.0, 1.])
p1 = np.array([0.5, 0.0])
p2 = np.array([-0.5, 0.0])
p3 = np.array([-0.5, 0.0])

# Potential Parameters
d0_12 = 2.5
De_12 = 15.0

d0_13 = 2.0
De_13 = 20.0

d0_23 = 4.0
De_23 = 30.0

paras = {"12":[d0_12, omega_12, De_12], 
         "13":[d0_13, omega_13, De_13], 
         "23":[d0_23, omega_23, De_23]

# Time for integration
t = np.linspace(0,tmax, 10000)            

# Derivative function with dampening
def dot_y(t,y):
    q = y[0]
    p = y[1]*np.exp(-t*np.log(2)*4/tmax)
    dot_q = p
    ra, rb, rc = q[0:2],q[2:4], q[4:6]
    nabla_ab = morse_potential_derivative(ra, rb, *paras["12"])
    nabla_ac = morse_potential_derivative(ra, rc, *paras["13"])
    nabla_bc = morse_potential_derivative(rb, rc, *paras["23"])
    nabla_a = nabla_ab[0] + nabla_ac[0]
    nabla_b = nabla_ab[1] + nabla_bc[0]
    nabla_c = nabla_ac[1] + nabla_bc[1]
    nabla_total = np.array([nabla_a, nabla_b, nabla_c]).flatten()
    dot_p = -nabla_total
    return np.array([dot_q, dot_p])

# Integration of the eom
# Create coordinate/momentum vector for the integration routine
R0 = np.array([r1,r2,r3]).flatten()
P0 = np.array([p1,p2,p3]).flatten()
y0 = [R0, P0]
yt = DoPri45integrate(dot_y, t, y0 )

# Unpack values for plot/animation
rt = yt[:,0] # rt[0] = xa0, ya0, xb0, yb0
rtx = rt[:,::2].flatten()
rty = rt[:,1::2].flatten()

fig, ax = plt.subplots()
colors = ["b", "r", "cyan"]
scat = ax.scatter(rtx[0:3], rty[0:3], c=colors)

# Initialize a text field that show the current time of the trajectory
text = ax.text(0.025, 1.0, "0", transform=ax.transAxes, fontsize=14,

# Set up a pause flag for the animation
pause = True
# Set up an iterator to step through the data during the animation.
N = len(t)
it = iter(range(0, N, it_step))

# Define function that toggles the pause flag
def onClick(event):
    global pause
    pause ^= True

# Define the function for the animation
def update_with_pause(l):
    global pause
    global it
    if pause == True:
            i = next(it)
            ri = rt[i]
            data_i = ri.reshape(-1,2)
            scat.set_offsets( data_i )
        except StopIteration:
            pause = True
            it = iter(range(0, N, it_step))

fig.canvas.mpl_connect('button_press_event', onClick)
# Start the animation
ani = animation.FuncAnimation(fig, update_with_pause, interval=50)
  • $\begingroup$ Great info. I would love to seem some example calculations. Reading on Morse potential is hard. $\endgroup$
    – John T
    Apr 25, 2021 at 11:38
  • 1
    $\begingroup$ @JohnT I have added an example written in Python. The code is not in good form but working and should illustrate what parts are required to model a simple system with a Morse potential. $\endgroup$
    – Hans Wurst
    Apr 25, 2021 at 16:43
  • $\begingroup$ What are 'Initial conditions' and 'Potential Parameters' ? $\endgroup$
    – John T
    Apr 25, 2021 at 18:25
  • $\begingroup$ The initial conditions are the initial positions and momenta of the particles that are required to integrate the equation of motion. The potential parameters are the parameters that appear in the definition of a Morse potential. You do understand that Morse potential is just a term to tell you which functional form is used ? An explicit instance of a potential needs numeric values for the parameters that define it so it can be evaluated numerically. $\endgroup$
    – Hans Wurst
    Apr 25, 2021 at 19:48

Not the answer you're looking for? Browse other questions tagged or ask your own question.