utils.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2021 S. Ohlmann
!!
!! 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 utils_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use iso_c_binding
  use loct_oct_m
  use messages_oct_m
  use mpi_oct_m
  use operate_f_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use string_oct_m
 
  implicit none
 
  private
  public ::                       &
    get_divisors,                 &
    index2axis,                   &
    index2axisbz,                 &
    output_tensor,                &
    output_dipole,                &
    print_header,                 &
    leading_dimension_is_known,   &
    lead_dim,                     &
    get_config_opts,              &
    get_optional_libraries,       &
    make_array_larger
 
  interface leading_dimension_is_known
    module procedure dleading_dimension_is_known, &
      zleading_dimension_is_known,  &
      ileading_dimension_is_known,  &
      lleading_dimension_is_known,  &
      dleading_dimension_is_known2, &
      zleading_dimension_is_known2, &
      ileading_dimension_is_known2, &
      lleading_dimension_is_known2
  end interface leading_dimension_is_known
 
  interface lead_dim
    module procedure dlead_dim, zlead_dim
    module procedure dlead_dim2, zlead_dim2
  end interface lead_dim
 
contains
 
  ! ---------------------------------------------------------
  subroutine get_divisors(nn, n_divisors, divisors)
    integer, intent(in)    :: nn
    integer, intent(inout) :: n_divisors
    integer, intent(out)   :: divisors(:)
 
    integer :: ii, max_d
 
    push_sub(get_divisors)
 
    assert(n_divisors > 1)
    max_d = n_divisors
 
    n_divisors = 1
    divisors(n_divisors) = 1
    do ii = 2, nn / 2
      if (mod(nn, ii) == 0) then
        n_divisors = n_divisors + 1
 
        if (n_divisors > max_d - 1) then
          message(1) = "Internal error in get_divisors. Please increase n_divisors"
          call messages_fatal(1)
        end if
 
        divisors(n_divisors) = ii
      end if
    end do
    n_divisors = n_divisors + 1
    divisors(n_divisors) = nn
 
    pop_sub(get_divisors)
  end subroutine get_divisors
 
 
  ! ---------------------------------------------------------
  character pure function index2axis(idir) result(ch)
    integer, intent(in) :: idir
 
    select case (idir)
    case (1)
      ch = 'x'
    case (2)
      ch = 'y'
    case (3)
      ch = 'z'
    case (4)
      ch = 'w'
    case default
      write(ch,'(i1)') idir
    end select
 
  end function index2axis
 
  ! ---------------------------------------------------------
  pure function index2axisbz(idir) result(ch)
    integer, intent(in) :: idir
    character(len=2) :: ch
 
    select case (idir)
    case (1)
      ch = "kx"
    case (2)
      ch = "ky"
    case (3)
      ch = "kz"
    case (4)
      ch = "kw"
    case default
      write(ch,'(i2)') idir
    end select
 
  end function index2axisbz
 
 
  ! ---------------------------------------------------------
  subroutine output_tensor(tensor, ndim, unit, write_average, iunit, namespace)
    real(real64),                intent(in) :: tensor(:,:)
    integer,                     intent(in) :: ndim
    type(unit_t),                intent(in) :: unit
    logical,           optional, intent(in) :: write_average
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    real(real64) :: trace
    integer :: jj, kk
    logical :: write_average_
 
    push_sub(output_tensor)
 
    write_average_ = optional_default(write_average, .true.)
 
    trace = m_zero
    message(1) = ""
    do jj = 1, ndim
      do kk = 1, ndim
        write(message(1), '(a,f20.6)') trim(message(1)), units_from_atomic(unit, tensor(jj, kk))
      end do
      trace = trace + tensor(jj, jj)
      call messages_info(1, iunit=iunit, namespace=namespace)
    end do
 
    if (write_average_) then
      write(message(1), '(a, f20.6)') 'Isotropic average', units_from_atomic(unit, trace/real(ndim, real64) )
      call messages_info(1, iunit=iunit, namespace=namespace)
    end if
 
    pop_sub(output_tensor)
  end subroutine output_tensor
 
 
  ! ---------------------------------------------------------
  subroutine output_dipole(dipole, ndim, iunit, namespace)
    real(real64),                intent(in) :: dipole(:)
    integer,                     intent(in) :: ndim
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: idir
 
    push_sub(output_dipole)
 
    write(message(1), '(a,a20,a17)') 'Dipole:', '[' // trim(units_abbrev(units_out%length)) // ']', &
      '[' // trim(units_abbrev(unit_debye)) // ']'
    do idir = 1, ndim
      write(message(1+idir), '(6x,3a,es14.5,3x,2es14.5)') '<', index2axis(idir), '> = ', &
        units_from_atomic(units_out%length, dipole(idir)), units_from_atomic(unit_debye, dipole(idir))
    end do
    call messages_info(1+ndim, iunit=iunit, namespace=namespace)
 
    pop_sub(output_dipole)
  end subroutine output_dipole
 
  subroutine print_header()
    use iso_fortran_env
 
    character(len=256) :: sys_name
 
    ! Let us print our logo
    call io_dump_file(stdout, trim(trim(conf%share) // '/logo'))
 
    ! Let us print the version
    message(1) = ""
    message(2) = str_center("Running octopus", 70)
    message(3) = ""
    call messages_info(3)
 
    message(1) = &
      "Version                : " // trim(conf%version)
    message(2) = &
      "Commit                 : "// trim(conf%git_commit)
    message(3) = &
      "Configuration time     : "// trim(conf%config_time)
    call messages_info(3)
 
    message(1) = 'Configuration options  : ' // trim(get_config_opts())
    message(2) = 'Optional libraries     :'  // trim(get_optional_libraries())
 
    message(3) = 'Architecture           : ' + tostring(oct_arch)
    call messages_info(3)
 
    call write_vectorization_level()
 
    message(1) = &
      "C compiler             : "//trim(conf%cc)
    message(2) = &
      "C compiler flags       : "//trim(conf%cflags)
    message(3) = &
      "C++ compiler           : "//trim(conf%cxx)
    message(4) = &
      "C++ compiler flags     : "//trim(conf%cxxflags)
#ifdef HAVE_FC_COMPILER_VERSION
    message(5) = "Fortran compiler       : "//trim(conf%fc) //" ("//compiler_version()//")"
#else
    message(5) = "Fortran compiler       : "//trim(conf%fc)
#endif
    message(6) = &
      "Fortran compiler flags : "//trim(conf%fcflags)
    call messages_info(6)
 
    message(1) = ""
    call messages_info(1)
 
    ! Let us print where we are running
    call loct_sysname(sys_name)
    write(message(1), '(a)') str_center("The octopus is swimming in " // trim(sys_name), 70)
    message(2) = ""
    call messages_info(2)
 
    call mpi_world%barrier()
 
    call print_date("Calculation started on ")
  end subroutine print_header
 
  character(len=256) function get_config_opts()
 
    get_config_opts = ''
#ifdef HAVE_OPENMP
    get_config_opts = trim(get_config_opts)//' openmp'
#endif
#ifdef HAVE_MPI
    get_config_opts = trim(get_config_opts)//' mpi'
#endif
#ifdef HAVE_OPENCL
    get_config_opts = trim(get_config_opts)//' opencl'
#endif
#ifdef HAVE_CUDA
    get_config_opts = trim(get_config_opts)//' cuda'
#endif
#ifdef HAVE_LIBXC_FXC
    get_config_opts = trim(get_config_opts)//' libxc_fxc'
#endif
#ifdef HAVE_LIBXC_KXC
    get_config_opts = trim(get_config_opts)//' libxc_kxc'
#endif
 
  end function get_config_opts
 
  character(len=256) function get_optional_libraries()
 
    ! keep in alphabetical order, for ease in seeing if something is listed
    get_optional_libraries = ''
#ifdef HAVE_BERKELEYGW
    get_optional_libraries = trim(get_optional_libraries)//' berkeleygw'
#endif
#ifdef HAVE_CGAL
    get_optional_libraries = trim(get_optional_libraries)//' cgal'
#endif
#ifdef HAVE_CLFFT
    get_optional_libraries = trim(get_optional_libraries)//' clamdfft'
#endif
#if (defined(HAVE_CLBLAS)) || (defined(HAVE_CLBLAST))
    get_optional_libraries = trim(get_optional_libraries)//' clblas'
#endif
#ifdef HAVE_DFTBPLUS
    get_optional_libraries = trim(get_optional_libraries)//' DFTBPlus'
#endif
#ifdef HAVE_ELPA
    get_optional_libraries = trim(get_optional_libraries)//' ELPA'
#endif
#ifdef HAVE_ETSF_IO
    get_optional_libraries = trim(get_optional_libraries)//' etsf_io'
#endif
#ifdef HAVE_GDLIB
    get_optional_libraries = trim(get_optional_libraries)//' gdlib'
#endif
#ifdef HAVE_PSOLVER
    get_optional_libraries = trim(get_optional_libraries)//' psolver'
#endif
#ifdef HAVE_LIBVDWXC
    get_optional_libraries = trim(get_optional_libraries)//' libvdwxc'
#endif
#ifdef HAVE_METIS
    get_optional_libraries = trim(get_optional_libraries)//' metis'
#endif
#ifdef HAVE_NETCDF
    get_optional_libraries = trim(get_optional_libraries)//' netcdf'
#endif
#ifdef HAVE_NFFT
    get_optional_libraries = trim(get_optional_libraries)//' nfft'
#endif
#ifdef HAVE_PARMETIS
    get_optional_libraries = trim(get_optional_libraries)//' parmetis'
#endif
#ifdef HAVE_PFFT
    get_optional_libraries = trim(get_optional_libraries)//' pfft'
#endif
#ifdef HAVE_PNFFT
    get_optional_libraries = trim(get_optional_libraries)//' pnfft'
#endif
#ifdef HAVE_PSPIO
    get_optional_libraries = trim(get_optional_libraries)//' pspio'
#endif
#ifdef HAVE_SCALAPACK
    get_optional_libraries = trim(get_optional_libraries)//' scalapack'
#endif
#ifdef HAVE_SPARSKIT
    get_optional_libraries = trim(get_optional_libraries)//' sparskit'
#endif
#ifdef HAVE_NLOPT
    get_optional_libraries = trim(get_optional_libraries)//' nlopt'
#endif
 
  end function get_optional_libraries
 
  ! ---------------------------------------------------------
 
  logical function dleading_dimension_is_known(array) result(known)
    real(real64), target, intent(in) :: array(:, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2)), 0_c_intptr_t)
      known = ubound(array, dim = 1) == (addr2-addr1)/c_sizeof(array(1, 1))
    end if
 
  end function dleading_dimension_is_known
 
 
  ! ---------------------------------------------------------
 
  logical function zleading_dimension_is_known(array) result(known)
    complex(real64), target, intent(in) :: array(:, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2)), 0_c_intptr_t)
      known = ubound(array, dim = 1) == (addr2-addr1)/c_sizeof(array(1, 1))
    end if
 
  end function zleading_dimension_is_known
 
  ! ---------------------------------------------------------
 
  logical function ileading_dimension_is_known(array) result(known)
    integer, target, intent(in) :: array(:, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2)), 0_c_intptr_t)
 
      known = ubound(array, dim = 1) == (addr2-addr1)/c_sizeof(array(1, 1))
    end if
 
  end function ileading_dimension_is_known
 
 
  logical function lleading_dimension_is_known(array) result(known)
    integer(int64), target, intent(in) :: array(:, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2)), 0_c_intptr_t)
 
      known = ubound(array, dim = 1) == (addr2-addr1)/c_sizeof(array(1, 1))
    end if
 
  end function lleading_dimension_is_known
 
 
  ! ---------------------------------------------------------
 
  logical function dleading_dimension_is_known2(array) result(known)
    real(real64), target, intent(in) :: array(:, :, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2,1)), 0_c_intptr_t)
      known = ubound(array, dim = 1) == (addr2 - addr1)/c_sizeof(array(1, 1, 1))
    end if
 
  end function dleading_dimension_is_known2
 
 
  ! ---------------------------------------------------------
 
  logical function zleading_dimension_is_known2(array) result(known)
    complex(real64), target, intent(in) :: array(:, :, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2,1)), 0_c_intptr_t)
      known = ubound(array, dim = 1) == (addr2 - addr1)/c_sizeof(array(1, 1, 1))
    end if
 
  end function zleading_dimension_is_known2
 
  ! ---------------------------------------------------------
 
  logical function ileading_dimension_is_known2(array) result(known)
    integer, target, intent(in) :: array(:, :, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2,1)), 0_c_intptr_t)
      known = ubound(array, dim = 1) == (addr2 - addr1)/c_sizeof(array(1, 1, 1))
    end if
 
  end function ileading_dimension_is_known2
 
 
  logical function lleading_dimension_is_known2(array) result(known)
    integer(int64), target, intent(in) :: array(:, :, :)
 
    integer(c_intptr_t) :: addr1, addr2
 
    known = .true.
 
    if (ubound(array, dim = 2) > 1) then
      addr1 = transfer(c_loc(array(1,1,1)), 0_c_intptr_t)
      addr2 = transfer(c_loc(array(1,2,1)), 0_c_intptr_t)
      known = ubound(array, dim = 1) == (addr2 - addr1)/c_sizeof(array(1, 1, 1))
    end if
 
  end function lleading_dimension_is_known2
 
  ! ---------------------------------------------------------
 
  integer function dlead_dim(array) result(lead_dim)
    real(real64), intent(in) :: array(:, :)
 
    assert(leading_dimension_is_known(array))
 
    lead_dim = ubound(array, dim = 1)
  end function dlead_dim
 
  ! ---------------------------------------------------------
 
  integer function zlead_dim(array) result(lead_dim)
    complex(real64), intent(in) :: array(:, :)
 
    assert(leading_dimension_is_known(array))
 
    lead_dim = ubound(array, dim = 1)
  end function zlead_dim
 
  ! ---------------------------------------------------------
 
  integer function dlead_dim2(array) result(lead_dim)
    real(real64), intent(in) :: array(:, :, :)
 
    assert(leading_dimension_is_known(array))
 
    lead_dim = ubound(array, dim = 1) * ubound(array, dim = 2)
  end function dlead_dim2
 
  ! ---------------------------------------------------------
 
  integer function zlead_dim2(array) result(lead_dim)
    complex(real64), intent(in) :: array(:, :, :)
 
    assert(leading_dimension_is_known(array))
 
    lead_dim = ubound(array, dim = 1) * ubound(array, dim = 2)
  end function zlead_dim2
 
  subroutine make_array_larger(array, new_size)
    integer(int64), allocatable, intent(inout) :: array(:)
    integer,                  intent(in)    :: new_size
 
    integer(int64), allocatable :: tmp(:)
    integer :: copy_size
 
    push_sub(make_array_larger)
 
    ! Use allocate here as move_alloc deallocate internally the from array
    allocate(tmp(1:new_size))
    copy_size = min(new_size, size(array))
    tmp(1:copy_size) = array(1:copy_size)
    safe_deallocate_a(array)
    call move_alloc(tmp, array)
 
    pop_sub(make_array_larger)
  end subroutine make_array_larger
 
  subroutine write_vectorization_level()
    character(len=32) :: vec
    character(kind=c_char) :: c_str(33)
 
    c_str = c_null_char
    call get_vectorization_level(c_str)
    call string_c_to_f(c_str, vec)
    message(1) = 'Vectorization level    : ' // trim(vec)
    call messages_info(1)
 
  end subroutine write_vectorization_level
 
end module utils_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: