#!/bin/bash

report=report.txt

print_header() {
   echo "* This is a fairly advanced application of the Hubbard model"
   echo "* The toy model is a Cr dimer surrounded by two H atoms"
   echo "* It is a spin polarized calculation with the spin is flipped on one of the Cr atoms"
   echo "* The d orbitals of the two Cr atoms get a Hubbard U correction"
   echo "* Three equivalent input setups are shown to achieve the same"
   echo "* In the output feedback is given after the header: H U B B A R D    S E T T I N G S"

   printf "\n%10s %10s %10s\n" "input" "Energy" "E(hubbard)"
}

myreport () {
   energy=`$AMSBIN/amsreport $AMS_JOBNAME.results -k "AMSResults%Energy#10.6f"`
   hubbard_energy=`$AMSBIN/amsreport $AMS_JOBNAME.results/band.rkf -k "Bond energy terms%Hubbard Energy#10.6f"`

   printf "%10s %10s %10s\n" $1 $energy $hubbard_energy
}


print_header > $report


# Method 1, set hubbard for all Cr atoms

export AMS_JOBNAME=inp=atom

$AMSBIN/ams --delete-old-results <<EOF

Task SinglePoint

System
   Atoms
      Cr     0.0000     0.0000     0.0000 
      H      10.0       0.0        0.0
      Cr     0.0000     0.0000     1.6790  region=SpinFlip
      H      -10.0       0.0       0.0
   End
End


Engine band
   Unrestricted yes

   Convergence  SpinFlipRegion=SpinFlip

   HubbardU
      Enabled yes
      PrintOccupations No
      Atom Element=Cr UValue=0.02 lValue=d
   End
EndEngine


EOF

myreport "atoms" >> $report

# Method 2, set hubbard for atoms in a region

export AMS_JOBNAME=inp=region

$AMSBIN/ams --delete-old-results <<EOF

Task SinglePoint

System
   Atoms
      Cr     0.0000     0.0000     0.0000  region=AllCr
      H      10.0       0.0        0.0
      Cr     0.0000     0.0000     1.6790  region=SpinFlip,AllCr
      H      -10.0       0.0       0.0
   End
End


Engine band
   Unrestricted yes

   Convergence  SpinFlipRegion=SpinFlip

   HubbardU
      Enabled yes
      PrintOccupations No
      Region Name=AllCr UValue=0.02 lValue=d
   End
EndEngine


EOF

myreport "regions" >> $report

# Method 3: the old and inconvenient way, relying on knowledge of Band's internal atom type system.

export AMS_JOBNAME=inp=obsolete

$AMSBIN/ams --delete-old-results <<EOF

Task SinglePoint

System
   Atoms
      Cr     0.0000     0.0000     0.0000 
      H      10.0       0.0        0.0
      Cr     0.0000     0.0000     1.6790  region=SpinFlip
      H      -10.0       0.0       0.0
   End
End


Engine band
   Unrestricted yes

   Convergence  SpinFlipRegion=SpinFlip

   HubbardU
      Enabled yes  # band sees three types, 1 and 3 are the Cr atoms
      PrintOccupations No
      UValue 0.02    0.0   0.02
      lValue 2      -1     2
   End
EndEngine


EOF

myreport "obsolete" >> $report

echo "begin report"
cat $report
echo "end report"