OK, after pulling out a significant fraction of the limited amount of hair remaining to me, 3dREMLfit will now accept all zero columns, if you use the -GOFORIT option.
Even with pure OLSQ, there will be slight differences between the 3dREMLfit and 3dDeconvolve results in the same case. These are
* 3dREMLfit adjusts the degrees of freedom so as not to count all zero columns as "real" regressors -- this has the effect of slightly boosting F and t statistics.
* The numerical algorithm used to patch the matrix singularity problem is somewhat different than used in 3dDeconvolve (since the numerical linear algebra in the 2 programs is considerably different at the basic algorithm level), and so the numbers will vary slightly from this cause -- which should in itself be a tiny effect.
* However, if you compare 3dREMLfit GLSQ output with the all zero column included, and then remove the column and repeat the analysis, you will also find a few voxels that are slightly different. These alterations are due to the (a,b) parameters being estimated slightly differently in these few voxels, where the log likelihood function doesn't depend strongly on (a,b). In my test case, this effect was minor. This effect can be further minimized at the cost of CPU time by using the '-Grid 5' option, which uses a finer grid for the estimation of the (a,b) parameters for each voxel. Personally, I wouldn't bother -- the run time about doubled in my test case when using '-Grid 5'.
The changes are now in the AFNI source code. All binaries will be recompiled tonight, and be available by St Patrick's Day. Black and Tans, anyone?