[brahms-md] r68 committed - Improved KA tool
brah...@googlecode.com
Author: orsimario
Date: Sat Feb 18 04:38:01 2012
Log: Improved KA tool
http://code.google.com/p/brahms-md/source/detail?r=68
Modified:
/trunk/analysisScripts/general/computeAreaCompressibilityModulus.py
=======================================
--- /trunk/analysisScripts/general/computeAreaCompressibilityModulus.py Thu
Aug 25 09:03:28 2011
+++ /trunk/analysisScripts/general/computeAreaCompressibilityModulus.py Sat
Feb 18 04:38:01 2012
@@ -1,120 +1,106 @@
-## Copyright (C) 2006, 2007, 2008, 2009, 2010, 2011 Mario Orsi
-## This file is part of Brahms.
-## Brahms is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
-## Brahms is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
-## You should have received a copy of the GNU General Public License
along with Brahms. If not, see <http://www.gnu.org/licenses/>.
-
-# This program opens the file "area_volume.dat", which is one of BRAHMS's
output files, and extract the area compressibility ("stretch") modulus
-# USAGE: cat area_volume.dat | python computeAreaCompressibilityModulus.py
-
+#!/usr/bin/env python
+
+# Script: computeAreaCompressibilityModulus.py
+# Author: Mario Orsi (orsimario at gmail.com, www.soton.ac.uk/~orsi)
+# Purpose: Reads a file containing the time evolution of the interfacial
+# area of a lipid bilayer (usually the xy area) and calculates
+# the corresponding area compressibility (stretch) modulus
+# Syntax: computeAreaCompressibilityModulus.py inputFile temperature
+# Note: input area data are in nm^2
+# Example: computeAreaCompressibilityModulus.py xyArea.dat 303.15
+# Reference: - Orsi et al, J Phys Condens Matter 22: 155106 (2010),
+# section 5.1.2
+# - Waheed & Edholm, Biophys J 97: 2754 (2009)
+
import sys, string
from math import sqrt
-lines = sys.stdin.readlines()
-
-#constants
-BOLTZMANN = 1.3806505e-23 # [J / K]
-
-# variables
-T = 303.15 # [K]
-nLipids = 128 # number of lipid molecules in the system
-
-print "\n***REMEMBER TO ADJUST THE NUMBER OF LIPID AND TEMPERATURE!!!\n"
-
-nLines = 0 # total number of lines in the file
+
+kB = 1.3806505e-23 # Boltzmann constant [J / K]
+J_nm2__in__dyn_cm = 1e21 # 1 J/nm^2 = 10^21 dyn/cm
+
+if len(sys.argv) != 3:
+ print "Syntax: calcAreaComp.py inputFile temperature"
+ sys.exit()
+
+inFileName = sys.argv[1]
+T = float(sys.argv[2])
+print "\nT = %.2f K\n" % T
+inFile = open(inFileName, "r")
+lines = inFile.readlines()
+inFile.close()
+
lineCounter = 0
-meanBoxAreaSum = 0
-squareBoxAreaSum = 0
-meanLipidAreaSum = 0
-squareLipidAreaSum = 0
-
-meanBoxAreaSum1 = 0
-squareBoxAreaSum1 = 0
-meanBoxAreaSum2 = 0
-squareBoxAreaSum2 = 0
-
-meanLipidAreaSum1 = 0
-squareLipidAreaSum1 = 0
-meanLipidAreaSum2 = 0
-squareLipidAreaSum2 = 0
-
-brahmsAreaVolFile = open( 'area_volume.dat', 'r' )
+areaSum = squareAreaSum = 0
+areaSum1 = squareAreaSum1 = 0
+areaSum2 = squareAreaSum2 = 0
+
+nData = nData1 = nData2 = 0 # counter for number of data (measurements)
for line in lines:
- nLines = nLines + 1
+ if line[0] != '#': # ignore comments
+ nData += 1
for line in lines:
- lineCounter = lineCounter + 1
- words = string.split( line )
- boxArea = string.atof(words[0]);
- meanBoxAreaSum = meanBoxAreaSum + boxArea
- squareBoxAreaSum = squareBoxAreaSum + boxArea**2
- if lineCounter < nLines/2 : # first half of lines
- meanBoxAreaSum1 = meanBoxAreaSum1 + boxArea
- squareBoxAreaSum1 = squareBoxAreaSum1 + boxArea**2
- else: # second half of lines
- meanBoxAreaSum2 = meanBoxAreaSum2 + boxArea
- squareBoxAreaSum2 = squareBoxAreaSum2 + boxArea**2
-
-brahmsAreaVolFile.close
-
-# computing statistical quantities
-
-meanBoxArea = meanBoxAreaSum / nLines;
-print "meanBoxArea = %4.2f A^2" % (meanBoxArea * 100)
-meanSquaredBoxAreaFluctuation = squareBoxAreaSum / nLines - meanBoxArea**2
-# print "meanSquaredAreaFluctuation = %6.3f A^4" %
(meanSquaredAreaFluctuation * 10000)
-
-meanLipidArea = meanBoxArea / (nLipids/2);
-print "meanLipidArea = %4.2f A^2" % (meanLipidArea * 100)
-
-# computing modulus
-KA_boxArea = BOLTZMANN * T * meanBoxArea / meanSquaredBoxAreaFluctuation #
[ J / nm^2 ]
-
+ if line[0] != '#': # ignore comments
+ lineCounter = lineCounter + 1
+ words = string.split( line )
+ area = string.atof(words[0]);
+ areaSum += area
+ squareAreaSum += area**2
+ if lineCounter <= nData/2 : # first half of lines
+ nData1 += 1
+ areaSum1 += area
+ squareAreaSum1 += area**2
+ else: # second half of lines
+ nData2 += 1
+ areaSum2 += area
+ squareAreaSum2 += area**2
+
+# check:
+if nData != nData1+nData2:
+ print "ERROR: wrong scanning of data - check script"
+
+# Calc modulus considering all data:
+meanArea = areaSum / nData;
+print "\nmeanArea = %4.2f nm^2" % meanArea
+meanSquaredAreaFluct = squareAreaSum / nData - meanArea**2
+print "meanSquaredAreaFluct = %6.3f nm^4" % meanSquaredAreaFluct
+# computing modulus:
+KA = kB*T * meanArea / meanSquaredAreaFluct # [ J / nm^2 ]
+# convert J/nm^2 -> 10^21 dyn/cm:
+KA *= J_nm2__in__dyn_cm
+print "KA_total = %6.1f dyn/cm\n" % KA
+
+# Calc modulus considering only first half of data:
+meanArea1 = areaSum1 / nData1;
+print "meanArea1 = %4.2f nm^2" % meanArea1
+print "squareAreaSum1 / nData1= %4.2f nm^2" % (squareAreaSum1/ nData1)
+print "meanArea1**2 = %4.2f nm^2" % (meanArea1**2)
+meanSquaredAreaFluct1 = squareAreaSum1 / nData1 - meanArea1**2
+print "meanSquaredAreaFluct1 = %6.3f nm^4" % meanSquaredAreaFluct1
+# computing modulus:
+KA1 = kB*T * meanArea1 / meanSquaredAreaFluct1 # [ J / nm^2 ]
# conversion considering that J/nm^2 = 10^21 dyn/cm
-KA_boxArea = KA_boxArea * pow(10,21)
-
-print "KA_total = %6.1f dyn/cm\n" % ( KA_boxArea )
-
-meanBoxArea1 = meanBoxAreaSum1 / (nLines/2);
-print "meanBoxAreaFirstHalf = %4.2f A^2" % (meanBoxArea1 * 100)
-meanSquaredBoxAreaFluctuation1 = squareBoxAreaSum1 / (nLines/2) -
meanBoxArea1**2
-# print "meanSquaredAreaFluctuation = %6.3f A^4" %
(meanSquaredAreaFluctuation * 10000)
-
-meanLipidArea1 = meanBoxArea1 / (nLipids/2);
-print "meanLipidAreaFirstHalf = %4.2f A^2" % (meanLipidArea1 * 100)
-
-# computing modulus
-KA_boxArea1 = BOLTZMANN * T * meanBoxArea1 /
meanSquaredBoxAreaFluctuation1 # [ J / nm^2 ]
-
+KA1 *= J_nm2__in__dyn_cm
+print "KA_1 = %6.1f dyn/cm\n" % ( KA1 )
+
+# Calc modulus considering only second half of data:
+meanArea2 = areaSum2 / nData2;
+print "meanAreaSecond2 = %4.2f nm^2" % meanArea2
+print "squareAreaSum2/ nData2 = %4.2f nm^2" % (squareAreaSum2/ nData2)
+print "meanArea2**2 = %4.2f nm^2" % (meanArea2**2)
+meanSquaredAreaFluct2 = squareAreaSum2 / nData2 - meanArea2**2
+print "meanSquaredAreaFluct2 = %6.3f nm^4" % meanSquaredAreaFluct2
+# computing modulus:
+KA2 = kB*T * meanArea2 / meanSquaredAreaFluct2 # [ J / nm^2 ]
# conversion considering that J/nm^2 = 10^21 dyn/cm
-KA_boxArea1 = KA_boxArea1 * pow(10,21)
-
-print "KA_firstHalf = %6.1f dyn/cm\n" % ( KA_boxArea1 )
-
-meanBoxArea2 = meanBoxAreaSum2 / (nLines/2);
-print "meanBoxAreaSecondHalf = %4.2f A^2" % (meanBoxArea2 * 100)
-meanSquaredBoxAreaFluctuation2 = squareBoxAreaSum2 / (nLines/2) -
meanBoxArea2**2
-# print "meanSquaredAreaFluctuation = %6.3f A^4" %
(meanSquaredAreaFluctuation * 10000)
-
-meanLipidArea2 = meanBoxArea2 / (nLipids/2);
-print "meanLipidAreaSecondHalf = %4.2f A^2" % (meanLipidArea2 * 100)
-
-# computing modulus
-KA_boxArea2 = BOLTZMANN * T * meanBoxArea2 /
meanSquaredBoxAreaFluctuation2 # [ J / nm^2 ]
-
-# conversion considering that J/nm^2 = 10^21 dyn/cm
-KA_boxArea2 = KA_boxArea2 * pow(10,21)
-
-print "KA_secondHalf = %6.1f dyn/cm\n" % ( KA_boxArea2 )
-
-KA_avg12 = 0.5*(KA_boxArea1 + KA_boxArea2)
-
-print "\nKA_avg12 = %6.1f dyn/cm\n" % ( KA_avg12 )
-
-standardDeviation = sqrt( (KA_boxArea1-KA_avg12)**2 +
(KA_boxArea2-KA_avg12)**2 )
+KA2 *= J_nm2__in__dyn_cm
+print "KA_2 = %6.1f dyn/cm\n" % ( KA2 )
+
+# Calc average:
+KA_avg12 = 0.5 * ( KA1 + KA2 )
+print "KA_avg12 = %6.1f dyn/cm" % KA_avg12
+standardDeviation = sqrt( ( KA1 - KA_avg12 )**2 + ( KA2 - KA_avg12 )**2 )
standardError = standardDeviation / sqrt(2)
-
-print "Standard error = %6.1f dyn/cm" % (standardError)
+print "Standard error = %6.1f dyn/cm\n" % standardError