#!/bin/bash

# ---------------------------
# set the name of the job
#$ -N   test

# request between 2 and 8 slots
#$ -pe mpich 16

# run in infini.q
#$ -q ether.q

export MPICH_PROCESS_GROUP=no


# ---------------------------
# set up the mpich version to use
# ---------------------------
# load the module
. /etc/profile.d/modules.sh
module purge
module add ofed
module add sge
module add intel/fce/10.1.015
module add intel/mkl/10.0.1.014
module add openmpi/intel/64/1.2.6
#----------------------------
# set up the parameters for qsub
# ---------------------------

#  Mail to user at beginning/end/abort/on suspension
# #$ -m beas
#  By default, mail is sent to the submitting user 
#  Use  $ -M username    to direct mail to another userid 

# Execute the job from the current working directory
# Job output will appear in this directory
#$ -cwd
#   can use -o dirname to redirect stdout 
#   can use -e dirname to redirect stderr

#  Export these environment variables
#$ -v MPI_HOME 
#$ -v LD_LIBRARY_PATH 
echo $LD_LIBRARY_PATH

# ---------------------------
# run the job
# ---------------------------


###############################################################
# calculation of 3D force field  - for AFM imaging calculations
###############################################################
# syntax:
#          run_AFM
###############################################################
# IMPORTNAT!  Before running the script:
#    (1) run SIESTA for the initial geometry
#    (2) prepare INPUT_FILE.fdf with the relaxed geometry
#        using e.g. option Ss in tetr
#    (3) Set SystemLabel in INPUT_FILE.fdf to INPUT_FILE
#    (4) Remove UseSaveData and DM.UseSaveDM from the file
#    (5) change entry_count to the count number at which to
#        restart the calculation if stopped with an ERROR
#        during the run of the script
###############################################################
# tetr refers to L. Kantorovitch's geometry creation/modification software
entry_count=0
#... tip atoms (to be moved at each step)
tip=( 10,15-17 )
zt0=(4.00000)
#... elementary displacement vectors within the surface
a1=(0.10 0.00) 
a2=(0.00 0.10)
#... templates for tetr
cd=(4  INPUT_FILE)
Sp=

#echo "" > tmp

echo > ttmp1
for word in  ${cd[0]} ${cd[1]} Ss ${cd[1]} Q 
do
   echo $word >> ttmp1
done
if [ $entry_count -eq 0 ]; then cp -p INPUT_FILE.fdf "INPUT_FILE_ini.fdf"; fi
echo > ttmp
for word in  ${cd[0]} ${cd[1]} $Sp M T $tip Mv 4 Displ S ${cd[0]} ${cd[1]}  Q Q
do
   echo $word >> ttmp
done

#[1]... make loops over displcements along 3 directions
count=0
for n1 in  0
do
   for n2 in    2
   do

#________ restore the initial geometry
     cp -p INPUT_FILE_ini.fdf "INPUT_FILE.fdf"
     echo '<============= n1= '$n1' and n2= '$n2' =========>'
     if [ ! -d $n1.$n2 ]; then mkdir -p "$n1.$n2"; fi

#________ calculate the actual displ. vector
     if [ $n1 -ne 0 -o $n2 -ne 0 ]; then
	 d1=`echo $n1 $n2 ${a1[0]} ${a1[1]} ${a2[0]} ${a2[1]} | awk '{printf("%.13f", $1*$3 + $2*$5)}'`
	 d2=`echo $n1 $n2 ${a1[0]} ${a1[1]} ${a2[0]} ${a2[1]} | awk '{printf("%.13f", $1*$4 + $2*$6)}'`
         displ=`printf "%.13f, %.13f, " $d1  $d2`
         echo '========== Displacement = ' $displ '=========='

#_____________ create the input file with the displaced tip
         sed -e 's/Displ/'${displ[0]}' '${displ[1]}' 0.0/g' ttmp > tinp
         cp -p tinp "$n1.$n2"/.
         tetr < tinp >& tout

#_____________ check for possible errors
         if [ `cat tout| grep -i err| wc -l` -ne 0 ]; then
             mv -f tout "$n1.$n2"/.
             echo 'STOP: ERROR in TETR when Mv; check file ' $n1.$n2/tout
             echo 'stopped at count = '$count
             exit
         fi
     fi
     echo 'DM.UseSaveDM F' >> INPUT_FILE.fdf
     cp -p INPUT_FILE.fdf "$n1.$n2"/

#[2]........ loop over vert. down displs (may be noneven)
     z=0
     zt=$zt0
     if [ ! -e $n1.$n2/tip_forces_z ]; then echo > $n1.$n2/tip_forces_z; fi
     for zd in  0.0  -1.00  -1.00  -1.00  -1.00  -1.00  -1.00  -0.5  -0.5 
     do
#____________ if restart: update input file to the nearest from the previous run
       count=$(($count + 1))
       if [ $entry_count -ne 0 -a $entry_count -eq $count]; then
	   if [ $z -gt 0 ]; then cp -p $n1.$n2/$z/INPUT_FILE-relaxed.fdf "INPUT_FILE.fdf"; fi
       fi
#____________ calculate down displacement and vertical coordinate
       z=$(($z + 1))
       echo '<== '$count' ............$z= '$z' --> zd= '$zd' .........==>'
       zt=`echo $zt $zd | awk '{printf("%.13f", $1+$2)}'`
       if [ ! -d $n1.$n2/$z ]; then mkdir -p $n1.$n2/$z; fi

#_______ ______ move the tip by the displcement
       if [ $z -ne 1 ]; then
	   sed -e 's/Displ/0.0 0.0 '$zd'/g' ttmp > tinp
	   cp -p tinp $n1.$n2/$z/.
	   tetr < tinp >& tout

#_____________ check for possible errors
	   if [ `cat tout| grep -i err| wc -l` -ne 0 ]; then
	       mv tout $n1.$n2/$z/.
	       echo 'STOP: ERROR in TETR when Mv; check file ' $n1.$n2/$z/tout
	       echo 'stopped at count = '$count
	       exit
	   fi
	   echo 'DM.UseSaveDM T' >> INPUT_FILE.fdf
       fi
       echo '============ Displacement = 0.0 0.0 ' $zd ' and height = ' $zt
       cp -p INPUT_FILE.fdf  $n1.$n2/$z/.

#______________ run the code
       if [ $entry_count -le $count ]; then
	   cvos-launcher $MPI_HOME/bin/mpirun  -np $NSLOTS -hostfile $TMPDIR/machines   ~/bin/siesta_2.0.1_vdW_ethern <INPUT_FILE.fdf | tee OUTCAR

#______________ create the final fdf file
#	   if( -e INPUT_FILE.XV ); then 
	   tetr < ttmp1 >& tout
	   if [ `cat tout| grep -i err| wc -l` -ne 0 ]; then
	       mv -f tout $n1.$n2/$z/.
	       echo 'STOP: ERROR in TETR when Ss; check file ' $n1.$n2/$z/tout
	       echo 'stopped at count = '$count
	       exit
	   fi
	   mv -f INPUT_FILE.XV $n1.$n2/$z/.
	   cp -p INPUT_FILE-relaxed.fdf $n1.$n2/$z/.
	   mv -f INPUT_FILE-relaxed.fdf "INPUT_FILE.fdf"
#               else
#                   ~/bin/cp INPUT_FILE.fdf $n1.$n2/$z/INPUT_FILE-relaxed.fdf
#               fi
#_______________ get forces, tip force and check if DFT finshed relaxation
	   $HOME_TETR/get_forces F OUTCAR 10,15-17 > $n1.$n2/$z/tip_force
	   mv -f OUTCAR $n1.$n2/$z/.
	   if [ `cat $n1.$n2/$z/tip_force| grep -i error| wc -l` -ne 0 ]; then
	       echo 'STOP: ERROR in get_forces; check file ' $n1.$n2/$z/tip_force
	       echo 'STOP: ERROR in get_forces; check file ' $n1.$n2/$z/OUTCAR
	       echo 'stopped at count = '$count
	       exit
	   fi
	   tf=`cat $n1.$n2/$z/tip_force`
	   echo $z $zt $tf[1] $tf[2] $tf[3] >> $n1.$n2/tip_forces_z
	   mv -f tmp.f $n1.$n2/$z/forces
       fi
     done
   done
done
