[gpaw-users] problem with LCAO coefficients - probably in spherical harmonics

[Ondřej Krejčí]

Dear GPAW developers!



I'm trying to calculate a STM current by means of the Chen rules and wave-
functions restored from LCAO coefficients (and some empirical exponential as
an radial function).  This works pretty well when the Pi orbital of 
molecules (pz orbitals of C, N etc.) is dominating in the current, however 
when I try to implement d-orbitals (for some metal organic molecules), I get
some really asymmetric signal, even though the wavefunction of the eigen-
state looks absolutely symmetric in the cube file. I looked into the files 
with radial functions and all the valence radial functions seems to be 
positive (especially in the far distance), therefore it is in agreement with
my STM code.

Thus, I'm thinking that there is some different sign convention in the  
spherical harmonics for the GPAW basis set. For me: p_y = Np*y/r, p_z=Np*z/
r,  p_x = Np*y/r, Nd*(x*y)/r^2 ... where N are just normalization constants 
> 0. Am I right, or the mistake is somewhere else? For other LCAO codes 
(Fireball, FHI-AIMS) I'm getting proper and symmetrical current maps.




If somebody would like to see the way I'm obtaining the LCAO coefficients, I
put it here:





n_max_ ... last eigen-state in the calculations

n_min_ ... first eigen-state in the calculations


num_at_ ... number of atoms in the calcultions

Ynum = 4 for 'sp' orbs, 9 for 'spd' orbs (only valence)




coef = np.zeros((n_max_-n_min_,num_at_,Ynum_))

if (orbs=='spd'):

from gpaw.setup_data import SetupData


chem_sym=slab.get_chemical_symbols()

d_orb=np.zeros((num_at_));

for i in range(num_at_):

if at_num[i]>2:

setup =  SetupData(chem_sym[i],'LDA','paw');l_j = setup.l_j;tmp=l_j[:l_j.
index(2)];a=[1,3];oo=0;

for j in range(len(tmp)):

oo +=a[tmp[j]];

d_orb[i]=oo;

for i in range(n_min_,n_max_):

h=0

for j in range(num_at_):

ii = i-n_min_

coef[ii,j,0] = calc.wfs.kpt_u[0].C_nM[i,h]

if (at_num[j]>2):

coef[ii,j,1] = calc.wfs.kpt_u[0].C_nM[i,h+1]

coef[ii,j,2] = calc.wfs.kpt_u[0].C_nM[i,h+2]

coef[ii,j,3] = calc.wfs.kpt_u[0].C_nM[i,h+3]

if ((orbs=='spd')and(d_orb[j]>1)):

coef[ii,j,4] = calc.wfs.kpt_u[0].C_nM[i,h+d_orb[j]]

coef[ii,j,5] = calc.wfs.kpt_u[0].C_nM[i,h+d_orb[j]+1]

coef[ii,j,6] = calc.wfs.kpt_u[0].C_nM[i,h+d_orb[j]+2]

coef[ii,j,7] = calc.wfs.kpt_u[0].C_nM[i,h+d_orb[j]+3]

coef[ii,j,8] = calc.wfs.kpt_u[0].C_nM[i,h+d_orb[j]+4]


h += calc.wfs.setups[j].nao





This procedure sometimes gives me a crazy wavefunction, if Hydorgens are not
the last  atoms of the molecule, it seems that - h += calc.wfs.setups[j].nao
- is sometimes shifted by one, but this can be easily overcome.




Thank you for any help and suggestions!

All the best,

Ondrej Krejci
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/gpaw-users/attachments/20170216/a780e496/attachment.html>