From: Ban Arn (ban.arn_at_gmail.com)
Date: Fri Apr 20 2012 - 05:30:38 CDT

Dear Andres

Thanks for the script.
I ran the the script by modifying the number of lipid molecules, since i
have 154 lipid molecules.
I got output as

Frame thickness

0 -0.30396162335781018
1 -0.1713354174906333
2 -0.30750430435066173
3 -0.26556339701586062
4 -0.34344032037145983
5 -0.25293993409356585
6 -0.14722334875272414
7 -0.14722335044838467
8 -0.1021871230540968
9 -0.21873117335423042
10 -0.227086795466171
11 -0.17481152609827177
12 -0.19106022334606568
13 -0.18053145265450471
14 -0.18053142557786195
15 -0.13761933177174124
16 0.056644464849847725
17 0.023749397786653331
18 -0.026576366662544731
19 -0.0044238167118522687
20 -0.060259717657212247

My system has POPC bilayer and i'm getting very low values in the output.
In the vmd terminal I have loaded only lipids pdb and trajectory.
Kindly advice.

Many Thanks
Balaji

2012/4/20 Andrés Morales <h.andres.m1986_at_gmail.com>

>
> 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
>>
>
>