voronoi_interface.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Interface for Voronoi Integration and output of BQB files
!> \par History
!>      11/2020 created [mbrehm]
!> \author Martin Brehm
! **************************************************************************************************
MODULE voronoi_interface
   USE input_section_types, ONLY: section_vals_type, &
                                  section_vals_val_get
   USE cp_log_handling, ONLY: cp_get_default_logger, &
                              cp_logger_get_default_io_unit, &
                              cp_logger_type
   USE bibliography, ONLY: rycroft2009, thomas2015, brehm2018, brehm2020, &
                           brehm2021, cite_reference
   USE kinds, ONLY: dp, default_path_length
   USE cell_types, ONLY: cell_type, pbc
   USE pw_types, ONLY: pw_r3d_rs_type
   USE physcon, ONLY: bohr, debye
   USE mathconstants, ONLY: fourpi
   USE orbital_pointers, ONLY: indco
   USE qs_environment_types, ONLY: get_qs_env, &
                                   qs_environment_type
   USE molecule_kind_types, ONLY: molecule_kind_type, &
                                  write_molecule_kind_set
   USE molecule_types, ONLY: molecule_type
   USE qs_rho_types, ONLY: qs_rho_get, &
                           qs_rho_type
   USE atomic_kind_types, ONLY: atomic_kind_type, &
                                get_atomic_kind
   USE particle_list_types, ONLY: particle_list_type
   USE particle_types, ONLY: particle_type
   USE cp_files, ONLY: file_exists, close_file, open_file
   USE qs_kind_types, ONLY: get_qs_kind, &
                            qs_kind_type
   USE message_passing, ONLY: mp_para_env_type
   USE qs_subsys_types, ONLY: qs_subsys_get, &
                              qs_subsys_type
   USE cp_control_types, ONLY: dft_control_type
   USE pw_grid_types, ONLY: pw_mode_local
   USE physcon, ONLY: angstrom, femtoseconds
   USE message_passing, ONLY: mp_comm_type
   USE input_constants, ONLY: &
      voro_radii_unity, voro_radii_vdw, voro_radii_cov, voro_radii_user, &
      bqb_opt_off, bqb_opt_quick, bqb_opt_normal, bqb_opt_patient, bqb_opt_exhaustive, do_method_gapw
   USE qs_rho0_types, ONLY: rho0_atom_type, &
                            rho0_mpole_type
   USE iso_c_binding, ONLY: c_int, c_double, c_char
 
#if defined(__HAS_IEEE_EXCEPTIONS)
   USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
                              ieee_set_halting_mode, &
                              ieee_all
#endif
 
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   ! Global parameters
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'voronoi_interface'
   PUBLIC :: entry_voronoi_or_bqb, finalize_libvori
   INTEGER :: step_count = 0
 
#if defined(__LIBVORI)
 
   ! The C interface to libvori
 
   INTERFACE
 
      INTEGER(C_INT) FUNCTION libvori_setbqbskipfirst(i) bind(C, NAME='libvori_setBQBSkipFirst')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setbqbskipfirst
 
      INTEGER(C_INT) FUNCTION libvori_setbqbstorestep(i) bind(C, NAME='libvori_setBQBStoreStep')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setbqbstorestep
 
      INTEGER(C_INT) FUNCTION libvori_setvoronoiskipfirst(i) bind(C, NAME='libvori_setVoronoiSkipFirst')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setvoronoiskipfirst
 
      INTEGER(C_INT) FUNCTION libvori_setbqbcheck(i) bind(C, NAME='libvori_setBQBCheck')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setbqbcheck
 
      INTEGER(C_INT) FUNCTION libvori_setbqbfilename(len, s) bind(C, NAME='libvori_setBQBFilename')
         USE iso_c_binding, ONLY: c_int, c_char
      INTEGER(C_INT), VALUE                              :: len
      CHARACTER(C_CHAR)                                  :: s(*)
 
      END FUNCTION libvori_setbqbfilename
 
      INTEGER(C_INT) FUNCTION libvori_setbqbparmstring(len, s) bind(C, NAME='libvori_setBQBParmString')
         USE iso_c_binding, ONLY: c_int, c_char
      INTEGER(C_INT), VALUE                              :: len
      CHARACTER(C_CHAR)                                  :: s(*)
 
      END FUNCTION libvori_setbqbparmstring
 
      INTEGER(C_INT) FUNCTION libvori_setbqbhistory(i) bind(C, NAME='libvori_setBQBHistory')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setbqbhistory
 
      INTEGER(C_INT) FUNCTION libvori_setbqboptimization(i) bind(C, NAME='libvori_setBQBOptimization')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setbqboptimization
 
      INTEGER(C_INT) FUNCTION libvori_processbqbframe(step, t) bind(C, NAME='libvori_processBQBFrame')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: step
      REAL(C_DOUBLE), VALUE                              :: t
 
      END FUNCTION libvori_processbqbframe
 
      INTEGER(C_INT) FUNCTION libvori_setprefix_voronoi() bind(C, NAME='libvori_setPrefix_Voronoi')
         USE iso_c_binding, ONLY: c_int
      END FUNCTION libvori_setprefix_voronoi
 
      INTEGER(C_INT) FUNCTION libvori_setprefix_bqb() bind(C, NAME='libvori_setPrefix_BQB')
         USE iso_c_binding, ONLY: c_int
      END FUNCTION libvori_setprefix_bqb
 
      INTEGER(C_INT) FUNCTION libvori_setrefinementfactor(i) bind(C, NAME='libvori_setRefinementFactor')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setrefinementfactor
 
      INTEGER(C_INT) FUNCTION libvori_setjitter(i) bind(C, NAME='libvori_setJitter')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setjitter
 
      INTEGER(C_INT) FUNCTION libvori_setjitteramplitude(amp) bind(C, NAME='libvori_setJitterAmplitude')
         USE iso_c_binding, ONLY: c_int, c_double
      REAL(C_DOUBLE), VALUE                              :: amp
 
      END FUNCTION libvori_setjitteramplitude
 
      INTEGER(C_INT) FUNCTION libvori_setjitterseed(seed) bind(C, NAME='libvori_setJitterSeed')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: seed
 
      END FUNCTION libvori_setjitterseed
 
      INTEGER(C_INT) FUNCTION libvori_setempoutput(i) bind(C, NAME='libvori_setEMPOutput')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setempoutput
 
      INTEGER(C_INT) FUNCTION libvori_setprintlevel_verbose() bind(C, NAME='libvori_setPrintLevel_Verbose')
         USE iso_c_binding, ONLY: c_int
      END FUNCTION libvori_setprintlevel_verbose
 
      INTEGER(C_INT) FUNCTION libvori_setradii_unity() bind(C, NAME='libvori_setRadii_Unity')
         USE iso_c_binding, ONLY: c_int
      END FUNCTION libvori_setradii_unity
 
      INTEGER(C_INT) FUNCTION libvori_setradii_covalent() bind(C, NAME='libvori_setRadii_Covalent')
         USE iso_c_binding, ONLY: c_int
      END FUNCTION libvori_setradii_covalent
 
      INTEGER(C_INT) FUNCTION libvori_setradii_user(factor, rad) bind(C, NAME='libvori_setRadii_User')
         USE iso_c_binding, ONLY: c_int, c_double
      REAL(C_DOUBLE), VALUE                              :: factor
      REAL(C_DOUBLE)                                     :: rad(*)
 
      END FUNCTION libvori_setradii_user
 
      INTEGER(C_INT) FUNCTION libvori_step(step, t) bind(C, NAME='libvori_step')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: step
      REAL(C_DOUBLE), VALUE                              :: t
 
      END FUNCTION libvori_step
 
      INTEGER(C_INT) FUNCTION libvori_sanitycheck(step, t) bind(C, NAME='libvori_sanitycheck')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: step
      REAL(C_DOUBLE), VALUE                              :: t
 
      END FUNCTION libvori_sanitycheck
 
      INTEGER(C_INT) FUNCTION libvori_setgrid(rx, ry, rz, ax, ay, az, bx, by, bz, cx, cy, cz, tax, tay, taz, tbx, tby, tbz, &
                                              tcx, tcy, tcz) bind(C, NAME='libvori_setGrid')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: rx, ry, rz
      REAL(C_DOUBLE), VALUE                              :: ax, ay, az, bx, by, bz, cx, cy, cz, tax, &
                                                            tay, taz, tbx, tby, tbz, tcx, tcy, tcz
 
      END FUNCTION libvori_setgrid
 
      INTEGER(C_INT) FUNCTION libvori_pushatoms(n, pord, pchg, posx, posy, posz) bind(C, NAME='libvori_pushAtoms')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: n
      INTEGER(C_INT)                                     :: pord(*)
      REAL(C_DOUBLE)                                     :: pchg(*), posx(*), posy(*), posz(*)
 
      END FUNCTION libvori_pushatoms
 
      INTEGER(C_INT) FUNCTION libvori_push_rho_zrow(ix, iy, length, buf) bind(C, NAME='libvori_push_rho_zrow')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: ix, iy, length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_push_rho_zrow
 
      INTEGER(C_INT) FUNCTION libvori_setbqboverwrite(i) bind(C, NAME='libvori_setBQBOverwrite')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setbqboverwrite
 
      INTEGER(C_INT) FUNCTION libvori_setvorioverwrite(i) bind(C, NAME='libvori_setVoriOverwrite')
         USE iso_c_binding, ONLY: c_int
      INTEGER(C_INT), VALUE                              :: i
 
      END FUNCTION libvori_setvorioverwrite
 
      INTEGER(C_INT) FUNCTION libvori_get_radius(length, buf) bind(C, NAME='libvori_get_radius')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_radius
 
      INTEGER(C_INT) FUNCTION libvori_get_volume(length, buf) bind(C, NAME='libvori_get_volume')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_volume
 
      INTEGER(C_INT) FUNCTION libvori_get_charge(length, buf) bind(C, NAME='libvori_get_charge')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_charge
 
      INTEGER(C_INT) FUNCTION libvori_get_dipole(component, length, buf) bind(C, NAME='libvori_get_dipole')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: component, length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_dipole
 
      INTEGER(C_INT) FUNCTION libvori_get_quadrupole(component, length, buf) bind(C, NAME='libvori_get_quadrupole')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: component, length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_quadrupole
 
      INTEGER(C_INT) FUNCTION libvori_get_wrapped_pos(component, length, buf) bind(C, NAME='libvori_get_wrapped_pos')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: component, length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_wrapped_pos
 
      INTEGER(C_INT) FUNCTION libvori_get_charge_center(component, length, buf) bind(C, NAME='libvori_get_charge_center')
         USE iso_c_binding, ONLY: c_int, c_double
      INTEGER(C_INT), VALUE                              :: component, length
      REAL(C_DOUBLE)                                     :: buf(*)
 
      END FUNCTION libvori_get_charge_center
 
      INTEGER(C_INT) FUNCTION libvori_finalize() bind(C, NAME='libvori_finalize')
         USE iso_c_binding, ONLY: c_int
      END FUNCTION libvori_finalize
 
   END INTERFACE
 
#endif
 
! **************************************************************************************************
 
CONTAINS
 
! **************************************************************************************************
!> \brief Does a Voronoi integration of density or stores the density to compressed BQB format
!> \param do_voro whether a Voronoi integration shall be performed
!> \param do_bqb whether the density shall be written to compressed BQB format
!> \param input_voro the input section for Voronoi integration
!> \param input_bqb the input section for the BQB compression
!> \param unit_voro the output unit number for the Voronoi integration results
!> \param qs_env the qs_env where to calculate the charges
!> \param rspace_pw the grid with the real-space electron density to integrate/compress
!> \author Martin Brehm
! **************************************************************************************************
   SUBROUTINE entry_voronoi_or_bqb(do_voro, do_bqb, input_voro, input_bqb, unit_voro, qs_env, rspace_pw)
      INTEGER                                            :: do_voro, do_bqb
      TYPE(section_vals_type), POINTER                   :: input_voro, input_bqb
      INTEGER, INTENT(IN)                                :: unit_voro
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(pw_r3d_rs_type)                             :: rspace_pw
 
#if defined(__LIBVORI)
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'entry_voronoi_or_bqb'
      INTEGER                                            :: handle, iounit, &
                                                            ret, i, tag, &
                                                            nkind, natom, ikind, iat, ord, source, dest, &
                                                            ip, i1, i2, reffac, radius_type, bqb_optimize, &
                                                            bqb_history, nspins, jitter_seed
      LOGICAL                                            :: outemp, bqb_skip_first, voro_skip_first, &
                                                            bqb_store_step, bqb_check, voro_sanity, &
                                                            bqb_overwrite, vori_overwrite, molprop, &
                                                            gapw, jitter
      REAL(kind=dp)                                      :: zeff, qa, fn0, fn1, jitter_amplitude
      TYPE(qs_rho_type), POINTER                         :: rho
      TYPE(cp_logger_type), POINTER                      :: logger
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: buf
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: particles_z
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: particles_r
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: particles_c
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: particles_radius
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_list_type), POINTER                  :: particles
      REAL(kind=dp)                                      :: r
      TYPE(qs_subsys_type), POINTER                      :: subsys
      INTEGER, DIMENSION(:), POINTER                     :: atom_list
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(cell_type), POINTER                           :: cell
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: voro_radii
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: voro_charge
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: voro_volume
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: voro_dipole
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: voro_quadrupole
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: voro_buffer
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: voro_wrapped_pos
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: voro_charge_center
      REAL(kind=dp), DIMENSION(:), POINTER               :: user_radii
      CHARACTER(len=default_path_length)                 :: bqb_file_name, mp_file_name
      CHARACTER(len=128)                                 :: bqb_parm_string
      TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
#if defined(__HAS_IEEE_EXCEPTIONS)
      LOGICAL, DIMENSION(5)                              :: halt
#endif
 
      CALL timeset(routinen, handle)
      NULLIFY (logger)
      logger => cp_get_default_logger()
      iounit = cp_logger_get_default_io_unit(logger)
 
      CALL get_qs_env(qs_env=qs_env, rho=rho, qs_kind_set=qs_kind_set, &
                      atomic_kind_set=atomic_kind_set, para_env=para_env, &
                      nkind=nkind, natom=natom, subsys=subsys, dft_control=dft_control, &
                      cell=cell)
 
      tag = 1
 
      IF (do_voro /= 0) THEN
         CALL section_vals_val_get(input_voro, "REFINEMENT_FACTOR", i_val=reffac)
         CALL section_vals_val_get(input_voro, "OUTPUT_EMP", l_val=outemp)
         CALL section_vals_val_get(input_voro, "JITTER", l_val=jitter)
         CALL section_vals_val_get(input_voro, "JITTER_AMPLITUDE", r_val=jitter_amplitude)
         CALL section_vals_val_get(input_voro, "JITTER_SEED", i_val=jitter_seed)
         CALL section_vals_val_get(input_voro, "VORONOI_RADII", i_val=radius_type)
         CALL section_vals_val_get(input_voro, "SKIP_FIRST", l_val=voro_skip_first)
         CALL section_vals_val_get(input_voro, "SANITY_CHECK", l_val=voro_sanity)
         CALL section_vals_val_get(input_voro, "OVERWRITE", l_val=vori_overwrite)
         IF (radius_type == voro_radii_user) THEN
            CALL section_vals_val_get(input_voro, "USER_RADII", r_vals=user_radii)
         END IF
         IF (qs_env%single_point_run) THEN
            voro_skip_first = .false.
         END IF
         CALL cite_reference(rycroft2009)
         CALL cite_reference(thomas2015)
         CALL cite_reference(brehm2018)
         CALL cite_reference(brehm2020)
         CALL cite_reference(brehm2021)
         !
         CALL section_vals_val_get(input_voro, "MOLECULAR_PROPERTIES", l_val=molprop)
         CALL section_vals_val_get(input_voro, "MOLPROP_FILE_NAME", c_val=mp_file_name)
      END IF
 
      IF (do_bqb /= 0) THEN
         CALL section_vals_val_get(input_bqb, "HISTORY", i_val=bqb_history)
         CALL section_vals_val_get(input_bqb, "OPTIMIZE", i_val=bqb_optimize)
         CALL section_vals_val_get(input_bqb, "FILENAME", c_val=bqb_file_name)
         CALL section_vals_val_get(input_bqb, "SKIP_FIRST", l_val=bqb_skip_first)
         CALL section_vals_val_get(input_bqb, "STORE_STEP_NUMBER", l_val=bqb_store_step)
         CALL section_vals_val_get(input_bqb, "CHECK", l_val=bqb_check)
         CALL section_vals_val_get(input_bqb, "OVERWRITE", l_val=bqb_overwrite)
         CALL section_vals_val_get(input_bqb, "PARAMETER_KEY", c_val=bqb_parm_string)
         IF (qs_env%single_point_run) THEN
            bqb_skip_first = .false.
            bqb_history = 1
         END IF
         IF (bqb_history < 1) THEN
            bqb_history = 1
         END IF
         CALL cite_reference(brehm2018)
      END IF
 
      CALL qs_subsys_get(subsys, particles=particles)
 
      ! Temporarily disable floating point traps because libvori raise IEEE754 exceptions.
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_get_halting_mode(ieee_all, halt)
      CALL ieee_set_halting_mode(ieee_all, .false.)
#endif
 
      associate(ionode => para_env%is_source(), my_rank => para_env%mepos, &
                 num_pe => para_env%num_pe, &
                 sim_step => qs_env%sim_step, sim_time => qs_env%sim_time, &
                 l1 => rspace_pw%pw_grid%bounds(1, 1), l2 => rspace_pw%pw_grid%bounds(1, 2), &
                 l3 => rspace_pw%pw_grid%bounds(1, 3), u1 => rspace_pw%pw_grid%bounds(2, 1), &
                 u2 => rspace_pw%pw_grid%bounds(2, 2), u3 => rspace_pw%pw_grid%bounds(2, 3))
      IF (ionode) THEN
 
         IF (iounit > 0) THEN
            WRITE (iounit, *) ""
         END IF
 
         IF (do_voro /= 0) THEN
            ret = libvori_setprefix_voronoi()
            ret = libvori_setrefinementfactor(reffac)
            ret = libvori_setjitter(merge(1, 0, jitter))
            ret = libvori_setjitteramplitude(jitter_amplitude*angstrom)
            ret = libvori_setjitterseed(jitter_seed)
            ret = libvori_setvoronoiskipfirst(merge(1, 0, voro_skip_first))
            ret = libvori_setvorioverwrite(merge(1, 0, vori_overwrite))
            ret = libvori_setempoutput(merge(1, 0, outemp))
         ELSE
            ret = libvori_setprefix_bqb()
         END IF
 
         IF (do_bqb /= 0) THEN
            ret = libvori_setbqbfilename(default_path_length, bqb_file_name)
            ret = libvori_setbqbparmstring(128, bqb_parm_string)
            SELECT CASE (bqb_optimize)
            CASE (bqb_opt_off)
               bqb_optimize = 0
            CASE (bqb_opt_quick)
               bqb_optimize = 1
            CASE (bqb_opt_normal)
               bqb_optimize = 2
            CASE (bqb_opt_patient)
               bqb_optimize = 3
            CASE (bqb_opt_exhaustive)
               bqb_optimize = 4
            END SELECT
            ret = libvori_setbqboptimization(bqb_optimize)
            ret = libvori_setbqbhistory(bqb_history)
            ret = libvori_setbqbskipfirst(merge(1, 0, bqb_skip_first))
            ret = libvori_setbqbcheck(merge(1, 0, bqb_check))
            ret = libvori_setbqboverwrite(merge(1, 0, bqb_overwrite))
            ret = libvori_setbqbstorestep(merge(1, 0, bqb_store_step))
         END IF
 
         ret = libvori_setgrid( &
               u1 - l1 + 1, &
               u2 - l2 + 1, &
               u3 - l3 + 1, &
               rspace_pw%pw_grid%dh(1, 1)*(u1 - l1 + 1), &
               rspace_pw%pw_grid%dh(2, 1)*(u1 - l1 + 1), &
               rspace_pw%pw_grid%dh(3, 1)*(u1 - l1 + 1), &
               rspace_pw%pw_grid%dh(1, 2)*(u2 - l2 + 1), &
               rspace_pw%pw_grid%dh(2, 2)*(u2 - l2 + 1), &
               rspace_pw%pw_grid%dh(3, 2)*(u2 - l2 + 1), &
               rspace_pw%pw_grid%dh(1, 3)*(u3 - l3 + 1), &
               rspace_pw%pw_grid%dh(2, 3)*(u3 - l3 + 1), &
               rspace_pw%pw_grid%dh(3, 3)*(u3 - l3 + 1), &
               cell%hmat(1, 1), &
               cell%hmat(2, 1), &
               cell%hmat(3, 1), &
               cell%hmat(1, 2), &
               cell%hmat(2, 2), &
               cell%hmat(3, 2), &
               cell%hmat(1, 3), &
               cell%hmat(2, 3), &
               cell%hmat(3, 3) &
               )
 
         IF (ret /= 0) THEN
            cpabort("The library returned an error. Aborting.")
         END IF
 
         ALLOCATE (particles_z(natom))
         ALLOCATE (particles_c(natom))
         ALLOCATE (particles_r(3, natom))
         ALLOCATE (particles_radius(natom))
 
         DO ikind = 1, nkind
            CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, vdw_radius=r)
            r = r*angstrom
            atomic_kind => atomic_kind_set(ikind)
            CALL get_atomic_kind(atomic_kind, atom_list=atom_list, z=ord)
            DO iat = 1, SIZE(atom_list)
               i = atom_list(iat)
               particles_c(i) = zeff
               particles_z(i) = ord
               particles_r(:, i) = particles%els(i)%r(:)
               particles_radius(i) = r
            END DO
         END DO
 
         ret = libvori_pushatoms(natom, particles_z, particles_c, particles_r(1, :), particles_r(2, :), particles_r(3, :))
 
         IF (ret /= 0) THEN
            cpabort("The library returned an error. Aborting.")
         END IF
 
      END IF
 
      IF (iounit > 0) THEN
         IF (do_voro /= 0) THEN
            WRITE (iounit, *) "VORONOI| Collecting electron density from MPI ranks and sending to library..."
         ELSE
            WRITE (iounit, *) "BQB| Collecting electron density from MPI ranks and sending to library..."
         END IF
      END IF
 
      ALLOCATE (buf(l3:u3))
 
      dest = 0
 
      DO i1 = l1, u1
         DO i2 = l2, u2
 
            ! cycling through the CPUs, check if the current ray (I1,I2) is local to that CPU
            IF (rspace_pw%pw_grid%para%mode .NE. pw_mode_local) THEN
               DO ip = 0, num_pe - 1
                  IF (rspace_pw%pw_grid%para%bo(1, 1, ip, 1) <= i1 - l1 + 1 &
                      .AND. rspace_pw%pw_grid%para%bo(2, 1, ip, 1) >= i1 - l1 + 1 .AND. &
                      rspace_pw%pw_grid%para%bo(1, 2, ip, 1) <= i2 - l2 + 1 &
                      .AND. rspace_pw%pw_grid%para%bo(2, 2, ip, 1) >= i2 - l2 + 1) THEN
                     source = ip
                  END IF
               END DO
            ELSE
               source = dest
            END IF
 
            IF (source == dest) THEN
               IF (my_rank == source) THEN
                  buf(:) = rspace_pw%array(i1, i2, :)
               END IF
            ELSE
               IF (my_rank == source) THEN
                  buf(:) = rspace_pw%array(i1, i2, :)
                  CALL para_env%send(buf, dest, tag)
               END IF
               IF (my_rank == dest) THEN
                  CALL para_env%recv(buf, source, tag)
               END IF
            END IF
 
            IF (my_rank == dest) THEN
               ret = libvori_push_rho_zrow(i1 - l1, i2 - l2, u3 - l3 + 1, buf)
               IF (ret /= 0) THEN
                  cpabort("The library returned an error. Aborting.")
               END IF
            END IF
 
            ! this double loop generates so many messages that it can overload
            ! the message passing system, e.g. on XT3
            ! we therefore put a barrier here that limits the amount of message
            ! that flies around at any given time.
            ! if ever this routine becomes a bottleneck, we should go for a
            ! more complicated rewrite
            CALL para_env%sync()
 
         END DO
      END DO
 
      DEALLOCATE (buf)
 
      IF (ionode) THEN
 
         gapw = .false.
         IF (dft_control%qs_control%method_id == do_method_gapw) gapw = .true.
 
         IF (do_voro /= 0) THEN
 
            IF (radius_type == voro_radii_unity) THEN
               ret = libvori_setradii_unity()
            ELSE IF (radius_type == voro_radii_cov) THEN
               ! Use the covalent radii from LIBVORI
               ret = libvori_setradii_covalent()
            ELSE IF (radius_type == voro_radii_vdw) THEN
               ! Use the van der Waals radii from CP2K
               ret = libvori_setradii_user(100.0_dp, particles_radius)
            ELSE IF (radius_type == voro_radii_user) THEN
               ! Use the user defined atomic radii
               IF (natom /= SIZE(user_radii)) THEN
                  CALL cp_abort(__location__, &
                                "Length of keyword VORONOI\USER_RADII does not "// &
                                "match number of atoms in the input coordinate file.")
               END IF
               ret = libvori_setradii_user(100.0_dp, user_radii)
            ELSE
               cpabort("No valid radius type was specified for VORONOI")
            END IF
 
            IF (voro_sanity) THEN
 
               ret = libvori_sanitycheck(sim_step, sim_time)
 
               IF (ret /= 0) THEN
                  cpabort("The library returned an error. Aborting.")
               END IF
 
            END IF
 
            ret = libvori_step(sim_step, sim_time)
 
            step_count = step_count + 1
 
            IF (ret /= 0) THEN
               cpabort("The library returned an error. Aborting.")
            END IF
 
            IF ((step_count > 1) .OR. (.NOT. voro_skip_first)) THEN
 
               ALLOCATE (voro_radii(natom))
               ALLOCATE (voro_charge(natom))
               ALLOCATE (voro_volume(natom))
               ALLOCATE (voro_dipole(natom, 3))
               ALLOCATE (voro_quadrupole(natom, 9))
               ALLOCATE (voro_buffer(natom))
               ALLOCATE (voro_wrapped_pos(natom, 3))
               ALLOCATE (voro_charge_center(natom, 3))
 
               ret = libvori_get_radius(natom, voro_radii)
 
               ret = libvori_get_charge(natom, voro_charge)
 
               ret = libvori_get_volume(natom, voro_volume)
 
               DO i1 = 1, 3
                  ret = libvori_get_dipole(i1, natom, voro_buffer)
                  voro_dipole(:, i1) = voro_buffer(:)
               END DO
 
               DO i1 = 1, 9
                  ret = libvori_get_quadrupole(i1, natom, voro_buffer)
                  voro_quadrupole(:, i1) = voro_buffer(:)
               END DO
 
               DO i1 = 1, 3
                  ret = libvori_get_wrapped_pos(i1, natom, voro_buffer)
                  voro_wrapped_pos(:, i1) = voro_buffer(:)
               END DO
 
               DO i1 = 1, 3
                  ret = libvori_get_charge_center(i1, natom, voro_buffer)
                  voro_charge_center(:, i1) = voro_buffer(:)
               END DO
 
               IF (gapw) THEN
                  CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
                  nspins = dft_control%nspins
                  DO i1 = 1, natom
                     voro_charge(i1) = voro_charge(i1) - sum(rho0_mpole%mp_rho(i1)%Q0(1:nspins))
                     fn0 = rho0_mpole%norm_g0l_h(1)/bohr*100._dp
                     voro_dipole(i1, 1:3) = voro_dipole(i1, 1:3) + rho0_mpole%mp_rho(i1)%Qlm_car(2:4)/fn0
                     qa = voro_charge(i1) - particles_c(i1)
                     voro_charge_center(i1, 1:3) = voro_dipole(i1, 1:3)/qa
                     fn1 = rho0_mpole%norm_g0l_h(2)/bohr/bohr*10000._dp
                     voro_quadrupole(i1, 1) = voro_quadrupole(i1, 1) + rho0_mpole%mp_rho(i1)%Qlm_car(5)/fn1
                     voro_quadrupole(i1, 2) = voro_quadrupole(i1, 2) + rho0_mpole%mp_rho(i1)%Qlm_car(6)/fn1
                     voro_quadrupole(i1, 3) = voro_quadrupole(i1, 3) + rho0_mpole%mp_rho(i1)%Qlm_car(7)/fn1
                     voro_quadrupole(i1, 4) = voro_quadrupole(i1, 4) + rho0_mpole%mp_rho(i1)%Qlm_car(6)/fn1
                     voro_quadrupole(i1, 5) = voro_quadrupole(i1, 5) + rho0_mpole%mp_rho(i1)%Qlm_car(8)/fn1
                     voro_quadrupole(i1, 6) = voro_quadrupole(i1, 6) + rho0_mpole%mp_rho(i1)%Qlm_car(9)/fn1
                     voro_quadrupole(i1, 7) = voro_quadrupole(i1, 7) + rho0_mpole%mp_rho(i1)%Qlm_car(7)/fn1
                     voro_quadrupole(i1, 8) = voro_quadrupole(i1, 8) + rho0_mpole%mp_rho(i1)%Qlm_car(9)/fn1
                     voro_quadrupole(i1, 9) = voro_quadrupole(i1, 9) + rho0_mpole%mp_rho(i1)%Qlm_car(10)/fn1
                  END DO
               END IF
 
               IF (unit_voro > 0) THEN
                  WRITE (unit_voro, fmt="(T2,I0)") natom
                  WRITE (unit_voro, fmt="(A,I8,A,F12.4,A)") "# Step ", sim_step, ",  Time ", &
                     sim_time*femtoseconds, " fs"
                  WRITE (unit_voro, fmt="(A,9F20.10)") "# Cell ", &
                     cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom, &
                     cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom, &
                     cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom
                  WRITE (unit_voro, fmt="(A,22A20)") "# Atom     Z", &
                     "Radius", "Position(X)", "Position(Y)", "Position(Z)", &
                     "Voronoi_Volume", "Z(eff)", "Charge", "Dipole(X)", "Dipole(Y)", "Dipole(Z)", &
                     "ChargeCenter(X)", "ChargeCenter(Y)", "ChargeCenter(Z)", &
                     "Quadrupole(XX)", "Quadrupole(XY)", "Quadrupole(XZ)", &
                     "Quadrupole(YX)", "Quadrupole(YY)", "Quadrupole(YZ)", &
                     "Quadrupole(ZX)", "Quadrupole(ZY)", "Quadrupole(ZZ)"
                  DO i1 = 1, natom
                     WRITE (unit_voro, fmt="(2I6,22F20.10)") &
                        i1, &
                        particles_z(i1), &
                        voro_radii(i1)/100.0_dp, &
                        particles_r(1:3, i1)*angstrom, &
                        voro_volume(i1)/1000000.0_dp, &
                        particles_c(i1), &
                        voro_charge(i1), &
                        voro_dipole(i1, 1:3), &
                        voro_charge_center(i1, 1:3)/100.0_dp, &
                        voro_quadrupole(i1, 1:9)
                  END DO
               END IF
 
               IF (molprop) THEN
                  CALL molecular_properties(subsys, cell, sim_step, sim_time, iounit, &
                                            particles_r, particles_c, &
                                            voro_charge, voro_charge_center, mp_file_name)
               END IF
 
               DEALLOCATE (voro_radii)
               DEALLOCATE (voro_charge)
               DEALLOCATE (voro_volume)
               DEALLOCATE (voro_dipole)
               DEALLOCATE (voro_quadrupole)
               DEALLOCATE (voro_buffer)
               DEALLOCATE (voro_wrapped_pos)
               DEALLOCATE (voro_charge_center)
 
            END IF ! not skip_first
 
            IF (iounit > 0) THEN
               WRITE (iounit, *) "VORONOI| Voronoi integration finished."
            END IF
 
         END IF ! do_voro
 
         IF (do_bqb /= 0) THEN
 
            ret = libvori_processbqbframe(sim_step, sim_time*femtoseconds)
 
            IF (ret /= 0) THEN
               cpabort("The library returned an error. Aborting.")
            END IF
 
            IF (do_voro /= 0) THEN
               IF (iounit > 0) THEN
                  WRITE (iounit, *) "VORONOI| BQB compression finished."
               END IF
            ELSE
               IF (iounit > 0) THEN
                  WRITE (iounit, *) "BQB| BQB compression finished."
               END IF
            END IF
 
         END IF  ! do_bqb
 
      END IF
 
      IF (ionode) THEN
         DEALLOCATE (particles_z)
         DEALLOCATE (particles_c)
         DEALLOCATE (particles_r)
         DEALLOCATE (particles_radius)
      END IF
      END associate
 
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_set_halting_mode(ieee_all, halt)
#endif
 
      CALL timestop(handle)
 
#else
 
      mark_used(do_voro)
      mark_used(do_bqb)
      mark_used(input_voro)
      mark_used(input_bqb)
      mark_used(unit_voro)
      mark_used(qs_env)
      mark_used(rspace_pw)
 
      CALL cp_warn(__location__, &
                   "Voronoi integration and BQB output require CP2k to be compiled"// &
                   " with the -D__LIBVORI preprocessor option.")
 
#endif
 
   END SUBROUTINE entry_voronoi_or_bqb
 
! **************************************************************************************************
!> \brief Call libvori's finalize if support is compiled in
! **************************************************************************************************
   SUBROUTINE finalize_libvori()
#if defined(__LIBVORI)
      INTEGER(KIND=C_INT) :: ret
      ret = libvori_finalize()
#endif
   END SUBROUTINE
 
! **************************************************************************************************
!> \brief ...
!> \param subsys ...
!> \param cell ...
!> \param sim_step ...
!> \param sim_time ...
!> \param iounit ...
!> \param particles_r ...
!> \param particles_c ...
!> \param voro_charge ...
!> \param voro_charge_center ...
!> \param mp_file_name ...
! **************************************************************************************************
   SUBROUTINE molecular_properties(subsys, cell, sim_step, sim_time, iounit, &
                                   particles_r, particles_c, voro_charge, &
                                   voro_charge_center, mp_file_name)
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(cell_type), POINTER                           :: cell
      INTEGER, INTENT(IN)                                :: sim_step
      REAL(kind=dp), INTENT(IN)                          :: sim_time
      INTEGER, INTENT(IN)                                :: iounit
      REAL(kind=dp), DIMENSION(:, :)                     :: particles_r
      REAL(kind=dp), DIMENSION(:)                        :: particles_c, voro_charge
      REAL(kind=dp), DIMENSION(:, :)                     :: voro_charge_center
      CHARACTER(len=default_path_length)                 :: mp_file_name
 
      CHARACTER(len=3)                                   :: fstatus
      CHARACTER(len=default_path_length)                 :: fname
      INTEGER                                            :: ia, imol, mk, mpunit, na, na1, na2, &
                                                            nmolecule
      REAL(kind=dp)                                      :: cm, ddip
      REAL(kind=dp), DIMENSION(3)                        :: dipm, posa, posc, ref
      TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
 
      IF (iounit > 0) THEN
         WRITE (iounit, *) "VORONOI| Start Calculation of Molecular Properties from Voronoi Integration"
      END IF
      CALL qs_subsys_get(subsys, molecule_set=molecule_set)
 
      IF (index(mp_file_name, "__STD_OUT__") /= 0) THEN
         mpunit = iounit
      ELSE
         fname = adjustl(mp_file_name)
         IF (fname(1:2) /= "./") THEN
            fname = trim(fname)//".molprop"
         END IF
         IF (file_exists(fname)) THEN
            fstatus = "old"
         ELSE
            fstatus = "new"
         END IF
         CALL open_file(file_name=fname, file_status=fstatus, file_action="write", &
                        file_position="append", unit_number=mpunit)
      END IF
      nmolecule = SIZE(molecule_set)
      WRITE (mpunit, fmt="(T2,I0)") nmolecule
      WRITE (mpunit, fmt="(A,I8,A,F12.4,A)") " # Step ", sim_step, ",  Time ", &
         sim_time*femtoseconds, " fs"
      WRITE (mpunit, fmt="(A,T25,A)") " #   Mol  Type    Charge", &
         "               Dipole[Debye]         Total Dipole[Debye]"
      DO imol = 1, nmolecule
         molecule_kind => molecule_set(imol)%molecule_kind
         mk = molecule_kind%kind_number
         na1 = molecule_set(imol)%first_atom
         na2 = molecule_set(imol)%last_atom
         na = na2 - na1 + 1
         ref(1:3) = 0.0_dp
         DO ia = na1, na2
            ref(1:3) = ref(1:3) + pbc(particles_r(1:3, ia), cell)
         END DO
         ref(1:3) = ref(1:3)/real(na, kind=dp)
         dipm = 0.0_dp
         DO ia = na1, na2
            posa(1:3) = particles_r(1:3, ia) - ref(1:3)
            posa(1:3) = pbc(posa, cell)
            posc(1:3) = posa(1:3) + bohr*voro_charge_center(ia, 1:3)/100.0_dp
            posc(1:3) = pbc(posc, cell)
            cm = -particles_c(ia) + voro_charge(ia)
            dipm(1:3) = dipm(1:3) + posa(1:3)*particles_c(ia) + posc(1:3)*cm
         END DO
         dipm(1:3) = dipm(1:3)*debye
         ddip = sqrt(sum(dipm**2))
         cm = sum(voro_charge(na1:na2))
         WRITE (mpunit, fmt="(I8,I6,F12.4,T25,3F12.4,8X,F12.4)") imol, mk, cm, dipm(1:3), ddip
      END DO
      IF (mpunit /= iounit) THEN
         CALL close_file(mpunit)
      END IF
 
   END SUBROUTINE molecular_properties
 
END MODULE voronoi_interface