Problem with using rolling resistance in ViscElCap physics

Hi all,
I am trying to add rolling resistance (mR) to my simulation . When I add mR to my script, the simulation does not run (I think it gets stuck somehow as it does not give me any output that it gives in the case I do not add mR). I was wondering what might be the problem.
I have attached my script below.

from yade import pack
from yade import utils, wrapper

###functions for exporting data
os.mkdir(os.getcwd()+'/VTK/')

def savePropData(O):
        from yade import export
        import numpy as np

        path = os.getcwd()+'/VTK/'
        vtkExporter = export.VTKExporter(path)
        vtkExporter.exportSpheres(numLabel = O.iter, what = dict( \
                                      dist = 'b.state.pos.norm()', \
                               linVelocity = 'b.state.vel', \
                               angVelocity = 'b.state.angVel', \
                               mass = 'b.state.mass', \
                               mat_rand = 'b.material.id', \
                               id='b.id' , \
                               numOfContacts = 'len(b.intrs())'))

######Material properties

fr=0.9

rho = 3000

r = 0.0005

Gamma = 0.0
Theta = 0.0
vB = 0.0

CapType="Rabinovich"

kkN= 2*(100)
kkS= 2*(30)
ccN= 2*(0.0005)
ccS= 2*(0.0005)
###en=0.3
#####time parameters
O.dt = 5*1e-6
it=math.floor(0.1/O.dt )
simT=10
#####drum
drumr=0.05
druml=0.03

mat=O.materials.append(
        ViscElCapMat(frictionAngle=fr, density=rho, mR=0.001, Vb=vB, gamma=Gamma, theta=Theta, Capillar=True, CapillarType=CapType, kn=kkN, ks=kkS, cn=ccN, cs=ccS )
)
mat2=O.materials.append(
        ViscElCapMat(frictionAngle=fr, density=rho, mR=0.001, Vb=vB, gamma=Gamma, theta=90, Capillar=True, CapillarType=CapType, kn=kkN, ks=kkS, cn=ccN, cs=ccS )
)

#defining the spheres

sp=pack.SpherePack()
sp.makeCloud((drumr-0.7*drumr,drumr-0.7*drumr,0.0015),(drumr+0.7*drumr,drumr+0.7*drumr,0.0285),rMean=r,num=45000)
sp.toSimulation(material=mat)

Nprtcl=len(O.bodies)
print(Nprtcl)

#liquidMigration
VV=0.0
Vmin=0.0

for s in O.bodies:
        if not type(s.shape)==wrapper.Sphere:
                continue
        s.state.Vf=VV * (4/3) * 3.14*(s.shape.radius)**3
        s.state.Vmin=Vmin

Drum=geom.facetCylinder(material=mat2,center=(0.05,0.05,0.015), segmentsNumber=32, wallMask=4, radius=drumr,height=20*druml,orientation=Quaternion(Vector3(0,0,1),(pi/2.0)))
walls = O.bodies.append(Drum)
O.periodic = True
O.cell.setBox(0.5,0.5,0.03)

##engine
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()], allowBiggerThanPeriod = True),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
                [Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
                [Law2_ScGeom_ViscElCapPhys_Basic()],

        ),

        NewtonIntegrator(gravity=[0, -9.8, 0]),
        RotationEngine(ids=walls,rotationAxis=[0,0,1],rotateAroundZero=True, zeroPoint=[0.05,0.05,0.015], angularVelocity=0.45),
        PyRunner(command='savePropData(O)', iterPeriod=it)
         ]

#Functions

import math

duration=simT/O.dt
O.run( 1 * math.floor(duration),True)

###saving for Restart
O.save('test.bz2')

import pandas as pd

intrState = pd.DataFrame(columns = ['id1','id2','Fn','Fv','sCrit','normalForce0','normalForce1','normalForce2', 'shearForce0','shearForce1','shearForce2'], dtype=object)

for ii in O.interactions:
    iiState = pd.DataFrame({'id1':[ii.id1],'id2':[ii.id2],'Fn':[ii.phys.Fn],'Fv':[ii.phys.Fv], 'sCrit':[ii.phys.sCrit], 'normalForce0':[ii.phys.normalForce[0]], 'normalForce1':[ii.phys.normalForce[1]], 'normalForce2':[ii.phys.normalForce[2]], 'shearForce0':[ii.phys.shearForce[0]], 'shearForce1':[ii.phys.shearForce[1]], 'shearForce2':[ii.phys.shearForce[2]]})
    intrState = intrState.append(iiState,ignore_index = True)

intrState.to_csv('tmpIntrState.csv')

###

Thanks in advance for your response!

Regards,
Roxana