Wavefunction plot

Dear Users,

I was trying to find the spatial localization of electrons trapped in an oxygen vacancy.

From task 53 (ELF 3D plot), plotting a selected plane with color maps, I found a nice blob at the vacancy site, with ELF values ranging from about 0 to about 1. The blob was absent in the same calculation with additional +2 charge, i.e. without the trapped electrons.

From task 63 (wavefunction 3D plot), for the vacancy eigenvalue, I've got a similar blob.
But, this time, the values ranged from 0.3605 to 0.3683.
I.e. the differences were very clear visually, and quite small numerically.

How should I understand that?
Should I interpret this values as the probability of finding the electron at particular point?
Or, should the values be normalized somehow?
Should the sum of the wavefunction over the cell equal 2, as the eigenvalue occupation is 2?

Thank you.
Best regards.
Andrew

Last edit: Andrew Shyichuk 2020-03-17

J. K. Dewhurst - 2020-03-17

Hi Andrew,

The ELF was developed for just this purpose, namely for making the bonding structure more visually obvious than just plotting |φ(r)|².

The wavefunctions used in task 63 are mod-squared and normalised to the unit cell, i.e. the function integrated over the unit cell volume will give 1.

Strictly speaking, they should be normalised to the Born-von Karmen cell but it's more useful quote an obseverable as per-unit-cell.

Regards,
Kay.

Andrew Shyichuk - 2020-03-17

Dear Kay,

Thus, can I interpret the wavefunction at point x,y,z as the probability of finding the electron at this point?
Please consider the attached image. My eves say, the electrons are localized at the vacancy site.
The numbers say, not quite.

And, how do I calculate the integral over cell from a WF3D.OUT file?
Because, with either 50x50x50 plot grid, or 100x100x100 plot grid, or 300x300x300 plot grid, I still get most of the values in the file roughly equal 0.36, i.e. they do not depend on the voxel size.

Thanks.
Andrew

- Last edit: Andrew Shyichuk 2020-03-17
  
  WF3D_plane.png

J. K. Dewhurst - 2020-03-17

Hi Andrew,

Yes: you can interpret the mod-square of the wavefunction as the probability of finding an electron at (x,y,z).

You can approximate the integral by summing the values in WF3D.OUT, dividing by the total number of points and multiplying by the unit cell volume (found in LATTICE.OUT).

Regards,
Kay.

Andrew Shyichuk - 2020-03-17

Dear Kay,

From your description of the approximate integral, the WF3D.OUT values must depend strongly on the number of points.
In my case, they barely do. Why is that so?

Thank you.
Andrew

J. K. Dewhurst - 2020-03-17

The values will be mostly independent of the number of points, but there will be N of them. Thus if you sum them up you also have to divide by N. This gives the average over the unit cell. To get the integral over the unit cell you then have to multiply by the volume.

Andrew Shyichuk - 2020-03-17

I did that, and the integral turned out to be the total number of electrons in the cell, 2648.

Thus, the average value (0.36068) is the cell-average probability to find one of the 2648 electrons in a given voxel.

My bigger concern, however, is the way this plot changes, how much delocalized it is.
I.e. I'd like to have some estimate of the localization character from the wavefunction plots.

For instance, for eigenvalue 1, which was Lu 5s, occupied, it's peak is 0.4750, while the background (in the yz plane containing the atom) is 0.3607 (see the attached image).
Would that be correct to say that 0.3604 (which, multiplied by the cell volume, gives 2646, total number of electrons minus 2) is the background, and thus the probility for the 5s electrons in question changes (in the presented plane) from 0.0003 to 0.1143?

- WF3D_Lu5s_plane.png

J. K. Dewhurst - 2020-03-17

Hi Andrew,

There is a bug in the wavefunction plotting routine. I'm not sure when it was introduced but it's because the charge density routine is used to calculate the mod-square of the wavefunction and then normalises it to the total charge.

You can apply a tempory fix in wfplot.f90 by changing the lines

! select a particular wavefunction using its occupancy occsv(:,:)=0.d0 occsv(ist,ik)=1.d0/wkpt(ik) ! no symmetrisation required nsymcrys=1

to

! select a particular wavefunction using its occupancy occsv(:,:)=0.d0 occsv(ist,ik)=1.d0/wkpt(ik) ! total charge should be 1 chgtot=1.d0 ! no symmetrisation required nsymcrys=1

The wavefunction will be normalised to 1.

Thanks for pointing this out.

Regards,
Kay.

Andrew Shyichuk - 2020-03-18

Dear Kay,

I introduced the fix. The values sum up to almost 1.
For the 5s example, I did get the range similar to (0.0003, 0.1143), a localized state as it should be.
That is, the wavefunction was not only normalized in a different way, it was also shifted, and probably changed in more ways.
In other words, as a comment for the future users, the old files cannot be simply shifted and renormalized in order to get the correct plots.

Thank you.
Andrew

J. K. Dewhurst - 2020-03-18

Thanks Andrew.

I'll have it permanently fixed for the next release.

Regards,
Kay.

Youzhao Lan - 2021-01-19

Dear Kay,
Had this bug been fixed in the ELK 6.8.4?
I have'n seen the fixed line in wfplot.f90.

Youzhao Lan
China

Dear Youzhao,

I guess it was, I did some plots with Elk 6.8.4 and they were fine. You can test it easily by calculating the plot of interest and then checing the values inside a 3D file. Also, you can integrate the file and see if the integral is about 2 for spin-unpolarized or about 1 for spin-polarized calculation.

Best regards.
Andrew

https://sourceforge.net/p/elk/discussion/897820/thread/793128d708/