ABINIT Discussion Forums

Dear Marc Torrent/positron spectroscopy researchers,

We are carrying out positron ground state calculations on some semiconductors (positron=1). We use LDA-PAW datasets. The calculations are OK when using ixcpositron=1, but fail with ixcpositron=3 (or 31). The calculations stop just before printing the lifetimes in the output file (or writing the mkdenpos output in the log file). We do not understand why this is. The only difference between the two calculations is the value of ixcpositron.

At first sight it seems that it is because the electronic GS density was obtained previously with an LDA exchange-correlation functional. With ixcpositron=3 one invokes the "GGA" electron-positron correlation functional (zero positron density limit parametrized by Arponen & Pajanne). But in principle the choice of electron-positron correlation functional does not really depend on the exchange-correlation functional used for the electronic ground state.

Could you please comment on this?

Thanks for your time!
R. Saniz
U. Antwerp

Dear R. Saniz,

In ABINIT, the positron is not treated as a separate particle.
For the sake of simplicity, we have schoosen to consider the positron as an electron with an opposite charge.
Doing this, we call each routine of the code once for N+1 particles (the last one having an opposite charge).
In particular, the standard XC routines are called and they do not allow the use of different XC functional types.
For that reason, it is only possible to use a GGA electron-positron correlation functional with a GGA electronic XC functional.
This is clearly a restriction; you're right when you say that we should, in principle, be able to perform GGA electron-positron correlation + LDA electron-electron XC.
This rectriction could be eliminated in ABINIT provided that you modify the XC routines (Plane-wave + PAW)... a relatively small work.

You could try to impose ixc=11 (PBE-GGA) in your ixcpositron=3 run... but, doing this the electron XC is done within GGA... this is not consistent with the pure electronic DFT run...

Marc Torrent

Dear M. Torrent,

Thanks for your reply, it's being useful.

I have another problem, though. In one case my calculations stopped with the "STOP xcpositron: problem, negative GGA espilon !" message.

I am a bit puzzled by this. In the xcpositron.F90 routine one can see that eps=grhoe2(ipt)/nqtf2. I found out that the code stopped because grhoe2(ipt) turns negative (for ipt=495). But grhoe2 is supposed to be the square of the gradient of the electronic density rhoe. How can its square be negative for some ipt's?

Thanks in advance for your comments!
R. Saniz

Deae R. Saniz,
If you look at the code:

you can see that esp is a ratio beween grhoe2 (necessarily positive) and rhoe which is the electronic density.
This one should be positive ; it is "positived" before each call to xcpositron...
If the electronic density becomes negative, it is likely because:
1- you have vaccuum in your cell
I suppose you are stuying defects, right?
2- you have a to small number of plane waves...

If you are using norm conserving pseudopotentials, you should increase the ecut input parameter; if you are using PAW, you should increase pawecutdg.
If you are using NCPP, there is only one call to the xcpositron routine : in the rhohxc routine. You could try to look at the value of the rho_b(:) or rhonow(:,1,1) variable at i=495.
If you are using PAW, there a possible "bypass"... try to put "pawxcdev 0" i nyour input file this makes your calculation run longer but this uses other XC routines...

In any case, try to change your input data, in order to make this negative density disappear...

I'm sorry not to be clearer... but thi implementation was done several years ago...

Marc Torrent

Dear Marc,

Thanks for your reply. It seems that there might be a bug in one of the positron-related routines. In the xcpositron.F90 routine I added a line, so that if eps<0, it would write out the variables used for its calculation, i.e.

write(98,*) ipt,kf,rhoe,nqtf2,grhoe2(ipt),eps,npt

There turned out to be many ipt points for which eps<0. But as an example, the output for ipt=495 is

495 0.72885999273403157 1.30770967550927555E-002 1.58699992065288744E-004 -0.23614862691367328 -1488.0191475782958 8000000

So clearly grhoe2 is negative for this ipt.

Could you please comment on this?

Thanks for your time!
Rolando