-   Programming (
-   -   BASH: IF statement not working? (

EdinburghLad 05-15-2012 12:49 PM

BASH: IF statement not working?

I'm trying to write a script which will generate the sum of the van der waal (VDW) radii for two atoms, then determine if a calculated bond distance is less than this sum and therefore whether the two atoms are bonded.

The script determines whether any atom is bonded to any other atom so runs as a loop within a loop (with outer loop counter FirstAtom and inner loop counter SecondAtom). In order to calculate the VDW sum the script needs to know which type of atom it is, so I wrote the fairly extensive IF statements you see in the code below to read from the AtomType array (which stores the atomic symbols) and then, depending on which atomic symbols they were, allocate a value to VDWa and VDWb, which are then added to get SumVDW. Finally the BondTest statement gives a boolean output for whether the atoms are bonded or not and writes it out to a file.

The problem is that when I run it, regardless of what the AtomType is, it always returns the "H" value for VDWa and VDWb. I've been staring at this for about two hours now and cant see where I've gone wrong, any help would be appreciated! (also if I haven't explained it well feel free to berate me and ask for more information!)

(ignore stuff in orange, it works fine)


while [[ $FirstAtom -le $NoAtoms ]]; do
        if [ ${AtomType[$FirstAtom]}="H" ]; then
        elif [ ${AtomType[$FirstAtom]}="C" ]; then
        elif [ ${AtomType[$FirstAtom]}="S" ]; then
        elif [ ${AtomType[$FirstAtom]}="N" ]; then
        elif [ ${AtomType[$FirstAtom]}="Pt" ]; then
        elif [ ${AtomType[$FirstAtom]}="As" ]; then
                while [[ $SecondAtom -le $NoAtoms ]]; do
                        Xd[$SecondAtom]=`echo "scale=4; (${Xc[$FirstAtom]}-${Xc[$SecondAtom]})^2" | bc -l`
                        Yd[$SecondAtom]=`echo "scale=4; (${Yc[$FirstAtom]}-${Yc[$SecondAtom]})^2" | bc -l`
                        Zd[$SecondAtom]=`echo "scale=4; (${Zc[$FirstAtom]}-${Zc[$SecondAtom]})^2" | bc -l`
                        Dist=`echo "scale=2; sqrt(${Xd[$SecondAtom]}+${Yd[$SecondAtom]}+${Zd[$SecondAtom]})" | bc -l`

                        echo ${AtomType[$SecondAtom]}
                        if [ ${AtomType[$SecondAtom]}="H" ]; then
                        elif [ ${AtomType[$SecondAtom]}="C" ]; then
                        elif [ ${AtomType[$SecondAtom]}="S" ]; then
                        elif [ ${AtomType[$SecondAtom]}="N" ]; then
                        elif [ ${AtomType[$SecondAtom]}="Pt" ]; then
                        elif [ ${AtomType[$SecondAtom]}="As" ]; then
                        SumVDW=`echo "$VDWa + $VDWb" | bc -l`
                        echo $VDWb
                        BondTest=`echo "$Dist <= $SumVDW" | bc -l`
                                echo -e $Bondtest "\t" >> Column[$FirstAtom].txt
                let "SecondAtom=$SecondAtom+1"
let "FirstAtom=$FirstAtom+1"

Snark1994 05-15-2012 01:21 PM

I believe the issue is that there's no whitespace around the '=', though there may be other issues as well. You may also want to consider a 'case' statement rather than lots of 'if/elif's, as it's easier to read.

(Also, bash may not be the best language to do arithmetic-heavy operations in - something like python might be a better choice)

pan64 05-15-2012 01:30 PM

yes, you need a lot of spaces to separate variables, equal sign and [ and ] . Also it was suggested to use [[ and ]] instead of [ and ].
From the other hand you can write a better code with case:

case ${AtomType[$FirstAtom]} in
  H) VDWb=1.2 ;;
  C) VDWb=1.7 ;;

(similar works for the second atom too)

Nominal Animal 05-15-2012 06:17 PM

It would be much, much easier to write in awk.

Here is a hopefully useful example:

#!/usr/bin/awk -f

    # Use any newline convention for record separators.
    RS = "(\r\n|\n\r|\r|\n)"

    # Fields are separated by whitespace.
    FS = "[\t\v\f ]+"

    # Table of Van der Waals radiis
    VdWr["H"] = 1.2
    VdWr["C"] = 1.7
    VdWr["S"] = 1.8
    VdWr["N"] = 1.55
    VdWr["Pt"] = 1.75
    VdWr["As"] = 1.85

    # Table of the squared distances to compare against
    for (t1 in VdWr)
        for (t2 in VdWr)
            VdWrr[t1 "-" t2] = (VdWr[t1] + VdWr[t2])^2

# This will be run once per record.
(NF >= 4) {
    # We have at least four fields in this record ($1 $2 $3 $4).
    # Assume it is "Type X Y Z", and remember the atom.
    atype[atoms]  = $1
    xcoord[atoms] = $2
    ycoord[atoms] = $3
    zcoord[atoms] = $4

# This is only run after all input is processed.
    # Find all bonded pairs, and output them as
    # Type X Y Z  Type X Y Z  distance
    for (i = 1; i < atoms; i++) {
        for (j = i + 1; j <= atoms; j++) {
            dx = xcoord[i] - xcoord[j]
            dy = ycoord[i] - ycoord[j]
            dz = zcoord[i] - zcoord[j]
            dd = dx^2 + dy^2 + dz^2
            if (dd <= VdWrr[atype[i] "-" atype[j]])
                printf("%s %s %s %s  %s %s %s %s  %.6f\n", atype[i], xcoord[i], ycoord[i], zcoord[i], atype[j], xcoord[j], ycoord[j], zcoord[j], sqrt(dd))

If you have a lot of data, see if you have mawk installed; if so, use #!/usr/bin/mawk -f for the first line instead. It seems to be the fastest awk variant for this kind of processing, much faster than GNU awk (gawk).

jlinkels 05-15-2012 07:47 PM

Dunno how you define array AtomType[]. Might be very relevant in this case. But I can't replay your script now.

Run your script with bash -x /path/to/your/script and see where it fails. Should be obvious. If not, post the output where we see the values of AtomType[$FirstAtom] and how the comparison fails. Indeed is a case statement better in terms of programming style, but the if/elif should show nicely why the comparison fails.


EdinburghLad 05-15-2012 08:46 PM

Solved! It was just a white space error, thanks guys, I hate tiny bits of annoying syntax. To those who suggested I code in a different language I would usually use fortran, which I'm a bit more familiar with, unfortunately I seem to be on the one system in the world without a decent compiler and no admin privileges to install one! Anyhow the script works beautifully and has already shown where my research is going wrong, so thanks guys!

Nominal Animal 05-15-2012 09:26 PM


Originally Posted by EdinburghLad (Post 4679428)
To those who suggested I code in a different language I would usually use fortran, which I'm a bit more familiar with, unfortunately I seem to be on the one system in the world without a decent compiler and no admin privileges to install one!

I suspected exactly that. That's why I suggested awk. It is always available (except perhaps on Windows systems).

If you ever generate test lattices, initial coordinates for nanoclusters, or duplicate molecules (with translation and rotation and optionally mirroring), to be fed to a simulation of some kind, you'll find that awk is superior tool for that. While it is nowhere as fast as compiled code, you can do in a single line or just a few lines something that would take you a couple of hours to write in Fortran.

(For rotations, I suggest using the quaternion approach to construct the rotation matrix: let the user specify some rotation axis, and an angle to rotate around that axis. Euler angles are difficult to generate uniformly ("random orientation"), and suffer from gimbal lock. The rotation axis and angle are trivial in comparison.)

Also, if you ran that shell script on anything more than a handful of atoms on a cluster where I was the administrator, I would be quite livid about the CPU wastage: calling bc to compute each interatomic distance.. I hope you at least run that code only on your own workstation, not on a cluster.

All times are GMT -5. The time now is 03:59 AM.