#!/bin/sh
#

for ls in 10.0 10.2 10.4 10.6 10.8 10.9 11.0 11.2 11.4 11.6 12.0
do

if [ ! -d results_Fe_sc8_ucalc.$ls ]; then
   mkdir  results_Fe_sc8_ucalc.$ls
fi
if [ ! -d temp ]; then
   mkdir  temp
fi
if [ ! -d temp1 ]; then
   mkdir  temp1
fi
rm -rf temp/* temp1/*

# unperturbed self-consistent calculation
rm -f fe.scf.in
cat > fe.scf.in << EOF
 &control
    pseudo_dir = '../pseudi'
    outdir='./temp'
    restart_mode='from_scratch'
    verbosity = 'high',
    prefix='fe',
 /
 &system
    ibrav=  1, celldm(1)=$ls, nat=  16, ntyp= 2,
    ecutwfc = 40.0, ecutrho = 420.0,
    occupations='smearing', smearing='mp', degauss=0.005,
    nbnd = 100,
    nspin=2, 
    starting_magnetization(1)= 0.6
    starting_magnetization(2)= 0.6
    lda_plus_u = .true.,
    U_projection_type = 'ortho-atomic',
    Hubbard_U(1)= 1.d-20
    Hubbard_U(2)= 1.d-20
 /
 &electrons
    mixing_beta = 0.5 
    conv_thr =  1.0d-9,
 /
ATOMIC_SPECIES
 Fe1 1.  Fe.pbe-nd-rrkjus.UPF
 Fe  1.  Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
Fe1  0.00000   0.00000   0.00000
Fe   0.50000   0.00000   0.00000
Fe   0.00000   0.50000   0.00000
Fe   0.00000   0.00000   0.50000
Fe   0.00000   0.50000   0.50000
Fe   0.50000   0.50000   0.00000
Fe   0.50000   0.00000   0.50000
Fe   0.50000   0.50000   0.50000
Fe   0.25000   0.25000   0.25000
Fe   0.75000   0.25000   0.25000
Fe   0.25000   0.75000   0.25000
Fe   0.25000   0.25000   0.75000
Fe   0.25000   0.75000   0.75000
Fe   0.75000   0.75000   0.25000
Fe   0.75000   0.25000   0.75000
Fe   0.75000   0.75000   0.75000
K_POINTS automatic
8 8 8 0 0 0
EOF
./pw.x < fe.scf.in  > results_Fe_sc8_ucalc.$ls/fe.scf.out

ethr=`grep ethr results_Fe_sc8_ucalc.$ls/fe.scf.out |tail -1 |awk '{print $3}'`

cp -r temp/* temp1/

for a in -0.1 0.1 -0.05 0.05 0.0
do

rm -rf temp/*
cp -r temp1/* temp/

# perturbed calculations (finite alpha)
rm -f fe.scf.in
cat > fe.scf.in << EOF
 &control
    pseudo_dir = '../pseudi',
    outdir='./temp'
    restart_mode='from_scratch'
    verbosity = 'high',
    prefix='fe',
 /
 &system
    ibrav=  1, celldm(1)=$ls, nat=  16, ntyp= 2,
    ecutwfc = 40.0, ecutrho = 420.0,
    occupations='smearing', smearing='mp', degauss=0.005,
    nbnd = 100,
    nspin=2, 
    starting_magnetization(1)= 0.6
    starting_magnetization(2)= 0.6
    lda_plus_u = .true.,
    U_projection_type = 'ortho-atomic',
    Hubbard_U(1)= 1.d-20
    Hubbard_U(2)= 1.d-20
    Hubbard_alpha(1)= $a
 /
 &electrons
    startingpot = 'file'
    startingwfc = 'file'
    diago_thr_init = $ethr
    mixing_beta = 0.3 
    conv_thr =  1.0d-9,
 /
ATOMIC_SPECIES
 Fe1 1.  Fe.pbe-nd-rrkjus.UPF
 Fe  1.  Fe.pbe-nd-rrkjus.UPF
ATOMIC_POSITIONS crystal
Fe1   0.00000   0.00000   0.00000
Fe   0.50000   0.00000   0.00000
Fe   0.00000   0.50000   0.00000
Fe   0.00000   0.00000   0.50000
Fe   0.00000   0.50000   0.50000
Fe   0.50000   0.50000   0.00000
Fe   0.50000   0.00000   0.50000
Fe   0.50000   0.50000   0.50000
Fe   0.25000   0.25000   0.25000
Fe   0.75000   0.25000   0.25000
Fe   0.25000   0.75000   0.25000
Fe   0.25000   0.25000   0.75000
Fe   0.25000   0.75000   0.75000
Fe   0.75000   0.75000   0.25000
Fe   0.75000   0.25000   0.75000
Fe   0.75000   0.75000   0.75000
K_POINTS automatic
8 8 8 0 0 0
EOF
./pw.x < fe.scf.in  > results_Fe_sc8_ucalc.$ls/fe.pert_$a.out

done

done