``` #!/bin/bash

# Check for input file
if [ -z "$1" ]; then
    echo "Usage: $0 <input.xyz>"
    exit 1
fi

INPUT="$1"
N=$((1+ $(grep -oPm1 'NAtoms=\s*\K\d+' $INPUT)))

mapfile -t -d ',' xyz < <(grep -Pzo '\d,\d\\([\s\S]*?)\\\\' $INPUT | tr -d '[:space:]' | tr '\\' ',' | tr -d '\000')

# Check if we got enough data points
if [ ${#xyz[@]} -lt 50 ]; then
    echo "Error: Not enough data points found in input file."
    exit 1
fi

# Parse atom 11 and 12 coordinates (indices are 0-based)
#      0  1
# Atom  0  X  Y  Z
#    2  3  4  5  6
#    7  8  9  10 11
# Atom 11: elements 42,43,44,45 (charge,x,y,z)
# Atom 12: elements 46,47,48,49 (charge,x,y,z)


vx=$(echo "${xyz[54]} - ${xyz[59]}" | bc -l)
vy=$(echo "${xyz[55]} - ${xyz[60]}" | bc -l)
vz=$(echo "${xyz[56]} - ${xyz[61]}" | bc -l)

# Normalize vector v
mag=$(echo "sqrt($vx^2 + $vy^2 + $vz^2)" | bc -l)
if (( $(echo "$mag < 1e-8" | bc -l) )); then
    echo "Error: atoms 11 and 12 are coincident; cannot determine direction."
    exit 1
fi

nx=$(echo "$vx / $mag" | bc -l)
ny=$(echo "$vy / $mag" | bc -l)
nz=$(echo "$vz / $mag" | bc -l)

# Target unit vector (Y-axis)
tx=0
ty=1
tz=0

# Compute dot product
dot=$(echo "$nx*$tx + $ny*$ty + $nz*$tz" | bc -l)


# Compute cross product to get rotation axis
ax=$(echo "$ny*$tz - $nz*$ty" | bc -l)
ay=$(echo "$nz*$tx - $nx*$tz" | bc -l)
az=$(echo "$nx*$ty - $ny*$tx" | bc -l)

# Normalize rotation axis
    mag_axis=$(echo "sqrt($ax^2 + $ay^2 + $az^2)" | bc -l)

    ax=$(echo "$ax / $mag_axis" | bc -l)
    ay=$(echo "$ay / $mag_axis" | bc -l)
    az=$(echo "$az / $mag_axis" | bc -l)

# Compute angle of rotation (using atan2 for better angle calculation)
    theta=$(echo "a(sqrt(1 - $dot^2)/$dot)" | bc -l)


    # Compute quaternion components
    qw=$(echo "c($theta / 2)" | bc -l)
    qx=$(echo "$ax * s($theta / 2)" | bc -l)
    qy=$(echo "$ay * s($theta / 2)" | bc -l)
    qz=$(echo "$az * s($theta / 2)" | bc -l)


    # Function to rotate a point using a quaternion
    function rotate_point {
        local x=$1
        local y=$2
        local z=$3

        # Compute conjugate of the quaternion
        local qc_x=$(echo "$qx * -1" | bc -l)
        local qc_y=$(echo "$qy * -1" | bc -l)
        local qc_z=$(echo "$qz * -1" | bc -l)
        local qc_w=$qw

        # First multiplication: q * p
        local w1=$(echo "- $qx*$x - $qy*$y - $qz*$z" | bc -l)
        local x1=$(echo "$qw*$x + $qy*$z - $qz*$y" | bc -l)
        local y1=$(echo "$qw*$y - $qx*$z + $qz*$x" | bc -l)
        local z1=$(echo "$qw*$z + $qx*$y - $qy*$x" | bc -l)


        # Second multiplication: result * qc

        local x_final=$(echo "$x1*$qc_w + $w1*$qc_x + $y1*$qc_z - $z1*$qc_y" | bc -l)
        local y_final=$(echo "$y1*$qc_w + $w1*$qc_y - $x1*$qc_z + $z1*$qc_x" | bc -l)
        local z_final=$(echo "$z1*$qc_w + $w1*$qc_z + $x1*$qc_y - $y1*$qc_x" | bc -l)

        echo "${x_final:0:10} ${y_final:0:10} ${z_final:0:10}"
    }

# Print header
echo '# nmr PBEQIDH/def2TZVPP' > $INPUT.gjf
echo "" >> $INPUT.gjf
echo "Test" >> $INPUT.gjf
echo "" >> $INPUT.gjf
echo "0 1" >> $INPUT.gjf



# Process each atom
for ((i=0; i<=$((5 * ($N - 1))); i+=5)); do
        echo -ne "${xyz[$i+2]} " >> $INPUT.gjf
        rotate_point ${xyz[$i+4]} ${xyz[$i+5]} ${xyz[$i+6]} >> $INPUT.gjf
done

echo "" >> $INPUT.gjf   ```