Cannot invert?!

Hi all,

I am getting a cannot invert error in attempting to run a CLPM in twins:

Error in runHelper(model, frontendStart, intervals, silent, suppressWarnings,  : 
 MZ.data: weight matrix is rank deficient; cannot invert

I suppose there is something wrong with my identity and beta matrices. Can someone please advise? The minimal example that gets you there:

library(umx)
iden <-mxMatrix(name = "I", type = "Iden", nrow = 6, ncol = 6)
psi_a <- umxMatrixFree("psi_a", type = "Symm", nrow = 6, 
    dimnames = list(c("X1", "X2","Y1",   "Y2","Z1",  "Z2"),
                   c("X1", "X2","Y1",   "Y2","Z1",  "Z2")),
   labels =   c("ax11",      # X1
                NA,     "ax22",      # X2
                "cAxy1",NA,     "ay11",      # Y1
                NA,     "cAxy2",NA,     "ay22",      # Y2
                "cAxz1",NA,     "cAyz2",NA,     "az11",     # Z1
                NA,     "cAxz2",NA,     NA,     NA,    "az22") # Z2
)
psi_e <- umxMatrixFree("psi_e", type = "Symm", nrow = 6, 
    dimnames = list(c("X1", "X2","Y1",   "Y2","Z1",  "Z2"),
                   c("X1", "X2","Y1",   "Y2","Z1",  "Z2")),
   labels =   c("ex11",      # X1
                NA,     "ex22",      # X2
                "cExy1",NA,     "ey11",      # Y1
                NA,     "cExy2",NA,     "ey22",      # Y2
                "cExz1",NA,     "cEyz2",NA,     "ez11",     # Z1
                NA,     "cExz2",NA,     NA,     NA,    "ez22") # Z2
)

betas <- umxMatrixFree("B", ncol = 6,
    dimnames = list(c("X1", "X2","Y1",   "Y2","Z1",  "Z2"),
                   c("X1", "X2","Y1",   "Y2","Z1",  "Z2")),
   labels =   c(NA,    NA,  NA,     NA, NA,    NA,     # X1
                "arx", NA,  "clyx", NA, "clzx",NA,     # X2
                NA,    NA,  NA,     NA, NA,    NA,     # Y1
                "clxy",NA,  "ary",  NA, "clzy",NA,     # Y2
                NA,    NA,  NA,     NA, NA,    NA,     # Z1
                "clxz",NA,  "clyz", "i","arz", NA)     # Z2
)  
solvBE <-  mxAlgebra('solvBE'  , expression =  solve(I - B) )
Av    <-        mxAlgebra("A", expression =  top.solvBE %&% top.psi_a)
Ev <-            mxAlgebra("E", expression =  top.solvBE %&% top.psi_e)
CovMZ <-    mxAlgebra("CovMZ", expression =  rbind(
                cbind(A +  E, A ),
                cbind(A , A  + E)
                ))
CovDZ <-    mxAlgebra("CovDZ", expression = rbind(
                cbind(A +  E, 0.5 %x% A ),
                cbind(0.5 %x% A , A  + E)
                ))
means <-    mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE,
                labels=c("mX1","mX2","mY1","mY2","mZ1","mZ1"),
                name="means" )
expMean <- mxAlgebra("expMean",
            expression = cbind(top.means ,
                top.means ))
expMZ = mxExpectationNormal("CovMZ",
            means = "top.expMean", dimnames = c("X1_T1","X2_T1","Y1_T1","Y2_T1","Z1_T1","Z2_T1",
        "X1_T2","X2_T2","Y1_T2","Y2_T2","Z1_T2","Z2_T2"))
expDZ =  mxExpectationNormal("CovDZ",
            means = "top.expMean", dimnames = c("X1_T1","X2_T1","Y1_T1","Y2_T1","Z1_T1","Z2_T1",
        "X1_T2","X2_T2","Y1_T2","Y2_T2","Z1_T2","Z2_T2"))
top <- mxModel("top",expMean, means,iden,solvBE, psi_a, psi_e,betas)
MZ = mxModel("MZ",  expMZ,  CovMZ, Av, Ev,mxFitFunction() )
DZ = mxModel("DZ",  expDZ,  CovDZ, Av, Ev,mxFitFunction())
model <- umxSuperModel("clpm3",top,MZ, DZ)
model <- umxModify(model, autoRun=FALSE, free = TRUE,
    update =  c("ax11",      
        "ax22",      
        "cAxy1",    "ay11",      
        "cAxy2",  "ay22",      
        "cAxz1",    "cAyz2",   "az11",     
        "cAxz2",   "az22","ex11",     
        "ex22",      
        "cExy1",    "ey11",      
        "cExy2",   "ey22",      
        "cExz1",  "cEyz2",     "ez11",     
        "cExz2",    "ez22","arx",    "clyx",   "clzx",      
        "clxy",   "ary",    "clzy",     
        "clxz",   "clyz", "i","arz"), 
    value = c(.7,      
        .7,     
        .6,    .6,      
        .6,  .6,      
        .6,    .6,   .65,     
        .6,   .65,.65,    
        .3,      
        .6,    .3,    
        .6,   .3,    
        .6,  .6,     .3,  
        .6,    .3,.8,    .2,   .2,  
        .2,  .8,    .2,      
        .2,   .2, .1,.8)
)
model <- mxGenerateData(model, nrows = 300, returnModel = T) |> mxAutoStart()
model <- umxModify(model, update = "i", autoRun = FALSE)
model <- mxTryHardOrdinal(model)
model$MZ$data$observed |> summary()
mxGetExpected(model, "covariances")$MZ |> eigen()