#!/bin/csh

if ( ! -e ~/.ezvasp.rc ) then
cat - >! ~/.ezvasp.rc <<EOF
#!/bin/csh
#enter name of vasp executable here
set VASPCMD="vasp"
#enter the directories containing the pseudopotentials here
set POTLDA="/usr/local/lib/pot_lda"
set POTGGA="/usr/local/lib/pot_gga"
set POTPAWLDA="/usr/local/lib/paw_lda"
set POTPAWGGA="/usr/local/lib/paw_gga"
set POTPAWPBE="/usr/local/lib/paw_pbe"
EOF
cat - <<EOF
I have created a default ~/.ezvasp.rc configuration file for you.
Please edit it to match your vasp installation before rerunning ezvasp.
EOF
exit 1
else
  source ~/.ezvasp.rc
endif

if ( $#argv == 0 ) then
  echo "Syntax: ezvasp [-options] input_file"
  echo "-h: help"
  echo "-c: cell patch"
  echo "-s: do static run afterward"
  echo "-n: do not run vasp"
  echo "-p: specifies the vasp command (default: vasp)"
  exit 1
endif

foreach option ( $* )
  if (x$1 == "x-p") then
    shift
    set VASPCMD="$1"
    shift
  endif
  if (x$1 == "x-h") then
cat - <<EOF
Syntax: ezvasp [-c] [-s] [-n] [-p "..."] input_file

-c: Cell patch. This option will convert the cell you entered into a cell that has the symmetry needed for
    vasp to work properly. The atom positions are altered correspondingly.

-s: Static run. After the relaxations are completed, this option makes vasp start again, turning off relaxations
    and turning on tetrahedron integration with bloech correction. Can acheive the same effect by
    including DOSTATIC in the INCAR portion of the input file.

-n: do Not run vasp. If you want to create the input files now and run them later, perhaps on another machine.

-p: name of the vasp command ("vasp" by default). You can include addition prefix such as:
    ezvasp -p "rsh remote.machine nice +5 /usr/local/bin/vasp" vasp.in

Here is an example of input file:

[INCAR]
SYSTEM = FeNi
ISPIN = 2
PREC = HIGH
ISMEAR = 1
NSW=41
SIGMA = 0.1
IBRION = 2
ISIF = 3
ENMAX = 400
EDIFF = 1e-6
EDIFFG = 1e-4
KPPRA = 1000
USEPOT = PAWGGA
DOSTATIC
SUBATOM = s/Fe\$/Fe_h/g
KSCHEME = GAMMA
MAGATOM =

[POSCAR]
 SOME TITLE
 3.55
  -.500000000     -.500000000    -1.000000000   
   .500000000    -1.000000000     -.500000000   
  -.500000000      .500000000     1.000000000   
 1 2
 Direct
  .33333     .33333     .00000     Fe+5
  .00000     .00000     .00000     Vac
  .66667     .66667     .00000     Fe-5

End of example.

Note that there 6 new tokens in the INCAR section.

1) KPPRA stands for "K-Point Per Reciprocal Atom". This is a way to automatically set the k-point mesh for
   a number of similar systems. Here is how it works. If you know that you need 500 k-points for a 2-atom
   structure, then you type in 1000. If you now try a similar structure with 4 atoms, the code automatically
   use 250 k-points. The mesh along the three axes is automatically chosen to make the mesh as uniform as possible.
   (For the algorithm, see atat/src/kmesh.cc).
   You can also specify KPPRA = UPDIR which indicates that the KPOINTS file
   found one directory above should be used.
   (This is useful for cell distorsions used in elastic constant calculations of constituent strain runs).
   If this token is ommited, you need to specify a [KPOINT] section in the input file that contains what
   you want the KPOINT file to be.

2) DOSTATIC indicates that a static run must be done after the relaxation run (NSW tag is removed and the ISIF and IBRION
   token are ajusted accordingly). The output files named *.relax contain the result of the (first)
   relaxation run while the files named *.static contain the results of the (second) static run.
   When DOSTATIC is not specified, all output files are named *.static , whether the VASP input parameters
   allow relaxations or not.

3) KSCHEME specifies the type mesh to use: Gamma shift or Monkorst-Pack. The default is Gamma if KSCHEME is ommited.

4) MAGATOM is similar to the usual token MAGMOM used to specify the moment on each atom. The novelty is that you specify the
   moments in the POSCAR section. (In version 3.15 and above, this token can be ommited. Having spins specified in the POSCAR
   section turns this on automatically. It also set ISPIN=2 automatically.)
   The usual MAGMOM token is still available (but you cannot use both MAGMOM and MAGATOM).

5) USEPOT 
   selects the potentials used: LDA GGA PAWLDA PAWGGA PAWPBE (or any other potential defined in .ezvasp.rc)
   By default, LDA (without PAW) is used.
   For backward compatibility, you can also use the tag DOGGA instead of USEPOT = GGA

6) The SUBATOM token lets you specify a sed-like substitution command to be
   applied to the name of the atomic species.
   For instance, SUBATOM= s/Fe\$/Fe_h/g  will cause ezvasp to use the Fe_h
   pseudopotential instead of the Fe one. Multiple SUBATOM token can be given,
   one for each substitution. Don't forget the s/ and /g and the \$ after the first element name
   (to avoid F matching Fe, or Y matching Yb).

There are also familiar VASP tokens that are robustified. For instance ISTART = 2 does not work unless you force
the FFT grid size to match the earlier run.
To help with this, ezvasp fetches the right NGX,NGY,NGZ from the previous run's OUTCAR or CHGCAR in the current directory.

There are new features in the POSCAR section as well.
1) The line that specify the number of atoms (here 1 2) can be ommited (it is ignored).
2) The atomic species are specified to the right of the atom positions.
   The corresponding POTCAR are automatically copied. You can specify pseudopotentials like Li_h or O_s.
3) The number following the species is the magnetic moment on the atom (if MAGMOM is in the INCAR section).
   Identical species with different moments will be considered symmetrically distinct to allow for
   symmetry-breaking relaxations.
4) Do not specify "selective dynamics". It will be written to the INCAR file
   automatically if you include T's and F's after the atom species, e.g.:
    0.5 0.5 0.5 Al T T F
5) The element "Vac" or "Va" stands for vacancy. It is ignored for everything but
   to compute the k-point mesh with the KPPRA option.
EOF
    exit 1
  endif
  if (x$1 == "x-c") then
    set latvectpatch
    shift
  endif
  if (x$1 == "x-s") then
    set dostatic
    shift
  endif
  if (x$1 == "x-n") then
    set norun
    shift
  endif
  if (x$1 == "x-d") then
    set debug
    shift
  endif
end

if (! -e $1) then
  echo Cannot open $1
  exit 1
endif

foreach mysection ( incar poscar )
  grep -q -i '\['$mysection $1
  if ( $status == 1 ) then
    echo "No \[${mysection}\] section in input file."
    exit 1
  endif
end

sspp < $1 >! vaspin.tmp

getlines -jaf "\[[Ii][Nn][Cc]" < vaspin.tmp | getlines -jbf "\\[" | grep -v '^[ 	]*$' | \
  grep -v -i -e subatom -e magatom -e kppra -e kscheme -e dostat -e dogga -e usepot -e 'istart[ 	]*=[ 	]*2' >! INCAR

grep -i subatom vaspin.tmp | sed 's/.*=//g' >! sed.tmp

(getlines -jaf "\[[Pp][Oo][Ss]" < vaspin.tmp | getlines -jbf "\\[") | grep -v '^[ 	]*$' >! preposcar.tmp
cat preposcar.tmp | sed -f sed.tmp >! poscar.tmp

cat preposcar.tmp | tail -n +7 | awk '{print $4}' >! atomlabel.tmp1
cat preposcar.tmp | tail -n +7 | awk '{print $4}' | sed -f sed.tmp >! atomlabel.tmp2
paste atomlabel.tmp2 atomlabel.tmp1 | sort -k 1,1 -s | awk '{n=split($0,a); print a[n]}' >! atomlabel.tmp

tail -n +3 poscar.tmp | head -3 >! oldcell.tmp

set posmode=`head -6 poscar.tmp | tail -n -1`

tail -n +7 poscar.tmp | awk '{print $5,$6,$7}' | grep -q '[TF]'
if ( $status == 0 ) then
  set seldyn
endif

if ($?latvectpatch) then
  cat oldcell.tmp | fixcell -c >! newcell.tmp
  cp -f newcell.tmp cell.tmp
else
  cp -f oldcell.tmp newcell.tmp
  cp -f oldcell.tmp cell.tmp
endif

tail -n +6 poscar.tmp | head -1 | grep -q -i '^d'
if ($status == 0 && $?latvectpatch) then
  tail -n +7 poscar.tmp | \
   awk 'BEGIN {for (i=1; i<=3; i++) {getline < "oldcell.tmp"; a[1,i]=$1; a[2,i]=$2; a[3,i]=$3; }} \
       {c[1]=0; c[2]=0; c[3]=0; b[1]=$1; b[2]=$2; b[3]=$3; \
        for (i=1; i<=3; i++) { \
          for (j=1; j<=3; j++) {c[i]=c[i]+a[i,j]*b[j]} \
        } \
        print $4,c[1],c[2],c[3],$5,$6,$7}' | sort -k 1,1 -s >! vacpos.tmp
    set posmode=Cartesian
else
    tail -n +7 poscar.tmp | awk '{print $4,$1,$2,$3,$5,$6,$7}' | sort -k 1,1 -s >! vacpos.tmp
endif

grep -v Va vacpos.tmp >! pos.tmp
#grep -v Vac vacpos.tmp >! pos.tmp

awk '{nb_atom[$1]++} END {for (atom in nb_atom) {print atom,nb_atom[atom]}}' pos.tmp | sort -k 1,1 -s >! atom.tmp
( head -2 poscar.tmp ; \
  cat cell.tmp ; \
  echo `awk '{print $2}' atom.tmp` ; \
  if ($?seldyn) echo Selective Dynamics ; \
  echo $posmode ; \
  awk '{print $2,$3,$4,$5,$6,$7,$1}' pos.tmp \
) >! POSCAR

grep -q -i dostat vaspin.tmp
if ($status == 0) then
  set dostatic
endif

grep -q -i kppra vaspin.tmp
if ($status == 0) then
  grep -q -i 'kppra *= *updir' vaspin.tmp
  if ($status == 0) then
    cp -f ../KPOINTS .
  else
    grep -q -i kscheme vaspin.tmp
    if ($status == 0) then
	set kscheme=`grep -i kscheme vaspin.tmp | sed 's/^.*=[ 	]*//g'`
    else
	set kscheme="Gamma"
    endif
    ( echo KPOINTS file ; echo 0 ; echo $kscheme ) >! KPOINTS
    echo $kscheme | grep -q "^[Mm]"
    if ( $status == 0 ) then  
	( echo `grep -i kppra vaspin.tmp | sed 's/^.*=//g'` `cat vacpos.tmp | wc -l` | awk '{print $1/$2}' ; cat cell.tmp ) | kmesh -q -e >>! KPOINTS
    else
	( echo `grep -i kppra vaspin.tmp | sed 's/^.*=//g'` `cat vacpos.tmp | wc -l` | awk '{print $1/$2}' ; cat cell.tmp ) | kmesh -q -r >>! KPOINTS
    endif
    echo 0 0 0 >>! KPOINTS
  endif
else
  (getlines -jaf "\[[Kk][Pp][Oo]" < vaspin.tmp | getlines -jbf "\\[") | grep -v '^[ 	]*$' >! KPOINTS
endif

grep -q -i dogga vaspin.tmp
if ($status == 0) then
  set POTPATH="$POTGGA"
else
  set POTPATH="$POTLDA"
endif

grep -q -i usepot vaspin.tmp
if ($status == 0) then
  echo set POTPATH=\$POT`grep -i usepot vaspin.tmp | tail -n -1 | sed -e 's/ *[^ =]*[ =]//g' -e 's/ *$//g'` > tmpez.sh
  source tmpez.sh
  rm -f tmpez.sh
endif

echo -n  >! POTCAR
foreach element ( `sed 's/[+-].*$//g' atom.tmp | awk '{print $1}'` )
  if ( -e $POTPATH/$element/POTCAR ) then
    cat $POTPATH/$element/POTCAR >>! POTCAR
  else
    zcat $POTPATH/$element/POTCAR.Z >>! POTCAR
  endif
  if ( $status == 1 ) then
    echo Cannot read $POTPATH/$element/POTCAR\*
    exit 1
  endif
end

awk '{print $1}' pos.tmp | sed 's/^[^-+]*//g' | awk '{print $1+0}' >! moment.tmp
grep -q -i '[1-9]' moment.tmp
if ( $status == 0 ) then
  echo "ISPIN= 2" >>! INCAR
  echo "MAGMOM= " `cat moment.tmp` >>! INCAR
endif

grep -q -i 'ISTART[ 	]*=[ 	]*2' vaspin.tmp
if ($status == 0) then
  if ( ( -e OUTCAR || -e CHGCAR ) && -e WAVECAR ) then
    echo ISTART=1 >>! INCAR
    if ( -e OUTCAR ) then
      grep "dimension x,y,z NGX =" OUTCAR | head -n -1 | sed -e 's/   dimension x,y,z //g'  -e 's/ NG/\nNG/g' >>! INCAR
    else
      getlines -jaf '^ *$' < CHGCAR | head -n 1 | awk '{print "NGX = " $1; print "NGY = " $2; print "NGZ = " $3}' >>! INCAR
    endif
  else
    echo ISTART=0 >>! INCAR
  endif
endif

if (! $?debug ) then
  rm -f sed.tmp cell.tmp oldcell.tmp newcell.tmp pos.tmp vacpos.tmp atom.tmp poscar.tmp moment.tmp cell.log atomlabel.tmp1 atomlabel.tmp2 preposcar.tmp vaspin.tmp >& /dev/null
endif

if (! $?norun ) then
  echo | $VASPCMD >&! vasp.out

  grep -q 'VERY BAD NEWS! internal error in subroutine IBZKPT' vasp.out
  if ( $status == 0 ) then
    cp INCAR INCAR.tmp
    ( grep -v SYMPREC INCAR.tmp ; echo SYMPREC = 1e-7 ) >! INCAR
    echo | $VASPCMD >&! vasp.out
    rm -f INCAR.tmp
  endif

  if ($?dostatic) then
    foreach file (INCAR POSCAR OSZICAR OUTCAR CONTCAR KPOINTS vasp.out)
      cp -f $file ${file}.relax
    end
    cp -f CONTCAR POSCAR
    cat INCAR.relax | \
        grep -v -i ismear | \
        grep -v -i isigma | \
        grep -v -i nsw | \
        grep -v -i ibrion | \
        grep -v  -i isif >! INCAR
    echo ISMEAR=-5 >>! INCAR
    echo IBRION=-1 >>! INCAR
    echo ISIF=2 >>! INCAR

    echo | $VASPCMD >&! vasp.out
  endif
  foreach file (INCAR POSCAR OSZICAR OUTCAR CONTCAR KPOINTS vasp.out DOSCAR)
    cp -f $file ${file}.static
  end
endif