10
\$\begingroup\$

I've been writing bash script on and off, with pretty good results in terms of getting the job done. However, I'm worried that my script might be very ugly, as I am a beginner. I'm looking for advice concerning this particular one.

I want to extract some part of a big text file (here a lammps .lmps file), and manipulate it to create another file (here .xyz).

The interesting parts of the big file are:

The beginning of the file:

300 atoms
300 bonds
450 angles
600 dihedrals
150 impropers

The part about atom tag (here either "1" or "2") and mass:

Masses

1         12.011150
2          1.007970

Pair Coeffs

The part concerning atom coordinates:

Atoms

1       1    1     -0.126800     20.511864     28.359121     11.290877
2       1    1     -0.126800     21.779636     28.644716     10.779171
3       1    1     -0.126800     20.381316     27.822484     12.573717
4       1    1     -0.126800     21.518471     27.571445     13.344853
5       1    1     -0.126800     22.786244     27.857074     12.833161
6       1    1     -0.126800     22.916794     28.393694     11.550321
7       1    2      0.126800     19.390282     27.599170     12.973874
8       1    2      0.126800     19.622826     28.555315     10.688110
9       1    2      0.126800     23.907808     28.617021     11.150121
10      1    2      0.126800     21.881943     29.064261      9.776262
11      1    2      0.126800     23.675251     27.660865     13.435963
12      1    2      0.126800     21.416213     27.151893     14.347761

Bonds

Now here's the script I wrote, with some comments to explain what I wanted to do:

#!/bin/bash/

echo "LAMMPS file name? (without .lmps)"
read -r filename

# Make sure file exists
if [ -r "$filename".lmps ]; then

# Appends number of atoms to xyz file
head "$filename".lmps | grep atoms | awk -F' ' '{print $1}' > "$filename".xyz
echo >> "$filename".xyz

# Extracts "atom coordinates" and "masses" sections from original lmps file
awk '/Atoms/,/Bonds/' "$filename".lmps | head -n -2 | tail -n +3 > coordinates.tmp
awk '/Masses/,/Pair Coeffs/' "$filename".lmps | head -n -2 | tail -n +3 > masses.tmp

# Iterates over all lines of atom coordinates
while read line_atoms
do
    # Iterates over all lines of masses
    while read line_masses
    do
            mass=`echo $line_masses | awk -F' ' '{print $2}'`
            tag=`echo $line_masses | awk -F' ' '{print $1}'`

            case $mass in
                    # If mass corresponds to carbon and current line has lammps "carbon" tag
                    # then append this line to the xyz file 
                    12.011150)
                    if [ "$tag" == `echo $line_atoms | awk -F' ' '{print $3}'` ]; then
                    echo -e "C\t`echo $line_atoms | awk -F' ' '{print $5,"\t",$6,"\t",$7}'`" >> "$filename".xyz
                    fi
                    ;;

                    # Same, checks for Hydrogen mass and lammps "hydrogen" tag
                    1.007970)
                    if [ "$tag" == `echo $line_atoms | awk -F' ' '{print $3}'` ]; then
                    echo -e "H\t`echo $line_atoms | awk -F' ' '{print $5,"\t",$6,"\t",$7}'`" >> "$filename".xyz
                    fi
                    ;;
            esac
    done < masses.tmp

done < coordinates.tmp

# In case the requested file doesn't exist
else echo "Error:"$filename".lmps doesn't exits"

fi

# Gets rid of temporary files
rm *.tmp

The script works properly, and yields something like:

300

C       20.511864        28.359121       11.290877
C       21.779636        28.644716       10.779171
C       20.381316        27.822484       12.573717
C       21.518471        27.571445       13.344853
C       22.786244        27.857074       12.833161
C       22.916794        28.393694       11.550321
H       19.390282        27.599170       12.973874
H       19.622826        28.555315       10.688110
H       23.907808        28.617021       11.150121
H       21.881943        29.064261       9.776262
H       23.675251        27.660865       13.435963
H       21.416213        27.151893       14.347761

I feel like it is very inefficient (the double while read line, the commands inside if [ ] conditions, and so on) so I'd very much like some input in order to improve my coding skills. Feel free to tell me if it looks good or awful to you, and why!

\$\endgroup\$
0

3 Answers 3

7
\$\begingroup\$

General remarks

The first line should be #!/bin/bash.

Don't use `cmd` style command substitution, it's deprecated. Use $(cmd) instead, it's better.

Tests like [ "$var" = "$other" ] are also deprecated. Use [[ $var = $other ]] instead. As you can see in this example, you can omit the double-quotes in this modern version.

The http://www.shellcheck.net/ site is great for checking your code for common mistakes.

If there's awk in the pipeline, use it well

When you have a pipeline like this:

head "$filename".lmps | grep atoms | awk -F' ' '{print $1}' > "$filename".xyz

That is, there is an awk in the pipeline, along with other operations that awk could do all by itself. This code is equivalent:

awk -F' ' '/atoms/ {print $1} NR == 10 { exit }' "$filename".lmps > "$filename".xyz

This is better, because instead of 3 processes (head + grep + awk), you've managed to do everything in just one process.

Note: some of your other pipelines with awk are not well-suited for this, for example:

awk '/Atoms/,/Bonds/' "$filename".lmps | head -n -2 | tail -n +3 > coordinates.tmp

This is different from the first case, because there's no easy way to do with awk the equivalent of head -n -2. Also, moving the tail -n +3 logic inside the awk would be possible, but in this example too complicated, so it's ok to leave this statement as it is. It's only executed once per run, so using 3 processes instead of 2 is not a big problem.

Reading multiple variables from a line

You can simplify this:

while read line_masses
do
        mass=`echo $line_masses | awk -F' ' '{print $2}'`
        tag=`echo $line_masses | awk -F' ' '{print $1}'`
        # ...

by writing like this:

while read tag mass
do
        # ...

This is much better, as you just got rid of 2 extra processes per iteration.

You can do similarly for the outer loop as well:

while read f1 f2 atag f4 f5 f6 f7

This will simplify your if statments in the case $mass in, like this:

12.011150)
if [[ $tag = $atag ]]; then
    echo -e "C\t$f5\t$f6\t$f7" >> "$filename".xyz
fi
;;

1.007970)
if [[ $tag = $atag ]]; then
    echo -e "H\t$f5\t$f6\t$f7" >> "$filename".xyz
fi
;;

Calculate once, save in variable to reuse

Be careful with code like this:

while read line_atoms
do
    while read line_masses
    do
        if [ "$tag" == `echo $line_atoms | awk -F' ' '{print $3}'` ]; then
            echo -e "C\t`echo $line_atoms | awk -F' ' '{print $5,"\t",$6,"\t",$7}'`"

A big problem here is repeated evaluation of those echo $line_atoms | awk commands for each mass line in the input, when it would have been more efficient to calculate these before starting the inner loop.

Reduce nesting

The main part of the script is wrapped inside this large if block:

if [ -r "$filename".lmps ]; then
    # do the main work
fi

It would be better to reverse this logic, like this:

if [ ! -r "$filename".lmps ]; then
    echo "Error:"$filename".lmps doesn't exits"
    exit 1
fi

# do the main work

Related to this, it's a good practice to exit 1 to indicate an error to the caller.

Finally, you did some cleanup at the end of the script:

# Gets rid of temporary files
rm *.tmp

However, this is pointless if the input file did not exist. One more reason to exit early, so there's nothing to clean up.

\$\endgroup\$
3
\$\begingroup\$

awk is a very heavy duty utility. Usually, if you are using awk you want its full power. Calling it just for line splitting seems like a huge overkill. Let's go over one at a time:

    while read line_masses
    do
        mass=`echo $line_masses | awk -F' ' '{print $2}'`
        tag=`echo $line_masses | awk -F' ' '{print $1}'`
        ...

A read buit-in splits lines itself. Consider instead:

    while read tag mass
    do
        ...

Similarly,

    while read index hunoz1 a_tag hunoz2 rest:

allows the core to be rewritten as

while read index hunoz1 a_tag hunoz2 rest; do
    while read tag mass; do
        case $mass in
            12.011150)
            if [ "$tag" == "$a_tag]; then echo -e "C\t" "$rest"; fi;;
             1.007970)
            if [ "$tag" == "$a_tag]; then echo -e "H\t" "$rest"; fi;;
        esac
    done < masses.tmp
done < coordinates.tmp

Next, I'd try to get rid of the inner loop. The masses file should be read just once, and its data stored in array. This is how the core should look like:

declare -A elements

while read tag mass; do
    elements[tag]=mass
done < masses.tmp

while read index hunoz1 tag hunoz2 rest; do
    mass = $elements[tag]
    case $mass in
        12.011150) symbol=C ;;
         1.007970) symbol=H ;;
    esac
    echo $symbol "\t" $rest
done < coordinates.tmp

The preparation of temporary files may also be greatly simplified. I'll update on that later.

\$\endgroup\$
3
\$\begingroup\$

Minor issues

head -n -2 is not portable. The head in GNU Coreutils supports a negative number, but it's not POSIX.

If you create temporary files, use mktemp(1), so that you won't accidentally overwrite an existing file that is coincidentally has the same name as what you chose for the temp file. Also, immediately set a trap handler to clean up properly, even if the program aborts or is interrupted by ControlC.

If you detect a fatal condition (input file doesn't exist), bail out immediately. Don't put the error handler at the end, far from where you detected the problem. You would also avoid a level of indentation for nearly the entire script.

The error message should be printed to standard error, to avoid contaminating standard output. The exit status of your script should be non-zero in case of an error.

Single-language approach

Mixing Bash, AWK, head, and tail is rarely pretty. Why not write the whole program in AWK alone?

#!/usr/bin/env awk -f

BEGIN {
    OFS = "\t";
    ELEMENT_FOR_ATOMIC_WEIGHT["12.011150"] = "C";
    ELEMENT_FOR_ATOMIC_WEIGHT["1.007970"] = "H";
}

/^$/ {
    next;               # Ignore all blank lines
}

$2 == "atoms" {
    print $1
}

/Masses/, /Pair Coeffs/ {
    ATOMIC_WEIGHT_FOR_ATOM_TYPE[$1] = $2
}

/Atoms/, /Bonds/ {
    print ELEMENT_FOR_ATOMIC_WEIGHT[ATOMIC_WEIGHT_FOR_ATOM_TYPE[$3]], $5, $6, $7
}
\$\endgroup\$
2
  • \$\begingroup\$ Cool thank you so much for the clever awk script. For the creation of temporary files and the cleanup, if I understood your advice, I should write something like TMP_COORD=$(mktemp); TMP_MASS=$(mktemp); trap "rm $TMP_COORD $TMP_MASS; exit" SIGHUP SIGINT SIGTERM;, is that correct? \$\endgroup\$
    – Pierre
    Aug 8, 2014 at 15:23
  • \$\begingroup\$ You're welcome. That trap handler looks right. \$\endgroup\$ Aug 8, 2014 at 15:38

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