[gpaw-users] Problems with LDOS calculations

[Caro Miguel]

Hello again,

As a follow up, if the WS LDOS is computed just with the assertion commented out, the complex part of the wavefunction is left out of the calculation of the weights, and the sum over local DOS does not yield the total DOS exactly: there is a small discrepancy between the two values (http://mcaroba.dyndns.org/misc/gpaw/LDOS-WS-float.png).

With this very simple modification to the code in the relevant block of gpaw/utilities/dos.py

            P_p = pack(np.outer(np.conjugate(P_i), P_i))

where the complex conjugate of projector P_i is introduced, the sum of local DOS now seems to yield the total DOS exactly (http://mcaroba.dyndns.org/misc/gpaw/LDOS-WS-complex.png).

Regards,
Miguel
________________________________________
From: gpaw-users-bounces at listserv.fysik.dtu.dk [gpaw-users-bounces at listserv.fysik.dtu.dk] on behalf of Caro Miguel [miguel.caro at aalto.fi]
Sent: 26 October 2013 14:54
To: gpaw-users at listserv.fysik.dtu.dk
Subject: [gpaw-users] Problems with LDOS calculations

Dear all,

I think I might have found some problems with the code regarding LDOS calculations. I am new to GPAW and Python, so I might have made some mistake below that lead to the problems with LDOS, but from the checks I did I think it comes from GPAW itself.

I have been trying to compute LDOS for a two-atom diamond unit cell using two of the three schemes available: Orbital-projected and Wigner-Seitz LDOS.

A minimal working example of my DFT calculation is the following (results and discussion at the end of this email):

from ase import Atoms
from gpaw import GPAW
import numpy as np
#
a=3.56972
h=0.15
diamond = Atoms('C2', positions=[(0, 0, 0), (a/4, a/4, a/4)], cell=[(0, a/2, a/2), (a/2, 0, a/2), (a/2, a/2, 0)], pbc=True)
calc = GPAW(nbands=8, xc='PBE', txt='diamond.out', kpts=(8, 8, 8), convergence={'energy': 0.00001}, h=h)
diamond.set_calculator(calc)
diamond.get_potential_energy()
calc.write('diamond.gpw', mode='all')

My DOS and LDOS calculation scripts are as follows:

from ase import Atoms
from gpaw import GPAW, restart
import numpy as np
#
diamond, calc = restart('diamond.gpw', txt=None)
#
# Get density of states
dos_energy, dos_states = calc.get_dos(spin=0, width=0.4)
dos_file = open('DOS', 'w')
ef = calc.get_fermi_level()
iterates = len(dos_energy)
for iteration in range(0, iterates):
  print >>dos_file, dos_energy[iteration] - ef, dos_states[iteration]
#
# Get local DOS (orbital)
ldos_file = open('LDOS', 'w')
for atom_number in range(0, 2):
  print >>ldos_file, '#LDOS for atom number', atom_number
  ldos_energy, ldos_states = calc.get_orbital_ldos(atom_number, spin=0, angular='sp', width=0.4)
  iterates = len(ldos_energy)
  for iteration in range(0, iterates):
    print >>ldos_file, ldos_energy[iteration] - ef, ldos_states[iteration]
  print >>ldos_file, ' '
#
# Get local DOS (Wigner-Seitz)
ldos_file = open('LDOS-WS', 'w')
for atom_number in range(0, 2):
  print >>ldos_file, '#LDOS for atom number', atom_number
  ldos_energy, ldos_states = calc.get_wigner_seitz_ldos(atom_number, spin=0, npts=201, width=0.4)
  iterates = len(ldos_energy)
  for iteration in range(0, iterates):
    print >>ldos_file, ldos_energy[iteration] - ef, ldos_states[iteration]
  print >>ldos_file, ' '

My DOS calculation works just fine. You can see a plot here: http://mcaroba.dyndns.org/misc/gpaw/DOS.png

My orbital LDOS works fine only if run in serial. If run in parallel (mpirun and gpaw-python), no matter the number of CPUs, GPAW will not keep the type of some variable and will print many instances of this message:

0 <type 'numpy.ndarray'>
1 <type 'NoneType'>

followed by this error:

IndexError: list index out of range

This can be traced back to file gpaw/utilities/dos.py and within that file, to the code block:

    for k, w in enumerate(w_k):
        eps = wfs.collect_eigenvalues(k=k, s=spin)
        print wfs.world.rank, type(eps)
        if eps is not None:
            energies[x:x + nb] = eps
        u = spin * nk + k
        P_ani = wfs.kpt_u[u].P_ani
        if a in P_ani:
            weights_xi[x:x + nb, :] = w * np.absolute(P_ani[a])**2
        x += nb

So apparently GPAW keeps the type for eps as 'numpy.ndarray' for the parent CPU (rank = 0) but does not assign a type for the other CPUs (rank = 1 in my example where I run the script with 2 CPUs). This should be very easy to fix in the code, I guess, by simply requiring that if the calculation is run in parallel only the parent CPU perform the orbital LDOS projections. In any case, when run in serial mode, the Orbital LDOS is computed correctly: http://mcaroba.dyndns.org/misc/gpaw/LDOS-orbital.png

When computing the Wigner-Seitz LDOS, the problem is that in my case it's not computed at all, either in serial or parallel, because the type of the wavefunction object is 'complex', while gpaw/utilities/dos.py requires it to be 'float':

assert wfs.dtype == float

which can be found below the definition of raw_wignerseitz_LDOS() in the same file. Just by commenting this out, GPAW will transform the WFs to float (message "ComplexWarning: Casting complex values to real discards the imaginary part") and the Wigner-Seitz LDOS will be computed: http://mcaroba.dyndns.org/misc/gpaw/LDOS-WS.png

In this case, both local contributions to the overall DOS should be the same. However they differ, I am guessing, because GPAW only uses the atoms within the unit cell of the calculation to check which atom is closest to any given point in space. In my opinion, it would be more intuitive if periodic replicas of the atoms where used in calculations with periodic boundary conditions to determine which atom is closest to any given point in space. If I move the two atoms along the diagonal of the cell, the two LDOS contributions get closer to each other. In this figure I placed the center of mass in the center of the unit cell: http://mcaroba.dyndns.org/misc/gpaw/LDOS-WS-moved.png My guess is that they are not exactly the same because of the finite discretization of space used, where the points close to the unit cell boundaries are placed preferentially on one of the two possible sides from each of the six faces of the unit cell.

Even if these issues are not critical, I think it would be helpful at least to know that they exist.

Regards,
Miguel
_______________________________________________
gpaw-users mailing list
gpaw-users at listserv.fysik.dtu.dk
https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users