1.1. Simple Molecular Device

Experimentally, monoatomic gold chains have been successfully synthesized using mechanically controllable break junctions or with a scanning tunneling microscope under elongation [YBBA98]. In this tutorial we’ll show how to use the NanoDCAL code to calculate the electronic transport through a carbon-carbon triple bond sandwiched between gold chain electrodes, as shown in Fig. 1.1.1. This simple model will be useful to explore all methodology steps.

Fig. 1.1.1 Structure of two-probe device.

The calculation steps are given as follows:

1. Structural Relaxation:

   • Build and relax the electrode unit cell;

   • Build and relax the scattering region.

2. Transport calculation:

   • Electrode SCF calculation;

   • Two-probe SCF-NEGF calculation.

   • Effective Potential analysis

   • Transmission spectrum

1.1.1. Structural relaxation

Build and relax the electrode unit cell

First, we have to build the electrode unit cell and generate the required input file. This could be done using Device Studio software, as described bellow:

1. Open Device Studio software ➟ Create a new Project ➟ on File name type Au-wire ➟ press Save;

2. In the menu bar, go to File ➟ New ➟ Crystal to pop up a new window;

   • Set a orthorhombic unit cell ➟ \(a\) = 8 , \(b\) = 8, \(c\) = 2.5 and \(\alpha = \beta = \gamma = 90^{\circ}\);

   • In the element tab, click on the plus button to add a Au atom ➟ set the fractional coordinates as [0.5 0.5 0.5], click Preview and Build.

3. In the menu bar of Device Studio, click Simulator ➟ NanoDCAL ➟ Relax to pop up the interface Relax Calculation for Crystal to set the parameters;

   • In the K-point Sampling set n3 = 100;

   • Go to the Relax ➟ Method and choose CG. Then deselect the Fix Central Cell Shape;

   • Go to Interaction control ➟ select TotalEnergy and type 1e-04

   • Click Generate files

The Device Studio will create a new folder containing the Relax.input file.

Note

A vacuum region of 8 Å was set along the a and b axis to avoid interaction between periodic images.

Now, we are ready to start the calculation. However, before doing that we’ll briefly present the simulation parameters specified in the generated file.

Electrode optimization file

%%What quantities should be calculated
calculation.name = relax
calculation.relaxation.method = 'CG'
calculation.relaxation.forceConvergenceCriterion = 0.03 eV/angstrom
calculation.relaxation.fixCentralCellShape = false

%Basic setting
calculation.occupationFunction.temperature = 300
calculation.realspacegrids.E_cutoff = 80 Hartree
calculation.xcFunctional.Type = LDA_PZ81

calculation.k_spacegrids.number = [ 1 1 100 ]'
system.centralCellVectors = [[8 0 0]' [0 8 0]' [0 0 2.5]']
system.spinType = NoSpin

%Iteration control
calculation.SCF.monitoredVariableName = {'rhoMatrix','hMatrix','totalEnergy','bandEnergy','gridCharge','orbitalCharge'}
calculation.SCF.convergenceCriteria = {1e-04,1e-04,1e-04,[],[],[]}
calculation.SCF.maximumSteps = 200
calculation.SCF.mixMethod = Pulay
calculation.SCF.mixRate = 0.01
calculation.SCF.mixingMode = H
calculation.SCF.startingMode = H
%calculation.SCF.donatorObject = NanodcalObject.mat

%Basic set
system.neutralAtomDataDirectory = '../'
system.atomBlock = 1
AtomType OrbitalType X Y Z
Au  LDA-DZP 4.00000000      4.00000000      0.00000000
end

The keywords specify the following:

• calculation.relaxation.method The method used for the structure relaxation. The NanoDCAL carries out structural relaxation by damped Newtonian dynamics and/or conjugate gradient (CG) schemes.

• calculation.relaxation.forceConvergenceCriterion The criterion used to determine relaxation convergence (maximal absolute force component).

• calculation.relaxation.fixCentralCellShape If true, the central cell shape will not be changed during the relaxation.

Note

The user should manually include the flags below to relax the cell volume while fixing the cell shape. I that case the structure possess the tetragonal unit cell a = b \(\neq\) c, and \(\alpha=\beta=\gamma=90^{\circ}\)

calculation.relaxation.stressConvergenceCriterion = 0.1 GPa
calculation.relaxation.history = true
calculation.relaxation.fixedCellScale = {1, 2, 3}

These keywords could be changed, depending on the system:

• calculation.relaxation.stressConvergenceCriterion The convergence criterion for the stresses (maximum component of the stress per volume). If the unit string is missing, it will be considered to use the unit implicitly defined by the parameters calculation.control.energyUnit and lengthUnit.

• calculation.relaxation.history If true, a file searchingHistoryOfTheStructures.txt is generated to record the searching history of the structures.

• calculation.relaxation.fixedCellScale Relax the unit cell volume under different external mechanical pressures. The three lattice vectors (\({i, j, k}\)) shall vary at the same rate during the relaxation.

To perform the calculation download the input file Relax-1Dgold-unitCell.input or use Device Studio as follow:

1. Go to the Project window, open the Crystal directory ➟ and right-click on scf.input to pop up the Run window.

   • In the Gateway location choose the machine;

   • In the Run parallel on select the number of cores to do the calculation.

2. Press Run.

After the calculation finished, the Atoms_relaxed.xyz and log.txt files will be generated. The optimized structure can be loaded in Device Studio.

Note that the vacuum distances slightly deviates from the initial value and the Au lattice parameter was found to be \(d_{Au-Au}\approx 2.54\) Å, in agreement with previous calculations [SDCS].

Build and relax the scattering region

Next, we have to build the scattering region. We’ll use the relaxed Au unit cell as the basis of the device, i.e., we’ll build a supercell from the relaxed one and include the \(-C\equiv C-\) between the Au leads, as we shown in Fig. 1.1.2:

1. Select the relaxed Au unit cell, go to the menu bar and click Build ➟ select Redefine Crystal:

   • In the Redefine Crystal window ➟ expand the unit cell 14 times along the c vector ➟ Click build.

2. Select the Au supercell, go to the toolbar and click Convert to Crystal:

   • In the Convert to Crystal window, go to the Primitive Vectors tab and add a vacuum region along the the z direction setting the c vector length to 70 Å ➟ press Build.

3. Let’s change the two Au atoms from the right end of the chain. Thus, select these Au atoms, go to the toolbar and click Replace Atom . In the PeriodicTable ➟ click on C atom.

From theoretical works we know that the Au \(-C\) and \(-C\equiv C-\) bonds are around 1.95 Å and 1.34 Å, respectively [SDCS]. So we’ll set manually the initial position of the \(-C\equiv C-\) bond from the electrode as follows:

• Click on the last but one C atom, go to the panel properties and set z value as 29.91 Å.

• Click on the last C atom, go to the panel properties and set z value as 31.34 Å.

Fig. 1.1.2 The coupling of the C-C bond at the Au chain.

So far there is only the left side of device, we have to mirror the Au chain to complete the scattering region, as shown in Fig. 1.1.3.

1. Select all Au atoms from the chain, go to the toolbar and click the Mirror Atom button;

   • In the Mirror Atom window, click on the last box of the Normal vector row and select the xy plane, go to the Distance box and set 2.545 Å, select copy and click Apply.

2. Go to the toolbar and click on Fit Cell Along c. Select and delete one of the Au atom from the cell border;

Fig. 1.1.3 Schematic diagram of mirror symmetry.

Once the scattering region is built we can generate the relaxation input file:

• In the menu bar of Device Studio, click Simulator ➟ NanoDCAL ➟ Relax to pop up the interface Relax Calculation for Crystal to set the parameters;

  • Go to the Relax ➟ Method and choose CG.;

  • Go to Interaction control ➟ select TotalEnergy and type 1e-04

  • Click Generate files

Important

In this example we’ll relax just few central atoms from the chain. This procedure is quite common since we have always to ensure that all structural bulk characteristics will be preserved next to the interface between the scattering and electrode regions.

The generated chain is comprised by two carbon atoms and nineteen Au atoms. The input file will be edited to relax the two central C atoms and two Au nearest neighbors (NB = 2).

The following keywords should be included in the file:

calculation.fixCentralCellShape = true
calculation.relaxation.movingAtoms = [12:14, 25]

The keywords specify the following:

• calculation.fixCentralCellShape If true, the central cell shape will not be changed during the relaxation;

• calculation.relaxation.movingAtoms The \(3\times n\) coordinates corresponding to the atoms listed in the movingAtoms will be allowed to change during the relaxation, and the others will be fixed.

Run the calculation using this input file (Fix-Relax-atoms.input).

1.1.2. Transport calculation

The required input files for transport calculation could be generated using Device Studio with NanoDCAL module simulator, as follows:

• In the menu bar of Device Studio, click Simulator ➟ NanoDCAL ➟ SCF Calculation to pop up a new window to set the calculation parameters:

  • In this example the default parameters will be used, except for the Electron temperature option, which will be set as 300 K.

  • Click Generate files:

The Device Studio will create new folders that contain the scf.input files of LeftElectrode, RightElectrode and Device.

1.1.3. Electrode SCF calculation

First, we have to perform an SCF calculation for the leads. Then, the physical properties from the leads can be used as boundary conditions applied to the central scattering region of the device. Before starting the calculation, we’ll briefly present the simulation parameters specified in the left and right electrode files.

Left electrode input file

%%What quantities should be calculated
calculation.name = scf

%Basic setting
calculation.occupationFunction.temperature = 300
calculation.realspacegrids.E_cutoff = 80 Hartree
calculation.xcFunctional.Type = LDA_PZ81
calculation.k_spacegrids.number = [ 1 1 100 ]'

system.centralCellVectors = [[8.14311 0 0]' [0 8.14311 0]' [0 0 7.6341]']
system.spinType = NoSpin

%Iteration control
calculation.SCF.monitoredVariableName = {'rhoMatrix','hMatrix','totalEnergy','bandEnergy','gridCharge','orbitalCharge'}
calculation.SCF.convergenceCriteria = {1e-04,1e-04,1e-04,[],[],[]}
calculation.SCF.maximumSteps = 200
calculation.SCF.mixMethod = Pulay
calculation.SCF.mixRate = 0.1
calculation.SCF.mixingMode = H
calculation.SCF.startingMode = H
%calculation.SCF.donatorObject = NanodcalObject.mat

%Basic set
system.neutralAtomDataDirectory = '../'
system.atomBlock = 3
AtomType OrbitalType X Y Z
Au  LDA-DZP 4.07155631      4.07155631      1.27236151
Au  LDA-DZP 4.07155631      4.07155631      3.81708449
Au  LDA-DZP 4.07155631      4.07155631      6.36180749
end

Right electrode input file

system.atomBlock = 3
AtomType OrbitalType X Y Z
Au  LDA-DZP 4.07155631      4.07155631      6.36173849
Au  LDA-DZP 4.07155631      4.07155631      3.81701551
Au  LDA-DZP 4.07155631      4.07155631      1.27229251
end

All of these keywords were already described here. Note that the left and right electrodes are equivalent. For example, we just have to perform the SCF calculation for the left electrode and repeat NanodcalObject.mat for the right one. Note that the electrode Hamiltonian is stored in this file and is required for the transport calculation. Use the left (or right) electrode input file and run the scf calculation.

1.1.4. Two-probe SCF-NEGF calculation

Having calculated the electrodes, we move forward to carry out the two-probe transport calculation.

Scattering region input file

%%What quantities should be calculated
calculation.name = scf

%Basic setting
calculation.occupationFunction.temperature = 300
calculation.realspacegrids.E_cutoff = 80 Hartree
calculation.xcFunctional.Type = LDA_PZ81
calculation.k_spacegrids.number = [ 1 1 1 ]'

%Description of electrode
system.numberOfLeads = 2
system.typeOfLead1 = left
system.voltageOfLead1 = 0
system.objectOfLead1 = ../LeftElectrode/NanodcalObject.mat
system.typeOfLead2 = right
system.voltageOfLead2 = 0
system.objectOfLead2 = ../RightElectrode/NanodcalObject.mat

%Contour integral
%calculation.complexEcontour.lowestEnergyPoint = 1.5 Hartree
calculation.complexEcontour.numberOfPoints = 40
calculation.realEcontour.interval = 0.0272114
calculation.realEcontour.eta = 0.0272114

system.centralCellVectors = [[8.14311 0 0]' [0 8.14311 0]' [0 0 38.1269]']
system.spinType = NoSpin

%Iteration control
calculation.SCF.monitoredVariableName = {'rhoMatrix','hMatrix','totalEnergy','bandEnergy','gridCharge','orbitalCharge'}
calculation.SCF.convergenceCriteria = {1e-04,1e-04,1e-04,[],[],[]}
calculation.SCF.maximumSteps = 200
calculation.SCF.mixMethod = Pulay
calculation.SCF.mixRate = 0.01
calculation.SCF.mixingMode = H
calculation.SCF.startingMode = H
%calculation.SCF.donatorObject = NanodcalObject.mat

%Basic set
system.neutralAtomDataDirectory = '../'
system.atomBlock = 16
AtomType OrbitalType X Y Z
Au  LDA-DZP 4.07155631      4.07155631      1.272430510
Au  LDA-DZP 4.07155631      4.07155631      3.817153500
Au  LDA-DZP 4.07155631      4.07155631      6.361876490
Au  LDA-DZP 4.07155631      4.07155631      8.906599490
Au  LDA-DZP 4.07155631      4.07155631      11.45132250
Au  LDA-DZP 4.07155631      4.07155631      13.99604550
Au  LDA-DZP 4.07155631      4.07155634      16.52975129
C   LDA-DZP 4.07155631      4.07155631      18.40918156
C   LDA-DZP 4.07155632      4.07155632      19.68854959
Au  LDA-DZP 4.07155631      4.07155631      36.85451207
Au  LDA-DZP 4.07155631      4.07155631      34.30978908
Au  LDA-DZP 4.07155631      4.07155631      31.76506609
Au  LDA-DZP 4.07155631      4.07155631      29.22034309
Au  LDA-DZP 4.07155631      4.07155631      26.67562008
Au  LDA-DZP 4.07155631      4.07155631      24.13089708
Au  LDA-DZP 4.07155631      4.07155634      21.59719129
end

The keywords for central region were already explained here.

Attention

It is very important to check the alignment of the leads and the device before moving forward to production runs. To do this, use the command:

nanodcal -structure scf.input

This will produce a file Atoms.xyz which can then be loaded into an atomic visualization program.

With the generated input files (Left-elec-Au-C-scf.input, Right-elec-Au-C-scf.input, Central-region-Au-C-scf.input), we can perform the calculation for the two-probe device.

Thus, in the Device Studio follow the steps:

1. Go to the Project window, open the Device directory ➟ and right-click on scf.input to pop up the Run window.

   • In the Gateway location choose the machine;

   • In the Run parallel on select the number of cores to do the calculation.

2. Press Run.

Important

In order to analyze of the results, we also perform the calculation for the pure gold chain. The generated input files for the pure gold example are available to download below:

• Left electrode scf.input (source) (Left-elec-Au-scf.input);

• Right electrode scf.input (drain) (Right-elec-Au-scf.input);

• Central region scf.input (Au \(-\) Au) (Central-region-Au-scf.input);

Once the self-consistent two-probe calculation is finished, it is time to compute fundamental quantum transport properties.

1.1.5. Effective potential analysis

For the convergence of transport properties, we have to be sure that the central region is a semi-infinite region with a large number of cells from the electrode. Thus, once the potential is matched at the boundary, charge density automatically goes to the bulk-electrode values at the boundaries. This could be an issue if these conditions are not satisfied and the user should increase the central and electrode regions to be sure that the effective potential plot has periodic features, typical of bulks systems. Next, we will plot the calculated effective potentials for both Au \(-C\equiv C-\) Au chain and pure gold chain.

In order to calculate the potential along the transport direction we use the analysis interface of Nanodcal simulator, also present in the Device Studio software:

1. Click on Simulator ➟ NanoDCAL ➟ Analysis to set physical quantity;

   • In the Analysis window, on the top left choose the Device option;

   • Navigate to the Analysis panel (left side) ➟ select Potential ➟ click on the Right arrow button to transfer Potential calculator to the Calculation Selected panel.

2. The Potential options are shown in the right side whatPotential window. Choose the Effective and click Generate files.

The generated file is showing below.

Effective potential input file

system.object = ../NanodcalObject.mat
calculation.name = potential
calculation.potential.whatPotential = 'Effective'
calculation.control.xml = true
calculation.control.dsf = 1

The keywords specify the following:

• calculation.potential.whatPotential If ‘Hartree’, ‘Coulomb’, ‘Vh’, ‘Vdh’, the \(delta\)-Hartree potential will be provided. If ‘Effective’, ‘Veff’, the effective potential, i.e. summation of \(delta\)-Hartree and exchange-correlation potential, will be provided. If ‘ExchangeCorrelation’, ‘Vxc’, the exchange-correlation potential will be provided. If ‘NeutralAtom’, ‘Vna’, the summation of all neutral atomic potential will be provided.

• calculation.control.dsf 3D visual analysis.

Run the calculation and visualize the results using the generated EffectivePotential.mat, EffectivePotential.xml and EffectivePotential.dsf files. Those results could be plotted with MATLAB (EffectivePotential.mat) or loaded in Device Studio software:

• In menu bar of Device Studio, click on Simulator ➟ NanoDCAL ➟ Analysis Plot ➟ import EffectivePotential.dsf file.

The cross-sectional view of the electrostatic potential in Fig. 1.1.4 shows that the pure Au nanowire gives rise to a periodic and smooth potential. However, for the Au \(-C\equiv C-\) Au device, the electrostatic potential is very low on the two carbon atoms. There are some important physical insights that we can track from these results.

• The perturbation that is clearly relates to the carbon atoms, can interpret as a perturbation that does not extend far into the gold wire, which is important to get a really bulk-like region to couple the electrode;

• The electrostatic potential of the junction is affected by the charge transfers across the Au \(-\) C interface;

• The electrostatic potential can be also used to explain the difference in the conductance between the devices. As the charge transfers across the Au-C interface reduce the locally electrostatic potential. The conductance might exhibit the same behavior.

Fig. 1.1.4 Effective electrostatic potential for both devices.

1.1.6. Transmission spectrum

The transmission of electrons can be understood as the probability of the electron with given energy transmits through the device. The transmission.input file can be generated following the steps below:

1. In menu bar of Device Studio, click on Simulator ➟ NanoDCAL ➟ Analysis to set physical quantity;

   • In the Analysis window, on the top left choose the Device option;

   • Navigate to the Analysis panel (left side) ➟ select Transmission ➟ click on the Right arrow button to transfer Transmission calculator to the Calculation Selected panel.

2. The TransmissionSpectrum options are shown in the right side of the Analysis window. Edit the k-point sampling set to n1 = 1, n2 = 1 and n3 = 1 and click Generate files.

Transmission input file

system.object = NanodcalObject.mat
calculation.name = transmission
calculation.transmission.kSpaceGridNumber = [ 1 1 1 ]'
calculation.transmission.energyPoints = -10:0.05:10
%calculation.transmission.plot = true
calculation.control.xml = true

The keywords specifies the following:

• calculation.transmission.kSpaceGridNumber the small k-space grid number in each direction which, together with kSpaceGridShift, are used to produce the parameter kSpacePoints.

• calculation.transmission.energyPoints The energy points at which the transmission will be calculated. Note that the energy values are measured from chemical potential of a lead having zero applied voltage.

• calculation.bandStructure.coordinatesOfTheSymmetryKPoints The k-space fractional coordinates of those high symmetry k-points listed in symmetryKPoints, each column corresponds to one of them.

Now we can run the simulation with the transmission input file Au-wire-Transmission.input.

Next, you should perform the analyses of the results from the produced transmission dataset (.xml and .mat) which can be plotted using MATLAB or Device Studio as follows:

• In menu bar of Device Studio, click on Simulator ➟ NanoDCAL ➟ Analysis Plot ➟ import Transmission.xml file.

The transmission spectrum for a perfect Au nanowire is shown in Fig. 1.1.5. In this system, there is a quantized transmission (step-like shape) and the height of each step corresponds to the number of available conductance channels. The user can also understand that in the coherent regime, the transmission function is proportional to the number of bands at a certain energy. On the other hand, in the presence of the \(-C\equiv C-\) bond we found a huge reduction in the transmission spectrum along the whole energy range (-4 eV to +4 eV), caused by the misalignment of electronic states between the conductor and the leads.

Fig. 1.1.5 Transmission spectrum for both devices

Addition relaxation and transmission spectrum

At this point we could also explore additional relaxation effects over our obtained transmission results. In fact, the number of atoms that should be free to relax in the scattering region will change depending on the system. For example, we could rerun the optimization procedure of the scattering region, allowing the relaxation of the two C atoms and four Au nearest neighbors (NB = 4). As one can see in Table 1.1.1, the Au \(-C\) and \(-C\equiv C-\) bond lengths present only slight deviations between the two relaxation procedures.

Table 1.1.1 Bond lengths (Å); NB: number of Au nearest neighbors free to relax.

NB	\(Au-C\)	\(C\equiv C\)
2	1.87	1.27
4	1.89	1.26

This result indicates that our first optimization was good enough to reproduce the geometrical parameter at the scattering region. In this context, the transmission spectrum have also the same features, as one can see in Fig. 1.1.6.

Fig. 1.1.6 Transmission spectrum when NB = 2 and 4.

[YBBA98]

1. Yanson, G. R. Bollinger, H. E. van den Brom, N. Agrait, and J. M. van Ruitenbeek. Formation and manipulation of a metallic wire of single gold atoms Nature 395 (1998), p. 783.

[SDCS] (1,2)

7. Sclauzero, A. Dal Corso, and A. Smogunov. Interaction of CO with an Au monatomic chain at different strains: Electronic structure and ballistic transport. Phys. Rev. B 85 (2012), p. 165411.

[TGJ01]

J. Taylor, H. Guo, and J. Wang. Ab initio modeling of quantum transport properties of molecular electronic devices. Phys. Rev. B 63, (2001), p. 245407.

[WHL06]

4. Waldron, P. Haney, B. Larade, A. MacDonald and H. Guo. Nonlinear Spin Current and Magnetoresistance of Molecular Tunnel Junctions. Phys. Rev. Lett. 96 (2006), p. 166804.

[WLG07]

4. Waldron, L. Liu and H. Guo. Ab initio simulation of magnetic tunnel junctions. Nanotechnology, 18 (2007), p. 424026; Derick Waldron, PhD thesis, McGill University (2007).