From: Andrés Morales (h.andres.m1986_at_gmail.com)
Date: Thu Oct 20 2011 - 01:12:14 CDT

Dear VMD users:

 Is there anyone here to know how to calculate the area per lipid and
thickness for a lipid bilayer in a membrane system?
 Or, could anyone share the scripts to calculate them?

I was looking at VMD and NAMD mailing lists, and I found some suggestions to
calculate a crude aproximation for area per lipid and thickness of a lipid
bilayer .
Acording with those suggestions I wrote some simple scripts to calculate
them. But I have some problems.

 First, the area per lipid calculated for each frame was around ~130 A^2 per
lipid. It is a too large value since I am workin with a DPPC bilayer (~62
A^2 per lipid).
Second, the thickness calcualed was good for the 50 first frames, but for
the other ones, it took too large values (~1000 A).

I calculated :
 Average area per lipid: "box size X * box size Y / no. of lipids" (lipids
are in the XY plane)
 Thickness: I brought the membrane center of mass at origin (0,0,0),
calculated average positive Z coordinates and average negative Z
coordinates, substracted the negative from the positive
average.
The scripts I use was (lipids 1 to 36 formed the monolayer above z=0 and
lipids 37 to 72 formed the monolayer beneath z=0 ):
**
*
###################Protocol to calculate average Z coordinates
##########################
proc prom_z {sel} {
 set sumz 0
 foreach atom [$sel get index] {
  set pos [lindex [[atomselect top "index $atom"] get {x y z}] 0]
  set z1 [lindex $pos 2]
   set sumz [expr abs($sumz )+ abs($z1)]
 }
 set promz [expr $sumz / 36]
 return $promz
}
##############################################################################################

set mol [mol new dppc.psf type psf]
mol addfile dppc_eq.dcd type dcd waitfor all

set outfile [open Analysis_DPPC.txt
w];
puts $outfile "Frame Area per lipid Thickness (N) Thickness (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 element N"]
set sel3 [atomselect top "lipids and resid 37 to 72 and element N"]
set sel4 [atomselect top "lipids and resid 1 to 36 and element P"]
set sel5 [atomselect top "lipids and resid 37 to 72 and element P"]

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

 $sel1 moveby [vecinvert [measure center $sel weight mass]]

 ## Area per lipid
 set min [lindex [measure minmax $sel] 0]
 set max [lindex [measure minmax $sel] 1]
 set minx "[lindex $min 0]"
 set miny "[lindex $min 1]"
 set maxx "[lindex $max 0]"
 set maxy "[lindex $max 1]"
 set lenght_x [expr $maxx - $minx]
 set lenght_y [expr $maxy - $miny]
 set area_por_lipido [expr ($lenght_x * $lenght_x) / 36]

 ## Thicknes
set z1 [prom_z $sel2]
 set z2 [prom_z $sel3]
set thinckness1 [expr abs($z1) + abs($z2)]

 set z3 [prom_z $sel4]
 set z4 [prom_z $sel5]
set thinckness2 [expr abs($z3) + abs($z4)]

puts $outfile "$i $area_por_lipido $thinckness1 $thinckness1"
}
close $outfile

*

Thanks for your suggestions

 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/