SE computations break when variances in bifactor model estimated

Code to reproduce issue:

library(mirt)

set.seed(12345)
nexaminee <- 2000
nitem <- 20
ntestlet <- 4 # 5 items per testlet
maxitemsc <- 3

## Generate item parameters
amat <- matrix(NA, nitem, ntestlet+1)
amat[,1] <- 1
amat[1:5,2] <- 1
amat[6:10,3] <- 1
amat[11:15,4] <- 1
amat[16:20,5] <- 1

bmat <- matrix(0, nitem, maxitemsc)
for (j in 1:nitem){
    bmat[j, ] <- sort(rnorm(maxitemsc)) # , decreasing=TRUE
}

## Transform item step parameters to intercept parameters
bmat_int <- matrix(0, nitem, maxitemsc)
for (j in 1:nitem){
    for (l in 1:maxitemsc){
        bmat_int[j, l] <- sum( bmat[j, 1:l] )
    }
}

mu <- rep(0, ntestlet+1)
vcov <- matrix(0, ntestlet+1, ntestlet+1)
# diag(vcov) <- c(1, seq(0.1, 1, len=ntestlet))
diag(vcov) <- 1

bmat_0int <- cbind(rep(0, nitem), bmat_int)
resp_sim_bf <- simdata(amat, bmat_0int, N=nexaminee, itemtype='gpcm', sigma=vcov)

testlet_loading = rep(1:ntestlet, each=(nitem/ntestlet))
model_spec <- mirt::mirt.model("G = 1-20
  START = (1-20, a1, 1.0), (1-5, a2, 1.0), (6-10, a3, 1.0), (11-15, a4, 1.0), (16-20, a5, 1.0)
  FIXED = (1-20, a1), (1-5, a2), (6-10, a3), (11-15, a4), (16-20, a5)
  COV = S1*S1, S2*S2, S3*S3, S4*S4")

mod_bf1 <- mirt::bfactor(as.data.frame(resp_sim_bf), model=testlet_loading, 
                         model2=model_spec, SE=TRUE, itemtype = 'gpcm') #
coef(mod_bf1, simplify=TRUE)

philchalmers / mirt

SE computations break when variances in bifactor model estimated #170