.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples_generated/tutorials/wannier.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_generated_tutorials_wannier.py: .. _wannier tutorial: ============================================ Partly occupied Wannier Functions ============================================ This tutorial walks through building **partly occupied Wannier functions** with the :mod:`ase.dft.wannier` module and the `GPAW `_ electronic structure code. For more information on the details of the method and the implementation, see | K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen | Partly occupied Wannier functions: Construction and applications | `Phys. Rev. B 72, 125119 (2005) `__ .. contents:: **Outline** :depth: 2 :local: Benzene molecule ================ Step 1 – Ground-state calculation --------------------------------- Run the script below to obtain the ground-state density and the Kohn–Sham (KS) orbitals. The result is stored in :file:`benzene.gpw`. .. GENERATED FROM PYTHON SOURCE LINES 32-57 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from gpaw import GPAW, restart from ase import Atoms from ase.build import molecule from ase.dft.kpoints import monkhorst_pack from ase.dft.wannier import Wannier atoms = molecule('C6H6') atoms.center(vacuum=3.5) calc = GPAW(mode='pw', xc='PBE', txt='benzene.txt', nbands=18) atoms.calc = calc atoms.get_potential_energy() calc = calc.fixed_density( txt='benzene-harris.txt', nbands=40, convergence={'bands': 35}, ) atoms.get_potential_energy() calc.write('benzene.gpw', mode='all') .. GENERATED FROM PYTHON SOURCE LINES 58-67 Step 2 – Maximally localized WFs for the occupied subspace (15 WFs) ------------------------------------------------------------------- There are 15 occupied bands in the benzene molecule. We construct one Wannier function per occupied band by setting ``nwannier = 15``. By calling ``wan.localize()``, the code attempts to minimize the spread functional using a gradient-descent algorithm. The resulting WFs are written to .cube files, which allows them to be inspected using e.g. VESTA. .. GENERATED FROM PYTHON SOURCE LINES 69-77 .. code-block:: Python atoms, calc = restart('benzene.gpw', txt=None) # Make wannier functions of occupied space only wan = Wannier(nwannier=15, calc=calc) wan.localize() for i in range(wan.nwannier): wan.write_cube(i, f'benzene15_{i}.cube') .. GENERATED FROM PYTHON SOURCE LINES 78-86 Step 3 – Adding three extra degrees of freedom (18 WFs) ------------------------------------------------------- To improve localization we augment the basis with three extra Wannier functions - so-called *extra degrees of freedom* (``nwannier = 18``, ``fixedstates = 15``). This will allow the Wannierization procedure to use the unoccupied states to minimize spread functional. .. GENERATED FROM PYTHON SOURCE LINES 88-97 .. code-block:: Python atoms, calc = restart('benzene.gpw', txt=None) # Make wannier functions using (three) extra degrees of freedom. wan = Wannier(nwannier=18, calc=calc, fixedstates=15) wan.localize() wan.save('wan18.json') for i in range(wan.nwannier): wan.write_cube(i, f'benzene18_{i}.cube') .. GENERATED FROM PYTHON SOURCE LINES 98-104 Step 4 – Spectral-weight analysis --------------------------------- The script below projects the WFs on the KS eigenstates. You should see the 15 lowest bands perfectly reconstructed (weight ≃ 1.0) while higher bands are only partially represented. .. GENERATED FROM PYTHON SOURCE LINES 106-131 .. code-block:: Python atoms, calc = restart('benzene.gpw', txt=None) wan = Wannier(nwannier=18, calc=calc, fixedstates=15, file='wan18.json') weight_n = np.sum(abs(wan.V_knw[0]) ** 2, 1) N = len(weight_n) F = wan.fixedstates_k[0] plt.figure(1, figsize=(12, 4)) plt.bar( range(1, N + 1), weight_n, width=0.65, bottom=0, color='k', edgecolor='k', linewidth=None, align='center', orientation='vertical', ) plt.plot([F + 0.5, F + 0.5], [0, 1], 'k--') plt.axis(xmin=0.32, xmax=N + 1.33, ymin=0, ymax=1) plt.xlabel('Eigenstate') plt.ylabel('Projection of wannier functions') plt.savefig('spectral_weight.png') plt.show() .. image-sg:: /examples_generated/tutorials/images/sphx_glr_wannier_001.png :alt: wannier :srcset: /examples_generated/tutorials/images/sphx_glr_wannier_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 132-143 Polyacetylene chain (1-D periodic) ================================== We now want to construct partially occupied Wannier functions to describe a polyacetylene chain. Step 1 – Structure & ground-state calculation --------------------------------------------- Polyacetylene is modelled as an infinite chain; we therefore enable periodic boundary conditions along *x*. .. GENERATED FROM PYTHON SOURCE LINES 145-172 .. code-block:: Python kpts = monkhorst_pack((13, 1, 1)) calc = GPAW( mode='pw', xc='PBE', kpts=kpts, nbands=12, txt='poly.txt', convergence={'bands': 9}, symmetry='off', ) CC = 1.38 CH = 1.094 a = 2.45 x = a / 2.0 y = np.sqrt(CC**2 - x**2) atoms = Atoms( 'C2H2', pbc=(True, False, False), cell=(a, 8.0, 6.0), calculator=calc, positions=[[0, 0, 0], [x, y, 0], [x, y + CH, 0], [0, -CH, 0]], ) atoms.center() atoms.get_potential_energy() calc.write('poly.gpw', mode='all') .. GENERATED FROM PYTHON SOURCE LINES 173-178 Step 2 – Wannierization ----------------------- We repeat the localization procedure, keeping the five lowest bands fixed and adding one extra degree of freedom to aid localization. .. GENERATED FROM PYTHON SOURCE LINES 180-215 .. code-block:: Python import numpy as np from gpaw import restart from ase.dft.wannier import Wannier atoms, calc = restart('poly.gpw', txt=None) # Make wannier functions using (one) extra degree of freedom wan = Wannier( nwannier=6, calc=calc, fixedenergy=1.5, initialwannier='orbitals', functional='var', ) wan.localize() wan.save('poly.json') wan.translate_all_to_cell((2, 0, 0)) for i in range(wan.nwannier): wan.write_cube(i, f'polyacetylene_{i}.cube') # Print Kohn-Sham bandstructure ef = calc.get_fermi_level() with open('KSbands.txt', 'w') as fd: for k, kpt_c in enumerate(calc.get_ibz_k_points()): for eps in calc.get_eigenvalues(kpt=k): print(kpt_c[0], eps - ef, file=fd) # Print Wannier bandstructure with open('WANbands.txt', 'w') as fd: for k in np.linspace(-0.5, 0.5, 100): ham = wan.get_hamiltonian_kpoint([k, 0, 0]) for eps in np.linalg.eigvalsh(ham).real: print(k, eps - ef, file=fd) .. GENERATED FROM PYTHON SOURCE LINES 216-221 Step 3 – High-resolution band structure --------------------------------------- Using the Wannier Hamiltonian we can interpolate the band structure on a fine 100-point *k* mesh and compare it to the original DFT result. .. GENERATED FROM PYTHON SOURCE LINES 223-247 .. code-block:: Python fig = plt.figure(dpi=80, figsize=(4.2, 6)) fig.subplots_adjust(left=0.16, right=0.97, top=0.97, bottom=0.05) # Plot KS bands k, eps = np.loadtxt('KSbands.txt', unpack=True) plt.plot(k, eps, 'ro', label='DFT', ms=9) # Plot Wannier bands k, eps = np.loadtxt('WANbands.txt', unpack=True) plt.plot(k, eps, 'k.', label='Wannier') plt.plot([-0.5, 0.5], [1, 1], 'k:', label='_nolegend_') plt.text(-0.5, 1, 'fixedenergy', ha='left', va='bottom') plt.axis('tight') plt.xticks( [-0.5, -0.25, 0, 0.25, 0.5], [r'$X$', r'$\Delta$', r'$\Gamma$', r'$\Delta$', r'$X$'], size=16, ) plt.ylabel(r'$E - E_F\ \rm{(eV)}$', size=16) plt.legend() plt.savefig('bands.png', dpi=80) plt.show() .. image-sg:: /examples_generated/tutorials/images/sphx_glr_wannier_002.png :alt: wannier :srcset: /examples_generated/tutorials/images/sphx_glr_wannier_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 248-253 Within the fixed-energy window—that is, for energies below the fixed-energy line—the Wannier-interpolated bands coincide perfectly with the DFT reference (red circles). Above this window the match is lost, because the degrees of freedom deliberately mix several Kohn–Sham states to achieve maximal real-space localisation. .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 0.266 seconds) .. _sphx_glr_download_examples_generated_tutorials_wannier.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: wannier.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: wannier.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: wannier.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_