Attractive force

Hello everybody,

I try to study the behavior of nanoparticles.
I initially coded the Brownian motion and now I'd like to model a simple force of attraction between the different particles. This force would appear only when the particles are quite closed and therefore depend on the distance between the particles.
Do you know a way to do this ?

I give you my script here. I'm new to Yade, if you have suggestions or remarks about my code, please, in particular in the last part which is obviously not perfect. However, the code works. I can use the results and trace the MSD, which looks good.
Thank you in advance.

//

import random, numpy
from yade import pack, qt
from yade import plot

# DATA COMPONENTS

a = 1e-6;
# rayon in meters
d = 2*a;
# diameter in meters
rho = 2.2626e+3;
#masse volumique
m = (4 * math.pi * a * a * a * rho)/3
#masse
mu = 1e-3;
# viscosity of water in SI units (Pascal-seconds) viscosite dynamique
nu =1e-6;
# viscosity of water in SI units (m au carre s moins un) viscosite cinematique
kB = 1.38e-23;
# Boltzmann constant
T = 293;
# Temperature in degrees Kelvin
D = kB * T / (6 * pi * mu * a);
#coef de diffusion
gamma = 6*pi*mu*a;
# Coefficient de friction
tau = m/gamma;
# Momentum relaxation time
Deltat = 10e-9;
#pas de temps que JE fixe

print("La valeur de m est {} ".format(m))
print("La valeur de d est {} ".format(d))
print("La valeur de D est {} ".format(D))
print("La valeur de gamma est {} ".format(gamma))
print("La valeur de Deltat est {} ".format(Deltat))
print("La valeur de tau est {} ".format(tau))

Mat=O.materials.append(FrictMat(density=rho))

O.bodies.append([utils.sphere(center =(0,0,0),radius = a, material=Mat, dynamic = True)])

sp=pack.SpherePack()
sp.makeCloud((-1e-5,-1e-5,-1e-5),(1e-5,1e-5,1e-5),rRelFuzz=0,rMean= a, num = 9)

sp.toSimulation()

def brownian():
  for b in O.bodies:
   Vpre=0
   Id = b.id
   f1 = - (gamma*Vpre) + (sqrt(2*kB*T*gamma/Deltat)*random.gauss(0,1))
   f2 = - (gamma*Vpre) + (sqrt(2*kB*T*gamma/Deltat)*random.gauss(0,1))
   f3 = - (gamma*Vpre) + (sqrt(2*kB*T*gamma/Deltat)*random.gauss(0,1))
   O.forces.addF(Id,(f1, f2,f3), permanent = True)
   Vpre = b.state.vel

# FUNCTIONAL COMPONENTS

# simulation loop -- see presentation for the explanation
O.engines=[

   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(

      [Ig2_Sphere_Sphere_L3Geom()], # collision geometry
      [Ip2_FrictMat_FrictMat_FrictPhys()], # collision "physics"
      [Law2_L3Geom_FrictPhys_ElPerfPl()] # contact law -- apply forces
   ),

NewtonIntegrator(damping=0.1),
PyRunner(command='addPlotData()', iterPeriod=1),
PyRunner(command='brownian()', iterPeriod=1)

]

O.dt = Deltat

O.saveTmp()

from yade import plot

#here i would like to do something simpler and save position in different files but i can not and i don't know whi

def addPlotData():
 t = O.time
 position = O.bodies[0].state.pos
 position1 = O.bodies[1].state.pos
 position2 = O.bodies[2].state.pos
 position3 = O.bodies[3].state.pos
 position4 = O.bodies[4].state.pos
 position5 = O.bodies[5].state.pos
 position6 = O.bodies[6].state.pos
 position7 = O.bodies[7].state.pos
 position8 = O.bodies[8].state.pos
 position9 = O.bodies[9].state.pos
 plot.addData(
    tempssurtau = t/tau,
    x=position[0],
    y=position[1],
    z=position[2],
    x1=position1[0],
    y1=position1[1],
    z1=position1[2],
    x2=position2[0],
    y2=position2[1],
    z2=position2[2],
    x3=position3[0],
    y3=position3[1],
    z3=position3[2],
    x4=position4[0],
    y4=position4[1],
    z4=position4[2],
    x5=position5[0],
    y5=position5[1],
    z5=position5[2],
    x6=position6[0],
    y6=position6[1],
    z6=position6[2],
    x7=position7[0],
    y7=position7[1],
    z7=position7[2],
    x8=position8[0],
    y8=position8[1],
    z8=position8[2],
    x9=position9[0],
    y9=position9[1],
    z9=position9[2],
    temps =t,
       )
O.run(100000, True)

plot.plots={
'x':('y'),
'tempssurtau':('x'),
'x1':('y1'),
'x2':('y2'),
#'x3':('y3'),
#'x4':('y4'),
#'x5':('y5'),
#'x6':('y6'),
#'x7':('y7'),
#'x8':('y8'),
#'x9':('y9')
}
plot.plot(subPlots=False)
plot.live=True
plot.autozoom=True
plot.saveDataTxt('Results.txt')

qt.View()
qt.center()