1 #!/usr/bin/env python
   2 
   3 
   4 # Localized Wannier orbitals for ferromagnetic iron
   5 # Generating the bandstructure 
   6 
   7 from Dacapo import Dacapo
   8 from ASE import Atom,ListOfAtoms
   9 from ASE.Trajectories import NetCDFTrajectory
  10 import Numeric as num
  11 import os
  12 
  13 # check if we have already calculated the nc file
  14 if not os.path.isfile('fe-bcc.nc'):
  15     
  16     # Pt wire
  17     lat = 2.87
  18     atoms = ListOfAtoms([
  19                          Atom('Fe', ([0,0,0]),magmom=2.32)])
  20                          
  21 
  22     cell = num.array([[-1/2.,1/2.,1/2.],[1/2.,-1/2.,1/2.],[1/2.,1/2.,-1/2.]])*2.87
  23     atoms.SetUnitCell(cell)
  24 
  25     # Dacapo calculator:
  26     calc = Dacapo(planewavecutoff=400,nbands=16,kpts=(5,5,5),xc='PBE',out='fe-bcc.nc', 
  27                   spinpol=True)
  28     atoms.SetCalculator(calc)
  29 
  30     # displace kpoints sligthly, so that the symmetry program in dacapo does
  31     # not use inversion symmetry to remove kpoints.
  32     kpoints = calc.GetBZKPoints()
  33     kpoints[:,0] += 2e-5
  34     calc.SetBZKPoints(kpoints)
  35     
  36     tot = atoms.GetPotentialEnergy()
  37 
  38 atoms = Dacapo.ReadAtoms('fe-bcc.nc')
  39 
  40 # Begin the Wannier part
  41 from Dacapo import Dacapo
  42 from ASE.Utilities.Wannier.Wannier import Wannier
  43 import Numeric as num
  44 
  45 atoms = Dacapo.ReadAtoms('fe-bcc.nc')
  46 calc = atoms.GetCalculator()
  47 
  48 # Use 5 d-orbitals (l=2,m=-2,..,2) centered on atom 0 with a radius of 0.4 A as start guess.
  49 # A random start guess will be used for the last (6th) WF.
  50 initialwannier = [[[0],2,0.4]]
  51 
  52 # Include all eigenstate below the Fermi level
  53 occenergy=0.0
  54 # Construct Wannier functions for spin down
  55 spin=1
  56 
  57 wannier = Wannier(numberofwannier=5+1,calculator=calc,
  58                   occupationenergy=occenergy, 
  59                   initialwannier=initialwannier, 
  60                   spin=spin) 
  61 
  62 # Store the localization matrix. 
  63 # wannier.ReadZIBlochMatrix('fe_bloch.pickle')
  64 wannier.SaveZIBlochMatrix('fe_bloch.pickle')
  65 # It can be read again in a later run by: wannier.ReadZIBlochMatrix('fe_bloch.pickle')
  66 
  67 # Perform the localization
  68 wannier.Localize(tolerance=1.0e-7)
  69 
  70 # Store the WFs. They can be read again in a later run by:
  71 wannier.SaveRotation('fe_wannier')
  72 # They can be read again in a later run by: wannier.ReadRotation('fe_wannier')
  73 
  74 # Translate all WFs to the unit cell (2,2,2)
  75 wannier.TranslateAllWannierFunctionsToCell((2,2,2))
  76 
  77 # Print the centers and radii of the WFs
  78 centers = wannier.GetCenters()
  79 for n in range(len(centers)):
  80     print n,centers[n]
  81 
  82 # Get the centers as an atom plot
  83 centers = wannier.GetCentersAsAtoms()
  84 traj1 = NetCDFTrajectory('centers-febcc-'+str(spin)+'.traj',centers)
  85 traj1.Update()
  86 traj1.Close()
  87 
  88 # Store a '.cube' file for each Wannier function
  89 for n in range(6): 
  90     wannier.WriteCube(n,'fe_'+str(n)+'_spin_'+str(spin)+'.cube')
  91 
  92 
  93 # band structure 
  94 fermilevel = calc.GetFermilevel()
  95 G = (0,0,0)
  96 P = (0.25,0.25,0.25)
  97 N = (0.0,0.0,0.5)
  98 H = (-0.5,0.5,0.5)
  99 
 100 from ASE.Utilities.Wannier import HoppingParameters
 101       
 102 npoints = 80
 103 cutoff = 12.0
 104       
 105 hop=HoppingParameters.HoppingParameters(wannier,cutoff)
 106 hop.WriteBandDiagramToNetCDFFile('GH'+str(cutoff)+'.nc',npoints,G,H,offset=fermilevel)
 107 hop.WriteBandDiagramToNetCDFFile('HP'+str(cutoff)+'.nc',npoints,H,P,offset=fermilevel)
 108 hop.WriteBandDiagramToNetCDFFile('PN'+str(cutoff)+'.nc',npoints,P,N,offset=fermilevel)
 109 hop.WriteBandDiagramToNetCDFFile('NG'+str(cutoff)+'.nc',npoints,N,G,offset=fermilevel)
 110 hop.WriteBandDiagramToNetCDFFile('GP'+str(cutoff)+'.nc',npoints,G,P,offset=fermilevel)
 111