!#/bin/ksh
# ---------------------------------------------------------------
#   batch file to run over several diplaced atoms
#   example for PWSCF; MgO crystal 
# ---------------------------------------------------------------
if test -f POSCAR.phon 
then 
  echo "POSCAR.phon exist remove before starting batchfile"
  exit 1
fi
cp POSCAR POSCAR.phon 
#
#  loop over all exited ions
#
j=0
for i in \
# here put the content of the DISP file
 "   1  0.00477204 -0.00477204  0.00477204 " \
 "   9  0.00477204 -0.00477204  0.00477204 "
do
j=`expr $j + 1`
echo "run number $j"

#
# use awk to displace atoms (somewhat tricky)
#
awk '
/SUBSIT/  { npos = $2 ; 
            x=$3; y=$4 ; z=$5; }
!/SUBSIT/ { line=line+1 
  if (line-7 == npos) {
    printf "%12.9f %12.9f %12.9f\n", $1+x, $2+y, $3+z
  } 
  else if (line > 7) {
    printf "%12.9f %12.9f %12.9f\n", $1, $2, $3
  } 
  else 
    print 
} ' >POSCAR <<!
SUBSIT $i
`cat POSCAR.phon`
!

natoms=`wc -l POSCAR |awk '{print $1-7}'`
alat=`echo 7.92 $natoms |awk '{print $1*($2/2)**.3333333333333}'`
awk -v n=$natoms '{if(NR>7&&NR<=n/2+7){print "Mg",$0}}{if(NR>n/2+7&&NR<=n+7){print "O",$0}}' POSCAR > tmp
cat > pw.d << **
MgO
NaCl
&control
    calculation='scf'
    restart_mode='from_scratch',
    prefix='mgo'
    pseudo_dir='../../pseudo/'
    tprnfor=.true.
/
&system
    ibrav = 0, celldm(1) = $alat
    nat= $natoms, ntyp= 2,
    ecutwfc = 30, ecutrho=240
/
&electrons
!    startingwfc='file'
!    startingpot='file'
    conv_thr = 1.0e-8
    mixing_beta = 0.7
    mixing_mode = 'plain'
/
ATOMIC_SPECIES
 Mg  1.0 Mg.pz-n-vbc.UPF
 O   1.0 O.pz-rrkjus.UPF
ATOMIC_POSITIONS { crystal }
**
cat tmp >> pw.d
cat >> pw.d << ***
K_POINTS automatic
4 4 4 0 0 0
CELL_PARAMETERS { cubic }
   0.500000000000000   0.500000000000000   0.000000000000000
   0.000000000000000   0.500000000000000   0.500000000000000
   0.500000000000000   0.000000000000000   0.500000000000000
***
# here the call to the PWSCF executable to compute forces 
#../../PW/pw.x < pw.d > out.$j 
#
# use awk to extract forces 
#
echo $i >DYNMAT.$j
awk -v n=$natoms 'BEGIN{start=1000000}/Forces acting/{start=NR}
  { if(NR>start+1&&NR<=start+n+1) { printf "%14.8f  %14.8f  %14.8f\n",$7*25.711,$8*25.711,$9*25.711  }}
' out.$j >>DYNMAT.$j

done

# restore POSCAR file
mv POSCAR.phon POSCAR

echo $j > FORCES
i=1
while test $i -le $j
do 
cat DYNMAT.$i >> FORCES
i=`expr $i + 1`
done