elf.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module elf_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use global_oct_m
  use grid_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
 
  implicit none
 
  private
  public ::   &
    elf_init, &
    elf_calc
 
  logical :: with_current_term = .true.
 
contains
 
  subroutine elf_init(namespace)
    type(namespace_t),    intent(in)    :: namespace
 
    push_sub(elf_init)
 
    !%Variable ELFWithCurrentTerm
    !%Type logical
    !%Default true
    !%Section Output
    !%Description
    !% The ELF, when calculated for complex wavefunctions, should contain
    !% a term dependent on the current. This term is properly calculated by
    !% default; however, for research purposes it may be useful not to add it.
    !% If this feature proves to be useless, this option should go away.
    !%End
    call parse_variable(namespace, 'ELFWithCurrentTerm', .true., with_current_term)
 
    pop_sub(elf_init)
  end subroutine elf_init
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine elf_calc(space, st, gr, kpoints, elf, de)
    class(space_t),        intent(in)    :: space
    type(states_elec_t),   intent(inout) :: st
    type(grid_t),          intent(in)    :: gr
    type(kpoints_t),       intent(in)    :: kpoints
    real(real64),     intent(inout) :: elf(:,:)
    real(real64),  optional, intent(inout):: de(:,:)
 
    real(real64) :: factor, D0, dens
    integer :: ip, is, nelfs
    real(real64), allocatable :: rho(:,:), grho(:,:,:), jj(:,:,:)
    real(real64), allocatable :: kappa(:,:)
    real(real64), parameter :: dmin = 1e-10_real64
 
    push_sub(elf_calc)
 
    assert(space%dim == 2 .or. space%dim == 3)
 
    ! We may or may not want the total ELF.
    ! If we want it, the argument elf should have three components.
    nelfs = size(elf, 2)
 
    safe_allocate(  rho(1:gr%np, 1:st%d%nspin))
    safe_allocate(kappa(1:gr%np, 1:st%d%nspin))
    rho = m_zero
    kappa = m_zero
    call density_calc(st, gr, rho)
 
    safe_allocate(grho(1:gr%np, 1:space%dim, 1:st%d%nspin))
    safe_allocate(  jj(1:gr%np, 1:space%dim, 1:st%d%nspin))
 
    call states_elec_calc_quantities(gr, st, kpoints, .false., kinetic_energy_density = kappa, &
      paramagnetic_current = jj, density_gradient = grho)
 
    ! spin-dependent quantities
    if (st%d%ispin == unpolarized) then
      rho = rho/m_two
      kappa = kappa/m_two
      jj = jj/m_two
      grho = grho/m_two
    end if
 
    if (.not. with_current_term) jj = m_zero
 
    ! kappa will contain rho * D
    do_is: do is = 1, st%d%nspin
      do ip = 1, gr%np
        kappa(ip, is) = kappa(ip, is)*rho(ip, is)        &    ! + tau * rho
          - m_fourth*sum(grho(ip, 1:space%dim, is)**2)      &    ! - | nabla rho |^2 / 4
          - sum(jj(ip, 1:space%dim, is)**2)                      ! - j^2
      end do
 
      ! pass this information to the caller if requested
      if (present(de)) de(1:gr%np,is) = kappa(1:gr%np,is)
 
    end do do_is
 
    safe_deallocate_a(grho)
    safe_deallocate_a(jj)   ! these are no longer needed
 
    select case (space%dim)
    case (3)
      factor = m_three / m_five * (6.0_real64 * m_pi**2)**m_twothird
    case (2)
      factor = m_two * m_pi
    end select
 
    select case (st%d%ispin)
    case (unpolarized)
      do ip = 1, gr%np
        if (rho(ip, 1) >= dmin) then
          select case (space%dim)
          case (3)
            d0 = factor * rho(ip, 1)**(8.0_real64 / m_three)
          case (2)
            d0 = factor * rho(ip, 1)**3
          end select
          elf(ip, 1) = d0 * d0 / (d0 * d0 + kappa(ip, 1)**2)
        else
          elf(ip, 1) = m_zero
        end if
      end do
 
    case (spin_polarized, spinors)
      if (nelfs == 3) then
        do ip = 1, gr%np
          dens = rho(ip, 1) + rho(ip, 2)
          if (dens >= dmin) then
            select case (space%dim)
            case (3)
              d0 = factor * dens ** (m_five/m_three) * rho(ip, 1) * rho(ip, 2)
            case (2)
              d0 = factor * dens ** 2 * rho(ip, 1) * rho(ip, 2)
            end select
            elf(ip, 3) = d0 * d0 / (d0 * d0 + (kappa(ip, 1) * rho(ip, 2) + kappa(ip,2) * rho(ip, 1))**2)
          else
            elf(ip, 3) = m_zero
          end if
        end do
      end if
      do ip = 1, gr%np
        do is = 1, st%d%spin_channels
          if (rho(ip, is) >= dmin) then
            select case (space%dim)
            case (3)
              d0 = factor * rho(ip, is)**(8.0_real64/m_three)
            case (2)
              d0 = factor * rho(ip, is)**3
            end select
            elf(ip, is) = d0 * d0 / (d0 * d0 + kappa(ip,is)**2)
          else
            elf(ip, is) = m_zero
          end if
        end do
      end do
    end select
 
    safe_deallocate_a(rho)
    safe_deallocate_a(kappa)
 
    pop_sub(elf_calc)
  end subroutine elf_calc
 
end module elf_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: