[gpaw-users] H_ks in realspace ?

[Ask Hjorth Larsen]

Hi

2012/7/5 Torsten Hahn <Torsten.Hahn at physik.tu-freiberg.de>:
> Hi,
>
> thx for that comment.
>
> However that is only part of what i need. Let me formulate the question a bit different.
>
> If i really want to calculate something like
>
> (1) x = < psi1 | H | psi2 >
>
> than i need psi1/psi2 and H (or some other operator) to be defined in the same basis.
>
> With psi1 = calc.get_pseudo_wavefunction() i have psi1 in a realspace representation. To calculate (1) i would therefore need H also in realspace representation.
>
> Or the other way around. If i take H in LCAO-Basis (which is the case i think if i use get_lcao_hamiltonian() ) i need psi also converted (back?) to LCAO basis.
>
> Dont get me wrong, i just do not want to reinvent the wheel. Im quite sure the required code is somewhere already present (and error checked) in GPAW. Would save me plenty of time if somebody could give me a hint.


Right, doing it on the real-space grid will be less efficient and less
precise.  You can obtain the potential matrix by calling
calc.wfs.basis_functions.calculate_potential_matrices(), which will
give you a (number of k-points X number of basis functions X number of
basis functions)-shaped array V.  Depending on your application, I am
guessing that you will need to add the nonlocal contribution from the
pseudopotential as well.  Maybe it's easier to just calculate the
whole Hamiltonian and subtract those parts that we do *not* want.
Let's do that.  In fact, please forget all of the above.

Call

  H = calc.wfs.eigensolver.calculate_hamiltonian_matrix(calc.hamiltonian,
calc.wfs, kpt)

for each kpt in calc.wfs.kpt_u.

Then you get the full Hamiltonian from which you subtract the kinetic matrix:

  H -= calc.wfs.T_qMM[kpt.q]

For each k-point you can then calculate C times H times C.T.conj(),
where C are the wavefunction coefficients calc.wfs.kpt_u[u].C_nM with
u being the index of that k-point.

This gives you the matrix elements of the Kohn-Sham potential,
including non-local atomic contributions, with respect to the
Kohn-Sham eigenstates.  You should be able to verify that C
diagonalizes H before the kinetic T_qMM is subtracted.

If you parallelize over states/orbitals, and/or you use scalapack,
this will require some parallel code as well to unravel the arrays.
If you are restarting from a file, you may need to call
calc.set_positions(atoms) to cause certain arrays to be deployed onto
grids.

Regards
Ask