[brahms-md] r47 committed - Cosmetic changes
brah...@googlecode.com
Author: orsimario
Date: Tue Jun 28 02:17:46 2011
Log: Cosmetic changes
http://code.google.com/p/brahms-md/source/detail?r=47
Modified:
/trunk/analysisScripts/general/computeAreaCompressibilityModulus.py
=======================================
--- /trunk/analysisScripts/general/computeAreaCompressibilityModulus.py Thu
Jan 6 07:00:01 2011
+++ /trunk/analysisScripts/general/computeAreaCompressibilityModulus.py Tue
Jun 28 02:17:46 2011
@@ -1,8 +1,12 @@
-# M. Orsi - Essex group - October 2006
-# modified December 2006, to consider total area + volume
-# modified december 2008, to deal with only two columns (area + volume) -
got rid of first column (time) to save on disk space
-# This program opens the file "areaVol.dat", which is one of BRAHMS's
output files, and extract the area compressibility (AKA stretch) modulus
-
+## 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
import sys, string
@@ -16,30 +20,51 @@
T = 303.15 # [K]
nLipids = 128 # number of lipid molecules in the system
-print "REMEMBER TO ADJUST THE NUMBER OF LIPID AND TEMPERATURE!!!"
-
+print "\n***REMEMBER TO ADJUST THE NUMBER OF LIPID AND TEMPERATURE!!!\n"
+
+nLines = 0 # total number of lines in the file
lineCounter = 0
-meanLipidAreaSum = 0
-squareLipidAreaSum = 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' )
-for line in lines: # only consider boxArea!
+for line in lines:
+ nLines = nLines + 1
+
+for line in lines:
lineCounter = lineCounter + 1
words = string.split( line )
boxArea = string.atof(words[0]);
- meanBoxAreaSum = meanBoxAreaSum + boxArea;
+ 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 / lineCounter;
+meanBoxArea = meanBoxAreaSum / nLines;
print "meanBoxArea = %4.2f A^2" % (meanBoxArea * 100)
-meanSquaredBoxAreaFluctuation = squareBoxAreaSum / lineCounter -
meanBoxArea**2
+meanSquaredBoxAreaFluctuation = squareBoxAreaSum / nLines - meanBoxArea**2
# print "meanSquaredAreaFluctuation = %6.3f A^4" %
(meanSquaredAreaFluctuation * 10000)
meanLipidArea = meanBoxArea / (nLipids/2);
@@ -51,4 +76,45 @@
# conversion considering that J/nm^2 = 10^21 dyn/cm
KA_boxArea = KA_boxArea * pow(10,21)
-print "KA = %6.1f dyn/cm" % ( KA_boxArea )
+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 ]
+
+# 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 "meanBoxAreaFirstHalf = %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 "meanLipidAreaFirstHalf = %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_firstHalf = %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 )
+standardError = standardDeviation / sqrt(2)
+
+print "Standard error = %6.1f dyn/cm" % (standardError)