Problem with forming bonds: 2

I tried RotBondPrms, but this doesn't work either -- ie, the bonds don't seem to form. After a few iterations the assembly flies apart.

Below is my source code

#import the appropriate ESyS-Particle modules:
import sys, math
from numpy import *
from scipy import *
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
import random
from POVsnaps import POVsnaps
from POVsnaps2 import POVsnaps2
from POVsnaps3 import POVsnaps3
from WallLoader import WallLoaderRunnable

# instantiate a simulation object
# and initialise the neighbour search algorithm:
sim = LsmMpi(numWorkerProcesses=1, mpiDimList=[1,1,1])
sim.initNeighbourSearch(particleType="RotSphere",gridSpacing=2.5,verletDist=0.5)

# some variable values
normalStiffness = 9.0e09
poisRatio = 0.25
shearStiffness = 3.6e09 # normalStiffness / (2*(1+poisRatio))
kinFriction = 0.4
statFriction = 0.6
compressiveStrength = 6.0e6
tensileStrength = 1.0e6
shearStrength = 3.5e6
tanPhi = 0.577

numsteps = 20
#stepsize = 0.2*particleRad*sqrt(particleRho/normalStiffness)
timestep = 5.0e-5
timeincrement = 1

#set the number of timesteps and timestep increment:
sim.setNumTimeSteps(numsteps)
sim.setTimeStepSize(timestep)

#specify the spatial domain for the simulation:
domain = BoundingBox(Vec3(-120,-120,-120), Vec3(120,120,120))
sim.setSpatialDomain(domain)

# create rectangular array, cubic packing
cube = CubicBlock(dimCount=[6,6,6], radius=1.0)
sim.createParticles(cube)

# create bonds between particles separated by less than the specified maxDist
sim.createConnections(ConnectionFinder(maxDist=0.05, bondTag=1, pList=cube))

# create rotational elastic-brittle bonds between particles
bondGrp = sim.createInteractionGroup(RotBondPrms(name = "pp_bonds", normalK=normalStiffness, shearK=shearStiffness, torsionK=shearStiffness, bendingK=normalStiffness, normalBrkForce=compressiveStrength, shearBrkForce=shearStrength, torsionBrkForce=shearStrength, bendingBrkForce=tensileStrength, scaling=True, tag=1))

# initialize frictional interactions for unbonded particles
sim.createInteractionGroup(FrictionPrms(name="friction", youngsModulus=normalStiffness, poissonsRatio=poisRatio, dynamicMu=kinFriction, staticMu=statFriction))

# create an exclusion between bonded and frictional interactions
sim.createExclusion(interactionName1="pp_bonds", interactionName2="friction")

# add a horizontal wall
sim.createWall(name = "icewall", posn=Vec3(0,-4.0,0), normal=Vec3(0,1,0))

# specify the type of interactions between wall and particles:
sim.createInteractionGroup(NRotElasticWallPrms(name = "elasticWall1", wallName = "icewall", normalK = normalStiffness))

#initialise gravity in the domain:
sim.createInteractionGroup(GravityPrms(name="earth-gravity", acceleration=Vec3(0,-9.81,0)))

#add local viscosity to simulate fluid resistance:
sim.createInteractionGroup(LinDampingPrms(name="linDamping",viscosity=0.95,maxIterations=100))

#add rotational damping
sim.createInteractionGroup(RotDampingPrms(name="rotDamping",viscosity=0.80,maxIterations=100))

# add a checkpointer
sim.createCheckPointer(CheckPointPrms(fileNamePrefix = 'snapshot', beginTimeStep = 0, endTimeStep = numsteps, timeStepIncr = timeincrement))

#add a POVsnaps Runnable:
povcam = POVsnaps(sim=sim, interval=timeincrement)
#povcam.configure(lookAt=Vec3(0,0,-25), camPosn=Vec3(20,0,-25), zoomFactor=0.2)
povcam.configure(lookAt=Vec3(-20,0,-10), camPosn=Vec3(-20,10,-10), zoomFactor=0.02)
sim.addPostTimeStepRunnable(povcam)

#add another POVsnaps Runnable:
povcam2 = POVsnaps2(sim=sim, interval=timeincrement)
#povcam.configure(lookAt=Vec3(0,0,-25), camPosn=Vec3(20,0,-25), zoomFactor=0.2)
povcam2.configure(lookAt=Vec3(0,0,-2), camPosn=Vec3(10,0,-2), zoomFactor=0.08)
sim.addPostTimeStepRunnable(povcam2)

#add another POVsnaps Runnable:
povcam3 = POVsnaps3(sim=sim, interval=timeincrement)
povcam3.configure(lookAt=Vec3(0,0,0), camPosn=Vec3(0,0,15), zoomFactor=0.06)
sim.addPostTimeStepRunnable(povcam3)

#execute the simulation
sim.run()

I found the source of the problem: the time step. When it's too large the solution is unstable and the particles fly apart.

I found the source of the problem: the time step. When it's too large the solution is unstable and the particles fly apart.