#!/bin/bash
atatdir=$(sed 's/.*=\(.*\)/\1/g' ~/.atat.rc)

if [[ $# -eq 0 || "x$1" == "x-h" ]]
then
  cat - <<-EOF
WYCKoff position EDitor, by Axel van de Walle
Syntax:
  wycked spcgrp [-r] list_of_wyckoff_position_letter_without_seperators [species1:species2:...]
The output is meant to first be edited to set the free parameters a,b,c,alpha,beta,gamma,x,y,z,etc
in each of the lines starting with #const 
and processed by sspp to create a valid ATAT structure file.
 If the -r option is specified, then:
   the structure is relaxed using a soft-sphere model (of radius 1.5 A) by invoking the GULP code
   (multiple times (10) with randomized initial positions).
   In that case, you can also specify a colon-separated list of species that will populate each wyckoff position.
   Also, the structure will be isotropically re-scaled so that the nearest neighbor distance are consistent
   with atom radii from the atat/data/radii.in file.

ACKNOWLEDGEMENTS
This code makes use of a Wyckoff position database downloaded from the Bilbao Crystallographic Server
at http://www.cryst.ehu.es/
Accordingly, if you use this code, you should cite, in addition to the usual ATAT citations,
the following papers:
  M. I. Aroyo, J. M. Perez-Mato, D. Orobengoa, E. Tasci, G. de la Flor, A. Kirov
      "Crystallography online: Bilbao Crystallographic Server" 
      Bulg. Chem. Commun. 43(2) 183-197 (2011).
  M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga,
  A. Kirov and H. Wondratschek,
      "Bilbao Crystallographic Server I: Databases and crystallographic computing programs" 
      Z. Krist. 221, 1, 15-27 (2006). doi:10.1524/zkri.2006.221.1.15 
  M. I. Aroyo, A. Kirov, C. Capillas, J. M. Perez-Mato and H. Wondratschek
      "Bilbao Crystallographic Server II: Representations of crystallographic point groups and
      space groups", Acta Cryst. A62, 115-128 (2006). doi:10.1107/S0108767305040286 
If you use the -r option, you should also acknowledge the GULP code:
  "GULP - a computer program for the symmetry adapted simulation of solids", J.D. Gale, JCS Faraday Trans., 93, 629 (1997)
  "The General Utility Lattice Program", J.D. Gale and A.L. Rohl, Mol. Simul., 29, 291-341 (2003)
  "GULP: Capabilities and prospects", J.D. Gale, Z. Krist., 220, 552-554 (2005)
EOF
exit
fi

dorelax=0
if [[ $1 == "-r" ]]
then
  dorelax=1
  shift
fi

sgonly=$(echo $1 | sed 's/s[12]//g')
sgset=$(echo $1 | sed 's/^[^s]*//g')

spcgrpdata=($(awk '{if ($1=='\"$sgonly\"' || $2=='\"$sgonly\"') {print $1,$3,$4;}}' ${atatdir}/data/wyckoff/spcgrp.txt) )

if [[ ${#spcgrpdata} == 0 ]]
then
  echo "Space Group $sgonly not found."
  exit 1
fi

sgnum=${spcgrpdata[0]}
sgaxes=${spcgrpdata[1]}
sgcell=${spcgrpdata[2]}

echo -n >| wycked${$}.tmp

cat ${atatdir}/data/wyckoff/axes_${sgaxes}.txt >> wycked${$}.tmp
cat ${atatdir}/data/wyckoff/cell_${sgcell}.txt >> wycked${$}.tmp

for ((i=0 ; i < ${#2} ; i++))
do
  site=${2:i:1}
  grep ${site}\$ ${atatdir}/data/wyckoff/wyckoff${sgnum}${sgset}.txt | awk '{print "{"$1"}","{"$2"}","{"$3"}","{"toupper($4)"}"}' >| wycktmp${$}.tmp
  echo -n "#const "
  for pvar in x y z
  do
     grep -q "$pvar" wycktmp${$}.tmp && echo -n "${pvar}=; "
  done
  #echo ${site^}=\"${site^}\"\;
  echo ${site} | awk '{print toupper($1)"=\"" toupper($1) "\";"}'
  cat wycktmp${$}.tmp
done >> wycked${$}.tmp
rm -f wycktmp${$}.tmp

if [[ $dorelax == 0 ]]
then
  cat wycked${$}.tmp
else

rtimes=10
gulp=gulp
which ${gulp} >& /dev/null
if [[ $? != 0 ]]
then
    echo Unable to find the ${gulp} command. Please install GULP from http://gulp.curtin.edu.au to use the -r option.
  exit 1
fi

tmpdir=tmpgulp_$$

mkdir $tmpdir
pushd $tmpdir > /dev/null

bestenergy=100

echo -n >| gulpenergy.out

for (( r=1; r<=rtimes; r++ ))
do

cat >| gulp.in <<EOF
opti conp gradient prop noelectrostatic
species
Al 0
morse 
Al Al 0.05 4.0 3.0 6.0
EOF

cat ../wycked${$}.tmp | awk 'BEGIN {
s=srand();
srand(s+'$r');
print "cell";
getline;
print  "#const a=1.+2.*"rand()"; b=1.+2.*"rand()"; c=1.+3.*"rand()"; alpha=60.+60.*"rand()"; beta=60.+60.*"rand()"; gamma=60.+60.*"rand();
getline; print $0;
for (i=1; i<=3; i++) {getline;}
print "spacegroup";
print '$1';
print "fractional";
}
{
if ($1=="#const") {
  for (i=1; i<=length($0); i++) {
    c=substr($0,i,1);
    if (c==";" && substr($0,i-1,1)=="=") {
      printf rand() ";";
    }
    else if (c=="\"") {
      printf "\"Al\"";
      i+=2;
    }
    else {
      printf c;
    }
  }
  print "";
  getline;
  print $0;
}
} END {print "dump gulpdump.out";}' | sspp -nb >> gulp.in

#cat gulp.in

${gulp} < gulp.in > gulp.out

energy=$(getvalue "Final energy" < gulp.out)

echo $energy >> gulpenergy.out

if [[ $(echo $energy $bestenergy | awk '{print ($1<$2)}') == 1 ]]
then
  bestenergy=$energy
  cp gulpdump.out gulpdumpkeep.out
fi

done

(echo nosymmetry ; grep -v -e opti -e iterations -e dump gulpdumpkeep.out ; echo "dump gulpdump2.out") | ${gulp} > gulp2.out

getlines -jbt fractional species < gulpdump2.out | grep -v 'totalenergy' | awk '{print "{"$3"}","{"$4"}","{"$5"}"}' > posfromgulp.out

if [[ $3 == "" ]]
then
  cat ../wycked${$}.tmp | tail -n +6 | awk 'BEGIN {last="";} /^#const/ {letter=substr($0,length($0)-2,1); if (letter!=last) {num=1;} else {num++;} last=letter;} ! /^#const/ {print letter num;}' > atomfromw.out
  ( getlines -jbt cell fractional < gulpdump2.out  ; echo 1 0 0 ; echo 0 1 0 ; echo 0 0 1 ; paste posfromgulp.out atomfromw.out | awk '{print $1,$2,$3,$4}' | sspp ) | cellcvrt -cc -abc
else
  cat ../wycked${$}.tmp | tail -n +6 | awk 'BEGIN {atomstr="'"$3"'"; split(atomstr,atoms,":"); s=0} /^#const/ {s++} ! /^#const/ {print atoms[s];}' > atomfromw.out
  ( getlines -jbt cell fractional < gulpdump2.out  ; echo 1 0 0 ; echo 0 1 0 ; echo 0 0 1 ; paste posfromgulp.out atomfromw.out | awk '{print $1,$2,$3,$4}' | sspp ) | cellcvrt -cc -abc | nntouch -l=- -r=${atatdir}/data/radii.in
fi

popd  > /dev/null

rm -r $tmpdir

fi

rm -f wycked${$}.tmp
