cp_units.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief unit conversion facility
!>
!>      Units are complex, this module does not try to be very smart, for
!>      example SI prefixes are not supported automatically, and
!>      which kinds are really basic can change depending on the system of
!>      units chosen, and equivalences are not always catched.
!>
!>      This is thought as a simple conversion facility for the input and output.
!>      If you need something more you are probably better off using the
!>      physcon module directly.
!> \note
!>      One design choice was not to use dynamically allocated elements to
!>      reduce the possibility of leaks.
!>      Needs to be extended (for example charge,...)
!>      I just added the units and kinds that I needed.
!>      Used by the parser
!>      Should keep an unsorted/uncompressed version for nicer labels?
!> \par History
!>      01.2005 created [fawzi]
!> \author fawzi
! **************************************************************************************************
MODULE cp_units
 
   USE cp_log_handling,                 ONLY: cp_to_string
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: radians,&
                                              twopi
   USE physcon,                         ONLY: &
        a_bohr, atm, bar, bohr, debye, e_charge, e_mass, evolt, femtoseconds, joule, kcalmol, &
        kelvin, kjmol, massunit, newton, pascal, picoseconds, seconds, wavenumbers
   USE string_utilities,                ONLY: compress,&
                                              s2a,&
                                              uppercase
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER          :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_units'
 
   INTEGER, PARAMETER, PUBLIC :: cp_ukind_none = 0, &
                                 cp_ukind_energy = 1, &
                                 cp_ukind_length = 2, &
                                 cp_ukind_temperature = 3, &
                                 cp_ukind_angle = 4, &
                                 cp_ukind_pressure = 5, &
                                 cp_ukind_time = 6, &
                                 cp_ukind_mass = 7, &
                                 cp_ukind_undef = 8, &
                                 cp_ukind_potential = 9, &
                                 cp_ukind_force = 10, &
                                 cp_ukind_efield = 11, &
                                 cp_ukind_dipole = 12, &
                                 cp_ukind_max = 12
 
   ! General
   INTEGER, PARAMETER, PUBLIC :: cp_units_none = 100, &
                                 cp_units_au = 101
   ! Mass
   INTEGER, PARAMETER, PUBLIC :: cp_units_m_e = 110, &
                                 cp_units_amu = 111, &
                                 cp_units_kg = 112
   ! Energy
   INTEGER, PARAMETER, PUBLIC :: cp_units_hartree = 130, &
                                 cp_units_wavenum = 131, &
                                 cp_units_joule = 132, &
                                 cp_units_kcalmol = 133, &
                                 cp_units_ry = 134, &
                                 cp_units_ev = 135, &
                                 cp_units_kjmol = 136, &
                                 cp_units_jmol = 137, &
                                 cp_units_kev = 138
 
   ! Length
   INTEGER, PARAMETER, PUBLIC :: cp_units_bohr = 140, &
                                 cp_units_angstrom = 141, &
                                 cp_units_m = 142, &
                                 cp_units_pm = 143, &
                                 cp_units_nm = 144
 
   ! Temperature
   INTEGER, PARAMETER, PUBLIC :: cp_units_k = 150
 
   ! Pressure
   INTEGER, PARAMETER, PUBLIC :: cp_units_bar = 161
   INTEGER, PARAMETER, PUBLIC :: cp_units_atm = 162
   INTEGER, PARAMETER, PUBLIC :: cp_units_kbar = 163
   INTEGER, PARAMETER, PUBLIC :: cp_units_pa = 164
   INTEGER, PARAMETER, PUBLIC :: cp_units_mpa = 165
   INTEGER, PARAMETER, PUBLIC :: cp_units_gpa = 166
 
   ! Angles
   INTEGER, PARAMETER, PUBLIC :: cp_units_rad = 170, &
                                 cp_units_deg = 171
 
   ! Time
   INTEGER, PARAMETER, PUBLIC :: cp_units_fs = 180, &
                                 cp_units_s = 181, &
                                 cp_units_wn = 182, &
                                 cp_units_ps = 183
 
   ! Potential
   INTEGER, PARAMETER, PUBLIC :: cp_units_volt = 190
 
   ! Force
   INTEGER, PARAMETER, PUBLIC :: cp_units_newton = 200, &
                                 cp_units_mnewton = 201
 
   ! Electric Field
   INTEGER, PARAMETER, PUBLIC :: cp_units_volt_per_m = 202, &
                                 cp_units_volt_per_nm = 203, &
                                 cp_units_volt_per_angstrom = 204
 
   ! Dipole moment
   INTEGER, PARAMETER, PUBLIC :: cp_units_debye = 210, &
                                 cp_units_e_angstrom = 211, &
                                 cp_units_c_m = 212
 
   INTEGER, PARAMETER, PUBLIC :: cp_unit_max_kinds = 8, cp_unit_basic_desc_length = 15, &
                                 cp_unit_desc_length = cp_unit_max_kinds*cp_unit_basic_desc_length
 
   PUBLIC :: cp_unit_type, cp_unit_set_type
   PUBLIC :: cp_unit_create, cp_unit_release, &
             cp_unit_to_cp2k, cp_unit_from_cp2k, cp_unit_desc, &
             cp_unit_set_create, cp_unit_set_release, &
             cp_unit_to_cp2k1, cp_unit_from_cp2k1, cp_unit_compatible, export_units_as_xml
 
! **************************************************************************************************
!> \brief stores a unit
!> \param kind the kind of unit (energy, length,...)
!> \param unit the actual unit (Joule, eV,...)
!> \author fawzi
! **************************************************************************************************
   TYPE cp_unit_type
      INTEGER :: n_kinds = -1
      INTEGER, DIMENSION(cp_unit_max_kinds):: kind_id = -1, unit_id = -1, power = -1
   END TYPE cp_unit_type
 
! **************************************************************************************************
!> \brief represent a pointer to a unit (to build arrays of pointers)
!> \param unit the pointer to the unit
!> \author fawzi
! **************************************************************************************************
   TYPE cp_unit_p_type
      TYPE(cp_unit_type), POINTER :: unit => null()
   END TYPE cp_unit_p_type
 
! **************************************************************************************************
!> \brief stores the default units to be used
!> \author fawzi
! **************************************************************************************************
   TYPE cp_unit_set_type
      TYPE(cp_unit_p_type), DIMENSION(cp_ukind_max) :: units = cp_unit_p_type()
   END TYPE cp_unit_set_type
 
CONTAINS
 
! **************************************************************************************************
!> \brief creates a unit parsing a string
!> \param unit the unit to initialize
!> \param string the string containing the description of the unit
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE cp_unit_create(unit, string)
      TYPE(cp_unit_type), INTENT(OUT)                    :: unit
      CHARACTER(len=*), INTENT(in)                       :: string
 
      CHARACTER(LEN=40)                                  :: formatstr
      CHARACTER(LEN=cp_unit_desc_length)                 :: desc
      CHARACTER(LEN=LEN(string))                         :: unit_string
      INTEGER                                            :: i_high, i_low, i_unit, len_string, &
                                                            next_power
      INTEGER, DIMENSION(cp_unit_max_kinds)              :: kind_id, power, unit_id
 
      unit_id = cp_units_none
      kind_id = cp_ukind_none
      power = 0
      i_low = 1
      i_high = 1
      len_string = len(string)
      i_unit = 0
      next_power = 1
      DO WHILE (i_low < len_string)
         IF (string(i_low:i_low) /= ' ') EXIT
         i_low = i_low + 1
      END DO
      i_high = i_low
      DO WHILE (i_high <= len_string)
         IF (string(i_high:i_high) == ' ' .OR. string(i_high:i_high) == '^' .OR. &
             string(i_high:i_high) == '*' .OR. string(i_high:i_high) == '/') EXIT
         i_high = i_high + 1
      END DO
      DO
         IF (i_high <= i_low .OR. i_low > len_string) EXIT
         i_unit = i_unit + 1
         IF (i_unit > cp_unit_max_kinds) THEN
            cpabort("Maximum number of combined units exceeded")
            EXIT
         END IF
         ! read unit
         unit_string = string(i_low:i_high - 1)
         CALL uppercase(unit_string)
         SELECT CASE (trim(unit_string))
         CASE ("INTERNAL_CP2K")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_undef
         CASE ("HARTREE")
            unit_id(i_unit) = cp_units_hartree
            kind_id(i_unit) = cp_ukind_energy
         CASE ("AU_E")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_energy
         CASE ("WAVENUMBER_E")
            unit_id(i_unit) = cp_units_wavenum
            kind_id(i_unit) = cp_ukind_energy
         CASE ("JOULE", "J")
            unit_id(i_unit) = cp_units_joule
            kind_id(i_unit) = cp_ukind_energy
         CASE ("KCALMOL")
            unit_id(i_unit) = cp_units_kcalmol
            kind_id(i_unit) = cp_ukind_energy
         CASE ("KJMOL")
            unit_id(i_unit) = cp_units_kjmol
            kind_id(i_unit) = cp_ukind_energy
         CASE ("JMOL")
            unit_id(i_unit) = cp_units_jmol
            kind_id(i_unit) = cp_ukind_energy
         CASE ("RY")
            unit_id(i_unit) = cp_units_ry
            kind_id(i_unit) = cp_ukind_energy
         CASE ("EV")
            unit_id(i_unit) = cp_units_ev
            kind_id(i_unit) = cp_ukind_energy
         CASE ("KEV")
            unit_id(i_unit) = cp_units_kev
            kind_id(i_unit) = cp_ukind_energy
         CASE ("K_E")
            unit_id(i_unit) = cp_units_k
            kind_id(i_unit) = cp_ukind_energy
         CASE ("ENERGY")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_energy
         CASE ("AU_L")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_length
         CASE ("BOHR")
            unit_id(i_unit) = cp_units_bohr
            kind_id(i_unit) = cp_ukind_length
         CASE ("M")
            unit_id(i_unit) = cp_units_m
            kind_id(i_unit) = cp_ukind_length
         CASE ("PM")
            unit_id(i_unit) = cp_units_pm
            kind_id(i_unit) = cp_ukind_length
         CASE ("NM")
            unit_id(i_unit) = cp_units_nm
            kind_id(i_unit) = cp_ukind_length
         CASE ("ANGSTROM")
            unit_id(i_unit) = cp_units_angstrom
            kind_id(i_unit) = cp_ukind_length
         CASE ("LENGTH")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_length
         CASE ("K", "K_TEMP")
            unit_id(i_unit) = cp_units_k
            kind_id(i_unit) = cp_ukind_temperature
         CASE ("AU_TEMP")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_temperature
         CASE ("TEMPERATURE")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_temperature
         CASE ("ATM")
            unit_id(i_unit) = cp_units_atm
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("BAR")
            unit_id(i_unit) = cp_units_bar
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("KBAR")
            unit_id(i_unit) = cp_units_kbar
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("PA")
            unit_id(i_unit) = cp_units_pa
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("MPA")
            unit_id(i_unit) = cp_units_mpa
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("GPA")
            unit_id(i_unit) = cp_units_gpa
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("AU_P")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("PRESSURE")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_pressure
         CASE ("RAD")
            unit_id(i_unit) = cp_units_rad
            kind_id(i_unit) = cp_ukind_angle
         CASE ("DEG")
            unit_id(i_unit) = cp_units_deg
            kind_id(i_unit) = cp_ukind_angle
         CASE ("ANGLE")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_angle
         CASE ("S")
            unit_id(i_unit) = cp_units_s
            kind_id(i_unit) = cp_ukind_time
         CASE ("FS")
            unit_id(i_unit) = cp_units_fs
            kind_id(i_unit) = cp_ukind_time
         CASE ("PS")
            unit_id(i_unit) = cp_units_ps
            kind_id(i_unit) = cp_ukind_time
         CASE ("WAVENUMBER_T")
            unit_id(i_unit) = cp_units_wn
            kind_id(i_unit) = cp_ukind_time
         CASE ("AU_T")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_time
         CASE ("TIME")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_time
         CASE ("KG")
            unit_id(i_unit) = cp_units_kg
            kind_id(i_unit) = cp_ukind_mass
         CASE ("AMU")
            unit_id(i_unit) = cp_units_amu
            kind_id(i_unit) = cp_ukind_mass
         CASE ("M_E")
            unit_id(i_unit) = cp_units_m_e
            kind_id(i_unit) = cp_ukind_mass
         CASE ("AU_M")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_mass
         CASE ("MASS")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_mass
         CASE ("VOLT")
            unit_id(i_unit) = cp_units_volt
            kind_id(i_unit) = cp_ukind_potential
         CASE ("AU_POT")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_potential
         CASE ("POTENTIAL")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_potential
         CASE ("N", "NEWTON")
            unit_id(i_unit) = cp_units_newton
            kind_id(i_unit) = cp_ukind_force
         CASE ("MN", "MNEWTON")
            unit_id(i_unit) = cp_units_mnewton
            kind_id(i_unit) = cp_ukind_force
         CASE ("AU_F")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_force
         CASE ("FORCE")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_force
         CASE ("VM-1", "VOLT_PER_M")
            unit_id(i_unit) = cp_units_volt_per_m
            kind_id(i_unit) = cp_ukind_efield
         CASE ("VNM-1", "VOLT_PER_NM")
            unit_id(i_unit) = cp_units_volt_per_nm
            kind_id(i_unit) = cp_ukind_efield
         CASE ("VA-1", "VOLT_PER_ANGSTROM")
            unit_id(i_unit) = cp_units_volt_per_angstrom
            kind_id(i_unit) = cp_ukind_efield
         CASE ("AU_EFIELD")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_efield
         CASE ("EFIELD")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_efield
         CASE ("DEBYE")
            unit_id(i_unit) = cp_units_debye
            kind_id(i_unit) = cp_ukind_dipole
         CASE ("E_ANGSTROM")
            unit_id(i_unit) = cp_units_e_angstrom
            kind_id(i_unit) = cp_ukind_dipole
         CASE ("C_M")
            unit_id(i_unit) = cp_units_c_m
            kind_id(i_unit) = cp_ukind_dipole
         CASE ("AU_DIP")
            unit_id(i_unit) = cp_units_au
            kind_id(i_unit) = cp_ukind_dipole
         CASE ("DIPOLE")
            unit_id(i_unit) = cp_units_none
            kind_id(i_unit) = cp_ukind_dipole
         CASE ("AU")
            CALL cp_abort(__location__, &
                          "au unit without specifying its kind not accepted, use "// &
                          "(au_e, au_f, au_t, au_temp, au_l, au_m, au_p, au_pot, "// &
                          "au_efield, au_dip)")
         CASE default
            cpabort("Unknown unit: "//string(i_low:i_high - 1))
         END SELECT
         power(i_unit) = next_power
         ! parse op
         i_low = i_high
         DO WHILE (i_low <= len_string)
            IF (string(i_low:i_low) /= ' ') EXIT
            i_low = i_low + 1
         END DO
         i_high = i_low
         DO WHILE (i_high <= len_string)
            IF (string(i_high:i_high) == ' ' .OR. string(i_high:i_high) == '^' .OR. &
                string(i_high:i_high) == '*' .OR. string(i_high:i_high) == '/') EXIT
            i_high = i_high + 1
         END DO
         IF (i_high < i_low .OR. i_low > len_string) EXIT
 
         IF (i_high <= len_string) THEN
            IF (string(i_low:i_high) == '^') THEN
               i_low = i_high + 1
               DO WHILE (i_low <= len_string)
                  IF (string(i_low:i_low) /= ' ') EXIT
                  i_low = i_low + 1
               END DO
               i_high = i_low
               DO WHILE (i_high <= len_string)
                  SELECT CASE (string(i_high:i_high))
                  CASE ('+', '-', '0', '1', '2', '3', '4', '5', '6', '7', '8', '9')
                     i_high = i_high + 1
                  CASE default
                     EXIT
                  END SELECT
               END DO
               IF (i_high <= i_low .OR. i_low > len_string) THEN
                  cpabort("an integer number is expected after a '^'")
                  EXIT
               END IF
               formatstr = "(i"//cp_to_string(i_high - i_low + 1)//")"
               READ (string(i_low:i_high - 1), formatstr) &
                  next_power
               power(i_unit) = power(i_unit)*next_power
               ! next op
               i_low = i_high
               DO WHILE (i_low < len_string)
                  IF (string(i_low:i_low) /= ' ') EXIT
                  i_low = i_low + 1
               END DO
               i_high = i_low
               DO WHILE (i_high <= len_string)
                  IF (string(i_high:i_high) == ' ' .OR. string(i_high:i_high) == '^' .OR. &
                      string(i_high:i_high) == '*' .OR. string(i_high:i_high) == '/') EXIT
                  i_high = i_high + 1
               END DO
            END IF
         END IF
         IF (i_low > len_string) EXIT
         next_power = 1
         IF (i_high <= len_string) THEN
            IF (string(i_low:i_high) == "*" .OR. string(i_low:i_high) == '/') THEN
               IF (string(i_low:i_high) == '/') next_power = -1
               i_low = i_high + 1
               DO WHILE (i_low <= len_string)
                  IF (string(i_low:i_low) /= ' ') EXIT
                  i_low = i_low + 1
               END DO
               i_high = i_low
               DO WHILE (i_high <= len_string)
                  IF (string(i_high:i_high) == ' ' .OR. string(i_high:i_high) == '^' .OR. &
                      string(i_high:i_high) == '*' .OR. string(i_high:i_high) == '/') EXIT
                  i_high = i_high + 1
               END DO
            END IF
         END IF
      END DO
      CALL cp_unit_create2(unit, kind_id=kind_id, unit_id=unit_id, &
                           power=power)
      desc = cp_unit_desc(unit)
   END SUBROUTINE cp_unit_create
 
! **************************************************************************************************
!> \brief creates and initializes the given unit of mesure (performs some error
!>      check)
!> \param unit the unit descriptor to be initialized
!> \param kind_id the kind of unit (length,energy,...), use the constants
!>        cp_ukind_*
!> \param unit_id the actual unit (use constants cp_units_*)
!> \param power ...
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE cp_unit_create2(unit, kind_id, unit_id, power)
      TYPE(cp_unit_type), INTENT(OUT)                    :: unit
      INTEGER, DIMENSION(:), INTENT(in)                  :: kind_id, unit_id
      INTEGER, DIMENSION(:), INTENT(in), OPTIONAL        :: power
 
      INTEGER                                            :: i, j, max_kind, max_pos
      LOGICAL                                            :: repeat
 
      cpassert(SIZE(kind_id) <= cp_unit_max_kinds)
      cpassert(SIZE(unit_id) <= cp_unit_max_kinds)
      unit%kind_id(1:SIZE(kind_id)) = kind_id
      unit%kind_id(SIZE(kind_id) + 1:) = cp_ukind_none
      unit%unit_id(1:SIZE(unit_id)) = unit_id
      unit%unit_id(SIZE(unit_id):) = cp_units_none
      IF (PRESENT(power)) THEN
         unit%power(1:SIZE(power)) = power
         unit%power(SIZE(power) + 1:) = 0
         DO i = 1, SIZE(unit%power)
            IF (unit%power(i) == 0) THEN
               unit%kind_id(i) = cp_ukind_none
               unit%unit_id(i) = cp_units_none
            END IF
         END DO
      ELSE
         DO i = 1, SIZE(unit%power)
            IF (unit%unit_id(i) /= 0) THEN
               unit%power(i) = 1
            ELSE
               unit%power(i) = 0
            END IF
         END DO
      END IF
 
      ! remove unnecessary units
      ! reorder & compress
      unit%n_kinds = 0
      DO i = 1, SIZE(unit%kind_id)
         ! find max and compress in the rest
         DO
            max_kind = unit%kind_id(i)
            max_pos = i
            repeat = .false.
            DO j = i + 1, SIZE(unit%kind_id)
               IF (unit%kind_id(j) >= max_kind) THEN
                  IF (unit%kind_id(j) /= 0 .AND. unit%kind_id(j) == max_kind .AND. &
                      unit%unit_id(j) == unit%unit_id(max_pos)) THEN
                     unit%power(max_pos) = unit%power(max_pos) + unit%power(j)
                     unit%kind_id(j) = cp_ukind_none
                     unit%unit_id(j) = cp_units_none
                     unit%power(j) = 0
                     IF (unit%power(max_pos) == 0) THEN
                        unit%kind_id(max_pos) = cp_ukind_none
                        unit%unit_id(max_pos) = cp_units_none
                        unit%power(max_pos) = 0
                        repeat = .true.
                        EXIT
                     END IF
                  ELSE IF (unit%kind_id(j) > max_kind .OR. &
                           (unit%kind_id(j) == max_kind .AND. &
                            unit%unit_id(j) > unit%unit_id(max_pos))) THEN
                     max_kind = unit%kind_id(j)
                     max_pos = j
                  END IF
               END IF
            END DO
            IF (.NOT. repeat) EXIT
         END DO
         IF (max_kind /= 0) unit%n_kinds = unit%n_kinds + 1
         ! put the max at pos i
         IF (max_pos /= i) THEN
            unit%kind_id(max_pos) = unit%kind_id(i)
            unit%kind_id(i) = max_kind
            max_kind = unit%unit_id(max_pos)
            unit%unit_id(max_pos) = unit%unit_id(i)
            unit%unit_id(i) = max_kind
            max_kind = unit%power(max_pos)
            unit%power(max_pos) = unit%power(i)
            unit%power(i) = max_kind
         END IF
         ! check unit
         CALL cp_basic_unit_check(basic_kind=unit%kind_id(i), &
                                  basic_unit=unit%unit_id(i))
      END DO
   END SUBROUTINE cp_unit_create2
 
! **************************************************************************************************
!> \brief releases the given unit
!> \param unit the unit to release
!> \author fawzi
!> \note
!>      at the moment not needed, there for completeness
! **************************************************************************************************
   ELEMENTAL SUBROUTINE cp_unit_release(unit)
      TYPE(cp_unit_type), INTENT(IN)                     :: unit
 
      mark_used(unit)
 
   END SUBROUTINE cp_unit_release
 
! **************************************************************************************************
!> \brief controls that the kind and contains meaningful information
!> \param basic_kind the kind of the unit
!> \param basic_unit the unit to check
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE cp_basic_unit_check(basic_kind, basic_unit)
      INTEGER, INTENT(in)                                :: basic_kind, basic_unit
 
      SELECT CASE (basic_kind)
      CASE (cp_ukind_undef)
         SELECT CASE (basic_unit)
         CASE (cp_units_none)
         CASE default
            cpabort("unknown undef unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_energy)
         SELECT CASE (basic_unit)
         CASE (cp_units_hartree, cp_units_wavenum, cp_units_joule, cp_units_kcalmol, &
               cp_units_kjmol, cp_units_ry, cp_units_ev, cp_units_kev, cp_units_au, cp_units_k, &
               cp_units_jmol, cp_units_none)
         CASE default
            cpabort("unknown energy unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_length)
         SELECT CASE (basic_unit)
         CASE (cp_units_bohr, cp_units_angstrom, cp_units_au, cp_units_none, cp_units_m, &
               cp_units_pm, cp_units_nm)
         CASE default
            cpabort("unknown length unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_temperature)
         SELECT CASE (basic_unit)
         CASE (cp_units_k, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown temperature unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_pressure)
         SELECT CASE (basic_unit)
         CASE (cp_units_bar, cp_units_atm, cp_units_kbar, cp_units_pa, cp_units_mpa, cp_units_gpa, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown pressure unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_angle)
         SELECT CASE (basic_unit)
         CASE (cp_units_rad, cp_units_deg, cp_units_none)
         CASE default
            cpabort("unknown angle unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_time)
         SELECT CASE (basic_unit)
         CASE (cp_units_s, cp_units_fs, cp_units_ps, cp_units_au, cp_units_wn, cp_units_none)
         CASE default
            cpabort("unknown time unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_mass)
         SELECT CASE (basic_unit)
         CASE (cp_units_kg, cp_units_amu, cp_units_m_e, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown mass unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_potential)
         SELECT CASE (basic_unit)
         CASE (cp_units_volt, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown potential unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_force)
         SELECT CASE (basic_unit)
         CASE (cp_units_newton, cp_units_mnewton, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown force unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_efield)
         SELECT CASE (basic_unit)
         CASE (cp_units_volt_per_m, cp_units_volt_per_nm, &
               cp_units_volt_per_angstrom, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown electric field unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_dipole)
         SELECT CASE (basic_unit)
         CASE (cp_units_debye, cp_units_e_angstrom, cp_units_c_m, cp_units_au, cp_units_none)
         CASE default
            cpabort("unknown dipole unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_none)
         IF (basic_unit /= cp_units_none) THEN
            CALL cp_abort(__location__, &
                          "if the kind of the unit is none also unit must be undefined,not:" &
                          //trim(cp_to_string(basic_unit)))
         END IF
      CASE default
         cpabort("unknown kind of unit:"//trim(cp_to_string(basic_kind)))
      END SELECT
   END SUBROUTINE cp_basic_unit_check
 
! **************************************************************************************************
!> \brief converts a value to the internal cp2k units
!> \param value the value to convert
!> \param basic_kind the kind of the unit of the value
!> \param basic_unit the unit of the value
!> \param power the power of the unit (defaults to 1)
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_basic_unit_to_cp2k(value, basic_kind, basic_unit, power) RESULT(res)
      REAL(kind=dp), INTENT(in)                          :: value
      INTEGER, INTENT(in)                                :: basic_kind, basic_unit
      INTEGER, INTENT(in), OPTIONAL                      :: power
      REAL(kind=dp)                                      :: res
 
      INTEGER                                            :: my_power
 
      my_power = 1
      IF (PRESENT(power)) my_power = power
      IF (basic_unit == cp_units_none .AND. basic_kind /= cp_ukind_undef) THEN
         IF (basic_kind /= cp_units_none) THEN
            CALL cp_abort(__location__, &
                          "unit not yet fully specified, unit of kind "// &
                          trim(cp_to_string(basic_unit)))
         END IF
      END IF
      SELECT CASE (basic_kind)
      CASE (cp_ukind_undef)
         SELECT CASE (basic_unit)
         CASE (cp_units_none)
            res = value
         CASE default
            cpabort("unknown energy unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_energy)
         SELECT CASE (basic_unit)
         CASE (cp_units_hartree, cp_units_au)
            res = value
         CASE (cp_units_wavenum)
            res = wavenumbers**(-my_power)*value
         CASE (cp_units_joule)
            res = joule**(-my_power)*value
         CASE (cp_units_kcalmol)
            res = kcalmol**(-my_power)*value
         CASE (cp_units_kjmol)
            res = kjmol**(-my_power)*value
         CASE (cp_units_jmol)
            res = (kjmol*1.0e+3_dp)**(-my_power)*value
         CASE (cp_units_ry)
            res = 0.5_dp**my_power*value
         CASE (cp_units_ev)
            res = evolt**(-my_power)*value
         CASE (cp_units_kev)
            res = (1.0e-3_dp*evolt)**(-my_power)*value
         CASE (cp_units_k)
            res = kelvin**(-my_power)*value
         CASE default
            cpabort("unknown energy unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_length)
         SELECT CASE (basic_unit)
         CASE (cp_units_bohr, cp_units_au)
            res = value
         CASE (cp_units_m)
            res = value*(1.0e10_dp*bohr)**my_power
         CASE (cp_units_pm)
            res = value*(0.01_dp*bohr)**my_power
         CASE (cp_units_nm)
            res = value*(10.0_dp*bohr)**my_power
         CASE (cp_units_angstrom)
            res = value*bohr**my_power
         CASE default
            cpabort("unknown length unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_temperature)
         SELECT CASE (basic_unit)
         CASE (cp_units_k)
            res = kelvin**(-my_power)*value
         CASE (cp_units_au)
            res = value
         CASE default
            cpabort("unknown temperature unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_pressure)
         SELECT CASE (basic_unit)
         CASE (cp_units_bar)
            res = bar**(-my_power)*value
         CASE (cp_units_atm)
            res = atm**(-my_power)*value
         CASE (cp_units_kbar)
            res = (1.0e-3_dp*bar)**(-my_power)*value
         CASE (cp_units_pa)
            res = pascal**(-my_power)*value
         CASE (cp_units_mpa)
            res = (1.0e-6_dp*pascal)**(-my_power)*value
         CASE (cp_units_gpa)
            res = (1.0e-9_dp*pascal)**(-my_power)*value
         CASE (cp_units_au)
            res = value
         CASE default
            cpabort("unknown pressure unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_angle)
         SELECT CASE (basic_unit)
         CASE (cp_units_rad)
            res = value
         CASE (cp_units_deg)
            res = value*(radians)**my_power
         CASE default
            cpabort("unknown angle unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_time)
         SELECT CASE (basic_unit)
         CASE (cp_units_s)
            res = value*seconds**(-my_power)
         CASE (cp_units_fs)
            res = value*femtoseconds**(-my_power)
         CASE (cp_units_ps)
            res = value*picoseconds**(-my_power)
         CASE (cp_units_au)
            res = value
         CASE (cp_units_wn)
            res = (twopi*wavenumbers)**(my_power)/value
         CASE default
            cpabort("unknown time unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_mass)
         SELECT CASE (basic_unit)
         CASE (cp_units_kg)
            res = e_mass**my_power*value
         CASE (cp_units_amu)
            res = massunit**my_power*value
         CASE (cp_units_m_e, cp_units_au)
            res = value
         CASE default
            cpabort("unknown mass unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_potential)
         SELECT CASE (basic_unit)
         CASE (cp_units_volt)
            res = evolt**(-my_power)*value
         CASE (cp_units_au)
            res = value
         CASE default
            cpabort("unknown potential unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_force)
         SELECT CASE (basic_unit)
         CASE (cp_units_newton)
            res = value*newton**(-my_power)
         CASE (cp_units_mnewton)
            res = value*(1.0e+3_dp*newton)**(-my_power)
         CASE (cp_units_au)
            res = value
         CASE default
            cpabort("unknown force unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_efield)
         SELECT CASE (basic_unit)
         CASE (cp_units_volt_per_m)
            res = (1.0e+10_dp*evolt*bohr)**(-my_power)*value
         CASE (cp_units_volt_per_nm)
            res = (10.0_dp*evolt*bohr)**(-my_power)*value
         CASE (cp_units_volt_per_angstrom)
            res = (evolt*bohr)**(-my_power)*value
         CASE (cp_units_au)
            res = value
         CASE default
            cpabort("unknown electric field unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_dipole)
         SELECT CASE (basic_unit)
         CASE (cp_units_debye)
            res = debye**(-my_power)*value
         CASE (cp_units_e_angstrom)
            res = value*bohr**my_power
         CASE (cp_units_c_m)
            res = (e_charge*a_bohr)**(-my_power)*value
         CASE (cp_units_au)
            res = value
         CASE default
            cpabort("unknown dipole unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_none)
         CALL cp_abort(__location__, &
                       "if the kind of the unit is none also unit must be undefined,not:" &
                       //trim(cp_to_string(basic_unit)))
      CASE default
         cpabort("unknown kind of unit:"//trim(cp_to_string(basic_kind)))
      END SELECT
   END FUNCTION cp_basic_unit_to_cp2k
 
! **************************************************************************************************
!> \brief returns the label of the current basic unit
!> \param basic_kind the kind of the unit of the value
!> \param basic_unit the unit of the value
!> \param power the power of the unit (defaults to 1)
!> \param accept_undefined ...
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_basic_unit_desc(basic_kind, basic_unit, power, accept_undefined) &
      result(res)
      INTEGER, INTENT(in)                                :: basic_kind, basic_unit
      INTEGER, INTENT(in), OPTIONAL                      :: power
      LOGICAL, INTENT(in), OPTIONAL                      :: accept_undefined
      CHARACTER(len=cp_unit_basic_desc_length)           :: res
 
      INTEGER                                            :: a, my_power
      LOGICAL                                            :: my_accept_undefined
 
      my_power = 1
      res = ""
      my_accept_undefined = .false.
      IF (accept_undefined) my_accept_undefined = accept_undefined
      IF (PRESENT(power)) my_power = power
      IF (basic_unit == cp_units_none) THEN
         IF (.NOT. my_accept_undefined .AND. basic_kind == cp_units_none) THEN
            CALL cp_abort(__location__, "unit not yet fully specified, unit of kind "// &
                          trim(cp_to_string(basic_kind)))
         END IF
      END IF
      SELECT CASE (basic_kind)
      CASE (cp_ukind_undef)
         SELECT CASE (basic_unit)
         CASE (cp_units_none)
            res = "internal_cp2k"
         CASE DEFAULT
            CALL cp_abort(__location__, &
                          "unit not yet fully specified, unit of kind "// &
                          trim(res))
         END SELECT
      CASE (cp_ukind_energy)
         SELECT CASE (basic_unit)
         CASE (cp_units_hartree, cp_units_au)
            res = "hartree"
         CASE (cp_units_wavenum)
            res = "wavenumber_e"
         CASE (cp_units_joule)
            res = "joule"
         CASE (cp_units_kcalmol)
            res = "kcalmol"
         CASE (cp_units_kjmol)
            res = "kjmol"
         CASE (cp_units_jmol)
            res = "jmol"
         CASE (cp_units_ry)
            res = "Ry"
         CASE (cp_units_ev)
            res = "eV"
         CASE (cp_units_kev)
            res = "keV"
         CASE (cp_units_k)
            res = "K_e"
         CASE (cp_units_none)
            res = "energy"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown energy unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_length)
         SELECT CASE (basic_unit)
         CASE (cp_units_bohr, cp_units_au)
            res = "bohr"
         CASE (cp_units_m)
            res = "m"
         CASE (cp_units_pm)
            res = "pm"
         CASE (cp_units_nm)
            res = "nm"
         CASE (cp_units_angstrom)
            res = "angstrom"
         CASE default
            res = "length"
            cpabort("unknown length unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_temperature)
         SELECT CASE (basic_unit)
         CASE (cp_units_k)
            res = "K"
         CASE (cp_units_au)
            res = "au_temp"
         CASE (cp_units_none)
            res = "temperature"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown temperature unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_pressure)
         SELECT CASE (basic_unit)
         CASE (cp_units_bar)
            res = "bar"
         CASE (cp_units_atm)
            res = "atm"
         CASE (cp_units_kbar)
            res = "kbar"
         CASE (cp_units_pa)
            res = "Pa"
         CASE (cp_units_mpa)
            res = "MPa"
         CASE (cp_units_gpa)
            res = "GPa"
         CASE (cp_units_au)
            res = "au_p"
         CASE (cp_units_none)
            res = "pressure"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown pressure unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_angle)
         SELECT CASE (basic_unit)
         CASE (cp_units_rad)
            res = "rad"
         CASE (cp_units_deg)
            res = "deg"
         CASE (cp_units_none)
            res = "angle"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown angle unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_time)
         SELECT CASE (basic_unit)
         CASE (cp_units_s)
            res = "s"
         CASE (cp_units_fs)
            res = "fs"
         CASE (cp_units_ps)
            res = "ps"
         CASE (cp_units_au)
            res = "au_t"
         CASE (cp_units_wn)
            res = "wavenumber_t"
         CASE (cp_units_none)
            res = "time"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown time unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_mass)
         SELECT CASE (basic_unit)
         CASE (cp_units_kg)
            res = "kg"
         CASE (cp_units_amu)
            res = "amu"
         CASE (cp_units_m_e, cp_units_au)
            res = "m_e"
         CASE (cp_units_none)
            res = "mass"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown mass unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_potential)
         SELECT CASE (basic_unit)
         CASE (cp_units_volt)
            res = "volt"
         CASE (cp_units_au)
            res = "au_pot"
         CASE (cp_units_none)
            res = "potential"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown potential unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_force)
         SELECT CASE (basic_unit)
         CASE (cp_units_newton)
            res = "N"
         CASE (cp_units_mnewton)
            res = "mN"
         CASE (cp_units_au)
            res = "au_f"
         CASE (cp_units_none)
            res = "force"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown potential unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_efield)
         SELECT CASE (basic_unit)
         CASE (cp_units_volt_per_m)
            res = "Vm-1"
         CASE (cp_units_volt_per_nm)
            res = "Vnm-1"
         CASE (cp_units_volt_per_angstrom)
            res = "Vangstrom-1"
         CASE (cp_units_au)
            res = "au_efield"
         CASE (cp_units_none)
            res = "electric field"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown efield unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_dipole)
         SELECT CASE (basic_unit)
         CASE (cp_units_debye)
            res = "debye"
         CASE (cp_units_e_angstrom)
            res = "e_angstrom"
         CASE (cp_units_c_m)
            res = "C_m"
         CASE (cp_units_au)
            res = "au_dip"
         CASE (cp_units_none)
            res = "dipole"
            IF (.NOT. my_accept_undefined) THEN
               CALL cp_abort(__location__, &
                             "unit not yet fully specified, unit of kind "// &
                             trim(res))
            END IF
         CASE default
            cpabort("unknown dipole unit:"//trim(cp_to_string(basic_unit)))
         END SELECT
      CASE (cp_ukind_none)
         CALL cp_abort(__location__, &
                       "if the kind of the unit is none also unit must be undefined,not:" &
                       //trim(cp_to_string(basic_unit)))
      CASE default
         cpabort("unknown kind of unit:"//trim(cp_to_string(basic_kind)))
      END SELECT
      IF (my_power /= 1) THEN
         a = len_trim(res)
         cpassert(len(res) - a >= 3)
         WRITE (res(a + 1:), "('^',i3)") my_power
         CALL compress(res, .true.)
      END IF
   END FUNCTION cp_basic_unit_desc
 
! **************************************************************************************************
!> \brief returns the "name" of the given unit
!> \param unit the unit to describe
!> \param defaults defaults for the undefined units, optional
!> \param accept_undefined if defaults is not present or is not associated
!>        whether undefined units should be accepted (defaults to false)
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_unit_desc(unit, defaults, accept_undefined) &
      result(res)
      TYPE(cp_unit_type), INTENT(IN)                     :: unit
      TYPE(cp_unit_set_type), INTENT(IN), OPTIONAL       :: defaults
      LOGICAL, INTENT(in), OPTIONAL                      :: accept_undefined
      CHARACTER(len=cp_unit_desc_length)                 :: res
 
      INTEGER                                            :: i, my_unit, pos
      LOGICAL                                            :: check, has_defaults, my_accept_undefined
 
      res = ""
      pos = 1
      my_accept_undefined = .false.
      IF (PRESENT(accept_undefined)) my_accept_undefined = accept_undefined
      DO i = 1, unit%n_kinds
         cpassert(unit%kind_id(i) /= 0)
         cpassert(pos < len(res))
         my_unit = unit%unit_id(i)
         has_defaults = .false.
         IF (PRESENT(defaults)) has_defaults = ASSOCIATED(defaults%units(1)%unit)
         IF (my_unit == 0) THEN
            IF (has_defaults) THEN
               my_unit = defaults%units(unit%kind_id(i))%unit%unit_id(1)
            ELSE
               check = my_accept_undefined .OR. unit%kind_id(i) /= 0
               cpassert(check)
            END IF
         END IF
         IF (i > 1) THEN
            res(pos:pos) = "*"
            pos = pos + 1
         END IF
         res(pos:) = trim(cp_basic_unit_desc(basic_kind=unit%kind_id(i), &
                                             basic_unit=my_unit, accept_undefined=my_accept_undefined, &
                                             power=unit%power(i)))
         pos = len_trim(res) + 1
      END DO
 
   END FUNCTION cp_unit_desc
 
! **************************************************************************************************
!> \brief transform a value to the internal cp2k units
!> \param value the value to convert
!> \param unit the unit of the result
!> \param defaults the defaults unit for those that are left free
!>        (cp_units_none)
!> \param power the power of the unit (defaults to 1)
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_unit_to_cp2k1(value, unit, defaults, power) RESULT(res)
      REAL(kind=dp), INTENT(in)                          :: value
      TYPE(cp_unit_type), INTENT(IN)                     :: unit
      TYPE(cp_unit_set_type), INTENT(IN), OPTIONAL       :: defaults
      INTEGER, INTENT(in), OPTIONAL                      :: power
      REAL(kind=dp)                                      :: res
 
      INTEGER                                            :: i_unit, my_basic_unit, my_power
 
      my_power = 1
      IF (PRESENT(power)) my_power = power
      res = value
      DO i_unit = 1, unit%n_kinds
         cpassert(unit%kind_id(i_unit) > 0)
         my_basic_unit = unit%unit_id(i_unit)
         IF (my_basic_unit == 0 .AND. unit%kind_id(i_unit) /= cp_ukind_undef) THEN
            cpassert(PRESENT(defaults))
            cpassert(ASSOCIATED(defaults%units(unit%kind_id(i_unit))%unit))
            my_basic_unit = defaults%units(unit%kind_id(i_unit))%unit%unit_id(1)
         END IF
         res = cp_basic_unit_to_cp2k(value=res, basic_unit=my_basic_unit, &
                                     basic_kind=unit%kind_id(i_unit), &
                                     power=my_power*unit%power(i_unit))
      END DO
   END FUNCTION cp_unit_to_cp2k1
 
! **************************************************************************************************
!> \brief converts from the internal cp2k units to the given unit
!> \param value the value to convert
!> \param unit the unit of the result
!> \param defaults the defaults unit for those that are left free
!>        (cp_units_none)
!> \param power the power of the unit (defaults to 1)
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_unit_from_cp2k1(value, unit, defaults, power) RESULT(res)
      REAL(kind=dp), INTENT(in)                          :: value
      TYPE(cp_unit_type), INTENT(IN)                     :: unit
      TYPE(cp_unit_set_type), INTENT(IN), OPTIONAL       :: defaults
      INTEGER, INTENT(in), OPTIONAL                      :: power
      REAL(kind=dp)                                      :: res
 
      INTEGER                                            :: my_power
 
      my_power = 1
      IF (PRESENT(power)) my_power = power
      IF (PRESENT(defaults)) THEN
         res = cp_unit_to_cp2k1(value=value, unit=unit, defaults=defaults, &
                                power=-my_power)
      ELSE
         res = cp_unit_to_cp2k1(value=value, unit=unit, power=-my_power)
      END IF
   END FUNCTION cp_unit_from_cp2k1
 
! **************************************************************************************************
!> \brief converts to the internal cp2k units to the given unit
!> \param value the value to convert
!> \param unit_str the unit of the result as string
!> \param defaults the defaults unit for those that are left free
!>        (cp_units_none)
!> \param power the power of the unit (defaults to 1)
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_unit_to_cp2k(value, unit_str, defaults, power) RESULT(res)
      REAL(kind=dp), INTENT(in)                          :: value
      CHARACTER(len=*), INTENT(in)                       :: unit_str
      TYPE(cp_unit_set_type), INTENT(IN), OPTIONAL       :: defaults
      INTEGER, INTENT(in), OPTIONAL                      :: power
      REAL(kind=dp)                                      :: res
 
      TYPE(cp_unit_type)                                 :: my_unit
 
      CALL cp_unit_create(my_unit, unit_str)
      IF (PRESENT(defaults)) THEN
         res = cp_unit_to_cp2k1(value=value, unit=my_unit, defaults=defaults, &
                                power=power)
      ELSE
         res = cp_unit_to_cp2k1(value=value, unit=my_unit, power=power)
      END IF
      CALL cp_unit_release(my_unit)
   END FUNCTION cp_unit_to_cp2k
 
! **************************************************************************************************
!> \brief converts from the internal cp2k units to the given unit
!> \param value the value to convert
!> \param unit_str the unit of the result as string
!> \param defaults the defaults unit for those that are left free
!>        (cp_units_none)
!> \param power the power of the unit (defaults to 1)
!> \return ...
!> \author fawzi
! **************************************************************************************************
   FUNCTION cp_unit_from_cp2k(value, unit_str, defaults, power) RESULT(res)
      REAL(kind=dp), INTENT(in)                          :: value
      CHARACTER(len=*), INTENT(in)                       :: unit_str
      TYPE(cp_unit_set_type), INTENT(IN), OPTIONAL       :: defaults
      INTEGER, INTENT(in), OPTIONAL                      :: power
      REAL(kind=dp)                                      :: res
 
      TYPE(cp_unit_type)                                 :: my_unit
 
      CALL cp_unit_create(my_unit, unit_str)
      IF (PRESENT(defaults)) THEN
         res = cp_unit_from_cp2k1(value=value, unit=my_unit, defaults=defaults, &
                                  power=power)
      ELSE
         res = cp_unit_from_cp2k1(value=value, unit=my_unit, power=power)
      END IF
      CALL cp_unit_release(my_unit)
   END FUNCTION cp_unit_from_cp2k
 
! **************************************************************************************************
!> \brief returs true if the two units are compatible
!> \param ref_unit ...
!> \param unit ...
!> \return ...
!> \author Teodoro Laino [tlaino] - 11.2007 - University of Zurich
! **************************************************************************************************
   FUNCTION cp_unit_compatible(ref_unit, unit) RESULT(res)
      TYPE(cp_unit_type), INTENT(IN)                     :: ref_unit, unit
      LOGICAL                                            :: res
 
      INTEGER                                            :: i
 
      res = .true.
      DO i = 1, SIZE(ref_unit%kind_id)
         IF (ref_unit%kind_id(i) == unit%kind_id(i)) cycle
         IF ((ref_unit%kind_id(1) == cp_ukind_undef) .AND. (all(ref_unit%kind_id(2:) == cp_ukind_none))) cycle
         res = .false.
         EXIT
      END DO
 
   END FUNCTION cp_unit_compatible
 
! **************************************************************************************************
!> \brief initializes the given unit set
!> \param unit_set the set to initialize
!> \param name the name of the set, used for the dafault initialization of
!>        the various units
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE cp_unit_set_create(unit_set, name)
      TYPE(cp_unit_set_type), INTENT(OUT)                :: unit_set
      CHARACTER(len=*), INTENT(in)                       :: name
 
      CHARACTER(len=cp_unit_desc_length)                 :: my_name
      INTEGER                                            :: i
 
      my_name = name
      CALL uppercase(my_name)
 
      DO i = 1, cp_ukind_max
         NULLIFY (unit_set%units(i)%unit)
         ALLOCATE (unit_set%units(i)%unit)
      END DO
      DO i = 1, cp_ukind_max
         SELECT CASE (name)
         CASE ('ATOM', 'ATOMIC', 'INTERNAL', 'CP2K')
            IF (i == cp_ukind_angle) THEN
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], &
                                    unit_id=[cp_units_rad], power=[1])
            ELSE
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], &
                                    unit_id=[cp_units_au], power=[1])
            END IF
         CASE ('OUTPUT')
            SELECT CASE (i)
            CASE (cp_ukind_undef)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_none], &
                                    power=[1])
            CASE (cp_ukind_energy)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_hartree], &
                                    power=[1])
            CASE (cp_ukind_length)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_angstrom], &
                                    power=[1])
            CASE (cp_ukind_temperature)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_k], &
                                    power=[1])
            CASE (cp_ukind_angle)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_deg], &
                                    power=[1])
            CASE (cp_ukind_pressure)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_bar], &
                                    power=[1])
            CASE (cp_ukind_time)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_fs], &
                                    power=[1])
            CASE (cp_ukind_mass)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_amu], &
                                    power=[1])
            CASE (cp_ukind_potential)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_volt], &
                                    power=[1])
            CASE (cp_ukind_force)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_newton], &
                                    power=[1])
            CASE (cp_ukind_efield)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_volt_per_m], &
                                    power=[1])
            CASE (cp_ukind_dipole)
               CALL cp_unit_create2(unit_set%units(i)%unit, kind_id=[i], unit_id=[cp_units_debye], &
                                    power=[1])
            CASE default
               cpabort("unhandled unit type "//trim(cp_to_string(i)))
               EXIT
            END SELECT
         CASE default
            cpabort('unknown parameter set name '//trim(name))
         END SELECT
      END DO
   END SUBROUTINE cp_unit_set_create
 
! **************************************************************************************************
!> \brief releases the given unit set
!> \param unit_set the unit set to release
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE cp_unit_set_release(unit_set)
      TYPE(cp_unit_set_type), INTENT(INOUT)              :: unit_set
 
      INTEGER                                            :: i
 
      DO i = 1, SIZE(unit_set%units)
         CALL cp_unit_release(unit_set%units(i)%unit)
         DEALLOCATE (unit_set%units(i)%unit)
      END DO
 
   END SUBROUTINE cp_unit_set_release
 
! **************************************************************************************************
!> \brief Exports all available units as XML.
!> \param iw ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE export_units_as_xml(iw)
      INTEGER, INTENT(IN)                                :: iw
 
      CALL format_units_as_xml("energy", s2a("hartree", "wavenumber_e", "joule", "kcalmol", &
                                             "kjmol", "Ry", "eV", "keV", "K_e"), iw)
      CALL format_units_as_xml("length", s2a("bohr", "m", "pm", "nm", "angstrom"), iw)
      CALL format_units_as_xml("temperature", s2a("K", "au_temp"), iw)
      CALL format_units_as_xml("pressure", s2a("bar", "atm", "kbar", "Pa", "MPa", "GPa", "au_p"), iw)
      CALL format_units_as_xml("angle", s2a("rad", "deg"), iw)
      CALL format_units_as_xml("time", s2a("s", "fs", "ps", "au_t", "wavenumber_t"), iw)
      CALL format_units_as_xml("mass", s2a("kg", "amu", "m_e"), iw)
      CALL format_units_as_xml("potential", s2a("volt", "au_pot"), iw)
      CALL format_units_as_xml("force", s2a("N", "Newton", "mN", "mNewton", "au_f"), iw)
      CALL format_units_as_xml("efield", s2a("Vm-1", "Vnm-1", "VA-1", "volt_per_m", "volt_per_nm", &
                                             "volt_per_angstrom", "au_efield"), iw)
      CALL format_units_as_xml("dipole", s2a("debye", "e_angstrom", "C_m", "au_dip"), iw)
 
   END SUBROUTINE export_units_as_xml
 
! **************************************************************************************************
!> \brief Format units as xml.
!> \param unit_kind ...
!> \param units_set ...
!> \param iw ...
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE format_units_as_xml(unit_kind, units_set, iw)
      CHARACTER(LEN=*), INTENT(IN)                       :: unit_kind
      CHARACTER(LEN=*), DIMENSION(:), INTENT(IN)         :: units_set
      INTEGER, INTENT(IN)                                :: iw
 
      INTEGER                                            :: i
 
      WRITE (iw, fmt='(T2,A)') '<UNIT_KIND name="'//trim(unit_kind)//'">'
      DO i = 1, SIZE(units_set)
         WRITE (iw, fmt='(T3,A)') '<UNIT>'//trim(units_set(i))//'</UNIT>'
      END DO
      WRITE (iw, fmt='(T3,A)') '<UNIT>'//trim(unit_kind)//'</UNIT>' ! internal unit
      WRITE (iw, fmt='(T2,A)') '</UNIT_KIND>'
   END SUBROUTINE format_units_as_xml
 
END MODULE cp_units