function v_DISTANCE_Z(cv, x) result(value)

   use nfe_lib_mod
   use parallel_dat_mod

   implicit none

   double precision :: value

   type(colvar_t), intent(in) :: cv

   double precision, intent(in) :: x(*)

   integer :: a1, a2

   nfe_assert(cv%type == COLVAR_DISTANCE_Z)

   nfe_assert(associated(cv%i))
   nfe_assert(size(cv%i) == 2)
   nfe_assert(size(cv%r) == 3)

   nfe_assert(cv%i(1) > 0 .and. cv%i(1) <= pmemd_natoms())
   a1 = 3*cv%i(1) - 2

   nfe_assert(cv%i(2) > 0 .and. cv%i(2) <= pmemd_natoms())
   a2 = 3*cv%i(2) - 2

   nfe_assert(a1 /= a2)

#ifdef MPI
   if (mytaskid.eq.0) then
#endif /* MPI */
      value = (x(a2)-x(a1))*cv%r(1)+(x(a2+1)-x(a1+1))*cv%r(2)+(x(a2+2)-x(a1+2))*cv%r(3)
#ifdef MPI
   else
      value = ZERO
   end if
#endif /* MPI */

end function v_DISTANCE_Z

!=============================================================================

subroutine f_DISTANCE_Z(cv, x, fcv, f)

   use nfe_lib_mod
   use parallel_dat_mod

   implicit none

   type(colvar_t), intent(in) :: cv

   double precision, intent(in) :: x(*), fcv

   double precision, intent(inout) :: f(*)

   double precision :: d1(3), d2(3)

   integer :: a1, a2
#ifdef MPI
   integer :: error
#endif /* MPI */

   nfe_assert(cv%type == COLVAR_DISTANCE_Z)

   nfe_assert(associated(cv%i))
   nfe_assert(size(cv%i) == 2)
   nfe_assert(size(cv%r) == 3)

   nfe_assert(cv%i(1) > 0 .and. cv%i(1) <= pmemd_natoms())
   a1 = 3*cv%i(1) - 2

   nfe_assert(cv%i(2) > 0 .and. cv%i(2) <= pmemd_natoms())
   a2 = 3*cv%i(2) - 2

   nfe_assert(a1 /= a2)

   NFE_MASTER_ONLY_BEGIN
   
   d1(1) = -cv%r(1); d1(2) = -cv%r(2); d1(3) = -cv%r(3)
   d2(1) =  cv%r(1); d2(2) =  cv%r(2); d2(3) =  cv%r(3)

   f(a1:a1 + 2) = f(a1:a1 + 2) + fcv*d1
   f(a2:a2 + 2) = f(a2:a2 + 2) + fcv*d2

   NFE_MASTER_ONLY_END
#ifdef MPI
   call mpi_bcast(f(1:3*pmemd_natoms()), 3*pmemd_natoms(), MPI_DOUBLE_PRECISION, 0, pmemd_comm, error)
   nfe_assert(error.eq.0)
#endif /* MPI */

end subroutine f_DISTANCE_Z

!=============================================================================

subroutine b_DISTANCE_Z(cv, cvno)

   use nfe_lib_mod
   use parallel_dat_mod

   implicit none

   type(colvar_t), intent(inout) :: cv
   integer,        intent(in)    :: cvno
   integer :: error

   nfe_assert(cv%type == COLVAR_DISTANCE_Z)
   call check_i(cv%i, cvno, 'DISTANCE_Z', 2)

   if ((size(cv%r).ne.3).or.(abs(norm3(cv%r)-1.0d0).gt.TINY)) then   
      NFE_MASTER_ONLY_BEGIN
         write (unit = ERR_UNIT, fmt = '(a,a,'//pfmt(cvno)//',a)') &
            NFE_ERROR, 'CV #', cvno, &
            ' (DISTANCE_Z) : unexpected unit vector for axis'
      NFE_MASTER_ONLY_END
      call terminate()
   end if

end subroutine b_DISTANCE_Z

!=============================================================================

subroutine p_DISTANCE_Z(cv, lun)

   use nfe_lib_mod

   implicit none

   type(colvar_t), intent(in) :: cv
   integer, intent(in) :: lun

   nfe_assert(is_master())
   nfe_assert(cv%type == COLVAR_DISTANCE_Z)
   nfe_assert(associated(cv%i))

   write (unit = lun, fmt = '(a,a,'//pfmt(cv%r(1),3)//','//pfmt(cv%r(2),3)//','//pfmt(cv%r(3),3)//')') &
      NFE_INFO, '      axis = ', cv%r(1), cv%r(2), cv%r(3)
   call print_i(cv%i, lun)

end subroutine p_DISTANCE_Z
