lattice_vectors.F90

Go to the documentation of this file.

!! Copyright (C) 2021 N. Tancogne-Dejean, 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 lattice_vectors_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
 
  public ::                   &
    lattice_vectors_t,        &
    lattice_iterator_t
 
  type lattice_vectors_t
    ! Components are public by default
    type(space_t), private :: space
    real(real64), allocatable :: rlattice_primitive(:, :)
    real(real64), allocatable :: rlattice          (:, :)
    real(real64), allocatable :: klattice_primitive(:, :)
    real(real64), allocatable :: klattice          (:, :)
    real(real64) :: alpha, beta, gamma
    real(real64) :: rcell_volume
    logical :: nonorthogonal = .false.
    ! Some notes:
    !  rlattice_primitive is the A matrix from Chelikowski PRB 78 075109 (2008)
    !  klattice_primitive is the transpose (!) of the B matrix, with no 2 pi factor included
    !  klattice is the proper reciprocal lattice vectors, with 2 pi factor, and in units of 1/bohr
    !  The F matrix of Chelikowski is matmul(transpose(latt%klattice_primitive), latt%klattice_primitive)
  contains
    procedure :: copy => lattice_vectors_copy
    generic   :: assignment(=) => copy
    procedure :: scale => lattice_vectors_scale
    procedure :: write_info => lattice_vectors_write_info
    procedure :: short_info => lattice_vectors_short_info
    procedure :: cart_to_red => lattice_vectors_cart_to_red
    procedure :: red_to_cart => lattice_vectors_red_to_cart
    procedure :: fold_into_cell => lattice_vectors_fold_into_cell
    procedure :: max_length => lattice_vectors_max_length
    procedure :: update => lattice_vectors_update
    final :: lattice_vectors_finalize
  end type lattice_vectors_t
 
  interface lattice_vectors_t
    module procedure lattice_vectors_constructor_from_input, lattice_vectors_constructor_from_rlattice
  end interface lattice_vectors_t
 
  type lattice_iterator_t
    private
    integer, public :: n_cells = 0
    integer, allocatable :: icell(:, :)
    type(lattice_vectors_t), pointer :: latt => null()
  contains
    procedure :: copy => lattice_iterator_copy
    generic   :: assignment(=) => copy
    procedure :: get => lattice_iterator_get
    final :: lattice_iterator_finalize
  end type lattice_iterator_t
 
  interface lattice_iterator_t
    module procedure lattice_iterator_constructor
  end interface lattice_iterator_t
 
  real(real64), parameter :: eps = 100.0_real64 * m_epsilon
 
contains
 
  !--------------------------------------------------------------
  type(lattice_vectors_t) function lattice_vectors_constructor_from_rlattice(namespace, space, rlattice) result(latt)
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    real(real64),                intent(in)    :: rlattice(:, :)
 
    integer :: idir, idir2
    real(real64) :: volume_element
 
    push_sub(lattice_vectors_constructor_from_rlattice)
 
    latt%space = space
 
    safe_allocate(latt%rlattice_primitive(1:space%dim, 1:space%dim))
    safe_allocate(latt%rlattice(1:space%dim, 1:space%dim))
    safe_allocate(latt%klattice_primitive(1:space%dim, 1:space%dim))
    safe_allocate(latt%klattice(1:space%dim, 1:space%dim))
 
    latt%rlattice(1:space%dim, 1:space%dim) = rlattice(1:space%dim, 1:space%dim)
    do idir = 1, space%dim
      latt%rlattice_primitive(:, idir) = latt%rlattice(:, idir) / norm2(latt%rlattice(:, idir))
    end do
 
    latt%nonorthogonal = .false.
    do idir = 1, space%dim
      do idir2 = 1, space%dim
        if (idir /= idir2 .and. abs(latt%rlattice_primitive(idir, idir2)) > m_epsilon) then
          latt%nonorthogonal = .true.
        end if
      end do
    end do
 
    call reciprocal_lattice(latt%rlattice, latt%klattice, latt%rcell_volume, space%dim, namespace)
    latt%klattice = latt%klattice * m_two*m_pi
 
    call reciprocal_lattice(latt%rlattice_primitive, latt%klattice_primitive, volume_element, space%dim, namespace)
 
    if (space%dim == 3) then
      call angles_from_rlattice_primitive(latt%rlattice_primitive, latt%alpha, latt%beta, latt%gamma)
    else
      ! Angles should not be used, so set them to zero
      latt%alpha = m_zero
      latt%beta  = m_zero
      latt%gamma = m_zero
    end if
 
    pop_sub(lattice_vectors_constructor_from_rlattice)
  end function lattice_vectors_constructor_from_rlattice
 
  !--------------------------------------------------------------
  type(lattice_vectors_t) function lattice_vectors_constructor_from_input(namespace, space, variable_prefix) result(latt)
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    character(len=*),  optional, intent(in)    :: variable_prefix
 
    type(block_t) :: blk
    real(real64) :: norm, lparams(space%dim), volume_element
    integer :: idim, jdim, ncols
    logical :: has_angles
    real(real64) :: angles(1:space%dim)
    character(len=:), allocatable :: prefix
 
    push_sub(lattice_vectors_constructor_from_input)
 
    if (present(variable_prefix)) then
      prefix = variable_prefix
    else
      prefix = ''
    end if
 
    latt%space = space
 
    safe_allocate(latt%rlattice_primitive(1:space%dim, 1:space%dim))
    safe_allocate(latt%rlattice(1:space%dim, 1:space%dim))
    safe_allocate(latt%klattice_primitive(1:space%dim, 1:space%dim))
    safe_allocate(latt%klattice(1:space%dim, 1:space%dim))
 
    latt%alpha = 90.0_real64
    latt%beta  = 90.0_real64
    latt%gamma = 90.0_real64
 
    has_angles = .false.
    angles = 90.0_real64
 
    if (space%is_periodic()) then
 
      !%Variable LatticeParameters
      !%Type block
      !%Section Mesh::Simulation Box
      !%Description
      !% The lattice parameters (a, b, c).
      !% This variable is mandatory for periodic systems and is ignored otherwise.
      !% When PeriodicDimensions = 3, a second optional line can be used to
      !% define the angles between the lattice vectors. If the angles are not
      !% provided, then the variable LatticeVectors must be set.
      !% The number of parameters specified in the block must be at least equal
      !% to the number of periodic dimensions, but it is not mandatory to
      !% specify parameters for the non-periodic dimensions (in that case they
      !% are set to 1).
      !%End
      if (parse_block(namespace, prefix//'LatticeParameters', blk) == 0) then
        ncols = parse_block_cols(blk, 0)
        if (ncols < space%periodic_dim) then
          call messages_input_error(namespace, prefix//'LatticeParameters', &
            'The number of columns must be at least PeriodicDimensions')
        end if
        do idim = 1, ncols
          call parse_block_float(blk, 0, idim - 1, lparams(idim))
        end do
 
        ! If some parameters for non-periodic dimensions are not set in the input file, with set them to 1.
        do idim = ncols + 1, space%dim
          lparams(idim) = m_one
        end do
 
        ! Parse angles, if available
        if (parse_block_n(blk) > 1) then
          if (space%dim /= 3) then
            call messages_input_error(namespace, prefix//'LatticeParameters', 'Angles can only be specified when Dimensions = 3')
          end if
 
          ncols = parse_block_cols(blk, 1)
          if (ncols /= space%dim) then
            call messages_input_error(namespace, prefix//'LatticeParameters', 'You must specify three angles')
          end if
          do idim = 1, space%dim
            call parse_block_float(blk, 1, idim - 1, angles(idim))
          end do
          has_angles = .true.
        end if
        call parse_block_end(blk)
 
      else
        call messages_input_error(namespace, prefix//'LatticeParameters', 'Variable is mandatory for periodic systems')
      end if
 
 
      if (has_angles) then
        latt%alpha = angles(1)
        latt%beta  = angles(2)
        latt%gamma = angles(3)
 
        if (parse_is_defined(namespace, prefix//'LatticeVectors')) then
          message(1) = 'LatticeParameters with angles is incompatible with LatticeVectors'
          call messages_print_var_info(prefix//"LatticeParameters", namespace=namespace)
          call messages_fatal(1, namespace=namespace)
        end if
 
        call build_metric_from_angles(latt, angles)
 
      else
 
        !%Variable LatticeVectors
        !%Type block
        !%Default simple cubic
        !%Section Mesh::Simulation Box
        !%Description
        !% Primitive lattice vectors. Vectors are stored in rows.
        !% LatticeVectors are normalized and then used, if both LatticeParameters and
        !% LatticeVectors are defined, norm of LatticeVectors are
        !% used as multipliers of lattice parameters.
        !% Default:
        !% <br><br><tt>%LatticeVectors
        !% <br>&nbsp;&nbsp;1.0 | 0.0 | 0.0
        !% <br>&nbsp;&nbsp;0.0 | 1.0 | 0.0
        !% <br>&nbsp;&nbsp;0.0 | 0.0 | 1.0
        !% <br>%<br></tt>
        !%End
        latt%rlattice_primitive(:, :) = diagonal_matrix(space%dim, m_one)
        latt%nonorthogonal = .false.
 
        if (parse_block(namespace, prefix//'LatticeVectors', blk) == 0) then
          do idim = 1, space%dim
            do jdim = 1, space%dim
              call parse_block_float(blk, idim - 1,  jdim - 1, latt%rlattice_primitive(jdim, idim))
              if (idim /= jdim .and. abs(latt%rlattice_primitive(jdim, idim)) > m_epsilon) then
                latt%nonorthogonal = .true.
              end if
            end do
          end do
          call parse_block_end(blk)
 
        end if
      end if
 
    else
      ! Non-periodic
      lparams = m_one
      latt%rlattice_primitive(:, :) = diagonal_matrix(space%dim, m_one)
    end if
 
    latt%rlattice = m_zero
    do idim = 1, space%dim
      norm = norm2(latt%rlattice_primitive(1:space%dim, idim))
      lparams(idim) = lparams(idim)*norm
      do jdim = 1, space%dim
        latt%rlattice_primitive(jdim, idim) = latt%rlattice_primitive(jdim, idim) / norm
        latt%rlattice(jdim, idim) = latt%rlattice_primitive(jdim, idim) * lparams(idim)
      end do
    end do
 
    call reciprocal_lattice(latt%rlattice, latt%klattice, latt%rcell_volume, space%dim, namespace)
    latt%klattice = latt%klattice * m_two*m_pi
 
    call reciprocal_lattice(latt%rlattice_primitive, latt%klattice_primitive, volume_element, space%dim, namespace)
 
    if (.not. has_angles .and. space%dim == 3) then
      !We compute the angles from the lattice vectors
      call angles_from_rlattice_primitive(latt%rlattice_primitive, latt%alpha, latt%beta, latt%gamma)
    end if
 
    pop_sub(lattice_vectors_constructor_from_input)
  end function lattice_vectors_constructor_from_input
 
  !--------------------------------------------------------------
  subroutine angles_from_rlattice_primitive(rlattice_primitive, alpha, beta, gamma)
    real(real64), intent(in)  :: rlattice_primitive(1:3, 1:3)
    real(real64), intent(out) :: alpha
    real(real64), intent(out) :: beta
    real(real64), intent(out) :: gamma
 
    real(real64)  :: rlatt(1:3, 1:3)
 
    push_sub(angles_from_rlattice_primitive)
 
    rlatt = matmul(transpose(rlattice_primitive), rlattice_primitive)
    alpha = acos(rlatt(2, 3)/sqrt(rlatt(2, 2)*rlatt(3, 3)))/m_pi*180.0_real64
    beta  = acos(rlatt(1, 3)/sqrt(rlatt(1, 1)*rlatt(3, 3)))/m_pi*180.0_real64
    gamma = acos(rlatt(1, 2)/sqrt(rlatt(1, 1)*rlatt(2, 2)))/m_pi*180.0_real64
 
    pop_sub(angles_from_rlattice_primitive)
  end subroutine angles_from_rlattice_primitive
 
  !--------------------------------------------------------------
  subroutine lattice_vectors_copy(this, source)
    class(lattice_vectors_t), intent(out) :: this
    class(lattice_vectors_t), intent(in)  :: source
 
    push_sub(lattice_vectors_copy)
 
    this%space = source%space
    safe_allocate_source_a(this%rlattice_primitive, source%rlattice_primitive)
    safe_allocate_source_a(this%rlattice, source%rlattice)
    safe_allocate_source_a(this%klattice_primitive, source%klattice_primitive)
    safe_allocate_source_a(this%klattice, source%klattice)
    this%alpha = source%alpha
    this%beta = source%beta
    this%gamma = source%gamma
    this%rcell_volume = source%rcell_volume
    this%nonorthogonal = source%nonorthogonal
 
    pop_sub(lattice_vectors_copy)
  end subroutine lattice_vectors_copy
 
  !--------------------------------------------------------------
  subroutine lattice_vectors_finalize(this)
    type(lattice_vectors_t), intent(inout) :: this
 
    push_sub(lattice_vectors_finalize)
 
    safe_deallocate_a(this%rlattice_primitive)
    safe_deallocate_a(this%rlattice)
    safe_deallocate_a(this%klattice_primitive)
    safe_deallocate_a(this%klattice)
    this%nonorthogonal = .false.
 
    pop_sub(lattice_vectors_finalize)
  end subroutine lattice_vectors_finalize
 
  !--------------------------------------------------------------
  subroutine lattice_vectors_scale(this, factor)
    class(lattice_vectors_t), intent(inout) :: this
    real(real64),             intent(in)    :: factor(this%space%dim)
 
    integer :: idir
 
    push_sub(lattice_vectors_scale)
 
    ! Scale the lattice in real space
    do idir = 1, this%space%dim
      this%rlattice(1:this%space%dim, idir) = this%rlattice(1:this%space%dim, idir)*factor(idir)
    end do
 
    ! Regenerate the lattice in reciprocal space
    call reciprocal_lattice(this%rlattice, this%klattice, this%rcell_volume, this%space%dim)
    this%klattice = this%klattice * m_two*m_pi
 
    pop_sub(lattice_vectors_scale)
  end subroutine lattice_vectors_scale
 
  !--------------------------------------------------------------
  subroutine lattice_vectors_update(this, rlattice)
    class(lattice_vectors_t), intent(inout) :: this
    real(real64),             intent(in)    :: rlattice(this%space%dim, this%space%dim)
 
    integer :: idir
    real(real64) :: volume_element, maxrlatt
 
    push_sub(lattice_vectors_update)
 
    this%rlattice(:, :) = rlattice(:, :)
    maxrlatt = maxval(abs(this%rlattice))
    where(abs(this%rlattice) < 1e-12_real64*maxrlatt) this%rlattice = m_zero
 
    ! Update the primitive lattice vectors
    do idir = 1, this%space%dim
      this%rlattice_primitive(:, idir) = this%rlattice(:, idir) / norm2(this%rlattice(:, idir))
    end do
 
    ! Regenerate the lattice in reciprocal space
    call reciprocal_lattice(this%rlattice, this%klattice, this%rcell_volume, this%space%dim)
    this%klattice = this%klattice * m_two*m_pi
 
    call reciprocal_lattice(this%rlattice_primitive, this%klattice_primitive, volume_element, this%space%dim)
 
    ! Recompute the angles
    if (this%space%dim == 3) then
      call angles_from_rlattice_primitive(this%rlattice_primitive, this%alpha, this%beta, this%gamma)
    else
      ! Angles should not be used, so set them to zero
      this%alpha = m_zero
      this%beta  = m_zero
      this%gamma = m_zero
    end if
 
    pop_sub(lattice_vectors_update)
  end subroutine lattice_vectors_update
 
  !--------------------------------------------------------------
  pure function lattice_vectors_cart_to_red(this, xx_cart) result(xx_red)
    class(lattice_vectors_t), intent(in) :: this
    real(real64),             intent(in) :: xx_cart(this%space%dim)
    real(real64) :: xx_red(this%space%dim)
 
    xx_red = matmul(xx_cart, this%klattice)/(m_two*m_pi)
 
  end function lattice_vectors_cart_to_red
 
  !--------------------------------------------------------------
  pure function lattice_vectors_red_to_cart(this, xx_red) result(xx_cart)
    class(lattice_vectors_t), intent(in) :: this
    real(real64),             intent(in) :: xx_red(this%space%dim)
    real(real64) :: xx_cart(this%space%dim)
 
    xx_cart = matmul(this%rlattice, xx_red)
 
  end function lattice_vectors_red_to_cart
 
  !--------------------------------------------------------------
  function lattice_vectors_fold_into_cell(this, xx) result(new_xx)
    class(lattice_vectors_t), intent(in) :: this
    real(real64),             intent(in) :: xx(this%space%dim)
    real(real64) :: new_xx(this%space%dim)
 
    integer :: idir
 
    if (this%space%is_periodic()) then
      ! Convert the position to reduced coordinates
      new_xx = this%cart_to_red(xx)
 
      do idir = 1, this%space%periodic_dim
        ! Change of origin
        new_xx(idir) = new_xx(idir) + m_half
 
        ! Fold into cell
        new_xx(idir) = new_xx(idir) - anint(new_xx(idir))
        if (new_xx(idir) < -1.0e-6_real64) then
          new_xx(idir) = new_xx(idir) + m_one
        end if
 
        ! Sanity checks
        assert(new_xx(idir) >= -1.0e-6_real64)
        assert(new_xx(idir) < m_one)
 
        ! Change origin back
        new_xx(idir) = new_xx(idir) - m_half
      end do
 
      ! Convert back to Cartesian coordinates
      new_xx = this%red_to_cart(new_xx)
    else
      new_xx = xx
    end if
 
  end function lattice_vectors_fold_into_cell
 
  !--------------------------------------------------------------
  subroutine lattice_vectors_write_info(this, namespace)
    class(lattice_vectors_t), intent(in) :: this
    type(namespace_t),        intent(in) :: namespace
 
    integer :: idir, idir2
 
    push_sub(lattice_vectors_write_info)
 
    assert(this%space%is_periodic())
 
    call messages_print_with_emphasis(msg="Lattice", namespace=namespace)
 
    write(message(1),'(a,3a,a)') '  Lattice Vectors [', trim(units_abbrev(units_out%length)), ']'
    do idir = 1, this%space%periodic_dim
      write(message(1+idir),'(9f12.6)') (units_from_atomic(units_out%length, this%rlattice(idir2, idir)), &
        idir2 = 1, this%space%dim)
    end do
    call messages_info(1+this%space%periodic_dim, namespace=namespace)
 
    select case (this%space%periodic_dim)
    case (1)
      write(message(1),'(a)') '  Cell length ='
    case (2)
      write(message(1),'(a)') '  Cell area   ='
    case default
      write(message(1),'(a)') '  Cell volume ='
    end select
    write(message(1),'(a,1x,f18.4,3a,i1.1,a)') &
      trim(message(1)), units_from_atomic(units_out%length**this%space%periodic_dim, this%rcell_volume), &
      ' [', trim(units_abbrev(units_out%length**this%space%periodic_dim)), ']'
    call messages_info(1, namespace=namespace)
    write(message(1),'(a,3a,a)') '  Reciprocal-Lattice Vectors [', trim(units_abbrev(units_out%length**(-1))), ']'
    do idir = 1, this%space%periodic_dim
      write(message(1+idir),'(3f12.6)') (units_from_atomic(unit_one / units_out%length, this%klattice(idir2, idir)), &
        idir2 = 1, this%space%dim)
    end do
    call messages_info(1+this%space%periodic_dim, namespace=namespace)
 
    if (this%space%dim == 3) then
      write(message(1),'(a)') '  Cell angles [degree]'
      write(message(2),'(a, f8.3)') '    alpha = ', this%alpha
      write(message(3),'(a, f8.3)') '    beta  = ', this%beta
      write(message(4),'(a, f8.3)') '    gamma = ', this%gamma
      call messages_info(4, namespace=namespace)
    end if
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(lattice_vectors_write_info)
  end subroutine lattice_vectors_write_info
 
  !--------------------------------------------------------------
  character(len=140) function lattice_vectors_short_info(this, unit_length) result(info)
    class(lattice_vectors_t), intent(in) :: this
    type(unit_t),             intent(in) :: unit_length
 
    integer :: idir1, idir2
 
    push_sub(lattice_vectors_short_info)
 
    write(info, '(a,a,a)') 'LatticeVectors [', trim(units_abbrev(unit_length)), '] = '
 
    do idir1 = 1, this%space%dim
      write(info, '(a,a)') trim(info), '['
      do idir2 = 1, this%space%dim
        write(info, '(a,1x,f11.6)') trim(info), units_from_atomic(unit_length, this%rlattice(idir2, idir1))
      end do
      write(info, '(a,a)') trim(info), ']'
      if (idir1 < this%space%dim) then
        write(info, '(a,a)') trim(info), ', '
      end if
    end do
 
    pop_sub(lattice_vectors_short_info)
  end function lattice_vectors_short_info
 
  !--------------------------------------------------------------
  real(real64) function lattice_vectors_max_length(this) result(length)
    class(lattice_vectors_t), intent(in) :: this
 
    push_sub(lattice_vectors_max_length)
 
    if (this%space%is_periodic()) then
      length = maxval(norm2(this%rlattice(:,1:this%space%periodic_dim), dim=1))
    else
      length = m_zero
    end if
 
    pop_sub(lattice_vectors_max_length)
  end function lattice_vectors_max_length
 
  !--------------------------------------------------------------
  subroutine build_metric_from_angles(this, angles)
    type(lattice_vectors_t), intent(inout) :: this
    real(real64),            intent(in)    :: angles(3)
 
    real(real64) :: cosang, a2, aa, cc
    real(real64), parameter :: tol_angle = 1.0e-6_real64
 
    push_sub(build_metric_from_angles)
 
    !Converting the angles to LatticeVectors
    !See 57_iovars/ingeo.F90 in Abinit for details
    if (abs(angles(1) - angles(2)) < tol_angle .and. abs(angles(2) - angles(3)) < tol_angle .and.  &
      (abs(angles(1) - 90.0_real64) + abs(angles(2) - 90.0_real64) + abs(angles(3) - 90.0_real64)) > tol_angle) then
 
      cosang = cos(m_pi*angles(1)/180.0_real64)
      a2 = m_two/m_three*(m_one - cosang)
      aa = sqrt(a2)
      cc = sqrt(m_one - a2)
      this%rlattice_primitive(1, 1) =  aa
      this%rlattice_primitive(2, 1) =  m_zero
      this%rlattice_primitive(3, 1) =  cc
      this%rlattice_primitive(1, 2) = -m_half*aa
      this%rlattice_primitive(2, 2) =  m_half*sqrt(m_three)*aa
      this%rlattice_primitive(3, 2) =  cc
      this%rlattice_primitive(1, 3) = -m_half*aa
      this%rlattice_primitive(2, 3) = -m_half*sqrt(m_three)*aa
      this%rlattice_primitive(3, 3) =  cc
    else
      this%rlattice_primitive(1, 1) = m_one
      this%rlattice_primitive(2, 1) = m_zero
      this%rlattice_primitive(3, 1) = m_zero
      this%rlattice_primitive(1, 2) = cos(m_pi*angles(3)/180.0_real64)
      this%rlattice_primitive(2, 2) = sin(m_pi*angles(3)/180.0_real64)
      this%rlattice_primitive(3, 2) = m_zero
      this%rlattice_primitive(1, 3) = cos(m_pi*angles(2)/180.0_real64)
      this%rlattice_primitive(2, 3) = (cos(m_pi*angles(1)/180.0_real64) - &
        this%rlattice_primitive(1, 2)*this%rlattice_primitive(1, 3))/this%rlattice_primitive(2,2)
      this%rlattice_primitive(3, 3) = sqrt(m_one - this%rlattice_primitive(1, 3)**2 - this%rlattice_primitive(2, 3)**2)
    end if
 
    this%nonorthogonal = any(abs(angles - 90.0_real64) > m_epsilon)
 
    pop_sub(build_metric_from_angles)
  end subroutine build_metric_from_angles
 
  !--------------------------------------------------------------
  subroutine reciprocal_lattice(rv, kv, volume, dim, namespace)
    real(real64),      intent(in)  :: rv(:, :)
    real(real64),      intent(out) :: kv(:, :)
    real(real64),      intent(out) :: volume
    integer,           intent(in)  :: dim
    type(namespace_t), optional, intent(in)  :: namespace
 
    integer :: ii
    real(real64) :: cross(1:3)
 
    push_sub(reciprocal_lattice)
 
    select case (dim)
    case (3)
      cross(1:3) = dcross_product(rv(1:3, 2), rv(1:3, 3))
      volume = dot_product(rv(1:3, 1), cross(1:3))
 
      kv(1:3, 1) = dcross_product(rv(:, 2), rv(:, 3))/volume
      kv(1:3, 2) = dcross_product(rv(:, 3), rv(:, 1))/volume
      kv(1:3, 3) = dcross_product(rv(:, 1), rv(:, 2))/volume
    case (2)
      volume = rv(1, 1)*rv(2, 2) - rv(2, 1)*rv(1, 2)
      kv(1, 1) =  rv(2, 2)/volume
      kv(2, 1) = -rv(1, 2)/volume
      kv(1, 2) = -rv(2, 1)/volume
      kv(2, 2) =  rv(1, 1)/volume
    case (1)
      volume = rv(1, 1)
      kv(1, 1) = m_one / rv(1, 1)
    case default ! dim > 3
      message(1) = "Reciprocal lattice for dim > 3 assumes no periodicity."
      call messages_warning(1, namespace=namespace)
      volume = m_one
      kv(:, :) = m_zero
      do ii = 1, dim
        kv(ii, ii) = m_one/rv(ii,ii)
        !  At least initialize the thing
        volume = volume * norm2(rv(:, ii))
      end do
    end select
 
    if (volume < m_zero) then
      message(1) = "Your lattice vectors form a left-handed system."
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(reciprocal_lattice)
  end subroutine reciprocal_lattice
 
  ! ---------------------------------------------------------
  !! The solution is based on the replies found in
  !! https://stackoverflow.com/questions/29538135/best-way-to-find-all-points-of-lattice-in-sphere
  !! https://math.stackexchange.com/questions/1230112/maximizing-a-coordinate-of-xt-at-a-x-leq-r2/1230160#1230160
  !!
  !! When the reduced coordinates of a center are given, the iterator uses
  !! asymmetric, non-zero centered periodic copies to minimize the number of
  !! replicas, and precomputes the cells that are actually relevant for that
  !! center. In both cases, %get(ii) returns the Cartesian translation of the ii-th cell
  function lattice_iterator_constructor(latt, range, center) result(iter)
    type(lattice_vectors_t),  target, intent(in)  :: latt
    real(real64),                     intent(in)  :: range
    real(real64), optional,           intent(in)  :: center(:)
    type(lattice_iterator_t) :: iter
 
    integer :: ii, jj, idir, n_valid
    integer :: n_size(latt%space%periodic_dim)
    real(real64)   :: temp(latt%space%dim)
    real(real64), allocatable :: delta(:)
    integer, allocatable :: n_count(:), icell_tmp(:, :)
    integer :: icell(latt%space%dim), n_min, n_max
    logical :: valid
 
    push_sub(lattice_iterator_constructor)
 
    iter%latt => latt
 
    ! Determine number of cells
    iter%n_cells = 1
 
    if (present(center)) then
      ! Allow for asymmetric, non-zero centered periodic copies, to minimize the number of replicas
      safe_allocate(delta(1:latt%space%dim))
      safe_allocate(n_count(1:latt%space%dim))
 
      delta(:) = m_zero
      n_count(:) = 1
      do idir = 1, latt%space%periodic_dim
        temp = range * matmul(transpose(latt%klattice), latt%klattice(:, idir)) &
          / norm2(latt%klattice(:, idir)) / (m_two * m_pi)
        delta(idir) = temp(idir)
        n_count(idir) = floor(m_one + m_two*delta(idir) + eps) + 1
        iter%n_cells = iter%n_cells*n_count(idir)
      end do
 
      ! Precompute the cells that fall inside the exact center-dependent bounds
      safe_allocate(icell_tmp(1:latt%space%dim, 1:iter%n_cells))
      n_valid = 0
      do ii = 1, iter%n_cells
        icell(:) = 0
        jj = ii - 1
        valid = .true.
        do idir = latt%space%periodic_dim, 1, -1
          n_min = ceiling(-delta(idir) - center(idir) - eps)
          n_max = floor(m_one + delta(idir) - center(idir) + eps)
          icell(idir) = n_min + mod(jj, n_count(idir))
          if (idir > 1) jj = jj/n_count(idir)
          if (icell(idir) < n_min .or. icell(idir) > n_max) then
            valid = .false.
            exit
          end if
        end do
        if (.not. valid) cycle
        n_valid = n_valid + 1
        icell_tmp(:, n_valid) = icell(:)
      end do
 
      iter%n_cells = n_valid
      safe_allocate(iter%icell(1:latt%space%dim, 1:n_valid))
      iter%icell(:, 1:n_valid) = icell_tmp(:, 1:n_valid)
 
      safe_deallocate_a(icell_tmp)
      safe_deallocate_a(delta)
      safe_deallocate_a(n_count)
 
    else
      do idir = 1, latt%space%periodic_dim
        temp = range * matmul(transpose(latt%klattice),latt%klattice(:, idir)) / norm2(latt%klattice(:, idir)) / (m_two * m_pi)
        n_size(idir) = ceiling(temp(idir))
        iter%n_cells = iter%n_cells*(2*n_size(idir) + 1)
      end do
 
      ! Indexes of the cells
      safe_allocate(iter%icell(1:latt%space%dim, 1:iter%n_cells))
 
      do ii = 1, iter%n_cells
        jj = ii - 1
        do idir = latt%space%periodic_dim, 1, -1
          iter%icell(idir, ii) = mod(jj, 2*n_size(idir) + 1) - n_size(idir)
          if (idir > 1) jj = jj/(2*n_size(idir) + 1)
        end do
        iter%icell(latt%space%periodic_dim + 1:latt%space%dim, ii) = 0
      end do
 
    end if
 
    pop_sub(lattice_iterator_constructor)
  end function lattice_iterator_constructor
 
  !--------------------------------------------------------------
  subroutine lattice_iterator_copy(this, source)
    class(lattice_iterator_t), intent(out) :: this
    class(lattice_iterator_t), intent(in)  :: source
 
    push_sub(lattice_iterator_copy)
 
    this%n_cells = source%n_cells
    safe_allocate_source_a(this%icell, source%icell)
    this%latt => source%latt
 
    pop_sub(lattice_iterator_copy)
  end subroutine lattice_iterator_copy
 
  function lattice_iterator_get(this, ii) result(coord)
    class(lattice_iterator_t), intent(in) :: this
    integer,                   intent(in) :: ii
    real(real64) :: coord(1:this%latt%space%dim)
 
    ! No push/pop, this is called too often
 
    assert(ii <= this%n_cells)
 
    coord(:) = matmul(this%latt%rlattice(:, :), this%icell(:, ii))
 
  end function lattice_iterator_get
 
  ! ---------------------------------------------------------
  subroutine lattice_iterator_finalize(this)
    type(lattice_iterator_t), intent(inout) :: this
 
    push_sub(lattice_iterator_finalize)
 
    safe_deallocate_a(this%icell)
    this%n_cells = 0
    nullify(this%latt)
 
    pop_sub(lattice_iterator_finalize)
  end subroutine lattice_iterator_finalize
 
end module lattice_vectors_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: