Raster calculator the difference between the true DEM and your interpolated one, square it, the run something like zonal stats as table to grab the mean, then take the root.
here's some code. I've not tested it. may have to save the raster objects as interim data. let me know if it doesn't work and i'll amend the code.
import arcpy
original_dem = r'c:\\myproject\mygeodatabase\originalDEM'
interp_dem1 = r'c:\\myproject\mygeodatabase\interpDEM1'
interp_dem2 = r'c:\\myproject\mygeodatabase\interpDEM2'
interp_dem3 = r'c:\\myproject\mygeodatabase\interpDEM3'
interp_dem4 = r'c:\\myproject\mygeodatabase\interpDEM4'
interp_dem5 = r'c:\\myproject\mygeodatabase\interpDEM5'
interp_dem_list = [interp_dem1, interp_dem2, interp_dem3, interp_dem4, interp_dem5]
original_dem_rasobj = arcpy.Raster(original_dem)
iteration_counter = 1
for interp_dem_path in interp_dem_list:
ras_object = arcpy.Raster(interp_dem_path)
ras_diff = (ras_object - original_dem_rasobj) ** 2
ras_stats = arcpy.CalculateStatistics_management(ras_diff)
ras_prop = arcpy.GetRasterProperties_management(ras_stats, "MEAN")
ras_mean = ras_prop.getOutput(0)
rmse = float(ras_mean) ** 0.5
print("RMSE of Interpolated DEM Raster " + str(iteration_counter) + " = " + str(rmse))
iteration_counter += 1
print("Processing Complete")