In this example, we will calculate crystalline orbitals of Cu₂O. The crystalline orbitals can be obtained either as Bloch functions (corresponding to the canonical molecular obitals in molecules) or as Wannier functions (localized crystalline orbitals). Wannier functions can only be obtained for systems with a band gap.

Prerequisites:

1. You should already have calculated and studied the electronic band structure as discussed in the section Electronic band structure and density of states for Cu2O (DFT-PBE0/TZVP level of theory).

1. Crystalline orbitals

We will run an ORBITALS calculation, which will give us botn occupied and virtual orbitals at all k-points that we request.

Let's start from the directory that contains the wavefunction file Cu2O_Pn-3m_PBE0_TZVP_opt.w from the geometry optimization. Create a new subdirectory orbitals and copy the wavefunction file there with a new name:

mkdir orbitals
cp Cu2O_Pn-3m_PBE0_TZVP_opt.w orbitals/Cu2O_Pn-3m_PBE0_TZVP_orbitals.w

Go to the orbitals directory and create the following property calculation input file Cu2O_Pn-3m_PBE0_TZVP_orbitals.d3 in the directory:

NEWK
0 8
8 8 8
1 0
ORBITALS
Cu2O
1
0
ENDORB
ENDPROP

Line 1, NEWK: Request CRYSTAL to calculate the Fock/Kohn-Sham eigenvectors at a number of k-points in the reciprocal space. The lines 0 8 and 8 8 8 are equal to the original SHRINK values (see the chapter on geometry optimization). You could also reduce the k-mesh here since you probably are not interested in seeing the orbitals in all k-points (it depends on what you are looking for - often just the Γ-point is enough). Line 1 0 requests CRYSTAL to calculate the Fermi energy using this k-mesh (1) and tells that no additional printing is required (0).

Line 5, ORBITALS: Request CRYSTAL to write out crystalline orbitals, naming them as Cu2O_number_kpoint_real/complex.molden. The number 1 asks CRYSTAL to use fractional coordinates (default for periodic systems) and the number 0 asks CRYSTAL to write the crystalline orbitals as Bloch functions. 

After this, you can submit the job:

jsub -np 8 -smp crystal Cu2O_Pn-3m_PBE0_TZVP_orbitals.d3

Running the calculation will produce 35 .molden files, corresponding to the number of k-points in the irreducible Brillouin zone (line NUMBER OF K POINTS IN THE IBZ 35 in output).

The orbital files are written in a custom Molden format that is supported only by relatively recent (post-2017) versions of Jmol. Start by opening the Γ-point orbital file Cu2O_00001_K000_real.molden in Jmol. After opening the file, execute the following console commands:

load "" {2 2 2}
connect (_Cu) (_Cu) delete
isosurface color green red cutoff 0.04 resolution 10 mo 66

This will:

1. Load a 2x2x2 supercell of Cu₂O
2. Delete the Cu-Cu bonds that Jmol draws by default
3. Plot the highest occupied crystalline orbital (mo 66) with a isodensity value of 0.04 a.u. (cutoff 0.04).

The image will look bit like below, especially if you also enter command color background white in console (interactive JSmol plot perhaps coming later):

Have a look at the other orbitals, too (mo 65, mo 64, etc.). And have a look at orbitals at different k-points. For example, R point (1/2, 1/2, 1/2) is in file Cu2O_00035_K444_real.molden (K444 = k-point (1/2, 1/2, 1/2) in terms of SHRINK 8 8 8. They represent the same crystalline orbitals at a different phase.

2. Wannier functions

Let's start from the directory that contains the wavefunction file Cu2O_Pn-3m_PBE0_TZVP_opt.w from the geometry optimization. Create a new subdirectory orbitals and copy the wavefunction file there with a new name:

mkdir -p orbitals
cp Cu2O_Pn-3m_PBE0_TZVP_opt.w orbitals/Cu2O_Pn-3m_PBE0_TZVP_wf.w

Go to the orbitals directory and create the following property calculation input file Cu2O_Pn-3m_PBE0_TZVP_wf.d3 in the directory:

NEWK
0 8
8 8 8
1 0
LOCALWF
VALENCE
ENDLOCAL
ORBITALS
Cu2O_WF
1
1
ENDORB
ENDPROP

Line 1, NEWK: Same as for crystalline orbitals above

Line 5: LOCALWF: Ask CRYSTAL to calculate localized crystalline orbitals (Wannier functions). Here we localize all valence orbitals, but you can also ask localization of a range of bands (INIFIBND) or just some bands (BANDLIST). See manual under keyword LOCALWF. If localizing a subset of occupied bands crashes, try localizing all valence bands. Note that only occupied bands can be localized.

Line 8, ORBITALS: Request CRYSTAL to write out crystalline orbitals, naming them as Cu2O_WF_number_kpoint_real/complex.molden. The first number 1 asks CRYSTAL to use fractional coordinates (default for periodic systems) and the second number 1 asks CRYSTAL to write out the lozalized orbitals. 

Note that LOCALWF cannot be run in parallel!. So, submit the job with one CPU:

jsub crystal Cu2O_Pn-3m_PBE0_TZVP_WF.d3

Running the calculation will produce file Cu2O_WF_00001_K000_real.molden that contains the localized orbitals.

Open the file in Jmol and execute the following console commands:

load "" {2 2 2}
connect (_Cu) (_Cu) delete
isosurface color green red cutoff 0.004 resolution 10 mo 66

Note how the isodensity value is 10 times smaller than for the crystalline orbitals above. This is because the localized orbitals are spatially much more constrained:

3. Example for α-silicon

Cu₂O has rather ionic bonds, let's see an example for a material with covalent bonding.

TODO.