proc phi {i1 i2 i3 i4} {
	set PI 3.14159265358979323846
	set atoms [atomselect top "index $i1 $i2 $i3 $i4"]
	set positions [$atoms get {x y z}]
	set rA [lindex $positions 0]
	set rB [lindex $positions 1]
	set rC [lindex $positions 2]
	set rD [lindex $positions 3]
	set phi [expr acos(-[vecdot [veccross [vecsub $rA $rB] [vecsub $rB $rC]] [veccross [vecsub $rB $rC] [vecsub $rC $rD]]] / ([veclength [veccross [vecsub $rA $rB] [vecsub $rB $rC]]] * [veclength [veccross [vecsub $rB $rC] [vecsub $rC $rD]]])) * 180 / $PI]
	return [list $phi]
}