from chimerax.atomic import Residues, Atoms, selected_residues
from math import sqrt
import numpy
models = selected_residues(session).unique_structures
from chimerax.core.commands import run

run(session, f'match {"|".join(["#"+m.id_string for m in models[1:]])} to #{models[0].id_string}')

for residues in zip(*[m.residues for m in models]):
    residues = Residues(residues)
    pas = Atoms(residues.principal_atoms)
    if not len(pas):
        residues.atoms.bfactors = 0
        continue
        
    coords = pas.scene_coords
    avg = numpy.mean(coords, axis=0)
    deviations = coords-avg
    distances = numpy.linalg.norm(deviations, axis=0)
    rms = sqrt(numpy.mean(distances**2))
    residues.atoms.bfactors = rms

run(session, f'color bfactor sel')