ATK 11.8 contains an extended set of pseudopotentials, with semicore potentials for several elements. In this mini-tutorial you will learn how to use one of these to predict correctly the half-metallic ferromagnetic properties of an organometallic molecular wire consisting of benzene and vanadium atoms. The aim is essentially to reproduce the results by Maslyuk et al., PRL 97, 097201 (2006).

For this system it is crucial to include the entire n=3 and n=4 shells in the valence. Vanadium has an electronic configuration of [1s²2s²2p⁶]3s²3p⁶4s²3d³; the semicore pseudopotential that will be used includes the orbitals in brackets (corresponding to Ne), while the bold-face orbitals will be valence orbitals and thus constitute the SingleZeta basis set. The basis set will be discussed a bit further below.

It will be assumed that you already have some basic experience of setting up and running calculations with ATK using the graphical interface VNL, for instance by first going through the basic VNL tutorial.

Setting up the geometry

The geometry will be built using a generic benzene molecule from the VNL database, and the vanadium atom will be added by hand. It is immediately obvious that the C-C bond lengths do not correspond to those quoted in the article, so we will have to optimize the geometry to make sure we have found the minimum total energy. In fact, as it turns out if you just use the original geometry, the result is not half-metallic.

1. Click the Database icon on the VNL toolbar, and go into the molecule database via the Databases menu. Locate "Benzene" and send it to the Builder by using the "Send to" button in the lower right-hand corner.
2. Inscribe the molecule in a unit cell by clicking the Bulk icon on the left-hand side toolbar.
3. Taking the cell size in the direction of the wire from the article, we set the Z-component of the unit cell vector C to 3.4 Å.
4. Use the Center function in the Transform menu to adjust the geometry after changing the cell size.
5. Now, right-click the Basis widget and choose "Insert atom". This will place a new hydrogen atom last in the list of atoms (index 12). Change this atom to Vanadium, and set its z-coordinate to 0 Å.
6. Again use the Center function to adjust the geometry. The geometry setup should now look like this:
   
7. Finally, send the structure to the Script Generator, again using the Send To button in the lower right-hand corner.

Setting up the calculation script

This calculation involves one small aspect that you will have to enter into the script by hand. The Script Generator can however be used to do most of the work.

1. Define the output file name, for instance V-C6H6.nc.
2. Double-click New Calculator in the Blocks panel to insert such a block in the script, then double-click the inserted script block in the Script panel to open it.
3. As you can test yourself, if the calculation is run with LDA, the results do not show half-metallicity, so you have to use SGGA (S for spin-polarization).
4. Choose a conservative guesstimate for the number of kC-points; the article used 24, so 31 seems like a good number.
   
5. Click "LCAO basis set" in the left-hand side panel in the New Calculator dialog, and change the basis set of Vanadium to DoubleZeta. This is not the basis set you want, but by changing the value from the default DoubleZetaPolarized VNL is forced to write out the basis set definition explicitly in the script, so that you can easily edit it later.
   
6. Close the New Calculator dialog. It is not enough to set the exchange-correlation functionality to spin-polarized GGA, you must also initialize the configuration in a spin-dependent state to actually perform a spin-polarized calculation. To do so, insert an Initial State block and open it. Select "User spin" from "Initial state type", and set the value of the initial spin to 1 for the vanadium atom (C and H atoms can be left at 0 initial spin), and close the dialog.
   
7. As mentioned above, the geometry need to be optimized. Therefore, double-click Optimization and insert an OptimizeGeometry block. Open it, and lower the force tolerance to 0.01 eV/Å. Also remove (uncheck) the cell constraint in Z, and lower the stress tolerance to 0.01 eV/Å^3.
   
8. Finally, insert some Analysis quantities, i.e. things to be calculated for the optimized geometry.
   • Bandstructure. Change the number of points per segment to 100 for more detailed curves)
   • DensityOfStates. To get accurate results, use 301 k-points in C (don't forget to uncheck "Sync"; you only want 1 point in A and B!).
     
   • Electron density. Twice - once for spin Up, once for Down.
     
   • MullikenPopulation
9. Make sure you have inserted all script blocks before proceeding to the next step. (Do you have 2 electron density blocks?)
   

Selecting the semicore basis set/pseudopotential

The final step in the setup of the script is to replace the "standard" ATK basis set with one based on a semicore pseudopotential. To do this, use the Send To icon in the Script Generator and transfer the script to the Editor .

Locate the definition of the basis set, and replace the vanadium basis set with HGHGGAPBEBasis.Vanadium_HGH_13_DoubleZetaPolarized, as shown in the figure below.

Save the script!

Running the calculation

To execute the script, send it to the Job Manager via the Send To button in the Editor, and click Run Queue. It will take 30-60 minutes depending on your computer. If you prefer you can also run this script in parallel, that will speed up the calculation quite nicely, as there are many k-points to parallelize over. When run on 16 MPI nodes on two double-socket quad-core machines, it took about 20 minutes. (We here used the cores as MPI nodes; this is possible and beneficial for speed only when the calculation doesn't use a lot of memory compared to the available RAM/core.)

Analyzing the results

Once the script finishes, select the output file V-C6H6.nc in VNL (it will be located in the directory where VNL was started, unless you specified an explicit path for it in the Script Generator). In the following sections the results will be compared to those in the article.

Band structure

To plot the band structure, select the BandStructure object in the Result Browser (in the upper right-hand corner of the main VNL window), and click "Plot>Show" in the planel below (Result Viewer). Zooming in around the Fermi level, you will see a nice correspondence to the published results, and indeed the system exhibits half-metallicity: there are no majority states (black lines) at the Fermi level. You can right-click the plot and tick "Anti aliasing" to improve the visual quality of the plot.

Magnetic moment

Next you will compute the magnetic moments of the atoms. In the Result Browser, select the Mulliken population object, and then click "Print>Show" in the Result Viewer.

From the Mulliken population report that is displayed, you can subtract the spin up populations from the spin down populations to find the following net magnetic moments, in excellent agreement with Maslyuk et al. (article values in parentheses):

• –0.0405μ_B per C atom (–0.0467μ_B)
• +1.21μ_B for V (+1.28μ_B)
• Total spin polarization of the cell: +0.98μ_B (+1μ_B)
  

You can also infer from this report that the population of the carbon and d-orbitals and hydrogen p-prbitals is very small, thus probably a DoubleZeta basis set would have been sufficient for these elements.

It is also instructive to inspect the symmetries of the vanadium orbitals, to understand the basis set notation a bit better. As mentioned in the introduction,the SingleZeta basis set for this semicore calculation is 3s²3p⁶4s²3d³; this corresponds to the first 4 sets of orbitals in the report above (the selected lines). After that follows the split orbitals, which creates a DoubleZeta basis set. The 3s orbitals are however not split, so there are only 3 more sets of orbitals (3p,4s,3d). Finally, the final set completes the DoubleZetaPolarized basis set by adding 4p orbitals. This is summarized in the picture below.

It can been seen that the Vanadium 3p shell is rather inert, the total population of the 3p basis orbitals is very close to 6. The 4s orbitals, on the other hand, donate their electrons to the 3d and 4p shells, which have total occupations 3.76 and 0.60, respectively.

Optimized geometry

It would be interesting to see how the the optimized geometry compares to that in the reference article. The NetCDF file contains two configurations; drag the optimized geometry (id=gID001) onto the Builder icon on the VNL toolbar. Close the Basis and Lattice panels, and instead open the Z-Matrix panel. Select 2 neighboring carbon atoms to find that the optimized C-C distance is 1.456 Å, in perfect agreement with the article. The original distance in the configuration we built was 1.40 Å, and this difference is crucial to obtaining the spin-gap at the Fermi level; if you skipped the optimization step the result would be purely metallic.

Projected density of states

The band structure already shows that there are no majority states at the Fermi level. This can be confirmed by plotting the density of states.

Select the Density of states object in the Result Browser, and click "PDOS>Show". The 3s semicore states are very deep, and extremely flat, so there will be a very large DOS peak at –67 eV, but if you zoom in around the Fermi level you can indeed see that there is a band gap for the majority spin (black curve).

You can project the DOS on the various atoms by clicking on them in the 3D view, and/or on angular momenta by ticking the corresponding boxes. (The data is actually computed on the fly when you select atoms or orbitals, so it can take a few seconds each time a new selection is made. Be patient, and avoid clicking several times before the plot updates!)

Spin-polarization density

Finally, how about a comparison with Fig 3. in the article, which shows the real-space spin-polarization density, i.e. the difference in electron density for spin up and down. The NetCDF file already contains the spin up and down densities, but you need to make a small script that subtracts them from each other.

Open the Editor, and paste the following 2 lines of code into it.

e_dens_list = nlread("V-C6H6.nc", ElectronDensity)
nlsave("V-C6H6.nc", e_dens_list[0]-e_dens_list[1], object_id="Spin Polarization Density")

Send this to the Job Manager and run the queue (insert the path to the NC file unless it's located where VNL was started).

When the script has finished (it takes a few seconds only), go back to the main VNL window and select another file, then the V-C6H6.nc file again (this is to re-read the file). You will see a new object added to it, which you can plot as a contour plot. By manipulating the view settings, the figure below can be produced.

Follow these instructions to obtain the precise view of the plot shown above.

• You will need to add 2 contour plots and the configuration to the same Viewer.

1.  
   1. To insert the first one, select the "Spin Polarization Density" object in the Result Browser and then click Contour>Show in the Result Viewer below.
   2. Now, to insert the second contour plot, select "Plots>Show/hide>Electron density>Contour" from the menu in the Viewer window that just opened. The second contour plot will not be seen at first, because by default it is placed at the same position as the first one, but you can see there are two color bars.
   3. Then, go back to the Result Browser and drag and drop the relaxed geometry (id=gID001) onto the open Viewer window.

• Next, the plot styles need to be configured.

1. Press choose "Plots>Properties" from the menu, and select the first contour plot. Change its Θ_y rotation angle to 90 degrees to align it with the benzene plane, and then set its origin to 0.85 Å to make it go through the vanadium atom. 
2. Set the color minimum/maximum to -0.06 and 0.06.
3. Select the second contour plot and set the same color minimum/maximum.
4. Open its "Color bar" subentry and hide it; since both contour plots use the same color scale we only need one colorbar.
5. You can change the colors of the atoms by selecting them, and then right-clicking and choosing Colors... from the context menu.
6. Rotate the plot into a suitable orientation.

Note: The red-white-blue colormap used in the image above will be available in ATK 12.2. For now use one of the other ones; the default "Rainbow" is quite suitable. Both contour plots should use the same colormap.