VMD-L Mailing List

From: Andrés Morales (h.andres.m1986_at_gmail.com)
Date: Tue Apr 10 2012 - 11:43:10 CDT

Dear VMD users:

I want to calculate the average angle between the lipid dipole vector (P-N
dipole vectors in DPPC headgroups ) and the bilayer normal for each frame.
 So I use the following script:

proc angle_ave {sel_1 sel_2 frame} {

global M_PI
set conv [expr 180.0/$M_PI]
set num1 [$sel_1 num]
set sum 0
set mol1 [$sel_1 molindex]
set mol2 [$sel_2 molindex]

set sel1 [atomselect $mol1 [$sel_1 text] frame $frame]
set sel2 [atomselect $mol2 [$sel_2 text] frame $frame]

foreach atom1 [$sel1 get index] atom2 [$sel2 get index] {
set pos1 [lindex [[atomselect top "index $atom1"] get {x y z}] 0]
set pos2 [lindex [[atomselect top "index $atom2"] get {x y z}] 0]
set v [vecsub $pos2 $pos1]
set z1 [lindex $v 2]
                set theta [expr $conv * asin($z1 / [veclength $v] )]
set tilt [expr 90 - $theta ]
  set sum [expr $sum + $tilt]
}
set prom [expr $sum / $num1]
return $prom
}

set outfile [open tilt angles.txt w];

puts $outfile "Frame Average tilt angle"
set dppc [mol load psf final.psf dcd final.dcd]
set nf [molinfo top get numframes]
set sel1 [atomselect $dppc "lipids and resid 1 to 36 and name P " ]
set sel2 [atomselect $dppc "lipids and resid 1 to 36 and name N " ]
for {set i 0 } {$i < $nf } { incr i } {
 $sel1 update
 $sel2 update
 set average [angle_ave $sel1 $sel2 $i]
puts $outfile "$i $average"
}

close $outfile

It calculates the average angle quite fine, but it only do it for the first
frame. The output looks like:

Frame Average tilt angle
0 72.3545824395983
1 72.3545824395983
2 72.3545824395983
... ....

Any suggestions?

Thanks a lot.

Hernán Andrés Morales Navarrete

Biophysics and Molecular Modeling Group
Physics Department
Escuela Politécnica Nacional, Quito - Ecuador
Ladrón de Guevara E11-253.
Casilla 17-01-1253
http://www.ciencias.epn.edu.ec/~biomod/