Longitudinal Cholesky decomposition with contrast effect

Thank you so so much Mike!
I included the additions you suggested and the model does run. However, I received the attached Mx script a few days after I posted my inquiry here (it can only handle 2 age groups, and I wasn't able to figure out how to change it to make it fit 4 ages, which is what I need), and the results for 2 age groups seem to differ between Mx and openMx. I'm probably missing something simple here, but I'm not sure what. (I attached the openMx output, the Mx script and the output from Mx.)

Would love to hear your thoughts on this and many thanks again for your time!!
This is the script with your suggestions:

Vars <- c("EXT3", "EXT5")
nvar <- length(Vars)               # number of variables
tnvar <- nvar*2                 # number of variables*max family size
nlower <-nvar*(nvar+1)/2           # number of free elements in a lower matrix nvar*nvar
 
selVars <- paste(Vars, c(rep(1,nvar), rep(2,nvar)), sep="_")
 
## using both DZ-S and DZ-O.
mzData <- as.matrix(subset(twin_data, zygDich==1, selVars))
dzData <- as.matrix(subset(twin_data, zygDich==0, selVars))
 
describe(mzData, skew=F)
describe(dzData, skew=F)
dim(mzData)
dim(dzData)
 
# Generate Descriptive Statistics
colMeans(mzData,na.rm=TRUE)
colMeans(dzData,na.rm=TRUE)
cov(mzData,use="complete")
cov(dzData,use="complete")
cor(mzData,use="complete")
cor(dzData,use="complete")
 
nv <- length(Vars)                 # number of variables
#tnvar <- nvar*2               # number of variables*max family size
#nlower <-nvar*(nvar+1)/2          # number of free elements in a lower matrix nvar*nvar
 
 
# Set Starting Values
meVals    <- c(0.5,0.5)          # start value for means
meanLabs  <- paste(Vars,"mean",sep="") # create labels for means
paVal     <- .6                        # start value for path coefficient
paLBo     <- .0001                     # start value for lower bounds
paValD    <- vech(diag(paVal,nv,nv))   # start values for diagonal of covariance matrix
paLBoD    <- diag(paLBo,nv,nv)         # lower bounds for diagonal of covariance matrix
paLBoD[lower.tri(paLBoD)] <- -10       # lower bounds for below diagonal elements
paLBoD[upper.tri(paLBoD)] <- NA        # lower bounds for above diagonal elements
 
# Create Labels for Lower Triangular Matrices
aLabs     <- paste("a",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")
cLabs     <- paste("c",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")
eLabs     <- paste("e",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_")
 
# Matrices declared to store a, c, and e Path Coefficients
pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
                       values=paValD, labels=aLabs, name="a" )
pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
                       values=paValD, labels=cLabs,  name="c" )
pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE,
                       values=paValD, labels=eLabs, name="e" )
 
# Matrices generated to hold A, C, and E computed Variance Components
covA      <- mxAlgebra( expression=a %*% t(a), name="A" )
covC      <- mxAlgebra( expression=c %*% t(c), name="C" ) 
covE      <- mxAlgebra( expression=e %*% t(e), name="E" )
 
# Algebra to compute total variances and standard deviations (diagonal only)
covP      <- mxAlgebra( expression=A+C+E, name="V" )
matI      <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD     <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")
 
# Algebra to obtain solve(I-B) for sib interaction
bMat <- mxMatrix(type="Full", nrow=2*nv, ncol=2*nv, free=c(F,T,T,F), labels=c(NA,"b","b",NA), name="bMat")
iMat <- mxMatrix(type="Iden", nrow=2*nv, ncol=2*nv, name="iMat")
sibInt <- mxAlgebra(solve(iMat-bMat), name="sibInt")
 
# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins
meanG     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE,
                       values=meVals, labels=meanLabs, name="Mean" )
meanT     <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" )
covMZ     <- mxAlgebra( expression= sibInt %&% (
                          rbind( cbind(A+C+E , A+C),
                                 cbind(A+C   , A+C+E))), name="expCovMZ" )
covDZ     <- mxAlgebra( expression= sibInt %&% (
                          rbind( cbind(A+C+E     , 0.5%x%A+C),
                                 cbind(0.5%x%A+C , A+C+E))), name="expCovDZ" )
 
# Standeraized parameters and CI (getting CIs for squared path estimates)
sta <- mxAlgebra((iSD %*% a)*(iSD %*% a), name="sta") 
stc <- mxAlgebra((iSD %*% c)*(iSD %*% c), name="stc") 
ste <- mxAlgebra((iSD %*% e)*(iSD %*% e), name="ste") 
StandardizedA <- mxAlgebra(A/V, name="StandardizedA") 
StandardizedC <- mxAlgebra(C/V, name="StandardizedC") 
StandardizedE <- mxAlgebra(E/V, name="StandardizedE")
 
CI <-  mxCI(c('sta', 'stc', 'ste'))
 
# Data objects for Multiple Groups
dataMZ    <- mxData( observed=mzData, type="raw" )
dataDZ    <- mxData( observed=dzData, type="raw" )
 
# Objective objects for Multiple Groups
objMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean",
 dimnames=selVars )
objDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean",
 dimnames=selVars )
funML     <- mxFitFunctionML()
 
# Combine Groups
pars      <- list( pathA, pathC, pathE, covA, covC, covE, covP,
                   matI, invSD, meanG, meanT, sta, stc, ste,
                   StandardizedA, StandardizedC, StandardizedE, bMat, iMat, sibInt)
modelMZ   <- mxModel( pars, covMZ, dataMZ, objMZ, funML, name="MZ" )
modelDZ   <- mxModel( pars, covDZ, dataDZ, objDZ, funML, name="DZ" )
fitML     <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction") )
CholAceModel  <- mxModel( "CholACE", pars, modelMZ, modelDZ, fitML, CI)
 
# Run Cholesky Decomposition ACE model
CholAceFit    <- mxRun(CholAceModel, intervals = TRUE)
CholAceSumm   <- summary(CholAceFit, verbose=T)
CholAceSumm

PS. Sorry for not providing the full script before! Wasn't sure what is needed.