#!/bin/sh

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

###############################################################################
#    Element and Structure-Specific Constants & Data (Meant To Be Edited)     #
#                              (!con cuidado!)                                #
#           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    #
#              06/01/08   Added flag to compute DOS link la2F                 #
#                         for e-p coupling                                    #
#              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'
mp_size='4'
#mp_size='8'
tot_charge='0.0'
#tot_charge='0.05'
#tot_charge='1.0'
#tot_charge='0.3'
#tot_charge='-0.3'
#tot_charge='0.5'
#tot_charge='0.6'
#tot_charge='-0.6'
#tot_charge='-1.0'
#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"`
ep_la2F='.false.'
#ep_la2F='.true.'
prefix=af-$TM$A

# 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'
#U_p_t='atomic'
#U_p_t='file'
U_p_t='norm-atomic'
#s_ns_121=	s_ns_221=	s_ns_321= 	s_ns_421=	s_ns_521=
#s_ns_112=	s_ns_212=	s_ns_312= 	s_ns_412=	s_ns_512=
s_ns_121=	s_ns_221=1	s_ns_321=1	s_ns_421=1	s_ns_521=1
s_ns_112=	s_ns_212=1	s_ns_312=1	s_ns_412=1	s_ns_512=1

# Run Control Paramenters
electron_maxstep='100'

###############################################################################
# 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
#echo $LAT_BZ_dir
#exit

# cat to $prefix.scf.data file
printf "$0 \n" > $prefix.scf.data
cat >> $prefix.scf.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
  U_Projection_Type = $U_p_t
  Starting_ns_eigenvalues:
	121	221	321	421	521
	$s_ns_121	$s_ns_221	$s_ns_321	$s_ns_421	$s_ns_521	
	112	212	312	412	512
	$s_ns_112	$s_ns_212	$s_ns_312	$s_ns_412	$s_ns_512
  Maximum Number of SCF Iterations = $electron_maxstep	

*******************************************************************************
Housekeeping
------------
  PWscf Distro = $pwscf_distro
  My Home      = $home_dir
  PP Home      = $home_pot
  PW Exec      = $pwscf_exec
EOF

###############################################################################
# Generate SCF Run Stuff                                                      #
###############################################################################
cat > $prefix.scf.in << EOF
 &control
  title = '$prefix.scf',		calculation = 'scf',
  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',
  la2F = $ep_la2F
  nspin = 2,			
  starting_magnetization(1) = $sm_TM,
  starting_magnetization(2) = -$sm_TM,
  starting_magnetization(3) = $sm_A,
  starting_magnetization(4) = $sm_A,
 /
 &electrons
  electron_maxstep      = $electron_maxstep
  mixing_mode		= 'plain',
  mixing_fixed_ns	= 0,
  diagonalization	= 'david',
  conv_thr		= 1.0d-6,
  mixing_beta		= 0.3,
 /
CELL_PARAMETERS
EOF
cat >> $prefix.scf.in < $LAT_BZ_dir/CP_RS.uvw
cat >> $prefix.scf.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
K_POINTS {automatic}
    $mp_size    $mp_size    $mp_size    0    0    0
EOF
#exit
###############################################################################
# Do SCF Run                                                                  #
###############################################################################
cat >> $prefix.scf.data << EOF
  Beginning $prefix SCF...
  Starting at `date`
*******************************************************************************
$prefix.scf.in follows...
----------------------------------
EOF

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

*******************************************************************************
  Running $prefix SCF...
EOF

# Exit Point for Testing
#exit

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

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

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

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

  d-orbital occupations
  ---------------------
`grep -n 'eigenvalues:' $prefix.scf.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,
  U_projection_type = '$U_p_t',
  starting_ns_eigenvalue(1,2,1) = $s_ns_121,
  starting_ns_eigenvalue(2,2,1) = $s_ns_221,
  starting_ns_eigenvalue(3,2,1) = $s_ns_321,
  starting_ns_eigenvalue(4,2,1) = $s_ns_421,
  starting_ns_eigenvalue(5,2,1) = $s_ns_521,
  starting_ns_eigenvalue(1,1,2) = $s_ns_112,
  starting_ns_eigenvalue(2,1,2) = $s_ns_212,
  starting_ns_eigenvalue(3,1,2) = $s_ns_312,
  starting_ns_eigenvalue(4,1,2) = $s_ns_412,
  starting_ns_eigenvalue(5,1,2) = $s_ns_512,