!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!     % SUBROUTINE CTRANS(...)
!     %
!     % M. Vossel 21.03.2023
!     %
!     % Routine to transform symmetryinput coordinates  to symmetrized
!     % coordinates. Distances Are discribet by Morse coordinates or
!     % TMC depending on Set Parameters in the Genetic Input.
!     %
!     % input variables
!     % q:
!     %       q(1):  H1x
!     %       q(2):     y
!     %       q(3):      z
!     %       q(4):  H2x
!     %       q(5):     y
!     %       q(6):      z
!     %       q(7): H3x
!     %       q(8):    y
!     %       q(9):     z
!
!
!
!     % Internal variables:
!     % t:    primitive coordinates (double[qn])
!     %       t(1):
!     %       t(2):
!     %       t(3):
!     %       t(4):
!     %       t(5):
!     %       t(6):
!     %       t(7):
!     %       t(8):
!     %       t(9):
!     % t:    dummy (double[qn])
!     % p:    parameter vector
!     % npar: length of parameter vector
!     %
!     % Output variables
!     % s: symmetrized coordinates (double[qn])
!     %    s(1): CH-symetric streatch
!     %    s(2): CH-asymetric streatch-ex
!     %    s(3): CH-asymetric streatch-ey
!     %    s(4): CH-bend-ex
!     %    s(5): CH-bend-ey
!     %    s(6): CH-umbrella
!     %    s(7): CH-umbrella**2
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
module ctrans_mod
   use accuracy_constants, only: dp, idp
   implicit none
!     precalculate  pi, 2*pi and angle to radian conversion
   real(dp), parameter :: pi = 4.0_dp*datan(1.0_dp)
   real(dp), parameter :: pi2 = 2.0_dp*pi
   real(dp), parameter :: ang2rad = pi/180.0_dp
!     precalculate roots
   real(dp), parameter:: sq2 = 1.0_dp/dsqrt(2.0_dp)
   real(dp), parameter:: sq3 = 1.0_dp/dsqrt(3.0_dp)
   real(dp), parameter:: sq6 = 1.0_dp/dsqrt(6.0_dp)
!     change distances for equilibrium
   !real(dp), parameter :: dchequi = 1.02289024_dp
   real(dp), parameter :: dchequi = 2.344419_dp ! NO3
   !real(dp), parameter :: dchequi = 2.34451900_dp
  
! see changes
   
contains
   subroutine ctrans(q, x1,x2, invariants)
      use dim_parameter, only: qn
      integer(idp) k                 !running indices
      real(dp), intent(in) :: q(qn)    !given coordinates
      real(dp), intent(out) :: x1(qn)    !output coordinates symmetry adapted and scaled
      real(dp), intent(out) :: x2(qn)    !output coordinates symmetry adapted but not scaled
!     ANN Variables
      real(dp), optional, intent(out) :: invariants(:)
      real(dp) :: s(qn),t(qn)
      !     kartesian coordianates copy from MeF+ so substitute c by n and removed f
      real(dp) ch1(3), ch2(3), ch3(3), c_atom(3)
      real(dp) nh1(3), nh2(3), nh3(3)
      real(dp) zaxis(3), xaxis(3), yaxis(3)
      real(dp) ph1(3), ph2(3), ph3(3)
!     primitive coordinates
      real(dp) dch1, dch2, dch3 !nh-distances
      real(dp) umb      !Umbrella Angle from xy-plane

!     Symmetry coordinates
      real(dp) aR !a1-modes H-Dist.,
      real(dp) exR, exAng !ex components  H-Dist., H-Ang.
      real(dp) eyR, eyAng !ey components  H-Dist., H-Ang.
!     debugging
      logical, parameter ::  dbg = .false.

!      initialize coordinate vectors
      s = 0.0_dp
      t = 0.0_dp

!     write kartesian coords for readability
      c_atom(1:3) = q(1:3)
      do k = 1, 3
         ch1(k) = q(k + 3)
         ch2(k) = q(k + 6)
         ch3(k) = q(k + 9)
      end do

!     construct z-axis
      nh1 = normalized(ch1)
      nh2 = normalized(ch2)
      nh3 = normalized(ch3)
      zaxis = create_plane(nh1, nh2, nh3)

      !     calculate bonding distance
      dch1 = norm(ch1)
      dch2 = norm(ch2)
      dch3 = norm(ch3)

      !     construct symmertic and antisymmetric strech
      aR = symmetrize(dch1 - dchequi, dch2 - dchequi, dch3 - dchequi, 'a')
      exR = symmetrize(dch1, dch2, dch3, 'x')
      eyR = symmetrize(dch1, dch2, dch3, 'y')

      !     construc x-axis and y axis
      ph1 = normalized(project_point_into_plane(nh1, zaxis, c_atom))
      xaxis = normalized(ph1)
      yaxis = xproduct(zaxis, xaxis) ! right hand side koordinates

!     project H atoms into C plane
      ph2 = normalized(project_point_into_plane(nh2, zaxis, c_atom))
      ph3 = normalized(project_point_into_plane(nh3, zaxis, c_atom))

      call construct_HBend(exAng, eyAng, ph1, ph2, ph3, xaxis, yaxis)
      umb = construct_umbrella(nh1, nh2, nh3, zaxis)

!     set symmetry coordinates and even powers of umbrella
      !s(1) = dch1- dchequi !aR
      !s(2) = dch2 - dchequi !exR
      !s(3) = dch3 - dchequi !eyR
      s(1) = aR
      s(2) = exR
      s(3) = eyR
      s(4) = exAng
      s(5) = eyAng
      s(6) = umb
      s(7) = umb**2
      s(8) = 0
      s(9) = 0
!     pairwise distances as second coordinate set
      t = 0._dp
      call pair_distance(q, t(1:6))

      if (dbg) write (6, '("sym coords s=",9f16.8)') s(1:qn)
      if (dbg) write (6, '("sym coords t=",9f16.8)') t(1:qn)
      if (present(invariants)) then
         call get_invariants(s, invariants)
      end if

      ! transform s and t to x1 and x2
      x1(1:qn)=s(1:qn)
      x1(5)=x1(5)
      ! set other  x coordinate to zero other than strech
      !X1(4:qn)=0.0d0
      x2(1:qn)=t(1:qn)
   end subroutine ctrans

   subroutine pair_distance(q, r)
      real(dp), intent(in)  :: q(9)
      real(dp), intent(out) :: r(6)
      real(dp) :: atom(3, 4)
      integer :: n, k, count

      !atom order: H1 H2 H3 N
      atom(:, 1:3) = reshape(q, [3, 3])
      atom(:, 4) = (/0.0_dp, 0.0_dp, 0.0_dp/)

      ! disntace order 12 13 14 23 24 34
      count = 0
      do n = 1, size(atom, 2)
         do k = n + 1, size(atom, 2)
            count = count + 1
            r(count) = sqrt(sum((atom(:, k) - atom(:, n))**2))
         end do
      end do
   end subroutine pair_distance

   function morse_and_symmetrize(x,p,pst) result(s)
      real(dp), intent(in),dimension(3) :: x
      real(dp), intent(in),dimension(11) :: p
      integer, intent(in),dimension(2) :: pst
      integer :: k
      real(dp), dimension(3) :: s
      real(dp), dimension(3) :: t

      !     Morse transform
      do k=1,3
         t(k) = morse_transform(x(k), p, pst)
      end do   
      s(1) = symmetrize(t(1), t(2), t(3), 'a')
      s(2) = symmetrize(t(1), t(2), t(3), 'x')
      s(3) = symmetrize(t(1), t(2), t(3), 'y')
   end function morse_and_symmetrize

   subroutine get_invariants(s, inv_out)
      use dim_parameter, only: qn
      use select_monom_mod, only: v_e_monom, v_ee_monom
      real(dp), intent(in) :: s(qn)
      real(dp), intent(out) :: inv_out(:)
      !  real(dp), parameter ::  ck = 1.0_dp, dk = 1.0_dp/ck ! scaling for higher order invariants
      real(dp) inv(24)
      integer, parameter :: inv_order(12) = & ! the order in which the invariants are selected
      & [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
      real(dp) Rch, umb, xR, yR, xAng, yAng
!     for readability
      Rch = s(1)
      xR = s(2)
      yR = s(3)
      xAng = s(4)
      yAng = s(5)
      umb = s(6)**2
!     invarianten
!     a moden
      inv(1) = Rch
      inv(2) = umb
!     invariante e pairs
      inv(3) = v_e_monom(xR, yR, 1)
      inv(4) = v_e_monom(xAng, yAng, 1)
!     third order e pairs
      inv(5) = v_e_monom(xR, yR, 2)
      inv(6) = v_e_monom(xAng, yAng, 2)
      !  invariant ee coupling
      inv(7) = v_ee_monom(xR, yR, xAng, yAng, 1)
      !     mode combinations
      inv(8) = Rch*umb

      inv(9) = Rch*v_e_monom(xR, yR, 1)
      inv(10) = umb*v_e_monom(xR, yR, 1)

      inv(11) = Rch*v_e_monom(xAng, yAng, 1)
      inv(12) = umb*v_e_monom(xAng, yAng, 1)

!     damp coordinates because of second order and higher invariants
      inv(3) = sign(sqrt(abs(inv(3))), inv(3))
      inv(4) = sign(sqrt(abs(inv(4))), inv(4))
      inv(5) = sign((abs(inv(5))**(1./3.)), inv(5))
      inv(6) = sign((abs(inv(6))**(1./3.)), inv(6))
      inv(7) = sign((abs(inv(7))**(1./3.)), inv(7))
      inv(8) = sign(sqrt(abs(inv(8))), inv(8))
      inv(9) = sign((abs(inv(9))**(1./3.)), inv(9))
      inv(10) = sign((abs(inv(10))**(1./3.)), inv(10))
      inv(11) = sign((abs(inv(11))**(1./3.)), inv(11))
      inv(12) = sign((abs(inv(12))**(1./3.)), inv(12))

      inv_out(:) = inv(inv_order(1:size(inv_out, 1)))

   end subroutine get_invariants

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!     % real part of spherical harmonics
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!     Ylm shifted to 0 for theta=0
   real(dp) function ylm(theta, phi, l, m)
      implicit none
      real(dp) theta, phi
      integer(idp) l, m
      ylm = plm2(dcos(theta), l, m)*cos(m*phi) - plm2(1.0_dp, l, m)
   end function ylm
!----------------------------------------------------------
   real(dp) function plm2(x, l, n)
      implicit none
      real(dp) x
      integer(idp) l, m, n

      real(dp) pmm, p_mp1m, pllm
      integer(idp) ll

!     negative m und bereich von x abfangen
      if ((l .lt. 0)&
      &.or. (abs(n) .gt. abs(l))&
      &.or. (abs(x) .gt. 1.)) then
         write (6, '(''bad arguments in legendre'')')
         stop
      end if

!     fix sign of m to compute the positiv m
      m = abs(n)

      pmm = (-1)**m*dsqrt(fac(2*m))*1./((2**m)*fac(m))&    !compute P(m,m) not P(l,l)
           &*(dsqrt(1.-x**2))**m

      if (l .eq. m) then
         plm2 = pmm                                           !P(l,m)=P(m,m)
      else
         p_mp1m = x*dsqrt(dble(2*m + 1))*pmm                      !compute P(m+1,m)
         if (l .eq. m + 1) then
            plm2 = p_mp1m                                       !P(l,m)=P(m+1,m)
         else
            do ll = m + 2, l
               pllm = x*(2*l - 1)/dsqrt(dble(l**2 - m**2))*p_mp1m& ! compute P(m+2,m) up to P(l,m) recursively
               &- dsqrt(dble((l - 1)**2 - m**2))&
               &/dsqrt(dble(l**2 - m**2))*pmm
!    schreibe m+2 und m+1 jeweils fuer die naechste iteration
               pmm = p_mp1m                                     !P(m,m) = P(m+1,m)
               p_mp1m = pllm                                  !P(m+1,m) = P(m+2,m)
            end do
            plm2 = pllm                                       !P(l,m)=P(m+k,m), k element N
         end if
      end if

!     sets the phase of -m term right (ignored to gurantee Ylm=(Yl-m)* for JT terms
!      if(n.lt.0) then
!         plm2 = (-1)**m * plm2 !* fac(l-m)/fac(l+m)
!      endif

   end function
!----------------------------------------------------------------------------------------------------
   real(dp) function fac(i)
      integer(idp) i
      select case (i)
      case (0)
         fac = 1.0_dp
      case (1)
         fac = 1.0_dp
      case (2)
         fac = 2.0_dp
      case (3)
         fac = 6.0_dp
      case (4)
         fac = 24.0_dp
      case (5)
         fac = 120.0_dp
      case (6)
         fac = 720.0_dp
      case (7)
         fac = 5040.0_dp
      case (8)
         fac = 40320.0_dp
      case (9)
         fac = 362880.0_dp
      case (10)
         fac = 3628800.0_dp
      case (11)
         fac = 39916800.0_dp
      case (12)
         fac = 479001600.0_dp
      case default
         write (*, *) 'ERROR: no case for given faculty, Max is 12!'
         stop
      end select
   end function fac

   ! Does the simplest morse transform possible
   ! one skaling factor + shift
   function morse_transform(x, p, pst) result(t)
      real(dp), intent(in) :: x
      real(dp), intent(in) :: p(11)
      integer, intent(in) :: pst(2)
      real(dp) :: t
      if (pst(2) == 11) then
         t = 1.0_dp - exp(-abs(p(2))*(x - p(1)))
      else
         error stop 'in morse_transform key required or wrong number of parameters'
      end if
   end function morse_transform

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!     % FUNCTION F(...) ! MAIK DEPRICATING OVER THE TOP MORSE FUNCTION FOR MYSELF
!     %
!     % Returns exponent of tunable Morse coordinate
!     % exponent is polynomial * gaussian (skewed)
!     % ilabel = 1 or 2 selects the parameters a and sfac to be used
!     %
!     % Background: better representation of the prefector in the
!     %             exponend of the morse function.
!     % Formular: f(r) = lest no3 paper
!     %
!     % Variables:
!     % x:  distance of atoms (double)
!     % p:  parameter vector (double[20])
!     % ii: 1 for CCl and 2 for CCH (int)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   pure function f(x, p, ii)
      integer(idp), intent(in) :: ii                !1 for CCL and 2 for CCH
      real(dp), intent(in) :: x        !coordinate
      real(dp), intent(in) :: p(11)    !parameter-vector

      integer(idp) i                 !running index

      real(dp) r        !equilibrium distance
      real(dp) gaus     !gaus part of f
      real(dp) poly     !polynom part of f
      real(dp) skew     !tanh part of f

      real(dp) f        !prefactor of exponent and returned value

      integer(idp) npoly(2)          !order of polynom

!     Maximum polynom order
      npoly(1) = 5
      npoly(2) = 5

!     p(1):        position of equilibrium
!     p(2):        constant of exponent
!     p(3):        constant for skewing the gaussian
!     p(4):        tuning for skewing the gaussian
!     p(5):        Gaussian exponent
!     p(6):        Shift of Gaussian maximum
!     p(7)...:     polynomial coefficients
!     p(8+n)...:   coefficients of Morse Power series

!     1-exp{[p(2)+exp{-p(5)[x-p(6)]^2}[Taylor{p(7+n)}(x-p(6))]][x-p(1)]}

!     Tunable Morse function
!     Power series in Tunable Morse coordinates of order m
!     exponent is polynomial of order npoly * gaussian + switching function

!     set r  r-r_e
      r = x
      r = r - p(1)

!     set up skewing function:
      skew = 0.5_dp*p(3)*(dtanh(dabs(p(4))*(r - p(6))) + 1.0_dp)

!     set up gaussian function:
      gaus = dexp(-dabs(p(5))*(r - p(6))**2)

!     set up power series:
      poly = 0.0_dp
      do i = 0, npoly(ii) - 1
         poly = poly + p(7 + i)*(r - p(6))**i
      end do
!     set up full exponent function:
      f = dabs(p(2)) + skew + gaus*poly

   end function
!----------------------------------------------------------------------------------------------------
   pure function xproduct(a, b) result(axb)
      real(dp), intent(in) :: a(3), b(3)
      real(dp) ::  axb(3) !crossproduct a x b
      axb(1) = a(2)*b(3) - a(3)*b(2)
      axb(2) = a(3)*b(1) - a(1)*b(3)
      axb(3) = a(1)*b(2) - a(2)*b(1)
   end function xproduct

   pure function normalized(v) result(r)
      real(dp), intent(in) :: v(:)
      real(dp) :: r(size(v))
      r = v/norm(v)
   end function normalized

   pure function norm(v) result(n)
      real(dp), intent(in) :: v(:)
      real(dp) n
      n = dsqrt(sum(v(:)**2))
   end function norm

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!     % FUNCTION Project_Point_Into_Plane(x,n,r0) result(p)
!     % return the to n orthogonal part of a vector x-r0
!     %  p: projected point in plane
!     %  x: point being projected
!     %  n: normalvector of plane
!     % r0: Point in plane
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   pure function project_point_into_plane(x, n, r0) result(p)
      real(dp), intent(in) :: x(:), n(:), r0(:)
      real(dp) :: p(size(x)), xs(size(x))
      xs = x - r0
      p = xs - plane_to_point(x, n, r0)
   end function project_point_into_plane

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!     % Function Plane_To_Point(x,n,r0) result(p)
!     %  p: part of n in x
!     %  x: point being projected
!     %  n: normalvector of plane
!     % r0: Point in plane
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   pure function plane_to_point(x, n, r0) result(p)
      real(dp), intent(in) :: x(:), n(:), r0(:)
      real(dp) p(size(x)), xs(size(x)), nn(size(n))
      nn = normalized(n)
      xs = x - r0
      p = dot_product(nn, xs)*nn
   end function plane_to_point

   subroutine check_coordinates(q)
!     check for faulty kartesain coordinates
      real(dp), intent(in) :: q(:)
      integer(idp) :: i
      if (all(abs(q) <= epsilon(0.0_dp))) then
         stop 'Error (ctrans): all kartesian coordinates are<=1d-8'
      end if
      do i = 1, 9, 3
         if (all(abs(q(i:i + 2)) <= epsilon(0.0_dp))) then
            write (*, *) q
            stop 'Error(ctrans):kartesian coordinates zero for one atom'
         end if
      end do
   end subroutine

   pure function rotor_a_to_z(a, z) result(r)
      real(dp), intent(in) :: a(3), z(3)
      real(dp) :: r(3, 3)
      real(dp) :: alpha
      real(dp) :: s1(3), s(3, 3), rotor(3, 3)
      s1 = xproduct(normalized(a), normalized(z))
      alpha = asin(norm(s1))
      s(:, 1) = normalized(s1)
      s(:, 2) = normalized(z)
      s(:, 3) = xproduct(s1, z)
      rotor = init_rotor(alpha, 0.0_dp, 0.0_dp)
      r = matmul(s, matmul(rotor, transpose(s)))
   end function

!     function returning Rz(gamma) * Ry(beta) * Rx(alpha) for basis order xyz
   pure function init_rotor(alpha, beta, gamma) result(rotor)
      real(dp), intent(in) :: alpha, beta, gamma
      real(dp) :: rotor(3, 3)
      rotor = 0.0_dp
      rotor(1, 1) = dcos(beta)*dcos(gamma)
      rotor(1, 2) = dsin(alpha)*dsin(beta)*dcos(gamma)&
      &- dcos(alpha)*dsin(gamma)
      rotor(1, 3) = dcos(alpha)*dsin(beta)*dcos(gamma)&
      &+ dsin(alpha)*dsin(gamma)

      rotor(2, 1) = dcos(beta)*dsin(gamma)
      rotor(2, 2) = dsin(alpha)*dsin(beta)*dsin(gamma)&
      &+ dcos(alpha)*dcos(gamma)
      rotor(2, 3) = dcos(alpha)*dsin(beta)*dsin(gamma)&
      &- dsin(alpha)*dcos(gamma)

      rotor(3, 1) = -dsin(beta)
      rotor(3, 2) = dsin(alpha)*dcos(beta)
      rotor(3, 3) = dcos(alpha)*dcos(beta)
   end function init_rotor

   pure function create_plane(a, b, c) result(n)
      real(dp), intent(in) :: a(3), b(3), c(3)
      real(dp) :: n(3)
      real(dp) :: axb(3), bxc(3), cxa(3)
      axb = xproduct(a, b)
      bxc = xproduct(b, c)
      cxa = xproduct(c, a)
      n = normalized(axb + bxc + cxa)
   end function create_plane

   function symmetrize(q1, q2, q3, sym) result(s)
      real(dp), intent(in) :: q1, q2, q3
      character, intent(in) :: sym
      real(dp) :: s
      select case (sym)
      case ('a')
         s = (q1 + q2 + q3)*sq3
      case ('x')
         s = sq6*(2.0_dp*q1 - q2 - q3)
      case ('y')
         s = sq2*(q2 - q3)
      case default
         write (*, *) 'ERROR: no rule for symmetrize with sym=', sym
         stop
      end select
   end function symmetrize

   subroutine construct_HBend(ex, ey, ph1, ph2, ph3, x_axis, y_axis)
      real(dp), intent(in) :: ph1(3), ph2(3), ph3(3)
      real(dp), intent(in) :: x_axis(3), y_axis(3)
      real(dp), intent(out) :: ex, ey
      real(dp) :: x1, y1, alpha1
      real(dp) :: x2, y2, alpha2
      real(dp) :: x3, y3, alpha3
!     get x and y components of projected points
      x1 = dot_product(ph1, x_axis)
      y1 = dot_product(ph1, y_axis)
      x2 = dot_product(ph2, x_axis)
      y2 = dot_product(ph2, y_axis)
      x3 = dot_product(ph3, x_axis)
      y3 = dot_product(ph3, y_axis)
!     -> calculate H deformation angles
      alpha3 = datan2(y2, x2)
      alpha2 = -datan2(y3, x3)     !-120*ang2rad
!      write(*,*)' atan2'
!      write(*,*) 'alpha2:' , alpha2/ang2rad
!      write(*,*) 'alpha3:' , alpha3/ang2rad
      if (alpha2 .lt. 0) alpha2 = alpha2 + pi2
      if (alpha3 .lt. 0) alpha3 = alpha3 + pi2
      alpha1 = (pi2 - alpha2 - alpha3)
!      write(*,*)' fixed break line'
!      write(*,*) 'alpha1:' , alpha1/ang2rad
!      write(*,*) 'alpha2:' , alpha2/ang2rad
!      write(*,*) 'alpha3:' , alpha3/ang2rad
      alpha1 = alpha1 !- 120.0_dp*ang2rad
      alpha2 = alpha2 !- 120.0_dp*ang2rad
      alpha3 = alpha3 !- 120.0_dp*ang2rad
!      write(*,*)' delta alpha'
!      write(*,*) 'alpha1:' , alpha1/ang2rad
!      write(*,*) 'alpha2:' , alpha2/ang2rad
!      write(*,*) 'alpha3:' , alpha3/ang2rad
!      write(*,*)

!     construct symmetric and antisymmetric H angles
      ex = symmetrize(alpha1, alpha2, alpha3, 'x')
      ey = symmetrize(alpha1, alpha2, alpha3, 'y')
   end subroutine construct_HBend

   pure function construct_umbrella(nh1, nh2, nh3, n)&
   &result(umb)
      real(dp), intent(in) :: nh1(3), nh2(3), nh3(3)
      real(dp), intent(in) :: n(3)
      real(dp) :: umb
      real(dp) :: theta(3)
!     calculate projections for umberella angle
      theta(1) = dacos(dot_product(n, nh1))
      theta(2) = dacos(dot_product(n, nh2))
      theta(3) = dacos(dot_product(n, nh3))
!     construct umberella angle
      umb = sum(theta(1:3))/3.0_dp - 90.0_dp*ang2rad
   end function construct_umbrella

   pure subroutine construct_sphericals&
   &(theta, phi, cf, xaxis, yaxis, zaxis)
      real(dp), intent(in) :: cf(3), xaxis(3), yaxis(3), zaxis(3)
      real(dp), intent(out) :: theta, phi
      real(dp) :: x, y, z, v(3)
      v = normalized(cf)
      x = dot_product(v, normalized(xaxis))
      y = dot_product(v, normalized(yaxis))
      z = dot_product(v, normalized(zaxis))
      theta = dacos(z)
      phi = -datan2(y, x)
   end subroutine construct_sphericals

   subroutine int2kart(internal, kart)
      real(dp), intent(in) :: internal(6)
      real(dp), intent(out) :: kart(9)
      real(dp) :: h1x, h1y, h1z
      real(dp) :: h2x, h2y, h2z
      real(dp) :: h3x, h3y, h3z
      real(dp) :: dch0, dch1, dch2, dch3
      real(dp) :: a1, a2, a3, wci

      kart = 0.0_dp
      dch1 = dchequi + sq3*internal(1) + 2*sq6*internal(2)
      dch2 = dchequi + sq3*internal(1) - sq6*internal(2) + sq2*internal(3)
      dch3 = dchequi + sq3*internal(1) - sq6*internal(2) - sq2*internal(3)
      a1 = 2*sq6*internal(4)
      a2 = -sq6*internal(4) + sq2*internal(5)
      a3 = -sq6*internal(4) - sq2*internal(5)
      wci = internal(6)

      ! Berechnung kartesische Koordinaten
      ! -----------------------
      h1x = dch1*cos(wci*ang2rad)
      h1y = 0.0
      h1z = -dch1*sin(wci*ang2rad)

      h3x = dch2*cos((a2 + 120)*ang2rad)*cos(wci*ang2rad)
      h3y = dch2*sin((a2 + 120)*ang2rad)*cos(wci*ang2rad)
      h3z = -dch2*sin(wci*ang2rad)

      h2x = dch3*cos((-a3 - 120)*ang2rad)*cos(wci*ang2rad)
      h2y = dch3*sin((-a3 - 120)*ang2rad)*cos(wci*ang2rad)
      h2z = -dch3*sin(wci*ang2rad)

      kart(1) = h1x
      kart(2) = h1y
      kart(3) = h1z
      kart(4) = h2x
      kart(5) = h2y
      kart(6) = h2z
      kart(7) = h3x
      kart(8) = h3y
      kart(9) = h3z
   end subroutine int2kart

end module ctrans_mod