Simulation of STM and AFM images for two short graphene nanoribbons with different chemical termination

In case you do not want to use the quantum-mobile VM, you will need to install the asetk and ProbeParticle packages: 
git clone https://github.com/ltalirz/asetk
pip install -e asetk

and 

git clone https://github.com/ProkopHapala/ProbeParticleModel.git
cd ProbeParticleModel/
git checkout dev

download from here the tar file exercise_10.tar, move the file to your exercise directory, and extract the content 

tar -xvf exercise_10.tar
cd exercise_10

We consider two possible chemical terminations for a finite size 7-AGNR. In TASK_1 the ribbon is terminated with a C-H2 bonding while in TASK_2 the termination is C-H The additional H atom present at the termini of the ribbon of TASK_1 will suppress the spin polarized edge states that are evident in the ribbon of TASK_2

TASK_1

cd TASK_1

Have a look to the cp2k input file cp2k.inp used to obtain quickly the optimized geometry of a ribbon adsorbed on a Au substrate. The ribbon is modelled within DFTB (similar to tight binding) while the substrate is modelled via Embedded Atom Model. An empirical potential in the form of C6/R^6 plus a Pauli repulsion term are added to couple the adsorbate/substrate systems.

Two geometry fiels are present: mol.xyz and all.xyz The cp2k program will need both of them.

Have a look at the geometry of the system using ASE or vmd for both all.xyz and mol.xyz: 

ipython
In [1]: from ase.io import read

In [2]: from ase.visualize import view

In [3]: s=read("all.xyz")

In [4]: view(s)

In [5]: exit()
submit the geometry optimisation run 
./run

Have a look at the final geometry optained (you can extract the last frame from the file PROJ-pos-1.xyz) After completion of the optimization you should extract the final coordinates of the molecule (first 80 atoms) and copy them in the STM directory (call them p.xyz) to compute the KS orbitals and to compute the STM images 

tail -442 PROJ-pos-1.xyz | head -82 > p.xyz
mv p.xyz STM

Now go to the directorySTM 

cd STM

and have a look to the input file cp2k.inp used to converge the wavefunction of the system and to plot the cube files for the KS orbitals. Execute the program 

cd STM
./run

The program will compute the 4 highest occupied and 4 lowest unoccupied KS orbitals. Visualize the orbitals with VMD (remember + and -)

To obtain the STM images you have to combine different KS orbitals (depending on the bias voltage applied) into a single cube file: 

./run_sumbias

you will then obtain a cube file for each desired bias voltage (see the script run_sumbias)

Now you can compute a constant current STM image running the script 

./run_stm

Please note that we are simulating a molecule, we do not include the electrons of the substrate thus we have a discrete spectrum of energies

why some of the STM images look empty?

Now we can simulate for the same ribbon a nc-AFM image:

Go the the AFM directory of TASK_1 (and have a look to the parameter file params.ini) copy there the p.xyz file that you havein the STM directory and execute:
./run_PP

It will take ~ 5 minutes, then you will find a dir containing the AFM simulated image.

TASK_2

Modify the geometry of TASK_1 removing one H atom from each C-H2 at the termini of the ribbon (remove two H atoms in total). Create the corresponding mol.xyz and all.xyz files, optimize the geometry, compute STM and nc-AFM images repeating all the instructions of TASK_1 for the scripts present in the dir TASK_2

Be careful: here we do a spin polarised simulation. When doing the STM simulation (ONLY for the STM) we have to distinguish the three C atoms of one terminus of the ribbon from the three of the opposite terminus calling them C1 and C2. For these atoms we will define a guess electronic configuration with spin up on one side and spin down on the opposite side. This is achieved defining a occupation unbalance in the alpha and beta orbitals (try to identify this section of the input and note that the calculation is performed for a spin multiplicity of 1)

The file p.xyz in teh STM directory should look similar to: 

     78 
 i =       49, E =      -140.2738100175
  H         4.2914607718       10.2130614763       21.2815435017
  H         4.2778729017        7.7509218987       21.2954986738
  H         4.2782723704       12.6751364096       21.2955091278
  C         7.4704758534        6.5236639413       21.2352922076
  .
  .
  .
  C1        5.3788157746        7.7465647443       21.2687198580
  .
  .
  C1        5.3936844253       10.2129317839       21.2797918647
  .
  .
  C1        5.3792136407       12.6792819903       21.2687263656
  .
  .
  .
  C2       21.1530397078        7.7456205579       21.2687376504
  .
  C2       21.1385072480       10.2118877383       21.2797955201
  .
  C2       21.1533012965       12.6781430186       21.2687326678
  .
  .
Look at the KS orbitals (especially HOMO and LUMO) for both spin UP and DOWN
Notice the difference between the images in TASK_2 and the images in TASK_1 In TASK_2 we have KS states localised at the termini of the ribbon. These states are suppressed by the addiitonal H atoms in TASK_1
why some STM images are remarkably asymmetric? Is this correct?