C     Fortran-90 program to convert Cartesian cooordinates of two water
C     molecules to two sets of Euler angles and an intermolecular
C     distance.
C
C     Program is for use with SAPT5s.f and assumes, just like SAPT5s.f,
C     that the standard water is in the xz-plane with the oxygen on the
C     positive z-axis and the center of mass in the origin.
C
C
C     Input:
C       - namelist/input/   (from standard I/O):
C           coord_file  = name of file containing Cartesian
C                         coordinates of water dimer (char*60).
C                         Default: 'coordinates.dat'
C
C           output_file = name of file containing output (is input file
C                         to sapt5s.f). Default: 'example.dat'
C
C           what_units  = length unit (string: 'angstrom' or 'bohr')
C                         Default: bohr
C
C Example namelist: (do not copy the exclamation marks)
C------------------------------------------------------
! &input coord_file  = 'flop.dat',
!        output_file = 'sample.dat',
!        what_units  = 'angstrom'     /
C-------------------------------------------------------
C
C       - coord_file contains free format lines:
C              atom, x, y, z
C                  ...
C              atom, x, y, z
C
C Example coord_file ('flop.dat'). Do not copy exclamation marks.
C --------------------------------------------------------------
!* Coordinates obtained from F. van Duijneveldt:
!* Monomer 1
! O1         0.873683    0.000000   -1.246754
! h1         0.287800    0.000000   -2.003704
! H2         0.287800    0.000000   -0.489804
!* Monomer 2
! o2        -0.778985    0.000000    1.131754
! H3        -0.666591    0.756950    1.706755
! h4        -0.666591   -0.756950    1.706754
C-----------------------------------------------------------------
C
C     Notes on structure of coord_file (ASCII)
C     ----------------------------------------
C     - Lines on coord_file beginning with '*' are skipped (allows for
C       comments).
C     - Empty lines are skipped
C     - Atom is a string starting with O, H, or D.
C     - First  three non-comment lines --> first  monomer
C     - Second three non-comment lines --> second monomer
C     - Specification of monomers must start with oxygen atoms.
C
C     Output:
C     - Monomer geometries, Euler angles, distance on standard I/O.
C     - Distance, Euler angles on file with name 'output_file'
C
C     Author: P.E.S. Wormer, February 2001.

      MODULE constants
      implicit      none
      character*60, public            ::  coord_file
      character*60, public            ::  output_file
      integer*4,    public, parameter ::  maxatoms = 6

      real*8, public            :: units
      real*8, public, parameter :: pi        =  3.1415926535897932d0
      real*8, public, parameter :: angle_HOH =104.69d0
      real*8, public, parameter :: OH_length =  1.836106337d0 ! bohr
      real*8, public, parameter :: toler_l   =  OH_length/10.d0
      real*8, public, parameter :: toler_a   =  angle_HOH/10.d0
      real*8, public, parameter :: bohr      =  0.5291772083d0
      real*8, public, parameter :: m_H       =  1.007825d0
      real*8, public, parameter :: m_D       =  2.0140d0
      real*8, public, parameter :: m_O       = 15.99491d0
      real*8, public, parameter :: smallest  = tiny(units)
      real*8, public, parameter :: mac_prec  = epsilon(units)

      end module constants

      MODULE vector_functions

      contains
      FUNCTION cart2polar(R)

C     From Cartesians to polars (radians). European definition of
C     polars.
C     Input 3-vector R.
C     Output vector (/ Length_R, theta, phi /).

      use constants, only: pi, smallest, mac_prec
      implicit none
      real*8 st, ct, sp, cp,
     .       theta, phi,
     .       Length_R,  Length_Rp, Length_diff,
     .       cart2polar(3),
     .       R(3),
     .       diff(3)
      intent(in) R

      Length_R = sqrt(dot_product(R,R))
      if ( abs(Length_R) < smallest )  then
        theta = 0.d0
        phi   = 0.d0
        return
      endif

      theta = acos( R(3)/Length_R )

      Length_Rp = sqrt( dot_product( R(1:2),R(1:2) ) )
      if (abs(Length_Rp) < smallest ) then
          phi = 0.d0
      else
          phi = acos(R(1)/Length_Rp)
          if ( R(2) <  -smallest ) then
             phi = 2.d0*pi - phi
          endif
      endif
      cart2polar = (/Length_R, theta, phi/)

C--Check
*     st   = sin(theta)
*     ct   = cos(theta)
*     sp   = sin(phi)
*     cp   = cos(phi)
*
*     diff        =  R - Length_R*(/st*cp, st*sp, ct/)
*     Length_diff =  sqrt(dot_product(diff, diff))
*
*     if (Length_diff > mac_prec*1.d7 ) then
*        write(*,*) ' Error in cart2polar, length differences: ',
*    .       Length_diff
*        stop 16
*     endif

      END FUNCTION cart2polar
      FUNCTION euler_rotmat(alpha, beta, gamma)

C     Euler rotation matrix  R(alpha, beta, gamma); active convention.

      implicit none
      real*8   alpha, beta, gamma,
     .         euler_rotmat(3,3),
     .         Z1(3,3), Y(3,3), Z(3,3)
      intent(in) alpha, beta, gamma

      Z = 0.d0; Y = 0.d0; Z1 = 0.d0;

      Z(1,1) =  cos(gamma)
      Z(2,1) =  sin(gamma)
      Z(1,2) = -sin(gamma)
      Z(2,2) =  cos(gamma)
      Z(3,3) =  1.d0

      Y(1,1) =  cos(beta)
      Y(3,1) = -sin(beta)
      Y(2,2) =  1.d0
      Y(1,3) =  sin(beta)
      Y(3,3) =  cos(beta)

      Z1(1,1) =  cos(alpha)
      Z1(2,1) =  sin(alpha)
      Z1(1,2) = -sin(alpha)
      Z1(2,2) =  cos(alpha)
      Z1(3,3) =  1.d0

      euler_rotmat = matmul(matmul(Z1,Y),Z)

      end function euler_rotmat
      FUNCTION cross(a,b)

C     Cross product a x b

      implicit none
      real*8 a(3), b(3), cross(3)
      intent(in)  a, b

      cross(1) = a(2)*b(3) - a(3)*b(2)
      cross(2) = a(3)*b(1) - a(1)*b(3)
      cross(3) = a(1)*b(2) - a(2)*b(1)

      end function cross

      end module vector_functions

      PROGRAM CART2EULER

      use constants, only: maxatoms
      implicit    none
      real*8      coords(3, maxatoms)
      character*4 atoms (maxatoms)
      integer     natoms

      Call rdstdio                         ! namelist from standard i/o
      Call getcoords(coords, atoms, natoms)! coordinates from coord_file

C--Overall rotation of system such that RAB coincides with lab-
C     fixed z-axis. Modify array coords. Print R_AB on standard I/O.
      Call  rot_overall(coords, atoms, natoms)

C--Get monomer geometries. Print Euler angles of monomers on standard
C     I/O.
      Call get_monomers(coords)

      END
      SUBROUTINE RDSTDIO

C     Read namelist from standard I/O

      use constants, only: coord_file, output_file, units, bohr
      implicit none
      character*8   what_units
      namelist/input/coord_file, output_file, what_units

C--Namelist defaults
      coord_file  = 'coordinates.dat'
      output_file = 'example.dat'
      what_units  = 'bohr'

      read(*,input)
      write(*,'(t10, a/t10,a/ 3(2a/) )')
     .  ' Namelist input ',
     .  ' ============== ',
     .  ' coord_file  =  ', trim(adjustl(coord_file)),
     .  ' output_file =  ', trim(adjustl(output_file)),
     .  ' what_units  =  ', what_units

C--Convert to bohr
      if     (what_units(1:1) .eq. 'a' .or. what_units(1:1) .eq. 'A' )
     .   then
          units = 1.d0/bohr
      elseif (what_units(1:1) .eq. 'b' .or. what_units(1:1) .eq. 'B' )
     .   then
          units = 1.d0
      else
          write(*,'(/a)')
     .    ' Wrong length units, must be angstrom or bohr '
          stop 16
      endif

      END
      SUBROUTINE getcoords(coords, atoms, natoms)

C     Get Cartesian coordinates from  'coord_file'.
C     Return in array atom(char*4) and coords(3, natoms)

      use constants, only: maxatoms, coord_file, units
      implicit      none
      integer       natoms, i
      character*1,  parameter   ::  comment = '*'
      real*8        coords(3, maxatoms)
      character*4   atoms(maxatoms)
      character*100 line
      intent(out)   coords, atoms, natoms

      open(unit=1, file=coord_file,form='formatted')

      natoms = 0
   10 read(unit= 1,fmt='(a100)',end=20) line
      if ( (line(1:1) .eq. comment)  .or. (len(trim(line)) .eq. 0)) then
         go to 10
      else
         if (natoms .eq. maxatoms) then
             write(*,'(//a/a)') ' *** Array COORDS too small. ',
     .           ' Increase maxatoms and recompile.'
             stop 16
         endif
         natoms = natoms + 1
         read(line,*)  atoms(natoms), coords(:,natoms)
         go to 10
      endif

   20 write(*,'(t15,a/t15,a)')
     .    ' Input coordinates ',
     .    ' ================= '
      do i = 1,natoms
         write(*,'(1x, a, 3(f12.6,1x) )')
     .           atoms(i), coords(:,i)
      enddo

C--Check the atoms
      if    (atoms(1)(1:1) .ne. 'o' .and. atoms(1)(1:1) .ne. 'O' ) then
          write(*,*) ' First atom must be oxygen '
          stop 16
      endif

      if    (atoms(2)(1:1) .ne. 'h' .and. atoms(2)(1:1) .ne. 'H'
     . .and. atoms(2)(1:1) .ne. 'd' .and. atoms(2)(1:1) .ne. 'D' ) then
          write(*,*) ' Second atom must be proton or deuteron '
          stop 16
      endif

      if    (atoms(3)(1:1) .ne. 'h' .and. atoms(3)(1:1) .ne. 'H'
     . .and. atoms(3)(1:1) .ne. 'd' .and. atoms(3)(1:1) .ne. 'D' ) then
          write(*,*) ' Third atom must be proton or deuteron '
          stop 16
      endif

      if    (atoms(4)(1:1) .ne. 'o' .and. atoms(4)(1:1) .ne. 'O' ) then
          write(*,*) ' Fourth atom must be oxygen '
          stop 16
      endif

      if    (atoms(5)(1:1) .ne. 'h' .and. atoms(5)(1:1) .ne. 'H'
     . .and. atoms(5)(1:1) .ne. 'd' .and. atoms(5)(1:1) .ne. 'D' ) then
          write(*,*) ' Fifth  atom must be proton or deuteron '
          stop 16
      endif

      if    (atoms(6)(1:1) .ne. 'h' .and. atoms(6)(1:1) .ne. 'H'
     . .and. atoms(6)(1:1) .ne. 'd' .and. atoms(6)(1:1) .ne. 'D' ) then
          write(*,*) ' Sixth atom must be proton or deuteron '
          stop 16
      endif

C--To bohr:
      coords = coords*units

      END
      SUBROUTINE  rot_overall(coords, atoms, natoms)

C     Overall rotation of dimer

      use constants, only: maxatoms, m_H, m_D, m_O, bohr, mac_prec,
     .                     output_file
      use vector_functions, only: cart2polar, euler_rotmat
      implicit none
      integer i,  natoms
      real*8      m_tot_A,
     .            m_tot_B,
     .            coords (3, maxatoms),
     .            masses (maxatoms),
     .            com_A  (3),
     .            com_B  (3),
     .            RAB    (3),
     .            RABnew (3),
     .            polarAB(3),
     .            rotR   (3,3)
      character*4 atoms(maxatoms)
      intent(in)    atoms, natoms
      intent(inout) coords

C--Set array masses
      do i = 1, natoms
         if     ( (atoms(i)(1:1).eq.'h') .or. (atoms(i)(1:1).eq.'H'))
     .     then
           masses(i) = m_H
         elseif ( (atoms(i)(1:1).eq.'d') .or. (atoms(i)(1:1).eq.'D'))
     .     then
           masses(i) = m_D
         elseif ( (atoms(i)(1:1).eq.'o') .or. (atoms(i)(1:1).eq.'O'))
     .     then
           masses(i) = m_O
         endif
      enddo

C--Set centers of mass
      m_tot_A = sum(masses(1:3))           ! total mass of A
      m_tot_B = sum(masses(4:6))           ! total mass of B
      com_A   = matmul(coords(:,1:3), masses(1:3))/m_tot_A
      com_B   = matmul(coords(:,4:6), masses(4:6))/m_tot_B

C--Matrix rotR rotates RAB such that it lies along the z-axis
      RAB      = com_B - com_A        ! points from A to B
      polarAB  = cart2polar(RAB)      ! returns  R, theta, phi
      rotR     = transpose(euler_rotmat( polarAB(3), polarAB(2), 0.d0))
      RABnew   = matmul(rotR, RAB)    ! Rotate RAB
      coords   = matmul(rotR, coords) ! Rotate total system

C--Check
      if ( sqrt(dot_product(RABnew(1:2), RABnew(1:2))) > mac_prec*1.d3)
     .    then
          write(*,*) ' R not along z-axis, RABnew = ', RABnew
          stop 16
      endif

      write(*,'(/a,f14.8,a,2x,f14.8,a)')
     .      ' Intermolecular distance: ', polarAB(1), ' [bohr]',
     .        polarAB(1)*bohr, ' [Angstrom]'

      open(unit=2, file=output_file,form='formatted')
      write(2,'(i1)') 1
      write(2,*) polarAB(1)

      END

      SUBROUTINE get_monomers(coords)

C     Get the geometry and Euler angles of the two monomers

      use constants, only: maxatoms, pi, OH_length, angle_HOH,
     .                     toler_a, toler_l
      use vector_functions, only: cross
      implicit none
      real*8   l_OH1, l_OH2, l_OH3,  l_OH4,
     .         angle1, angle2,
     .         alphaA, betaA, gammaA,
     .         alphaB, betaB, gammaB,
     .         coords(3, maxatoms),
     .         OH1   (3),
     .         OH2   (3),
     .         OH3   (3),
     .         OH4   (3),
     .         monA  (3,3),
     .         monB  (3,3)
      integer  i
      intent(in) coords

C--Get the geometry of the monomers
      OH1   = (coords(:,2) - coords(:,1))
      OH2   = (coords(:,3) - coords(:,1))
      OH3   = (coords(:,5) - coords(:,4))
      OH4   = (coords(:,6) - coords(:,4))
      l_OH1 = sqrt(dot_product(OH1, OH1))
      l_OH2 = sqrt(dot_product(OH2, OH2))
      l_OH3 = sqrt(dot_product(OH3, OH3))
      l_OH4 = sqrt(dot_product(OH4, OH4))

      write(*,'(/a/a/a,a/ (a,f12.7,a,f12.7) )')
     .    '       Monomer Geometries [bohr/degrees]',
     .    '       =================================',
     .    '        A           ',   '           B           ',
     .    ' OH1 = ',l_OH1,     '      OH3 = ', l_OH3,
     .    ' OH2 = ',l_OH2,     '      OH4 = ', l_OH4

      OH1    = OH1/l_OH1
      OH2    = OH2/l_OH2
      OH3    = OH3/l_OH3
      OH4    = OH4/l_OH4
      angle1 = acos(dot_product(OH1,OH2))*180.d0/pi
      angle2 = acos(dot_product(OH3,OH4))*180.d0/pi

C--Continue writing ...
      write(*,'( a, f9.4, a, f9.4 )')
     .    ' angle1 = ', angle1,  '      angle2 = ', angle2

C--Error checks
      if  ( abs(OH_length-l_OH1) > toler_l ) then
         write(*,'(/a/a,f12.7,a,f12.7)')
     .    ' *** WARNING: O-H1 does not have the right length ',
     .    '               must be: ', OH_length, ' is ',  l_OH1
      endif

      if  ( abs(OH_length-l_OH2) > toler_l ) then
         write(*,'(/a/a,f12.7,a,f12.7)')
     .    ' *** WARNING: O-H2 does not have the right length ',
     .    '               must be: ', OH_length, ' is ',  l_OH2
      endif

      if  ( abs(OH_length-l_OH3) > toler_l ) then
         write(*,'(/a/a,f12.7,a,f12.7)')
     .    ' *** WARNING: O-H3 does not have the right length ',
     .    '              must be: ', OH_length, ' is ',  l_OH3
      endif

      if  ( abs(OH_length-l_OH4) > toler_l ) then
         write(*,'(/a/a,f12.7,a,f12.7)')
     .    ' *** WARNING: O-H4 does not have the right length ',
     .    '              must be: ', OH_length, ' is ',  l_OH4
      endif

      if  ( abs(angle_HOH-angle1) > toler_a ) then
         write(*,'(/a/a,f12.7,a,f12.7)')
     .    ' *** WARNING: First  H-O-H angle not correct ',
     .    '              must be: ',angle_HOH , ' is ',  angle1
      endif

      if  ( abs(angle_HOH-angle2) > toler_a ) then
         write(*,'(/a/a,f12.7,a,f12.7)')
     .    ' *** WARNING: Second H-O-H angle not correct ',
     .    '              must be: ', angle_HOH , ' is ',  angle2
      endif

      monA(:,3) = -(OH1 + OH2)
      monA(:,1) =  (OH1 - OH2)
      monA(:,3) =  monA(:,3)/sqrt(dot_product(monA(:,3), monA(:,3)))
      monA(:,1) =  monA(:,1)/sqrt(dot_product(monA(:,1), monA(:,1)))
      monA(:,2) =  cross( monA(:,3), monA(:,1) )

      write(*,'(/a/a)')
     .    '     Local system of axes on A',
     .    '     ========================='
      do i = 1, 3
        write(*, '( 1x,i1,a,1x, 3(f12.7,2x) )' ) i,':', monA(:,i)
      enddo

      monB(:,3) = -(OH3 + OH4)
      monB(:,1) =  (OH3 - OH4)
      monB(:,3) =  monB(:,3)/sqrt(dot_product(monB(:,3), monB(:,3)))
      monB(:,1) =  monB(:,1)/sqrt(dot_product(monB(:,1), monB(:,1)))
      monB(:,2) =  cross( monB(:,3), monB(:,1) )

      write(*,'(/a/a)')
     .    '     Local system of axes on B',
     .    '     ========================='
      do i = 1, 3
        write(*, '( 1x,i1,a,1x, 3(f12.7,2x) )' ) i,':', monB(:,i)
      enddo

      Call Euler(monA, alphaA,betaA, gammaA)
      Call Euler(monB, alphaB,betaB, gammaB)
      alphaB = alphaB - alphaA
      alphaA = 0.d0

      alphaA  = alphaA*180.d0/pi
      betaA   = betaA *180.d0/pi
      gammaA  = gammaA*180.d0/pi

      alphaB  = alphaB*180.d0/pi
      betaB   = betaB *180.d0/pi
      gammaB  = gammaB*180.d0/pi

      write(*,'(/a/a/ (a,f13.7) )')
     .    ' Euler angles of A: ',
     .    ' ================== ',
     .    ' alpha = ',  alphaA,
     .    ' beta  = ',  betaA,
     .    ' gamma = ',  gammaA
      write(2,'(3 (f14.10,2x) )') alphaA, betaA, gammaA

      write(*,'(/a/a/ (a,f13.7) )')
     .    ' Euler angles of B: ',
     .    ' ================== ',
     .    ' alpha = ',  alphaB,
     .    ' beta  = ',  betaB,
     .    ' gamma = ',  gammaB
      write(2,'(3 (f14.10,2x) )')  alphaB, betaB,  gammaB
      close(2)

      END
      Subroutine Euler(A, alpha,beta, gamma)

C     Determine Euler angles of the proper rotation matrix A.
C     Convention according to Biedenharn & Louck.

      use vector_functions, only: cart2polar
      implicit none
      intent(in)  A
      intent(out) alpha, beta, gamma
      real*8  alpha, beta, gamma,
     .        A(3,3),
     .        polar3(3)

C--Check if matrix is orthogonal with unit determinant.
      Call Check(A)

C--Get polar angles of third column
      polar3  = cart2polar( A(:,3) )    ! returns  R, theta, phi
      alpha   = polar3(3)
      beta    = polar3(2)

C--Compute gamma
      Call Comgamma(A, alpha, beta, gamma)

      end
      Subroutine Check(A)

C     Check if matrix is orthogonal with unit determinant.
      use constants, only: mac_prec
      implicit none
      real*8  t, det,
     .       A(3,3),
     .       unit(3,3)
      intent(in) A

      unit      = 0.d0
      unit(1,1) = 1.d0
      unit(2,2) = 1.d0
      unit(3,3) = 1.d0
      if (abs(sum(matmul(transpose(A),A)-unit)) .gt. mac_prec*1.d2 )
     .  then
         write(*,'(/a/a/3(d10.3) )')
     .       ' Non-orthogonal matrix in subroutine check,',
     .       ' overlap matrix: ',
     .         matmul(transpose(A),A)
      endif

      t = det(A)
      if ( abs(t-1.d0) .gt. mac_prec*1.d2 ) then
         write(*,'(/a,1x,d14.6)')
     .     '  Non-unit determinant:', t
         stop 16
      endif

      end
      real*8 function det(A)

C     Determinant of A

      implicit none
      real*8 A(3,3), B(3)
      integer i
      intent(in) A

      B(1) = A(2,2)*A(3,3) - A(3,2)*A(2,3)
      B(2) = A(3,2)*A(1,3) - A(1,2)*A(3,3)
      B(3) = A(1,2)*A(2,3) - A(2,2)*A(1,3)
      det = 0.d0
      do i=1,3
         det = det + A(i,1)*B(i)
      enddo
      end

      Subroutine Comgamma(A, alpha, beta, gamma)

C     Compute third Euler angle gamma.

      use constants, only: pi, mac_prec
      implicit none
      intent(in)  A, alpha, beta
      intent(out) gamma
      real*8 alpha, beta, gamma,
     .       cb, sb, ca, sa, cg, sg,
     .       A(3,3),
     .       B1(3),
     .       B2(3)
      integer i

      cb    =  cos(beta)
      sb    =  sin(beta)
      ca    =  cos(alpha)
      sa    =  sin(alpha)

      B1(1) =  ca*cb
      B1(2) =  sa*cb
      B1(3) = -sb

      B2(1) = -sa
      B2(2) =  ca
      B2(3) =  0.d0

      cg    =  0.d0
      sg    =  0.d0

      do i=1,3
        cg = cg + B1(i)*A(i,1)
        sg = sg + B2(i)*A(i,1)
      enddo

      if ( cg .ge. 1.d0 ) then
         gamma = 0.d0
      elseif ( cg .le. -1.d0 ) then
         gamma = pi
      else
         gamma = acos(cg)
      endif

      if ( sg  < - mac_prec )  gamma = 2.d0*pi - gamma

      end
