The results in random vibration and response spectrum are based on squaring mode shapes. So the signs of displacement/stresses/strains etc are lost. So taking the cut boundary displacement from a random or response spectrum and using that to do submodeling will not be correct. To get the submodel results correct you have to get the submodel mode shapes and combine these modes. This is not straight forward at this point since when you do submodel of each mode shape it will go in as a separate load step in the results file. Hence you cannot use the file.mcom directly (the file.mcom assumes there is only one load step with each mode being a substep). So you may need to modify the file.mcom. For PSD you can get the modal covariance terms by *GET (only diagonal terms can be done this way) and then combine the modes.