From: Andrés Morales (h.andres.m1986_at_gmail.com)
Date: Thu Apr 19 2012 - 20:58:35 CDT

Dear Balaji
The script you mention does not work well, so I made another one that has
worked quite fine for me.

proc prom_z {selection frame} {

set mol [$selection molindex]
set sel1 [atomselect $mol [$selection text] frame $frame]
set coords1 [$sel1 get {x y z}]
set sumz 0
foreach coord1 $coords1 {
  set z1 [lindex $coord1 2]
   set sumz [expr $sumz + $z1]
 }
return [expr $sumz / ([$selection num] + 0.0)]

}
set outfile [open thickness.txt
w];
puts $outfile "Frame thickness"

set nf [molinfo top get numframes]
set sel1 [atomselect top "lipids and resid 1 to 36" ]
set sel2 [atomselect top "lipids and resid 37 to 72"]

for {set i 0 } {$i < $nf } { incr i } {

set z1 [prom_z $sel1 $i]
set z2 [prom_z $sel2 $i]
set thickness [expr $z1 - $z2]

puts $outfile "$i $thickness"
}
close $outfil
I wait it is useful for you. If you have any question let me know.

Kind regards

Andres

2012/4/16 Ban Arn <ban.arn_at_gmail.com>

> Dear Andreas
>
> I got your script in vmd mailing list for calculating thickness of bilayer.
>
> set outfile [open Espesor.txt w];
> puts $outfile "Frame Espesor(P)"
>
> set nf [molinfo top get numframes]
> set sel [atomselect top "lipids"]
> set sel1 [atomselect top all]
> set sel2 [atomselect top "lipids and resid 1 to 36 and name P"]
> set sel3 [atomselect top "lipids and resid 37 to 72 and name P"]
>
> for {set i 0 } {$i < $nf } { incr i } {
> $sel frame $i
> $sel1 frame $i
> $sel2 frame $i
> $sel3 frame $i
>
> $sel1 moveby [vecinvert [measure center $sel weight mass]]
> set sumz 0
> foreach atom [$sel2 get index] {
> set pos [lindex [[atomselect top "index $atom"] get {x y z}] 0]
> set z1 [lindex $pos 0]
> set sumz [expr abs($sumz)+ abs($z1)]
> }
> set promz1 [expr $sumz / 36]
>
> set sumz2 0
> foreach atom [$sel3 get index] {
> set pos2 [lindex [[atomselect top "index $atom"] get {x y z}] 0]
> set z2 [lindex $pos2 0]
> set sumz2 [expr abs($sumz2)+ abs($z2)]
> }
> set promz2 [expr $sumz2 / 36]
>
> set thinckness1 [expr abs($promz1) + abs($promz2)]
> puts $outfile "$i $thinckness1"
> }
> close $outfile
>
> However, the script works for 50 frames and later shows area more than
> ~1000A.
>
> I tried with modifying the lines as mentioned, it doesn't work for me.
>
> Kindly advice.
>
> Many Thanks
> Balaji
>