VMD-L Mailing List

From: Justin Gullingsrud (justin_at_ks.uiuc.edu)
Date: Wed Jun 04 2003 - 12:29:43 CDT

On Wed, Jun 04, 2003 at 06:37:41PM +0200, Jan Saam wrote:
> Hi Thomas,
> this script does the job.

True, although the dihedral calculation done in the label command will be
dozens of times faster, mainly because it doesn't have to create atom
selections to fetch the coordinates.

Here's a little script in the same format as Jan's script that does the
same thing using the label command. It uses the current timestep by
default, instead of the first timestep. It also blows away any existing
Dihedral labels; this is because of a deficiency in the scripting
interface.

proc my_dihed_angle { a1 a2 a3 a4 { frame now }} {
  # Delete all existing Dihdral labels so that the one we add has index 0.
  label delete Dihedrals

  # Use the top molecule
  set molid [molinfo top]

  # Override the current frame
  set curframe [molinfo top get frame]
  molinfo top set frame $frame

  # Add the label
  label add Dihedrals $molid/$a1 $molid/$a2 $molid/$a3 $molid/$a4

  # Get its value
  set curval [lindex [lindex [label list Dihedrals] 0] 4]

  # Restore the old frame
  molinfo top set frame $curframe

  return $curval
}

Note that if you just want the value of the dihedral for all timesteps,
there's a much faster way to do it: the "label graph" command. Here's
a script that returns the values for all timesteps:

proc all_dihed_angle { a1 a2 a3 a4 } {
  # Delete all existing Dihdral labels so that the one we add has index 0.
  label delete Dihedrals

  # Use the top molecule
  set molid [molinfo top]

  # Add the dihedral monitor
  label add Dihedrals $molid/$a1 $molid/$a2 $molid/$a3 $molid/$a4

  # Get all values
  return [label graph Dihedrals 0]
}

Cheers,
Justin

>
> Jan
>
> # Computes the dihedral angle of a bond
> # (takes the indices of the four dihedral atoms as input)
> proc dihed_angle { a1 a2 a3 a4 {frame 0}} {
> set coord1 [lindex [[atomselect top "index $a1" frame $frame] get {x y z}] 0]
> set coord2 [lindex [[atomselect top "index $a2" frame $frame] get {x y z}] 0]
> set coord3 [lindex [[atomselect top "index $a3" frame $frame] get {x y z}] 0]
> set coord4 [lindex [[atomselect top "index $a4" frame $frame] get {x y z}] 0]
> set v1 [vecsub $coord1 $coord2]
> set v2 [vecsub $coord3 $coord2]
> set v3 [vecsub $coord4 $coord3]
>
> set cross1 [vecnorm [veccross $v2 $v1]]
> set cross2 [vecnorm [veccross $v2 $v3]]
> set dot [vecdot $cross1 $cross2]
> if {$dot>1.0 && $dot<1.0001} {set dot 1.0}
> if {$dot<-1.0 && $dot>-1.0001} {set dot -1.0}
> if {$dot>1.0 || $dot<-1.0} {
> puts "dihed_angle: dot>1.0, cannot compute acos($dot)"
> }
> set angle [rad2deg [expr acos($dot)]]
> return $angle
>
>
> Am Mittwoch, 4. Juni 2003 17:58 schrieb Thomas Hedegaard Pedersen:
> > Hi,
> >
> > Does anyone know if there is an easy way to obtain dihedral angles
> > from within a tcl-script. I'm aware that for a protein it is easy to
> > obtain the phi and psi angles by using the: get phi (or get psi),
> > calls, but are there someway that the dihedral for a given selection
> > of 4 atoms can be obtained (by a simple function call)????
> >
> > In addition, is it possible to rotate some parts of a molecule by
> > specifying a new dihedral, as it is possible to use the move command
> > to move atomselections by some given distance x, y, z????
> >
> > Thank you
> >
> > Stud. Polyt.
> > Thomas Hedegaard Pedersen
> > Department of Chemistry
> > Technical University of Denmark
> 

-- 
  Justin Gullingsrud        3111 Beckman Institute        217-244-8946
  I been dropping the new science, and I be kicking the new knowledge,
  and I'm seeing to a degree that you can't get in college.  -- b.boys