Failing to dope a system by setting Geometry.Hartree

You are effectively only setting the Hartree potential in the specified region.

I.e. it isn't self-consistent Hartree solution with a boundary condition.

Basically Geometry.Hartree does:

V(r) = V(r) + Geometry.Hartree(r)

Where V(r) is solved using Poisson solution.

Thank you Nick.
Do you mean that V(r) is changed after SCF steps? Does Geometry.Charge follow the same procedure? Why would .Charge cause Fermi level shift in bandstructure?

I plotted Rho.gird.nc for both case (using Geometry.Hartree and Geometry.Charge), they seem unchanged compared with pristine case. It is understandable for the .Hartree case because you said it is not changed in the SCF calculation. But is the charge density also unchanged when setting Geometry.Charge? From the bandstructure, Geometry.Charge setting should have affected the charge in the slab since Fermi level is shifted from Dirac point.

Can I attract some charge from the atoms or give some additional charge to the atom plane using Geometry.Charge? How can I run SCF calculation with the additional potential?

Thank you very much for helping!

Best,
Yuefei

Sorry, I meant "Can I attract some charge from the atoms or give some additional charge to the atom plane using Geometry.Hartree? ".
Thank you again!
Yuefei

The equation

V(r) = V(r) + Geometry.Hartree(r)

is for every SCF step. So charges may attract but probably only if rho(r) and Geometry.Hartree(r) overlaps in real space.
So yes in that sense it is in the SCF loop. But remark, not as a boundary condition to the Poisson solution!

As for Geometry.Charge. That is *exactly* the same as NetCharge in the sense that the charge addition/removal is also present when solving the Poisson equation.
In this case the Rho.grid.nc is different from the non Geometry.Charge calculation. Please check this (it may not be visible so differences are required).
By using Geometry.Charge you can attract charges by adding planes. You may find this article https://pubs.rsc.org/en/content/articlelanding/2016/cp/c5cp04613k/unauth#!divAbstract interesting since it was implemented for exactly this purpose.

Thank you Nick!
So is the process like:
1. Rho(r) = Rho(r)+Geometry.Charge(r)
2. Rho(r) -> V(r)
3. V(r) = V(r)+Geometry.Hartree(r)
4. V(r)->Rho(r)
And step 2,3,4 is repeated?

So is a jellium background of charge is applied to the cell to maintain charge neutrality? Is the background jellium distribution determined by the Geometry.Charge()?

In TranSIESTA, if I set Geometry.Charge for two electrode and device differently at the space that they are supposed to match, i.e. electrode atoms matches with left/right ends of the device but Geometry.Charge does not match. How does transiesta handle such case?

Thank you very much!

Sorry I miss-clicked 'Problem is solved'. Thank you!

No, step 1, 2, 3 and 4 are all repeated. :)

Using Geometry.Charge will retain a charge neutral cell since the charge is moved (NOT removed) from your system to the defined region.

As for transiesta. It is YOUR responsibility that the charged gate is equivalent in the two simulations (TS will not check this). I.e. for a pristine graphene sample you would define a charge-gate in the electrode, then you would calculate the charge/area and then define a gate in the device with the same charge/area. However, sometimes you need to take more care depending on your geometry etc.
I.e. it is your responsibility to ensure that the charge-density and electrostatics of the electrode regions are the same in the device region as in your electrode calculation.

Thank you Nick!

For charge density and electrostatics, do you mean the output of Rho.grid.nc, ElectrostaticPotential.grid.nc and TotalPotential.grid.nc in the L/R of L-C-R should match the pure electrode outputs point by point?

I am confused with this, because when I do a calculation with two different electrodes, and plot the above three .grid.nc files from the Scattering Region calculation, it seems like they are periodic in transport (z) direction. Are these files from the initial SIESTA calculation actually? Because for TranSIESTA calculaton I thought the left and right boundary should not match each other. How can I output these grids from the TranSIESTA step?

When you do calculation with two different electrodes with different work function. Say contact of two different metals, under 0 bias, will the electrodes Hamiltonians be rigidly shifted w.r.t each other to match the work function?

Thank you very much!

1) No not point by point. It is just one way of defining the proper boundary conditions. How you actually check it is your choice, for instance the DOS should match in those regions.

2) The rho.grid.nc is output after transiesta, so that *is* calculated for the right boundary. However, transiesta does not really remove the boundaries as the simulation cell is still periodic. Hence you will still handle a periodic cell.
For cases with different L/R electrodes it is VERY important to either a) add buffer atoms, or b) add vacuum between the L/R (not in the device obviously).

3) The chemical potentials of each electrode will be shifted so that they coincide with the fermi-level of the device. As for different work-functions, this is where you need to do testing with my suggestion above (a or b option).

Thank you Nick!
I added a left electrode cell and a right electrode cell as buffer layer on the two ends. It did increase some computation time in the SIESTA step, but resulted in better convergence in TranSIESTA step. It is very helpful.

I guess there is a typo in the SIESTA manual S I E S T A 4.1-b3 July 04, 2017 about buffer atoms:
'%block TBT.Atoms.Buffer' should be '%block TS.Atoms.Buffer'

I met another problem when I was calculating transport of lateral junction formed by 2 different 2D materials at different doping concentration. My system is a slab in yz plane at x=10Ang, and two identical delta Geometry.Charge plane at x=5Ang and x=15Ang (to maintain symmetry in x direction). Left and Right are two different 2D materials.

In some doping concentration range, from about -2*10^-13 e /cm^2 to -0, I get converged scf calculations but 0 conductance. It was quite strange because when I check Rho and orbital occupation of my whole system, the electrons in device region are totally drained. Rho plot is zero in Device region and Atomic Population is exactly zero in Device region in the output file. In Electrode and Buffer region Rho and Atomic Population looks fine.

At higher n doping, result was fine (conductance was finite and Rho is reasonable). At p doping, result was fine. Problem appears at 0V bias and higher bias as well.

I tried several ways to fix the problem:
1. Change PAO.BasisSize from DZP to SZP. Turns out good. Here is some output compare of DZP and SZP:
DZP:
Before TranSIESTA cycle:
transiesta: Charge distribution, target = 649.13320
ts-q: D E1 C1 E2 C2 B dQ
ts-q: 397.446 92.589 2.681 124.336 -405.185 220.693 -2.091E+2
ts-Vha: 0.56442E+03 eV
siesta: FreeEng = 209004.472253

After 1st TranSIESTA step:
ts-q: D E1 C1 E2 C2 B dQ
ts-q: 0.028 92.589 -0.059 124.336 -0.009 220.693 -2.041E+2
ts-Vha: 0.37213E+03 eV

SZP:
Before TranSIESTA cycle:
transiesta: Charge distribution, target = 649.13320
ts-q: D E1 C1 E2 C2 B dQ
ts-q: 198.808 94.783 1.264 126.581 1.372 222.905 -4.167E-1
ts-Vha: 0.72416E+00 eV
siesta: FreeEng = -24528.920797

After 1st TranSIESTA step:
ts-q: D E1 C1 E2 C2 B dQ
ts-q: 212.629 94.783 3.554 126.581 3.669 222.905 1.799E+1
ts-Vha: -0.24832E+02 eV

2. Run SIESTA steps only, set 'TS.SIESTA.Only .true.', Rho output is good. So the electrons are drained in the TranSIESTA cycle.

3. Added 'TS.Elecs.DM.Init bulk’ and ‘TS.Elecs.DM.Update none’ . This doesn't help.

4. Added 'TS.ChargeCorrection fermi' and 'TS.ChargeCorrection.Factor 0.75'. Doesn't help.

Have you ever seen such problem? Can you give me some idea on why all valence electrons would leak out of the device region? In SZP basis and other doping concentrations this would not happen, which is more strange. How can I solve this problem.?

Thank you very much!

Best,
Yuefei

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

1) "I guess there is a typo", this has already been fixed. Please do not use b3, immediately bump your used version to the latest beta (if so desired). The beta releases are quickly changing and they are not "stable" in the sense that one should stick with one version.

2) If you *ever* get a zero charge in your SCF device region then for sure you have a non-converged system. In fact you are doing something wrong. Also, when you have a Geometry.Charge region for transiesta, it is VERY important that your electrodes have the same charge density Geometry.Charge in the electrode calculation as in the device calculation. Otherwise you are forcefully attaching two systems at different chemical potentials which is quite wrong.

3) The first step in transiesta is unimportant since it always gives a kick. What is important is the final SCF step and whether you have a dQ which is relatively low.

4) If the electrons leak out it is a good indication that your system is erroneously setup since for 0 bias siesta and transiesta should be very similar in SCF iterations.