1 #!/usr/bin/env python
   2 
   3 # Localized Wannier orbitals for a Pt wire
   4 
   5 from Dacapo import Dacapo
   6 from ASE import Atom,ListOfAtoms
   7 from ASE.Visualization.RasMol import RasMol
   8 import os
   9 
  10 # check if we have already calculated the nc file
  11 if not os.path.isfile('ptwire.nc'):
  12     lat = 2.41
  13     atoms = ListOfAtoms([
  14           Atom('Pt', ([0*lat,2*lat,2*lat]),tag=0)],
  15           cell=([lat,4*lat,4*lat]), periodic=True)
  16 
  17     # Dacapo calculator:
  18     calc = Dacapo(planewavecutoff=350, nbands=8,kpts=(15,1,1),xc='PBE',out='ptwire.nc')
  19     atoms.SetCalculator(calc)
  20 
  21     # Displace kpoints sligthly, so that the symmetry program in dacapo does
  22     # not use inversion symmetry to reduce kpoints.
  23     kpoints = calc.GetBZKPoints()
  24     kpoints[:,0] += 2e-5
  25     calc.SetBZKPoints(kpoints)
  26 
  27 energy = atoms.GetPotentialEnergy()
  28 
  29 
  30 # Wannier part
  31 from ASE.Utilities.Wannier.Wannier import Wannier
  32 
  33 atoms = Dacapo.ReadAtoms('ptwire.nc')
  34 calc = atoms.GetCalculator()
  35 
  36 # Use 5 d-orbitals centered at the atom and 1 s-orbital centered at the Pt-Pt bond
  37 # The radius of both types of orbitals is set to 0.4 Angstroms  
  38 initialwannier=[
  39 [[0],2,0.4],
  40 [[0.5,0.5,0.5],0,0.4]
  41 ]
  42 
  43 # Initialize the Wannier class
  44 wannier = Wannier(numberofwannier=6,calculator=calc,occupationenergy=2.0,initialwannier=initialwannier)
  45 
  46 # Perform the localization
  47 wannier.Localize()
  48 
  49 # Translate all the WFs to cell (1,0,0)
  50 wannier.TranslateAllWannierFunctionsToCell([1,0,0])
  51 
  52 # Band diagram
  53 from ASE.Utilities.Wannier.HoppingParameters import HoppingParameters
  54 
  55 hop=HoppingParameters(wannier,couplingradius=9.0)
  56 
  57 # Get a band with 50 k-points
  58 k,band=hop.GetBandDiagram(50,[0,0,0],[0.5,0,0])
  59 
  60 # Plot the band using Gnuplot
  61 import Gnuplot
  62 plot=Gnuplot.Gnuplot()
  63 plot('set data style lines')
  64 plot.plot(Gnuplot.Data(k,band[0]),Gnuplot.Data(k,band[1]),Gnuplot.Data(k,band[2]),Gnuplot.Data(k,band[3]),Gnuplot.Data(k,band[4]),Gnuplot.Data(k,band[5]))
  65 
  66 # Write the band diagram as a NetCDF file.
  67 hop.WriteBandDiagramToNetCDFFile('ptband.nc',50,[0,0,0],[0.5,0,0])
  68 
  69 # make a 3D isosurface plot using VTK
  70 from ASE.Visualization.VTK import VTKPlotElectronicState,VTKPlotAtoms
  71 state = wannier.GetElectronicState(0)
  72 
  73 plot = VTKPlotElectronicState(state)
  74 plot.SetRepresentationToIsoSurface2([0.1])
  75 atomsplot = VTKPlotAtoms(atoms,parent=plot) 
  76 plot.Update()
  77 atomsplot.RemoveAvatar(atomsplot.GetAvatars()[0])
  78 plot.Render()