You are not logged in. [Log In] CHARMM Development Project » Forums » CHARMM Community » Script Archive » Is this possible in CHARMM (adding a routine) ? Register User    Forum List        Calendar     Active Topics    Search     FAQ
 Page 2 of 2 < 1 2
 Topic Options
 #7806 - 09/12/05 06:27 PM Angle between two least-square-planes with LSQP [Re: cola] ilja Forum Member Registered: 05/19/04 Posts: 34 Loc: La Jolla, California 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.14159265358979define plane1 sele ... enddefine plane2 sele ... enddefine nc1 sele atom1_from_plane1 enddefine nc2 sele atom2_from_plane1 enddefine nc3 sele atom3_from_plane1 enddefine pc1 sele atom1_from_plane2 enddefine pc2 sele atom2_from_plane2 enddefine pc3 sele atom3_from_plane3 endopen write card unit 22 name out_file.lsqp!Your standard way of reading trajectoriestraj query unit @istartcalc nframes = ?nfile * @pcstraj iread @istart nunits @pcs nfile @nframes skip 1set f 1label tlooptraj readformat (f16.8)coor axis sele nc1 end sele nc2 endset vx1 = ?xaxiset vy1 = ?yaxiset vz1 = ?zaxicoor axis sele nc1 end sele nc3 endset vx2 = ?xaxiset vy2 = ?yaxiset vz2 = ?zaxicalc cpx = @vy1*@vz2 - @vy2*@vz1calc cpy = - @vx1*@vz2 + @vx2*@vz1calc cpz = @vx1*@vy2 - @vx2*@vy1coor lsqp norm mass sele plane1 endset x1 = ?xaxiset y1 = ?yaxiset z1 = ?zaxiset r1 = ?raxi!now let's check if the vectors are in phasecalc dp = @cpx*@x1 + @cpy*@y1 + @cpz*@z1if @dp .gt. 0 then goto skipinv1set x1 = -@x1set y1 = -@y1set z1 = -@z1label skipinv1coor axis sele pc1 end sele pc2 endset vx1 = ?xaxiset vy1 = ?yaxiset vz1 = ?zaxicoor axis sele pc1 end sele pc3 endset vx2 = ?xaxiset vy2 = ?yaxiset vz2 = ?zaxicalc cpx = @vy1*@vz2 - @vy2*@vz1calc cpy = - @vx1*@vz2 + @vx2*@vz1calc cpz = @vx1*@vy2 - @vx2*@vy1coor lsqp norm mass sele plane2 endset x2 = ?xaxiset y2 = ?yaxiset z2 = ?zaxiset r2 = ?raxi!now let's check if the vectors are in phasecalc dp = @cpx*@x2 + @cpy*@y2 + @cpz*@z2if @dp .gt. 0 then goto skipinv2set x2 = -@x2set y2 = -@y2set z2 = -@z2label skipinv2calc angl = @x1*@x2 + @y1*@y2 + @z1*@z2calc angl = @angl/@r1calc angl = @angl/@r2calc angl = acos(@angl)*180.0/@picalc angl2 = 180.0 - @anglwrite title unit 22* @f @angl @angl2*formatincr f by 1if @f .le. @nframes then goto tloopstop Edited by ilja (09/13/05 12:13 PM) Top
 Page 2 of 2 < 1 2

 Hop to: User Discussion & Questions ------   Setup, I/O, and Basic questions   Energy terms, Constraints, Restraints, and Solvation   Parameter Set Discussion   Molecular Dynamics   Minimization, Normal modes, Monte Carlo,...   Installation and Testing   Performance and Benchmarking   Other User Discussion and Questions   General Chemistry DiscussionsCHARMM Interfaces ------   QM/MM Discussion and Questions   CHARMMing   CHARMM-GUI   MMTSB   AccelrysCHARMM Community ------   Script Archive   Bug Reports & Fixes   Meetings & Events   Position Openings   CHARMM Course   Suggestions and Requests   Forum News and Announcements
Moderator:  chmgr, John Legato, petrella