Indeed, this was a topology that we 'failed to expect' when we wrote the Pythia colour tracing for epsilons. It won't be totally trivial to implement, so I can't promise how fast I can do it, but will be working on it as soon as possible. Nishita, the topology is basically a resonant squark production via BNV, but with an ISR gluon already added by MG, so it looks e.g. like (using negative numbers for anti colour tags and positive for colours)

  \  qbar(-101)
   \ 
    X - - - - ~t(103) - - - --> decay to top etc 
   / 
  /~~~ g(102,-104) 
 /
/ qbar(-104)

Except that the virtual quark propagator carrying the intermediate -102 colour tag is not explicit in MadGraph, so what we get in the actual LHE event looks like

  \  qbar(-101)
   \ 
    XXX - - - - ~t(103) - - - --> decay to top etc 
    XXX
    XXX ~~~ g(102,-104) 
   /
  / qbar(-104)

where the big block of Xs is just ASCII for a blob, the MG hard interaction. This colour flow kills Pythia now, and somehow we need to make it able to recognise that it just corresponds to the one I drew above. One possible (albeit fragile) help is that MG seems to really strictly have adopted the proposed convention that the colour tags in an RPV vertex should be following each other in sequence, so we know there has to be an (101,102,103) junction. The tricky bit is realising that although the 102 tag here is in the final state, it actually corresponds to an *incoming* colour line from the point of view of the junction.