symmetrizer.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 symmetrizer_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use global_oct_m
  use index_oct_m
  use lalg_adv_oct_m
  use math_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mpi_oct_m
  use par_vec_oct_m
  use profiling_oct_m
  use space_oct_m
  use symm_op_oct_m
  use symmetries_oct_m
  use types_oct_m
 
  implicit none
 
  private
  public ::                             &
    symmetrizer_t,                      &
    symmetrizer_init,                   &
    symmetrizer_end,                    &
    symmetrizer_build_map_buffer,       &
    dsymmetrizer_apply,                 &
    zsymmetrizer_apply,                 &
    dsymmetrizer_apply_single,          &
    zsymmetrizer_apply_single,          &
    dsymmetrizer_apply_batch,          &
    zsymmetrizer_apply_batch,          &
    dsymmetrize_tensor_cart,            &
    zsymmetrize_tensor_cart,            &
    dsymmetrize_magneto_optics_cart,    &
    zsymmetrize_magneto_optics_cart,    &
    symmetrize_lattice_vectors
 
  type symmetrizer_t
    private
    type(symmetries_t), pointer :: symm
    integer(int64), allocatable :: map(:,:)
    integer(int64), allocatable :: map_inv(:,:)
    type(accel_mem_t) :: map_buff
  contains
    procedure symmetrize_lattice_vectors
  end type symmetrizer_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine symmetrizer_init(this, mesh, symm)
    type(symmetrizer_t),         intent(out) :: this
    class(mesh_t),               intent(in)  :: mesh
    type(symmetries_t),  target, intent(in)  :: symm
 
    integer :: nops, ip, iop, idir, idx(3)
    real(real64) :: destpoint(3), srcpoint(3), srcpoint_inv(3), lsize(3), offset(3)
 
    push_sub(symmetrizer_init)
 
    assert(mesh%box%dim <= 3)
 
    this%symm => symm
 
    !For each operation, we create a mapping between the grid point and the symmetric point
    nops = symmetries_number(symm)
 
    safe_allocate(this%map(1:mesh%np, 1:nops))
    safe_allocate(this%map_inv(1:mesh%np, 1:nops))
 
    call profiling_in("SYMMETRIZER_INIT")
 
    lsize = real(mesh%idx%ll, real64)
    offset = real(mesh%idx%nr(1, :) + mesh%idx%enlarge, real64)
 
    do ip = 1, mesh%np
      call mesh_local_index_to_coords(mesh, ip, idx)
      destpoint = real(idx, real64)  - offset
      ! offset moves corner of cell to origin, in integer mesh coordinates
 
      assert(all(nint(destpoint) >= 0))
      assert(all(nint(destpoint) < lsize))
 
      ! move to center of cell in real coordinates
      destpoint = destpoint + offset
 
      !convert to proper reduced coordinates
      destpoint = destpoint/lsize
 
      ! iterate over all points that go to this point by a symmetry operation
      do iop = 1, nops
        srcpoint = symm_op_apply_red(symm%ops(iop), destpoint)
        srcpoint_inv = symm_op_apply_inv_red(symm%ops(iop), destpoint)
 
        !We now come back to what should be an integer, if the symmetric point beloings to the grid
        !At this point, this is already checked
        srcpoint = srcpoint*lsize
        srcpoint_inv = srcpoint_inv*lsize
 
        ! move back to reference to origin at corner of cell
        srcpoint = srcpoint - offset
        srcpoint_inv = srcpoint_inv - offset
        ! apply periodic boundary conditions in periodic directions
        do idir = 1, symm%periodic_dim
          if (nint(srcpoint(idir)) < 0 .or. nint(srcpoint(idir)) >= mesh%idx%ll(idir)) then
            srcpoint(idir) = real(modulo(nint(srcpoint(idir)), mesh%idx%ll(idir)), real64)
          end if
          if (nint(srcpoint_inv(idir)) < 0 .or. nint(srcpoint_inv(idir)) >= mesh%idx%ll(idir)) then
            srcpoint_inv(idir) = real(modulo(nint(srcpoint_inv(idir)), mesh%idx%ll(idir)), real64)
          end if
        end do
        assert(all(nint(srcpoint) >= 0))
        assert(all(nint(srcpoint) < mesh%idx%ll))
        srcpoint = srcpoint + offset
 
        assert(all(nint(srcpoint_inv) >= 0))
        assert(all(nint(srcpoint_inv) < mesh%idx%ll))
        srcpoint_inv = srcpoint_inv + offset
 
        this%map(ip, iop) = mesh_global_index_from_coords(mesh, nint(srcpoint))
        assert(this%map(ip, iop) <= mesh%np_global)
        this%map_inv(ip, iop) = mesh_global_index_from_coords(mesh, nint(srcpoint_inv))
        assert(this%map_inv(ip, iop) <= mesh%np_global)
      end do
    end do
 
    call symmetrizer_build_map_buffer(this, mesh)
 
    call profiling_out("SYMMETRIZER_INIT")
 
    pop_sub(symmetrizer_init)
  end subroutine symmetrizer_init
 
  ! ---------------------------------------------------------
 
  subroutine symmetrizer_build_map_buffer(this, mesh)
    type(symmetrizer_t), intent(inout) :: this
    class(mesh_t),       intent(in)    :: mesh
 
    integer :: nops, ip, iop
    integer, allocatable :: map_tmp(:,:)
 
    push_sub(symmetrizer_build_map_buffer)
 
    if (accel_is_enabled() .and. .not. mesh%parallel_in_domains .and. .not. this%map_buff%allocated) then
      assert(allocated(this%map))
      nops = symmetries_number(this%symm)
      safe_allocate(map_tmp(1:mesh%np, 1:nops))
      ! Shift the map to C 0-based indexing here on the host, so the device kernel in
      ! symmetrize_batch.cu can index the field buffer directly
      do iop = 1, nops
        do ip = 1, mesh%np
          map_tmp(ip, iop) = int(this%map(ip, iop) - 1, int32)
        end do
      end do
      call accel_create_buffer(this%map_buff, accel_mem_read_only, type_integer, int(mesh%np, int64)*nops)
      call accel_write_buffer(this%map_buff, mesh%np, nops, map_tmp)
      safe_deallocate_a(map_tmp)
    end if
 
    pop_sub(symmetrizer_build_map_buffer)
  end subroutine symmetrizer_build_map_buffer
 
  ! ---------------------------------------------------------
 
  subroutine symmetrizer_end(this)
    type(symmetrizer_t), intent(inout) :: this
 
    push_sub(symmetrizer_end)
    nullify(this%symm)
 
    safe_deallocate_a(this%map)
    safe_deallocate_a(this%map_inv)
 
    if (this%map_buff%allocated) call accel_free_buffer(this%map_buff)
 
    pop_sub(symmetrizer_end)
  end subroutine symmetrizer_end
 
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine symmetrize_lattice_vectors(this, size, initial_rlattice, rlattice, symmetrize)
    class(symmetrizer_t), intent(in)    :: this
    integer,              intent(in)    :: size
    real(real64),         intent(in)    :: initial_rlattice(size,size)
    real(real64),         intent(inout) :: rlattice(size,size)
    logical,              intent(in)    :: symmetrize
 
    real(real64) :: strain(size,size), inv_initial_rlattice(size,size)
    real(real64), parameter :: tol_small = 1.0e-14_real64
 
    push_sub(symmetrize_lattice_vectors)
 
    ! Remove too small elements
    call dzero_small_elements_matrix(rlattice, tol_small)
 
    inv_initial_rlattice = initial_rlattice
    call lalg_inverse(size, inv_initial_rlattice, 'dir')
 
    ! Compute strain as rlattice * initial_rlattice^{-1}
    strain = matmul(rlattice, inv_initial_rlattice)
 
    if (symmetrize) then
      ! Symmetrize the strain tensor
      call dsymmetrize_tensor_cart(this%symm, strain, use_non_symmorphic=.true.)
    else ! The tensor should be at least symmetric, as we forbid rotations
      strain = m_half * (strain + transpose(strain))
    end if
 
    ! Remove too small elements
    call dzero_small_elements_matrix(strain, tol_small)
 
    ! Get the symmetrized lattice vectors
    rlattice = matmul(strain, initial_rlattice)
 
    ! Remove too small elements
    call dzero_small_elements_matrix(rlattice, tol_small)
 
    pop_sub(symmetrize_lattice_vectors)
  end subroutine symmetrize_lattice_vectors
 
 
#include "undef.F90"
#include "real.F90"
#include "symmetrizer_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "symmetrizer_inc.F90"
 
end module symmetrizer_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: