VMD-L Mailing List

From: Andrés Morales (h.andres.m1986_at_gmail.com)
Date: Mon Apr 09 2012 - 12:55:17 CDT

---------- Mensaje reenviado ----------
De: Andrés Morales <h.andres.m1986_at_gmail.com>
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. So I use
the following script (this is the correct one):

proc angle_ave {sel1 sel2} {

global M_PI
set conv [expr 180.0/$M_PI]
set num1 [$sel1 num]
set sum 0
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 abs($sum)+ abs($theta)]
return $sum
}
set prom [expr $sum / $num1]
return $prom
}

set sel1 [atomselect top "lipids and resid 1 to 36 and name P "]
set sel2 [atomselect top "lipids and resid 1 to 36 and name N "]
set average [angle_ave $sel1 $sel2]

But it only returns a wrong value for the average tilt angle.

Any suggestions?

Thanks a lot

Andres