spin-polarized transport calculation with magnetic electrodes

[wilson]

Dear all:
Does anyone has ever successfuly did spin-polarized transport calculation using Nickel or other magnetic metal as electrodes? I did some test calculation these days and I find that if I use spin-polarized Nickel or iron as electrodes, then no matter what kind of molecule I choose and what kinds of calculation parameters I used, this error will always appear. It seems that this is a bug of the program. This error bothers me a lot and I don’t know how to avoid this error right now. Can anyone help me?a

[Nordland]

I have performed a calculation with Nickel electrodes with a organic molecule in between.

Spin is simply harder. If you could share the output, I can properly give an hint or two.

[wilson]

Hi, Nordland, thanks for your reply!
I did another test calculation with a simplest molecule, that is , a benzenedithion between two Nickel 4*4(111) surface, the initial spin of the Nickel is set to 0.1 while others set to 0. The error also occurs. It seems a really trickey problem. The input and output files are both enclosed.

[zh]

One comment for your input file:
For a metallic system,  a very small value (e.g., 4K in your input file) of electron temperature in the Fermi-Dirac distribution may lead to a slow convergence of self-consistent calculation.

For a spin polarized calculation, the initial spin is also very important. The feasible procedure to do a spin-polarized calculations of a two-probe system:
i) do the spin polarized calculations for a bulk system of the material that constitutes the electrode. The aim of this step is to obtain the local magnetic moment of each atom in the electrode. These values will be good initial spins for the electrode atom in the spin-polarized calculations of the two-probe system.
ii) similarly, do the spin polarized calculations for a molecular system of the molecule that constitutes the scattering region.
One can also obtain the local magnetic moment of each atom in such molecule.
iii) After these two separated calculations, one starts to do the spin polarized calculations for the two-probe system, meanwhile the initial spins are set as  the local magnetic moment obtained in the above two steps.

Other tip:
It is well known that the magnetic ground state of iron (Fe) can not be described well by the LDA. That is to say, one had better choose the GGA for iron.

For the error message in your output file, i.e., "ATKError: inverse(DZMatrix const &) : Could not LU factorize!", maybe you should check the geometry structure of your two-probe system, e.g., are some of two atoms  too close to each other?

[zh]

I checked the geometry structure of your two-probe system. It is wrong. The followings are a xyz file for the center region of your two-probe system. Please visualize it to check the geometry structure:

66
##             
              Ni     0.00     0.00     6.10
              Ni    -1.25     2.16     6.10
              Ni    -2.49     4.32     6.10
              Ni    -3.74     6.47     6.10
              Ni     2.49     0.00     6.10
              Ni     1.25     2.16     6.10
              Ni     0.00     4.32     6.10
              Ni    -1.25     6.47     6.10
              Ni     4.98     0.00     6.10
              Ni     3.74     2.16     6.10
              Ni     2.49     4.32     6.10
              Ni     1.25     6.47     6.10
              Ni     7.48     0.00     6.10
              Ni     6.23     2.16     6.10
              Ni     4.98     4.32     6.10
              Ni     3.74     6.47     6.10
              Ni    -3.74     7.91     8.14
              Ni     0.00     1.44     8.14
              Ni    -1.25     3.60     8.14
              Ni    -2.49     5.75     8.14
              Ni    -1.25     7.91     8.14
              Ni     2.49     1.44     8.14
              Ni     1.25     3.60     8.14
              Ni     0.00     5.75     8.14
              Ni     1.25     7.91     8.14
              Ni     4.98     1.44     8.14
              Ni     3.74     3.60     8.14
              Ni     2.49     5.75     8.14
              Ni     3.74     7.91     8.14
              Ni     7.48     1.44     8.14
              Ni     6.23     3.60     8.14
              Ni     4.98     5.75     8.14
              H      2.27     3.08     9.41
              H      0.15     5.61    14.25
              H      0.15     5.61    11.75
              H      2.93     2.29    11.75
              H      2.93     2.29    14.25
              C      1.54     3.95    14.40
              C      0.76     4.88    13.70
              C      0.76     4.88    12.30
              C      1.54     3.95    11.60
              C      2.32     3.02    12.30
              C      2.32     3.02    13.70
              H      2.27     3.08    16.59
              S      1.43     4.08    16.17
              S      1.43     4.08     9.83
              Ni    -3.74     7.91    18.14
              Ni     0.00     1.44    18.14
              Ni    -1.25     3.60    18.14
              Ni    -2.49     5.75    18.14
              Ni    -1.25     7.91    18.14
              Ni     2.49     1.44    18.14
              Ni     1.25     3.60    18.14
              Ni     0.00     5.75    18.14
              Ni     1.25     7.91    18.14
              Ni     4.98     1.44    18.14
              Ni     3.74     3.60    18.14
              Ni     2.49     5.75    18.14
              Ni     3.74     7.91    18.14
              Ni     7.48     1.44    18.14
              Ni     6.23     3.60    18.14
              Ni     4.98     5.75    18.14
              Ni    -2.49     7.19    20.17
              Ni     1.25     0.72    20.17
              Ni     0.00     2.88    20.17
              Ni    -1.25     5.04    20.17

[wilson]

Thanks for your advice, zh.
For the electron temperature, the value I previously used was the default value (300K), after I came across this error, I change it to 1000K and then 4K under the advise of Anders Blom, but the error still not been eliminated.
Your advice for spin-polarized calculation is very important and valuable, I will follow it in my future calculation.
In anther topic "Error about runtime_parameters" of the forum, I have already talk about this error with Anders Blom and Nordland. Actually, I had already do a bulk calculation for the spin-polarized Ni (follow the instruction of the atk manual on spin polarized iron), the obtained spin polarization is about 0.33(I made a mistake in previous post). I have also tried many different initial spin for Ni, but the error still always unavoidable.

[wilson]

Hi,zh, I checked the xyz files, it seems that some Nickel atoms of the right surface is missing. My input file is all generated by VNL. I have just checked it again and didn't find significant structural errors, moreover, the calculation can complete correctly for exactly the same system without spin, so I guess the error should not due to the geometry structure.

[zh]

I am quite sure that the problem you met is caused by the wrong geometry structure. As clearly defined in your script file,  the center region of your two-probe system consists of 32 Ni atom in the left surface layers, 14 atoms for the molecule, and 32 Ni atoms in the right surface layers. So the printed xyz should include the coordinates of 78 atoms, rather than 66 atoms.

The problem you met is not related to the initial spin for Ni atoms defined your script file. Actually, it is caused by the improper input for the atomic coordinates in the center region. All the z-coordinates of these atoms should be shifted by 6.10374689, so that the z coordinates of the atoms in the first layer of the left side of center region are adjusted to be zero.

For example:

# Scattering coordinates
left_surface_coordinates  = [[  0.00000000e+00,   0.00000000e+00,   6.10374689e+00-6.10374689],
                             [ -1.24592221e+00,   2.15800047e+00,   6.10374689e+00-6.10374689],
                             [ -2.49184442e+00,   4.31600094e+00,   6.10374689e+00-6.10374689],
                             [ -3.73776650e+00,   6.47400141e+00,   6.10374689e+00-6.10374689],
                             [  2.49184442e+00,   0.00000000e+00,   6.10374689e+00-6.10374689],
                               .............

central_atom_coordinates  = [[  2.2682999 ,   3.08442769,   9.40920067-6.10374689],
                             [  0.14764471,   5.61160397,  14.25049067-6.10374689],
                             [  0.14755506,   5.61171081,  11.74980688-6.10374689],
                             [  2.93169612,   2.2937296 ,  11.74950635-6.10374689]
                             ......

right_surface_coordinates = [[ -3.73776650e+00,   7.91266823e+00,   1.81383300e+01-6.10374689],
                             [  1.55431223e-15,   1.43866694e+00,   1.81383300e+01-6.10374689],
                             [ -1.24592221e+00,   3.59666753e+00,   1.81383300e+01-6.10374689],
                             [ -2.49184442e+00,   5.75466776e+00,   1.81383300e+01-6.10374689],
                             [ -1.24592221e+00,   7.91266823e+00,   1.81383300e+01-6.10374689],
                             ......

After the above modification, the script file will work well.

[Anders Blom]

Zh, I do agree the geometry might be causing the problems, but where did you get your xyz list from? It's not what ATK itself produces; in the user's output file we see 78 atoms listed with ok coordinates.

Also, a pure shift of all z coordinates has no effect, "thanks" to the alignment atoms.

[wilson]

Hi, zh, Anders, bad news. In accordance with zh's suggestion, I modified the input file and re-do the calculation on this system. But the error still came out when the equivalentbulk calculation finished and the two-probe calculation starts. Following is the error message.

# sc 50 : Fermi Energy =   -0.29203 Ry  Etot = -13430.82223 Ry  dRho =  2.1265E-05  dEtot =  1.5297E-05 Ry
# sc 51 : Fermi Energy =   -0.29203 Ry  Etot = -13430.82224 Ry  dRho =  1.1465E-04  dEtot = -8.7508E-06 Ry
# ----------------------------------------------------------------
# TwoProbe Calculation
# ----------------------------------------------------------------
# sc  0 : q =  682.01332 e
Traceback (most recent call last):
  File "111niben2.py", line 312, in ?
    runtime_parameters = runtime_parameters
ATKError: inverse(DZMatrix const &) : Could not LU factorize!

[nori]

I think electrode_constraint is the main cause of unstable SCF in your calculation.
If you select "ElectrodeConstraints.RealSpaceDensity" instead of "ElectrodeConstraints.Off", the SCF convergence will be dramatically improved.
I've verified it with "111niben_2.py".
In "111niben_2.py", several other parameters are changed but probably these changes are trivial but only "ElectrodeConstraints.RealSpaceDensity" is indispensable for solving this issue.

[wilson]

Thank you, nori, your proposal is indeed effective.
I select "ElectrodeConstraints.RealSpaceDensity" and keep other parameters unchanged then re-do the calculation, the error message doesn't appears again.
But, just as Anders Blom said in another post,
"I generally advise against using electrodeConstraints.RealSpaceDensity. It's legacy method that used to be the old default,
but we have discovered a lot of things with it that work less than satisfactorily.",
I wonder whether the results obtained with this option (ElectrodeConstraints.RealSpaceDensity) is reliable? or this option will have any significant effect on the calculated results?

[Anders Blom]

The key is to consider the screening.

Actually all constraints are equivalent in the limit of infinitely many screening layers, but the Off constraint allows for an accurate calculation of the transmission with fewer screening layers. The voltage drop, however, is wrong.

RealSpaceDensity gives a good voltage drop AND a good transmission but only if you have a lot of screening layers. This method is also sensitive to how the atoms are placed in the electrode unit cell; most notably, the way VNL places them is bad for this method (that was the underlying reason for advising against it, because most people will not want to edit their scripts).

[wilson]

Hi, Anders, I did some test calculation these days, it seems that use "ElectrodeConstraints.RealSpaceDensity" is the most effective way to avoid this error.
Since the way VNL places them is bad for this method as you said, if I still want to use RealSpaceDensity, how should I arrange the atoms in the electrode unit cell?
Is there any handbook to tell me how to do this? Or I can do it follow the ATK manual?

[Anders Blom]

The trick to getting good results with RealSpaceDensity is to center the atoms in the unit cell, in the z direction. So, instead of having coordinates like

z=0 Ang
z=1 Ang
z=2 Ang

in a 3 Ang cell, one should use 0.5, 1.5, and 2.5 Ang.

I still think DensityMatrix is better, if you can get it to work...