Just thought that this might be useful for some who would want to compute angles between two mean planes, say for two stacking aromatic rings.

Because the plane equation as solved in CHARMM can have two solutions (defined up to inversion of the normal vector) straightforward computation of angles may not give consistent results.

This can be trivially addressed with the following script chunk that monitors the normal direction w.r.t. selected cross products generated by arbitrary three atoms from each of the two planes:

set pi = 3.14159265358979

define plane1 sele ... end

define plane2 sele ... end

define nc1 sele atom1_from_plane1 end

define nc2 sele atom2_from_plane1 end

define nc3 sele atom3_from_plane1 end

define pc1 sele atom1_from_plane2 end

define pc2 sele atom2_from_plane2 end

define pc3 sele atom3_from_plane3 end

open write card unit 22 name out_file.lsqp

!Your standard way of reading trajectories

traj query unit @istart

calc nframes = ?nfile * @pcs

traj iread @istart nunits @pcs nfile @nframes skip 1

set f 1

label tloop

traj read

format (f16.8)

coor axis sele nc1 end sele nc2 end

set vx1 = ?xaxi

set vy1 = ?yaxi

set vz1 = ?zaxi

coor axis sele nc1 end sele nc3 end

set vx2 = ?xaxi

set vy2 = ?yaxi

set vz2 = ?zaxi

calc cpx = @vy1*@vz2 - @vy2*@vz1

calc cpy = - @vx1*@vz2 + @vx2*@vz1

calc cpz = @vx1*@vy2 - @vx2*@vy1

coor lsqp norm mass sele plane1 end

set x1 = ?xaxi

set y1 = ?yaxi

set z1 = ?zaxi

set r1 = ?raxi

!now let's check if the vectors are in phase

calc dp = @cpx*@x1 + @cpy*@y1 + @cpz*@z1

if @dp .gt. 0 then goto skipinv1

set x1 =

-@x1set y1 =

-@y1set z1 =

-@z1label skipinv1

coor axis sele pc1 end sele pc2 end

set vx1 = ?xaxi

set vy1 = ?yaxi

set vz1 = ?zaxi

coor axis sele pc1 end sele pc3 end

set vx2 = ?xaxi

set vy2 = ?yaxi

set vz2 = ?zaxi

calc cpx = @vy1*@vz2 - @vy2*@vz1

calc cpy = - @vx1*@vz2 + @vx2*@vz1

calc cpz = @vx1*@vy2 - @vx2*@vy1

coor lsqp norm mass sele plane2 end

set x2 = ?xaxi

set y2 = ?yaxi

set z2 = ?zaxi

set r2 = ?raxi

!now let's check if the vectors are in phase

calc dp = @cpx*@x2 + @cpy*@y2 + @cpz*@z2

if @dp .gt. 0 then goto skipinv2

set x2 =

-@x2set y2 =

-@y2set z2 =

-@z2label skipinv2

calc angl = @x1*@x2 + @y1*@y2 + @z1*@z2

calc angl = @angl/@r1

calc angl = @angl/@r2

calc angl = acos(@angl)*180.0/@pi

calc angl2 = 180.0 - @angl

write title unit 22

* @f @angl @angl2

*

format

incr f by 1

if @f .le. @nframes then goto tloop

stop