#!/bin/ksh


# Script to convert the CTRL file for the LMTO code
# to an ABINIT compatible input file for the purposes
# of visualizing unit cells. No guarantee that it works
# in general and the analysis follows no particular 
# logic.
# J M Sullivan NRL 2/21/02

# JMS NRL 2/25/02
# Modified to relabel atom types in convenient format disregarding
# symmetry labels assigned by the ground state program

# JMS NRL 2/26/02
# Added spinat token to input file.
ctrl=CTRL

out=tmp
out=input.$out
rm -f $out

# Set up an array of atomic labels for atoms 1<=Z<=109
set -A atom_labels \
E H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe Cs Ba La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Db Sg Bh Hs Mt sullivanium
#echo ${atom_labels[*]}
#exit 0

# xmolprep parameters
echo "# xmolnext nx ny nz makes periodic images of the unit cell with " >| $out
echo "# nx translations in the a1 direction, ny translations in the a2 direction, etc." >> $out
echo "echo 0 xmolcent 0 0 0 xmolnext 1 1 1" >> $out
echo >> $out
echo "# xmolside ax ay az produces atoms in a box of size ax*ay*az" >> $out
echo "#echo 0 xmolcent 0 0 0 xmolside 3*20 " >> $out
echo >> $out

# Lattice constant
acell=$(grep ALAT $ctrl | sed 's/=/\ /g' | awk '{printf "%12.8f \n", $3}')
echo "Lattice constant:" $acell
echo "acell" $acell $acell $acell >> $out

# Spin polarization
nsppol=$(grep NSPIN $ctrl | sed 's/=/\ /g' | awk '{printf "%s \n", $3}')
echo "Spin polarization:", $nsppol
echo "nsppol" $nsppol >> $out

# Rprim
ns=$(grep -n PLAT $ctrl | sed 's/:/\ /'g | head -1 | awk '{print $1}')
a1=$(sed -n ''$ns' p' $ctrl | sed 's/=/\ /g' | awk '{print $2, $3, $4}')
ns=$((ns+1))
a2=$(sed -n ''$ns' p' $ctrl | awk '{print $1, $2, $3}')
ns=$((ns+1))
a3=$(sed -n ''$ns' p' $ctrl | awk '{print $1, $2, $3}')
echo "Rprim"
echo $a1
echo $a2
echo $a3
echo "rprim" >> $out
echo $a1 >> $out
echo $a2 >> $out
echo $a3 >> $out

# Default Epw
ecut=10.0
echo "ecut" $ecut >> $out


# Atomic types
set -A atoms $(grep "ATOM=" CTRL | grep "Z=" | sed 's/=/\ /g' | sed 's/CLASS//g' | awk '$4>0 {print $2}')
set -A zatoms $(grep "ATOM=" CTRL | grep "Z=" | sed 's/=/\ /g' | sed 's/CLASS//g' | awk '$4>0 {print $4}')

# Minimize the atom count according to the atomic numbers
rm -f tmp.dump
i=0
nt=${#zatoms[*]}
while [ $i -le $((nt-1)) ]; do
    echo ${atoms[$i]} ${zatoms[$i]} >> tmp.dump
    i=$((i+1))
done
mv -f tmp.dump tmp.dump2
cat tmp.dump2 | sort -n -k 2 >| tmp.dump
#exit 0
old=-1
j=0
set -A zt
set -A lab
set -A nc
while read trash1 trash
do
    if [ $trash -ne $old ]; then
	j=$((j+1))
	zt[$j]=$trash
	lab[$j]=${atom_labels[$trash]}
	old=$trash
	nc[$j]=1
#    else
#	echo "#"
    fi
done < tmp.dump
echo "Number of distinict atomic numbers:" ${#zt[*]}
ntype=${#zt[*]}
#exit 0

echo ${zt[*]}
set -A atoms ${lab[*]}
echo ${atoms[*]}
set -A zatoms ${zt[*]}
#exit 0



# Number of atoms
echo "ntype" $ntype >> $out
echo "Atom types:" ${atoms[*]}
echo "Atomic numbers:" ${zatoms[*]}
echo "zatnum" ${zatoms[*]} >> $out
#echo "Number of atomic numbers:" ${#zatoms[*]}

#exit 0
# Atomic coordinates
rm -f tmp.dump
i=0
natoms=0
nt2=1
rm -f tmp.dump3
for a in ${atoms[*]}
do
#    echo $a
    grep "ATOM=$a" CTRL | grep "POS=" | sed 's/SITE//g' | sed 's/=/\ /g' | awk '{printf "%12.8f %12.8f %12.8f \n", $4, $5, $6}' >| tmp.dump
    nt=$(wc -l tmp.dump | awk '{print $1}')
    nc[$((i+1))]=$nt
#    echo "Number of atoms of this type:" $nt
    natoms=$((natoms+nt))
    while read x y z
    do
#	echo $a $x $y $z
#    echo $i ${zatoms[$i]}
#    exit 0
#	    echo $acell $x $y $z $nt2 ${zatoms[$i]} | awk '{printf "%12.8f %12.8f %12.8f %s %s %s \n", $1*$2, $1*$3, $1*$4, "#", $5, $6}' >> $out    
	    echo $acell $x $y $z $nt2 ${zatoms[$i]} | awk '{printf "%12.8f %12.8f %12.8f %s %s %s \n", $1*$2, $1*$3, $1*$4, "#", $5, $6}' >> tmp.dump3   
	nt2=$((nt2+1))
    done < tmp.dump
    rm -f tmp.dump
    i=$((i+1))
done

echo "natom" $natoms >> $out
echo "xcart" >> $out
sort -n -k 6 tmp.dump3 >> $out

i=1
echo "type" >> $out
while [ $i -le $ntype ]; do
    echo ${nc[$i]} $i | awk '{printf "%s \n", $1"*"$2}' >> $out
    i=$((i+1))
done

rm -f tmp.dump2
rm -f tmp.dump3

# Generate default spinat tokens
# Should be generalized to grab the spin moments
# from the atomic charges.
if [ $nsppol -eq 2 ]; then
#    echo "spinat" >> $out
    set -A spinat
    i=1
    while [ $i -le $natoms ]; do
	spinat[$i]="0.0 0.0 0.0"
#	echo ${spinat[$i]} >> $out
	i=$((i+1))
    done
fi

echo 
echo "########################################"
blah=$(echo $out | sed 's/input.//g')
echo "Run xmp" $blah " to generate an xmol (jmol) file."

exit 0