Negative Current obtained from positive bias Voltage

[Anirban Basak]

Hi,

       Recently I've tried to simulate I-V characteristics of Platinum contacted (10, 0) Zigzag CNT. I've used Huckel calculator with cedra basis sets and (3, 3, 100) k-point sampling. And the attached pictures show my result. The transmission spectra for each bias did not have any negative value (but gaps where the curves did not appear above Tr=0). The I-V analyzer interpolated those to some negative value. So why do I get negative current for some positive bias voltages?

Thank you.

[nori]

The transmission spectra for each bias did not have any negative value

Apparently there are negative values in the transmission spectra.
They are probably due to inaccuracy of Krylov method (and of course have no physical meaning).
That's why your problem would be solved by using Recursion method or Direct method instead.

[Anirban Basak]

Is there ANY way to get rid of negative transmission spectra in Huckel calculation? I did as suggested in above but did not work.

[Anders Blom]

Your electrodes are clearly too narrow in Z. For DFT they might be ok technically/numerically, although it's bare-bottom minimum. For Huckel, however, the basis set range is often considerably longer and you run into the problem described at http://quantumwise.com/publications/tutorials/mini-tutorials/99 about too short electrodes. Basically you will need to double the electrode in Z, otherwise you will continue to have convergence problems and/or negative current - the approximations you make with this electrodes are just too crude.

This, unfortunately, probably outweighs any performance benefit of the Huckel method for this system, and I would actually suggest you focus on DFT and add 1-2 layers of Pt in the central region to improve the screening (you can probably keep the electrodes as they are).

[Anirban Basak]

But Anders, Huckel requiring more basis set length goes against the electrode length checker script. Thus I'm puzzled. 

The DFT is proving to be way too much time consuming and a resource hog in my system. It fills up all my RAM and swap space and asks for more 

Also does nonSCF Huckel tries to converge at best possible state in one step while SCF Huckel just tries to match the tolerance in a number of steps?

EDIT: I mean are the nonSCF calculation and the first iteration of SCF calculation exactly same?

[Anders Blom]

This strongly depends on which Huckel basis you use. With the Muller/Hoffmann basis sets, the 3-layer electrode is indeed long enough, but have you tested what kind of band structure these parameters give for Pt? If instead you use the Cerda Pt basis set, then you get the same message as for DFT, i.e. that the electrodes are too short (even first-order interactions are truncated).

Now, you don't use the Cerda basis set so in principle this is not a (numerical) concern, but I'm not sure what kind of quality you can expect from the results with the Hoffmann/Muller Pt basis set.

I did a quick calculation with 9x9x9 k-points comparing the DOS and band structure between Hoffmann and Cerda. They are similar, but not identical. In particular for Hoffmann, the Fermi level lies on the band edge.

Yes, the non-SCF is pretty much just the first SCF step. There is a fundamental difference between models which are inherently self-consistent and those which are not. Many Slater-Koster models are for instance non-self-consistent and do not describe charge transfer properly when you run self-consistently. For Huckel device calculations, you should however use the self-consistent approach.

For reduce the size of the DFT calculation you can try a SingleZetaPolarized basis set for Pt, that should probably be fine enough for accuracy. It may also be possible to go to a 4x4 surface instead of 5x5. But for DFT you clearly need 6 electrode layers, so that will be a problem...

So, to conclude, I note that you do use the Cerda fcc Pt basis set. But when you check the electrode length, you probably used the default, which is Hoffmann - thus you didn't notice that you will 6 layers also in Huckel (with the Cerda basis set). Maybe you can try the Hoffmann basis set and hope that the Pt is decently described (good enough to capture the physics you are interested in, anyway). But for DFT or Cerda Pt, you will need 6 layers (or go to 110 or 100 surfaces, where you can do 4 layers).

[Anirban Basak]

Thank you Anders. I wrongly assumed for SE all basis has same cut-off. So I did not use to attach calculator while checking.

I see both Hoffman and Cedra DOS for Pt shows generally metallic characteristics. So I think I'd go for the faster one (Hoffman). BTW, what is Muller basis good for?

Also how lagre can I keep the convergence tolerance for a (well if not reasonably good) moderate accuracy in I-V analysis ? (It would be good if I can limit number of iteration for each bias around 10)
I used Hoffman and after around 10 iterations dH = 4 Hartree. The rate of dH decrement was also 1 Hartree/iteration. However, dE was comparatively lower like around 1e-2.

[Anders Blom]

We are actually changing the tolerance to 1e-4 in ATK 12.2, in our experience all physical quantities are very well converged already at that level, and it's just extra time to reach 4e-5 without any real advantages in terms of accuracy.

To reduce each calculation to 10 iterations is a bit unrealistic (nice, but not really possible), especially at finite bias, but another way to reduce the iteration count is increase the number of mixing steps; if you set it to 10, 15, even 20, it can have a large effect on the number of steps.

The Hoffmann and Muller basis sets are as far as I know primarily designed for molecules and isolated atoms, but you can check the original publications and those referring to them for hints of application areas.

[Anirban Basak]

Thank you Anders. You are amazing. 

BTW, say I have converged my calculation with 3x3x100 k-pts. While I calculate transmission spectra should I reduce them to 1x1x100 or 3x3x10 or keep 3x3x100 ? (the goal is again to speed up  without loosing much accuracy of course)

[Anders Blom]

This is a very important point!!!

There is actually no relation at all between the k-point sampling for the SCF loop and that for the transmission. (Except in the case of a 1D or 2D system, of course, where you have 1x1 or 1xN by pure symmetry.)

In the self-consistent loop, the k-point sampling is used to convert an integral of the real-space density into a sum. The density often doesn't vary very rapidly with k, so 3x3 up to 9x9 for "difficult" system is usually sufficient (although it needs to be checked, from case to case).

For the transmission, on the other hand, the original integral that you discretize is over self-energies and Green's functions, and these can have a very complicated k-point dependence, reflecting the nature of the tunneling mechanism. It may be that you have simple Gamma-point dominated "standard" direct tunneling through a quantum barrier, in which case perhaps 3x3 or 5x5 is enough, but it may also be that you see resonant tunneling, as in the case of FeMgO magnetic tunnel junction, where the main contributions to the total tunneling probability comes from isolated and extremely narrow peaks far out in the Brillouin zone. As it turns out you need 200x200 or 400x400 k-points in such a system to get an accurate current. In other cases you have a mix of these behaviors, and you will find convergence at 21x21 or 50x50 k-points.

The point is, that there is no immediate way to say beforehand how many points are sufficient, until you have analyzed the tunneling mechanism. So, for each system it is crucial to do a convergence study, and keep increasing the k-points until the result (the current, or the transmission spectrum itself) does not change significantly any more.

The tricky part here is, that what is sufficient at zero bias might not be enough at finite bias, so it's a good idea to check the convergence at say 0.5 V bias.

The good thing about the calculation of the transmission coefficients is that it scales pretty much linearly with the number of parallel MPI nodes, even up to hundred of nodes. So you can easily gain a factor of 10 or 20 if you have a parallel cluster to run the calculations on - in fact, even if you have a single machine (at least if you have 2 or more sockets) you may want to parallelize precisely the transmission part since it doesn't use so much memory.

[Anirban Basak]

Thank you very much for the enlightment on the number of k-points.

I'd like to state a trend I've seen in bias point calculation. A system which is difficult to converge, converges easily if the mesh cutoff is increased without changing k-points. However this thing does not scale linearly with bias points. Sometimes a high bias point may be easier to converge than low bias point.

I also have a question. Say I've converged upto some bias points with 50 Hartree cutoff. Now I face difficulty in convergence and try to increase the cutoff to 100 or such Hartree for further bias points. Does the wavefunction data that gets saved becomes bulky onwards (like it is saved for more closely spaced real axis points)? Any problem if I do this?

Another thing is that I had inspected Pt and ZZ CNT contacts and it converged for all biases (in one step! even when bias is changed) where the mesh cutoff was 10 Hartree. After calculating current I got them around 1E-20 A. I think it is non-physical. Should I increase the mesh-points and redo my calculation?

[Anders Blom]

The mesh cut-off is also used for rendering 3D grids like eigenfunctions and the real-space electron density. Internally, during the calculations, it is however only used for solving the Poisson equation; the actual basis orbitals are defined on a much finer grid. For this reason you can restart a calculation with a different mesh cut-off, both higher and lower.

Therefore, your question about "wavefunction data" is not relevant for the calculations as such, only if you produce some plot data (3D grids). Of course, the calculations with higher mesh cut-off will be more accurate and this can possibly create a kink in the data, and it might make sense to go back and rerun the lower bias points with the higher mesh cut-off. As mentioned above, you can restart from the converged calculation so hopefully it will converge in just a few steps.

As for your last question, more than likely it has converged to a state of zero charge, all electrons have left the system, and indeed the solution is unphysical. I don't think tuing the mesh cut-off alone will be enough to solve it, you may need to look at details like how long hte electrode is, if you have enough k-points, and even the self-energy evaluations method (Direct is more stable than Krylov).

[esp]

helpful discussion here, i too have negative currents in a graphene system .. i will have to read through this carefully