#!/bin/sh

# run from directory where this script is
cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
EXAMPLE_DIR=`pwd`

# check whether echo has the -e option
if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi

$ECHO
$ECHO "$EXAMPLE_DIR : starting"
$ECHO
$ECHO "This example shows how to use cp.x to calculate Wannier functions and"
$ECHO "to perform dynamics with an external electric field."

# set the needed environment variables
. ../../../environment_variables

# required executables and pseudopotentials
BIN_LIST="cp.x"
PSEUDO_LIST="H_US.van O_US.van"

$ECHO
$ECHO "  executables directory: $BIN_DIR"
$ECHO "  pseudo directory:      $PSEUDO_DIR"
$ECHO "  temporary directory:   $TMP_DIR"
$ECHO
$ECHO "  checking that needed directories and files exist...\c"

# check for directories
for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
    if test ! -d $DIR ; then
        $ECHO
        $ECHO "ERROR: $DIR not existent or not a directory"
        $ECHO "Aborting"
        exit 1
    fi
done
for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
    if test ! -d $DIR ; then
        mkdir $DIR
    fi
done
cd $EXAMPLE_DIR/results

# check for executables
for FILE in $BIN_LIST ; do
    if test ! -x $BIN_DIR/$FILE ; then
        $ECHO
        $ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
        $ECHO "Aborting"
        exit 1
    fi
done

# check for pseudopotentials
for FILE in $PSEUDO_LIST ; do
    if test ! -r $PSEUDO_DIR/$FILE ; then
       $ECHO
       $ECHO "Downloading $FILE to $PSEUDO_DIR...\c"
            $WGET $PSEUDO_DIR/$FILE $NETWORK_PSEUDO/$FILE 2> /dev/null         
    fi
    if test $? != 0; then
        $ECHO
        $ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
        $ECHO "Aborting"
        exit 1
    fi
done
$ECHO " done"

# how to run executables
CP_COMMAND="$PARA_PREFIX $BIN_DIR/cp.x $PARA_POSTFIX"
$ECHO
$ECHO "  running cp.x as: $CP_COMMAND"
$ECHO

# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/cp*
$ECHO " done"

cat > h2o.wannier.in << EOF
&CONTROL
  calculation   = "cp-wf",
  restart_mode  = "from_scratch",
  nstep         = 100,
  iprint        = 100,
  isave         = 100,
  dt            = 4.D0,
  etot_conv_thr = 1.D-16,
  ekin_conv_thr = 1.D-16,
  prefix        = "h2o_mol",
  pseudo_dir    = "$PSEUDO_DIR/",
  outdir        = "$TMP_DIR/",
/
&SYSTEM
  ibrav     = 1,
  celldm(1) = 20.0,
  nat       = 3,
  ntyp      = 2,
  ecutwfc   = 25.D0,
  nr1b = 10, nr2b = 10, nr3b = 10,
/
&ELECTRONS
  emass             = 350.D0,
  emass_cutoff      = 3.D0,
  ortho_eps         = 5.D-8,
  ortho_max         = 250,
  electron_dynamics = "damp",
  electron_damping  = 0.15,
/
&IONS
  ion_dynamics = "none",
/
&WANNIER
  adapt  = .true.,
  nsteps = 10,
/
ATOMIC_SPECIES
O  16.D0  O_US.van
H   2.D0  H_US.van
ATOMIC_POSITIONS (bohr)
O     10.0000    10.0000    10.000
H     11.7325     9.6757    10.000
H      9.6757    11.7325    10.000
EOF
$ECHO "  running the Wannier functions calculation...\c"
$CP_COMMAND < h2o.wannier.in > h2o.wannier.out
check_failure $?
$ECHO "  Wannier functions spreads...\c">> h2o.wannier.out
$ECHO "">> h2o.wannier.out
  tail -4 $TMP_DIR/h2o_mol.spr  >> h2o.wannier.out
$ECHO "  Wannier functions centers...\c">> h2o.wannier.out
$ECHO "">> h2o.wannier.out
  tail -4 $TMP_DIR/h2o_mol.wfc  >> h2o.wannier.out
$ECHO " done"

# calculation with an electric field
cat > h2o.efield.in << EOF
&CONTROL
  calculation   = "cp-wf",
  restart_mode  = "reset_counters",
  ndw           = 51
  nstep         = 200,
  iprint        = 100,
  isave         = 100,
  dt            = 6.D0,
  prefix        = "h2o_mol",
  pseudo_dir    = "$PSEUDO_DIR/",
  outdir        = "$TMP_DIR/",
/
&SYSTEM
  ibrav     = 1,
  celldm(1) = 20.0,
  nat       = 3,
  ntyp      = 2,
  ecutwfc   = 25.D0,
  nr1b = 10, nr2b = 10, nr3b = 10,
/
&ELECTRONS
  emass             = 350.D0,
  emass_cutoff      = 3.D0,
  ortho_eps         = 5.D-8,
  ortho_max         = 250,
  electron_dynamics = "verlet",
/
&IONS
  ion_dynamics = "verlet",
/
&WANNIER
  wf_efield = .TRUE.,
  wf_switch = .TRUE.,
  sw_len    = 100,
  efz1      = 0.5D-2,
/
ATOMIC_SPECIES
O  16.D0  O_US.van
H   2.D0  H_US.van
ATOMIC_POSITIONS (bohr)
O     10.0000    10.0000    10.000
H     11.7325     9.6757    10.000
H      9.6757    11.7325    10.000
EOF
$ECHO "  running the calculation with electric field...\c"
$CP_COMMAND < h2o.efield.in > h2o.efield.out
check_failure $?

$ECHO
$ECHO "$EXAMPLE_DIR: done"