Does anybody know of an existing implementation of a harmonic restraint on rotational degrees of freedom? 

Namely, I have two reference frame of two particles p1 and p2. For simplicity, assume xyz(p1) = xyz(p2). I would like to add a harmonic restraint over the delta in their reference frame, for some definition of that delta.

Thank you
Barak