4
\$\begingroup\$

I am a computational chemist working with the program Gaussian 09. After I manually check the output(s) I want to create a summary for easier processing of the obtained values. Also avoid opening all the files again and again. The script searches for the last line of the energy statement. A portion of the output will be at the end of the post.

The outputs can become quite long and I am not completely satisfied with the performance of the script. I does its job, but if there are many big files, I can get a coffee in between. I do know that it is still faster than doing it by hand, but if it could be improved I would be quite happy.

I am pretty certain the problem comes from actually finding the string, i.e. the grep command. I am using tac here, since there could be multiple occurrences, but I am only interested in the last one. I have tried some of these solutions, too, but the tac|grep was the fastest. Depending on the steps of the optimisation it is therefore easier to read from the back. Since I already checked the file I also know the last value is the one I want.
Unfortunately this is slow if it also computes properties, where there is a huge block of properties at the end of the file. I don't know how I could skip it.

#!/bin/bash

# Find energy statement from a Gaussian 09 calculation
# Find energy statement from all G09 log files in working directory

findEnergy ()
{
    # Initiate variables necessary for parsing output
    local readWholeLine pattern functional energy cycles
    # Find match from the end of the file
    # Ref: https://unix.stackexchange.com/q/112159/160000
    readWholeLine=$(tac $1 | grep -m1 'SCF Done')
    # Gaussian output has following format, trap important information
    pattern="(E\(.+\)) = (.+) A\.U\.[^0-9]+([0-9]+) cycles"
    if [[ $readWholeLine =~ $pattern ]]
    then
        functional="${BASH_REMATCH[1]}"
        energy="${BASH_REMATCH[2]}"
        cycles="${BASH_REMATCH[3]}"
    fi

    # Print the line, format it for table like structure
    printf "%-25s %-15s = %20s ( %6s )\n" ${1%.*} $functional $energy $cycles

}

getAll ()
{
    # run over all commandfiles
    # ToDo: specify file suffixes
    local commandfile logfile
    printf "%-25s %s\n" "Summary for " ${PWD#\/*\/*\/}
    printf "%-25s %s\n\n" "Created " "$(date +"%Y/%m/%d %k:%M:%S")"
    # Print a header
    printf "%-25s %-15s   %20s ( %6s )\n" "Command file" "Functional" "Energy / Hartree" "cycles"
    for commandfile in *com; do
        logfile=${commandfile/com/log}
        if [[ -e $logfile ]]; then
            findEnergy $logfile
        else
            printf "%-25s No log file found.\n" ${commandfile%.*}
        fi
    done
}

getOnly ()
{
    # run only for commandline arguments
    local logfile=${1/com/log}
    if [[ -e $logfile ]]; then
        findEnergy $logfile
    else
        printf "%-25s No log file found.\n" ${1%.*}
    fi
}

if [[ -z $1 ]]; then getAll; fi
while [[ ! -z $1 ]]; do
    getOnly $1
    shift
done

A cut down sample of a Gaussian09, the relevant line is the one containing SCF Done.

Entering Gaussian System, Link 0=g09
 Input=bp86tzvpp-head.gjf
 Output=bp86tzvpp-head.log
 Initial command:
 /home/g09d01/g09/l1.exe "/scr/martin/g09job4401/Gau-29919.inp" -scrdir="/scr/martin/g09job4401/"
 Entering Link 1 = /home/g09d01/g09/l1.exe PID=     29920.

 Copyright (c) 1988,1990,1992,1993,1995,1998,2003,2009,2013,
            Gaussian, Inc.  All Rights Reserved.

 This is part of the Gaussian(R) 09 program.  It is based on
 the Gaussian(R) 03 system (copyright 2003, Gaussian, Inc.),
 the Gaussian(R) 98 system (copyright 1998, Gaussian, Inc.),
 the Gaussian(R) 94 system (copyright 1995, Gaussian, Inc.),
 the Gaussian 92(TM) system (copyright 1992, Gaussian, Inc.),
 the Gaussian 90(TM) system (copyright 1990, Gaussian, Inc.),
 the Gaussian 88(TM) system (copyright 1988, Gaussian, Inc.),
 the Gaussian 86(TM) system (copyright 1986, Carnegie Mellon
 University), and the Gaussian 82(TM) system (copyright 1983,
 Carnegie Mellon University). Gaussian is a federally registered
 trademark of Gaussian, Inc.

 This software contains proprietary and confidential information,
 including trade secrets, belonging to Gaussian, Inc.

 This software is provided under written license and may be
 used, copied, transmitted, or stored only in accord with that
 written license.

 The following legend is applicable only to US Government
 contracts under FAR:

                    RESTRICTED RIGHTS LEGEND

[...]
[...]

 Coeff:      0.166D-02 0.409D-02 0.195D-01 0.543D-01 0.157D+00 0.287D+00
 Coeff:      0.477D+00
 Gap=     0.077 Goal=   None    Shift=    0.000
 RMSDP=2.09D-09 MaxDP=3.55D-07 DE=-1.08D-08 OVMax= 1.13D-06

 SCF Done:  E(RB-P86) =  -4629.25655340     A.U. after   19 cycles
            NFock= 19  Conv=0.21D-08     -V/T= 2.0040
 KE= 4.610985852390D+03 PE=-3.763151266964D+04 EE= 1.501854952249D+04
 Leave Link  502 at Mon May  2 02:32:29 2016, MaxMem=   209715200 cpu:    175765.9
 (Enter /home/g09d01/g09/l508.exe)
 QCSCF skips out because SCF is already converged.
 Leave Link  508 at Mon May  2 02:32:29 2016, MaxMem=   209715200 cpu:         0.0
 (Enter /home/g09d01/g09/l601.exe)
 Copying SCF densities to generalized density rwf, IOpCl= 0 IROHF=0.

[...]
[...]

 IF YOU'RE NOT PART OF THE SOLUTION,
 THEN YOU'RE PART OF THE PRECIPITATE.
 Job cpu time:       2 days  9 hours 56 minutes 32.0 seconds.
 File lengths (MBytes):  RWF=   5328 Int=      0 D2E=      0 Chk=    359 Scr=      1
 Normal termination of Gaussian 09 at Mon May  2 05:26:15 2016.
\$\endgroup\$
  • 2
    \$\begingroup\$ IF YOU'RE NOT PART OF THE SOLUTION, THEN YOU'RE PART OF THE PRECIPITATE. =) \$\endgroup\$ – I'll add comments tomorrow Jun 1 '16 at 14:34
4
\$\begingroup\$

I don't see obvious signs why this script should be slow. There are no unnecessary sub-processes, no compute-heavy operations or nested loops, and so I don't know how to help you speed this up. The tac | grep -m1 combo are well utilized for their intended purpose, I don't think a specialized custom implementation for your purpose would make a significant difference.

I have a few tips only in terms of technique. In getAll, part of the code is identical to what you have in getOnly. As such, you can reuse that, and call getOnly from getAll:

getAll() {
    # run over all commandfiles
    # ToDo: specify file suffixes
    local commandfile logfile
    printf "%-25s %s\n" "Summary for " ${PWD#\/*\/*\/}
    printf "%-25s %s\n\n" "Created " "$(date +"%Y/%m/%d %k:%M:%S")"
    # Print a header
    printf "%-25s %-15s   %20s ( %6s )\n" "Command file" "Functional" "Energy / Hartree" "cycles"
    for commandfile in *com; do
        getOnly "$commandfile"
    done
}

This code is a bit sloppy, because if $1 is empty, then the while doesn't need to be executed:

if [[ -z $1 ]]; then getAll; fi
while [[ ! -z $1 ]]; do
    getOnly $1
    shift
done

It would be more appropriate to move the loop inside the else branch:

if [[ -z $1 ]]; then
    getAll
else
    while [[ ! -z $1 ]]; do
        getOnly $1
        shift
    done
fi

I'm guessing that the intention here is to call getAll if there are no arguments. But strictly speaking [[ -z $1 ]] doesn't mean there are no arguments, it just means that the first argument is empty. The correct way to check that there are no arguments:

if [[ $# == 0 ]]; then

And instead of a while loop, it would be more natural to use a for loop here:

    for commandfile in "$@"; do
        getOnly "$commandfile"
    done

Lastly, at many places you did not quote variables that are paths. I'm guessing you did that because you are certain they will never contain spaces. Even so, it's a good habit to double-quote such variables.

\$\endgroup\$
4
\$\begingroup\$

The performance issue hinges on

readWholeLine=$(tac $1 | grep -m1 'SCF Done')

The only suggestion I have would be to try $(grep 'SCF Done' | tail -n 1) instead, so that it doesn't have to seek backwards. If, as you say, you have already tried that alternative, I don't have much else to suggest. I doubt that even reimplementing it in another language would help much, since GNU grep is quite fast.


On the whole, this Bash script is well organized. You've split the work into functions and localized the variables.

For readability, I suggest adopting a convention of immediately unpacking the parameters at the top of every function. It's a pain to puzzle out what the functions' inputs are by scanning the entire code for a hidden $1, possibly disguised as $(1%.*}, and possibly a $2 as well. So, do this consistently:

findEnergy() {
    local logfile="$1"
    local readWholeLine pattern functional energy cycles
    …
}

In getOnly, this string substitution looks sloppy:

for commandfile in *com; do
    local logfile=${commandfile/com/log}
    …
done

I infer from ${commandfile%.*} that you expect it to end in .com. But what if the input filename is computation-20160603.com? You would end up with logputation-20160603.com. You should write ${1/%.com/.log} to ensure that the substitution only occurs for the extension.

As @janos has already pointed out, double-quoting variables is an important habit in shell scripts. Otherwise, your scripts can break in nasty ways just because an unexpected filename was passed in.

\$\endgroup\$
  • \$\begingroup\$ Thank you very much for these helpful comments. I have indeed never thought about that a file name would contain com except for the ending. | The parsing is certainly what takes longest, the problem is, that I don't know if it is an optimisation or not before parsing the file. I guess since there is like 50/50 of the files I handle it'll average out over the long time use. Experience will tell. And after all we are talking seconds, not hours. When I have gathered more information I'll post the updated code. \$\endgroup\$ – Martin - マーチン Jun 3 '16 at 11:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.