#!/bin/sh

###############################################################################
#  A Sort-Of Generalized PWscf NSCF Script for AF Transition Metal Monoxides  #
#                             P. M. Grant                                     #
#                           February, 2008                                    #
###############################################################################

###############################################################################
#    Element and Structure-Specific Constants & Data (Meant To Be Edited)     #
#                              (!con cuidado!)                                #
#           Be especially careful to link to the right PP folder              #
#     Run this with proper 'kpoints_array' prior to each of the below:        #
#          1) Density of States and Projected DOS                             #
#          2) Band Structure                                                  #
#          3) Fermi Surface                                                   #
#           Change Log:                                                       #
#              04/20/08 - Include ability to select Cell Parameter file       #
#              05/11/08 - Change above to refer to LAT_BZ directory           #
#                         Naming Convention:                                  #
#                          LAT_BZ_ab_c                                        #
#                          where ab is ab lattice constant and c=c/a ratio    #
#              09/16/10   Added explicitly naming ab_c to append to LAT_BZ_   #
#              06/21/11   Added a "Hubbard Bank" to temporarily store Hubbard #
#                         parameters from the &system namelist in order to    #
#                         effect a true U=0 calculation.  This "bank"         #
#                         is placed below an exit command at the end of this  #
#                         script.                                             #
#                                                                             #
#         THIS SCRIPT IS ONLY FOR ANTIFERROMAGNETIC COPPER MONOXIDE!          #
###############################################################################
#
# Conversion Factors (If Needed)
Ryd_eV=13.6057
Bohr_Ang=0.529
e_Charge=1.60217646e-19
e_Mass=9.10938188e-31
Pi=3.14159265

# Constants, Variables & Symbols Setup
TM='Cu';	A='O'
TM_m='63.456';	A_m='15.9994'
TM1=$TM'1'
TM2=$TM'2'
A1=$A'1'
A2=$A'2'
tot_charge='0.0'
#tot_charge='0.05'
#tot_charge='1.0'
#tot_charge='0.3'
#tot_charge='0.5'
#tot_charge='0.6'
#tot_charge='-0.6'
#tot_charge='-1.0'
#tot_charge='-0.3'
#tot_charge='-0.5'
#ecutwfc_size='60'  
ecutwfc_size='40'  
#ecutrho_size=`echo "scale=5; 4*$ecutwfc_size" | "bc"`
ecutrho_size=`echo "scale=5; 8*$ecutwfc_size" | "bc"`
prefix=af-$TM$A
nbnd='20'

# PseudoPots List
# Copper
TM_pseudopot='Cu.pz-3d9_4s2-rrkjus.UPF';	
#TM_pseudopot='029-Cu-gpw-n-campos.uspp.UPF';	
#TM_pseudopot='CuLDA.RRKJ3.PAW'		
#TM_pseudopot='Cu.pz-d-rrkjus.UPF';
#TM_pseudopot='Cu.pz-n-van_ak.UPF';
#TM_pseudopot='Cu.blyp-n-van_ak.UPF';
#TM_pseudopot='Cu.pw91-n-van_ak.UPF';
# Oxygen
#A_pseudopot='008-O-ca--vgrp.uspp.UPF'
#A_pseudopot='OLDA.RRKJ3.PAW'
#A_pseudopot='O.pz-rrkjus.UPF'
#A_pseudopot='O.LDA.US.RRKJ3.UPF'
#A_pseudopot='O.vdb.UPF'
A_pseudopot='O.pz-van_ak.UPF'
#A_pseudopot='O.pw91-van_ak.UPF'

# Magnetism Related Parameters
sm_A='0.0'
sm_TM='0.5'
H_U='1.d-7'
#H_U='1.0'
#H_U='2.0'
#H_U='3.0'
#H_U='4.0'
#H_U='5.0'
#H_U='6.0'
#H_U='7.0'
#H_U='8.0'
#H_U='9.0'
#H_U='10.0'



###############################################################################
# Housekeeping and Various Setup Stuff                                        #
###############################################################################
sysinfo=`uname '-a'`
kubuntu=`cat '/etc/issue'` 
now=`date`
cd `echo $0 | sed 's/\(.*\)\/.*/\1/'`
home_dir=`pwd`
cd ..
home_dir_up1=`pwd`
cd $home_dir
#pwscf_distro='espresso-3.2.3'
#pwscf_distro='espresso4.0cvs'
#pwscf_distro='espresso-4.0'
pwscf_distro='espresso-4.2.1'
#pwscf_distro='O-sesame'
pwscf_exec='/home/pmpgrant/'$pwscf_distro'/PW/pw.x'
home_pot=$home_dir_up1'/PSEUDOPOTS'
home_bz_kvects=$home_dir_up1'/BZ_&_K-VECTS'
#ab_c='3905_60'
#ab_c='3905_136'
ab_c='3905_1'
LAT_BZ_dir=$home_bz_kvects/'LAT_BZ_'$ab_c

# k-points Arrays
#kpoints_array=$LAT_BZ_dir/'k-paths_100.pwscf'
#kpoints_array='k-MP_4'42
#kpoints_array=$home_bz_kvects/'k-MP_16'
kpoints_array=$LAT_BZ_dir/kvecs_$prefix-fs
#kpoints_array=kvecs_$prefix-fs

# cat to $prefix.nscf.data file
printf "$0 \n" > $prefix.nscf.data
cat >> $prefix.nscf.data << EOF
  A PWscf LDA+U script to calculate the electronic properties of $prefix
  Starting at $now
*******************************************************************************
System & OS Info
---------------------------------------------------------------------
  $sysinfo
  $kubuntu
*******************************************************************************
Element, Structure and Run Parameter Editable Options (!con cuidado!)
---------------------------------------------------------------------
  TM = $TM		A = $A
  TM_mass = $TM_m	A_mass = $A_m
  LAT_BZ_dir = $LAT_BZ_dir	
  MP_mesh = $mp_size
  AB_C = $ab_c
  TM_PP = $TM_pseudopot	A_PP = $A_pseudopot
  Starting Magnetizations:
  A = $sm_A		TM = $sm_TM
  Hubbard_U = $H_U
 
*******************************************************************************
Housekeeping
------------
  PWscf Distro = $pwscf_distro
  My Home      = $home_dir
  PP Home      = $home_pot
  PW Exec      = $pwscf_exec

  K-VECTS      = $home_bz_kvects
  AB_C         = $ab_c
  LAT_BZ       = $LAT_BZ_dir
  K-POINTS     = $kpoints_array
EOF

###############################################################################
# Generate NSCF Run Stuff                                                      #
###############################################################################
cat > $prefix.nscf.in << EOF
 &control
  title = '$prefix.nscf',		calculation = 'nscf',
  restart_mode = 'from_scratch',	prefix = '$prefix',
  pseudo_dir = '$home_pot',
  outdir = '$home_dir',
  tprnfor = .false.,			tstress = .false.,	
  verbosity = 'minimal',
 /
 &system
  ibrav = 0,			celldm(1)=0,
  nat = 4,			ntyp = 4,
  ecutwfc = $ecutwfc_size		ecutrho = $ecutrho_size,
  tot_charge = $tot_charge
  smearing = 'methfessel-paxton',	degauss = 0.02,	occupations = 'smearing',
  nspin = 2,			nbnd=$nbnd			
  starting_magnetization(1) = $sm_TM,
  starting_magnetization(2) = -$sm_TM,
  starting_magnetization(3) = $sm_A,
  starting_magnetization(4) = $sm_A,
 /
 &electrons
 /
CELL_PARAMETERS
EOF
cat >> $prefix.nscf.in < $LAT_BZ_dir/CP_RS.uvw
cat >> $prefix.nscf.in << EOF
ATOMIC_SPECIES
    $TM1     $TM_m         $TM_pseudopot
    $TM2     $TM_m         $TM_pseudopot
    $A1      $A_m          $A_pseudopot
    $A2      $A_m          $A_pseudopot
ATOMIC_POSITIONS {crystal}
    $TM1     0.0     0.0     0.0
    $TM2     0.5     0.5     0.5 
    $A1      0.25    0.25    0.25
    $A2      0.75    0.75    0.75
EOF
cat $kpoints_array >> $prefix.nscf.in 
# NB...The k-points input for FS calcs must be specified as below and
# KPOINTS inserted above the EOF
#cat $kpoints_array >> $prefix.nscf.in 
#exit
###############################################################################
# Do NSCF Run                                                                 #
###############################################################################
cat >> $prefix.nscf.data << EOF
  Beginning $prefix NSCF...
  Starting at `date`
*******************************************************************************
$prefix.nscf.in follows...
----------------------------------
EOF

cat >> $prefix.nscf.data < $prefix.nscf.in
cat >> $prefix.nscf.data << EOF

*******************************************************************************
  Running $prefix NSCF...
EOF

# Exit Point for Testing
#exit

$pwscf_exec < $prefix.nscf.in  | tee $prefix.nscf.out
cat >> $prefix.nscf.data << EOF
  Finished...at `date`
  Collecting key data
  `grep 'kinetic-energy cutoff     ' $prefix.nscf.out`
  `grep 'charge density cutoff     ' $prefix.nscf.out`
  `grep 'convergence threshold     ' $prefix.nscf.out`
  `grep 'mixing beta               ' $prefix.nscf.out`
  `grep 'number of electrons       ' $prefix.nscf.out`
  `grep 'number of Kohn-Sham states' $prefix.nscf.out`
  `grep 'number of iterations used ' $prefix.nscf.out`
  `grep 'Exchange-correlation      ' $prefix.nscf.out`

  # k-points = `grep 'number of k points=' $prefix.nscf.out | cut -c 25-30`
  Fermi Energy = `grep 'Fermi' $prefix.nscf.out | cut -c 28-36`
  Total Energy = `grep '!' $prefix.nscf.out | cut -c 35-48`

  Total CPU Time = `grep 'CPU time' $prefix.nscf.out | cut -c 20-30`
  Total Wall Time = `grep 'WALL time' $prefix.nscf.out | cut -c 40-53`

  Iteration Summary
  -----------------
`grep 'estimated scf accuracy  ' $prefix.nscf.out`
`grep 'total magnetization     ' $prefix.nscf.out`
`grep 'absolute magnetization  ' $prefix.nscf.out`
`grep 'Tr' $prefix.nscf.out`
`grep 'nsum =' $prefix.nscf.out`

  d-orbital occupations
  ---------------------
`grep -n 'eigenvalues:' $prefix.nscf.out`
EOF
exit

# Hubbard parameter &system "bank."  To use, copy below last "starting mag     #
 lda_plus_u =.true.,
  Hubbard_U(1) = $H_U,	Hubbard_U(2) = $H_U,