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 = -@x1
set y1 = -@y1
set z1 = -@z1

label 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 = -@x2
set y2 = -@y2
set z2 = -@z2

label 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


Edited by ilja (09/13/05 12:13 PM)