# 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.

The calculation steps are given as follows:

Structural Relaxation:

Build and relax the electrode unit cell;

Build and relax the scattering region.

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:

Open Device Studio software ➟

**Create a new Project**➟ on File name type**Au-wire**➟ press**Save**;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**.

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:

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.

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:

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**.

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**.

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 Å.

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.

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.

Go to the

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

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:

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.

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:

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.

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.

## 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:

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.

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.

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.

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.

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.

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.

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.

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.

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).