#!/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 both pw.x and cp.x to calculate the"
$ECHO "minimum energy path (MEP) of the collinear proton transfer reaction:"
$ECHO "  2H2O2 <==> O2+H2O."

# set the needed environment variables
. ../environment_variables

# required executables and pseudopotentials
BIN_LIST="pw.x cp.x"
PSEUDO_LIST="H.vbc.UPF"

$ECHO
$ECHO "  executables directory: $BIN_DIR"
$ECHO "  pseudo directory:      $PSEUDO_DIR"
$ECHO "  temporary directory:   $TMP_DIR"
$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 "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
        $ECHO "Aborting"
        exit 1
    fi
done
$ECHO " done"

# how to run executables
PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
CP_COMMAND="$PARA_PREFIX $BIN_DIR/cp.x $PARA_POSTFIX"
$ECHO
$ECHO "  running pw.x as: $PW_COMMAND"
$ECHO "  running cp.x as: $CP_COMMAND"
$ECHO
 
# clean TMP_DIR
$ECHO "  cleaning $TMP_DIR...\c"
rm -rf $TMP_DIR/*
$ECHO " done"

# NEB calculation. Automatic choice of the climbing image.
cat > H2O2.in << EOF
&CONTROL
  calculation    = "neb",
  restart_mode   = "from_scratch",
  prefix         = "H2O2"
  outdir         = "$TMP_DIR",
  pseudo_dir     = "$PSEUDO_DIR",
/
&SYSTEM
  ibrav                  = 0,
  celldm(1)              = 1.D0,
  nat                    = 3,
  ntyp                   = 1,
  ecutwfc                = 20.0D0,
  ecutrho                = 100.0D0,
  nspin                  = 2,
  starting_magnetization = 0.5D0,
  occupations            = "smearing",
  degauss                = 0.03D0,
/
&ELECTRONS
  conv_thr    = 1.D-8,
  mixing_beta = 0.3D0,
/
&IONS
  ds                = 2.D0,
  opt_scheme        = "broyden",
  num_of_images     = 7,
  k_max             = 0.3D0,
  k_min             = 0.2D0,
  CI_scheme         = "auto",
  pot_extrapolation = "second_order",
  wfc_extrapolation = "second_order",
  path_thr          = 0.1D0,  
/
ATOMIC_SPECIES
H  1.00794  HUSPBE.RRKJ3
O  15.9994  O.LDA.US.RRKJ3.UPF
ATOMIC_POSITIONS { bohr }
first_image
O    		0.0000         0.0000    	0.0000
O   	     -0.5600         1.3200    	0.0000
H          -0.5914        -0.5985        -0.4622
H           0.0314         1.9185         0.4622
O    		5.0000         0.0000    	0.0000
O   	      4.5600         1.3200    	0.0000
H           4.5914        -0.5985        -0.4622
H           5.0314         1.9185         0.4622
last_image
O           0.00000        0.00000        0.00000
H           0.96000        0.00000        0.00000
H          -0.32000       -0.87950       -0.21400
O           5.00000        0.00000        0.00000
H           5.96000        0.00000        0.00000
H          -4.32000       -0.87950       -0.21400
O           6.05000       11.10700       10.62100
O           5.26200       10.18300       10.83900
K_POINTS { gamma }
CELL_PARAMETERS { cubic }
  12.00000  0.00000  0.00000
   0.00000 12.00000  0.00000
   0.00000  0.00000 12.00000
EOF
$ECHO "  running the NEB calculation with pw.x for 2H2O2 => O2+H2O...\c"
$PW_COMMAND < H2O2.in > H2O2.out
$ECHO " done"

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

# NEB calculation. Climbing image is not used
cat > symmetric_H2O2.in << EOF
&CONTROL
  calculation    = "neb",
  restart_mode   = "from_scratch",
  prefix         = "symmetric_H2O2"
  outdir         = "$TMP_DIR",
  pseudo_dir     = "$PSEUDO_DIR",
/
&SYSTEM
  ibrav                  = 0,
  celldm(1)              = 1.D0,
  nat                    = 3,
  ntyp                   = 1,
  ecutwfc                = 20.0D0,
  ecutrho                = 100.0D0,
  nspin                  = 2,
  starting_magnetization = 0.5D0,
  occupations            = "smearing",
  degauss                = 0.03D0,
/
&ELECTRONS
  conv_thr    = 1.D-8,
  mixing_beta = 0.3D0,
/
&IONS
  ds                = 2.D0,
  opt_scheme        = "broyden",
  num_of_images     = 8,
  k_max             = 0.3D0,
  k_min             = 0.2D0,
  pot_extrapolation = "second_order",
  wfc_extrapolation = "second_order",
  path_thr          = 0.2D0,
/
ATOMIC_SPECIES
H  1.00794  HUSPBE.RRKJ3
O  15.9994  O_US.van
ATOMIC_POSITIONS { bohr }
first_image
O    		0.0000         0.0000    	0.0000
O   	     -0.5600         1.3200    	0.0000
H          -0.5914        -0.5985        -0.4622
H           0.0314         1.9185         0.4622
O    		5.0000         0.0000    	0.0000
O   	      4.5600         1.3200    	0.0000
H           4.5914        -0.5985        -0.4622
H           5.0314         1.9185         0.4622
last_image
O           0.00000        0.00000        0.00000
H           0.96000        0.00000        0.00000
H          -0.32000       -0.87950       -0.21400
O           5.00000        0.00000        0.00000
H           5.96000        0.00000        0.00000
H          -4.32000       -0.87950       -0.21400
O           6.05000       11.10700       10.62100
O           5.26200       10.18300       10.83900
K_POINTS { gamma }
CELL_PARAMETERS { cubic }
  12.00000  0.00000  0.00000
   0.00000 12.00000  0.00000
   0.00000  0.00000 12.00000
EOF
$ECHO "  running the NEB calculation with pw.x for symmetric 2H2O2 => O2+H2O...\c"
$PW_COMMAND < symmetric_H2+H.in > symmetric_H2O2.out
$ECHO " done"

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

# the name of the restart file is changed in order to conform to the
# prefix of the new run 
# the restart file asymmetric_H2+H.neb is modified (second row)
# since a new simulation (from the old path) is started
cat symmetric_H2O2.path | \
    awk '{if(NR==2){printf"   0\n"}; if(NR!=2){print}}' > asymmetric_H2O2.path

# NEB calculation. The image that has to climb is manually chosen
cat > asymmetric_H2O2.in << EOF
&CONTROL
  calculation    = "neb",
  restart_mode   = "restart",
  prefix         = "asymmetric_H2O2"
  outdir         = "$TMP_DIR",
  pseudo_dir     = "$PSEUDO_DIR",
/
&SYSTEM
  ibrav                  = 0,
  celldm(1)              = 1.D0,
  nat                    = 3,
  ntyp                   = 1,
  ecutwfc                = 20.0D0,
  ecutrho                = 100.0D0,
  nspin                  = 2,
  starting_magnetization = 0.5D0,
  occupations            = "smearing",
  degauss                = 0.03D0,
/
&ELECTRONS
  conv_thr    = 1.D-8,
  mixing_beta = 0.3D0,
/
&IONS
  ds                = 2.D0,
  opt_scheme        = "broyden",
  num_of_images     = 8,
  k_max             = 0.3D0,
  k_min             = 0.2D0,
  CI_scheme         = "manual",
  pot_extrapolation = "second_order",
  wfc_extrapolation = "second_order",
  path_thr          = 0.05D0,
/
CLIMBING_IMAGES
 5
ATOMIC_SPECIES
H  1.00794  HUSPBE.RRKJ3
O  15.9994  
ATOMIC_POSITIONS { bohr }
first_image
O    		0.0000         0.0000    	0.0000
O   	     -0.5600         1.3200    	0.0000
H          -0.5914        -0.5985        -0.4622
H           0.0314         1.9185         0.4622
O    		5.0000         0.0000    	0.0000
O   	      4.5600         1.3200    	0.0000
H           4.5914        -0.5985        -0.4622
H           5.0314         1.9185         0.4622
last_image
O           0.00000        0.00000        0.00000
H           0.96000        0.00000        0.00000
H          -0.32000       -0.87950       -0.21400
O           5.00000        0.00000        0.00000
H           5.96000        0.00000        0.00000
H          -4.32000       -0.87950       -0.21400
O           6.05000       11.10700       10.62100
O           5.26200       10.18300       10.83900
K_POINTS { gamma }
CELL_PARAMETERS { cubic }
  12.00000  0.00000  0.00000
   0.00000 12.00000  0.00000
   0.00000  0.00000 12.00000
EOF
$ECHO "  running the NEB calculation with pw.x for asymmetric 2H2O2 => O2+H2O...\c"
$PW_COMMAND < asymmetric_H2+H.in > asymmetric_H2O2.out
$ECHO " done"

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

# NEB calculation. Automatic choice of the climbing image.
cat > H2O2-cp.in << EOF
&CONTROL
  dt             = 10.D0,
  calculation    = "neb",
  restart_mode   = "from_scratch",
  prefix         = "H2O2-cp"
  outdir         = "$TMP_DIR/",
  pseudo_dir     = "$PSEUDO_DIR/",
  iprint         = 50,
/
&SYSTEM
  ibrav     = 0,
  celldm(1) = 1.D0,
  nat       = 3,
  ntyp      = 1,
  nspin     = 2,
  nelec     = 3,
  nelup     = 2,
  neldw     = 1,
  ecutwfc   = 20.D0,
  ecutrho   = 100.D0,
  nr1b      = 12,
  nr2b      = 12,
  nr3b      = 12,
/
&ELECTRONS
  electron_dynamics = "damp",
  electron_damping  = 0.2D0,
  emass             = 250.D0,
  ortho_max         = 250,
/
&IONS
  ds                = 2.D0,
  opt_scheme        = "broyden",
  num_of_images     = 7,
  k_max             = 0.3D0,
  k_min             = 0.2D0,
  CI_scheme         = "auto",
  path_thr          = 0.1D0,  
/
ATOMIC_SPECIES
H  1.00794  HUSPBE.RRKJ3
O  15.9994  
ATOMIC_POSITIONS { bohr }
first_image
O    		0.0000         0.0000    	0.0000
O   	     -0.5600         1.3200    	0.0000
H          -0.5914        -0.5985        -0.4622
H           0.0314         1.9185         0.4622
O    		5.0000         0.0000    	0.0000
O   	      4.5600         1.3200    	0.0000
H           4.5914        -0.5985        -0.4622
H           5.0314         1.9185         0.4622
last_image
O           0.00000        0.00000        0.00000
H           0.96000        0.00000        0.00000
H          -0.32000       -0.87950       -0.21400
O           5.00000        0.00000        0.00000
H           5.96000        0.00000        0.00000
H          -4.32000       -0.87950       -0.21400
O           6.05000       11.10700       10.62100
O           5.26200       10.18300       10.83900
K_POINTS { gamma }
CELL_PARAMETERS { cubic }
  12.00000  0.00000  0.00000
   0.00000 12.00000  0.00000
   0.00000  0.00000 12.00000
EOF
$ECHO "  running the NEB calculation with cp.x for H2O2 => O2+H2O...\c"
$CP_COMMAND < H2O2-cp.in > H2O2-cp.out
$ECHO " done"

$ECHO
$ECHO "$EXAMPLE_DIR: done"