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.
EDIT:
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 = ",
humanize_time_hours(total_wall_time),
" 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
omega_12=2.5
De_12 = 15.0
d0_13 = 2.0
omega_13=2.0
De_13 = 20.0
d0_23 = 4.0
omega_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
tmax=50.0
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()
plt.close("all")
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,
verticalalignment='top')
# Set up a pause flag for the animation
pause = True
# Set up an iterator to step through the data during the animation.
it_step=50
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:
return
else:
try:
i = next(it)
ri = rt[i]
data_i = ri.reshape(-1,2)
scat.set_offsets( data_i )
text.set_text("{:6.1f}".format(t[i]))
except StopIteration:
pause = True
it = iter(range(0, N, it_step))
return
fig.canvas.mpl_connect('button_press_event', onClick)
# Start the animation
ani = animation.FuncAnimation(fig, update_with_pause, interval=50)
ax.set_xlim(-25,25)
ax.set_ylim(-25,25)
plt.show()
