From a706dda880c5a37b0dd3f812170e6c94a47f3338 Mon Sep 17 00:00:00 2001 From: "James W. Barnett" Date: Sat, 23 May 2015 12:56:36 -0500 Subject: [PATCH 1/2] Read in temperature from Gromacs files. Add a function to read in temperature from Gromacs dhdl.xvg files, so the user does not need to specify it directly. The temperature is read in from one of the files, whichever is first in the pre-sorted list (as they should be the same anyway). The temperature read in overrides the default temp and anything specified on the command line. If there is a problem reading in the temperature from dhdl.xvg (I can't see how at this point) it will fall back to the default or the one specified on the command line. --- alchemical_analysis.py | 5 +++++ parsers/parser_gromacs.py | 29 +++++++++++++++++++++++++---- 2 files changed, 30 insertions(+), 4 deletions(-) diff --git a/alchemical_analysis.py b/alchemical_analysis.py index 507cad6..74e60bb 100644 --- a/alchemical_analysis.py +++ b/alchemical_analysis.py @@ -101,6 +101,11 @@ def addRemove(string): def checkUnitsAndMore(units): + # If using Gromacs we can read in the temperature directly from dhdl.xvg + if P.software.title() == 'Gromacs': + import parsers.parser_gromacs + P.temperature = parsers.parser_gromacs.readTempGromacs(P) + kB = 1.3806488*6.02214129/1000.0 # Boltzmann's constant (kJ/mol/K). beta = 1./(kB*P.temperature) diff --git a/parsers/parser_gromacs.py b/parsers/parser_gromacs.py index 6160e94..82c55d1 100644 --- a/parsers/parser_gromacs.py +++ b/parsers/parser_gromacs.py @@ -6,10 +6,6 @@ import unixlike # some implemented unixlike commands -#=================================================================================================== -# FUNCTIONS: This is the Gromacs dhdl.xvg file parser. -#=================================================================================================== - def readDataGromacs(P): """Read in .xvg files; return nsnapshots, lv, dhdlt, and u_klt.""" @@ -280,3 +276,28 @@ def parseLog(self): f.iter_loadtxt(nf) return nsnapshots, lv, dhdlt, u_klt + +#=================================================================================================== +# FUNCTIONS: This reads in just the temperature from one of the dhdl.xvg files +#=================================================================================================== + +def readTempGromacs(P): + + datafile_tuple = P.datafile_directory, P.prefix, P.suffix + filename = glob( '%s/%s*%s' % datafile_tuple )[1] + + skip_lines = 0 # Number of lines from the top that are to be skipped. + + print "Reading temperature from %s..." % filename + with open(filename,'r') as infile: + for line in infile: + if line.startswith('@'): + elements = unixlike.trPy(line).split() + if 'subtitle' in elements: + temp = float(elements[4]) + print "Temperature is %g K." % temp + return temp + skip_lines += 1 + + print "Error reading in temperature. Using default." + return P.temperature From 7e27bda9b44644780458a8112f6539dd3def157b Mon Sep 17 00:00:00 2001 From: "James W. Barnett" Date: Tue, 26 May 2015 16:17:01 -0500 Subject: [PATCH 2/2] readTempGromacs: read in temps from all xvg files Update readTempGromacs to read in temperature from all .xvg files for the set of simulations. Now checks to ensure that the temperature is the same in every file, else stop the program. If the temperature is not present in all files, then use the default / command-line value. --- parsers/parser_gromacs.py | 51 +++++++++++++++++++++++++-------------- 1 file changed, 33 insertions(+), 18 deletions(-) diff --git a/parsers/parser_gromacs.py b/parsers/parser_gromacs.py index 2d78e83..cea4169 100644 --- a/parsers/parser_gromacs.py +++ b/parsers/parser_gromacs.py @@ -26,9 +26,11 @@ import unixlike # some implemented unixlike commands -#=================================================================================================== -# FUNCTIONS: This is the Gromacs dhdl.xvg file parser. -#=================================================================================================== +import sys + +#=================================================================================================== +# FUNCTIONS: This is the Gromacs dhdl.xvg file parser. +#=================================================================================================== def readDataGromacs(P): """Read in .xvg files; return nsnapshots, lv, dhdlt, and u_klt.""" @@ -302,26 +304,39 @@ def parseLog(self): return nsnapshots, lv, dhdlt, u_klt #=================================================================================================== -# FUNCTIONS: This reads in just the temperature from one of the dhdl.xvg files +# FUNCTIONS: This reads in just the temperature from the dhdl.xvg files #=================================================================================================== def readTempGromacs(P): datafile_tuple = P.datafile_directory, P.prefix, P.suffix - filename = glob( '%s/%s*%s' % datafile_tuple )[1] - skip_lines = 0 # Number of lines from the top that are to be skipped. + temperature = [] + print "Reading temperature from all .xvg files..." - print "Reading temperature from %s..." % filename - with open(filename,'r') as infile: - for line in infile: - if line.startswith('@'): - elements = unixlike.trPy(line).split() - if 'subtitle' in elements: - temp = float(elements[4]) - print "Temperature is %g K." % temp - return temp - skip_lines += 1 + # Cycles through all .xvg files (not currently sorted). Opens them, then + # gets the temperature from the file. + for filename in glob( '%s/%s*%s' % datafile_tuple ): + skip_lines = 0 # Number of lines from the top that are to be skipped. + with open(filename,'r') as infile: + for line in infile: + if line.startswith('@'): + elements = unixlike.trPy(line).split() + if "T" in elements: + temperature.append(float(elements[4])) + skip_lines += 1 + + # If there were no temperatures in any of the files, use the default + if len(temperature) == 0: + print "Temperature not present in xvg files. Using %g K." % P.temperature + return P.temperature + + # Check that all of the temperatures are the same. If not, script should + # not run. + if all(x==temperature[0] for x in temperature): + print "Temperature is %g K." % float(temperature[0]) + return temperature[0] + else: + print "ERROR: Temperature is not the same in all .xvg files." + sys.exit(0) - print "Error reading in temperature. Using default." - return P.temperature