! <compile=optimized>

#include "ncsu-utils.h"
#include "ncsu-config.h"

! vbabin-at-ncsu-dot-edu, 09/08/2010
!
! COM_TORSION -- dihedral angle formed by the centers
!                of mass of four atom groups
!
! input: i = (a1, ..., aN, 0, b1, ..., bM, 0, c1, ..., cK, 0, d1, ..., dL, 0)
! last zero is optional; a?/b?/c?/d? -- indices of the participating atoms

module ncsu_cv_COM_TORSION

implicit none

private

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

public :: colvar_value
public :: colvar_force

public :: colvar_bootstrap
public :: print_details

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

contains

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

function colvar_value(cv, x) result(value)

   NCSU_USE_AFAILED

   use ncsu_constants, only : ZERO
   use ncsu_colvar_type
   use ncsu_colvar_math
   use ncsu_colvar_utils

   implicit none

   NCSU_REAL :: value

   type(colvar_t), intent(in) :: cv
   NCSU_REAL, intent(in) :: x(*)

#  include "ncsu-mpi.h"

   integer :: n

   NCSU_REAL :: cm1(3), cm2(3), cm3(3), cm4(3)

   ncsu_assert(cv%type == COLVAR_COM_TORSION)

   ncsu_assert(associated(cv%i))
   ncsu_assert(size(cv%i).ge.7)
   ncsu_assert(associated(cv%r))
   ncsu_assert(size(cv%r).eq.size(cv%i))

#ifdef MPI
   if (sanderrank.eq.0) then
#endif /* MPI */
   n = 1
   call group_com(cv, x, n, cm1)
   ncsu_assert(n.le.size(cv%i))
   call group_com(cv, x, n, cm2)
   ncsu_assert(n.le.size(cv%i))
   call group_com(cv, x, n, cm3)
   ncsu_assert(n.le.size(cv%i))
   call group_com(cv, x, n, cm4)

   value = torsion(cm1, cm2, cm3, cm4)
#ifdef MPI
   else
      value = ZERO
   end if
#endif /* MPI */

end function colvar_value

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

subroutine colvar_force(cv, x, fcv, f)

   NCSU_USE_AFAILED

   use ncsu_constants, only : ERR_UNIT
   use ncsu_colvar_type
   use ncsu_colvar_math
   use ncsu_colvar_utils
   use ncsu_sander_proxy

   implicit none

   type(colvar_t), intent(in) :: cv
   NCSU_REAL, intent(in) :: x(*), fcv
   NCSU_REAL, intent(inout) :: f(*)

#  include "ncsu-mpi.h"

   NCSU_REAL :: d1(3), d2(3), d3(3), d4(3)
   NCSU_REAL :: cm1(3), cm2(3), cm3(3), cm4(3)

   NCSU_REAL :: d12, d23, d34
   real(8), parameter :: tiny = 1.0d-8

   integer :: n

   ncsu_assert(cv%type == COLVAR_COM_TORSION)

   ncsu_assert(associated(cv%i))
   ncsu_assert(size(cv%i).ge.7)
   ncsu_assert(associated(cv%r))
   ncsu_assert(size(cv%r).eq.size(cv%i))

   NCSU_MASTER_ONLY_BEGIN

   n = 1
   call group_com(cv, x, n, cm1)
   call group_com(cv, x, n, cm2)
   call group_com(cv, x, n, cm3)
   call group_com(cv, x, n, cm4)

   d12 = distance(cm1, cm2)
   d23 = distance(cm2, cm3)
   d34 = distance(cm3, cm4)

   if (d12.lt.tiny.or.d23.lt.tiny.or.d34.lt.tiny) then
      write (unit = ERR_UNIT, fmt = '(/a,a/)') NCSU_ERROR, &
         'COM_TORSION : centers of mass got too close to each other'
      call terminate()
   end if ! too close

   call torsion_d(cm1, cm2, cm3, cm4, d1, d2, d3,  d4)

   d1 = fcv*d1
   d2 = fcv*d2
   d3 = fcv*d3
   d4 = fcv*d4

   n = 1
   call group_com_d(cv, f, d1, n)
   call group_com_d(cv, f, d2, n)
   call group_com_d(cv, f, d3, n)
   call group_com_d(cv, f, d4, n)

   NCSU_MASTER_ONLY_END

end subroutine colvar_force

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

subroutine colvar_bootstrap(cv, cvno, amass)

   use ncsu_utils
   use ncsu_constants
   use ncsu_colvar_type
   use ncsu_colvar_utils
   use ncsu_sander_proxy

   implicit none

   type(colvar_t), intent(inout) :: cv
   integer,        intent(in)    :: cvno
   NCSU_REAL,      intent(in)    :: amass(*)

#  include "ncsu-mpi.h"

   integer :: error, ngroups

   ncsu_assert(cv%type == COLVAR_COM_TORSION)

   call com_check_i(cv%i, cvno, 'COM_TORSION', ngroups)
   if (ngroups.ne.4) then
      NCSU_MASTER_ONLY_BEGIN
         write (unit = ERR_UNIT, fmt = '(/a,a,'//pfmt(cvno)//',a/)') &
               NCSU_ERROR, 'CV #', cvno, &
               ' (COM_TORSION) : number of atom groups is not four'
      NCSU_MASTER_ONLY_END
      call terminate()
   end if ! ngroups.ne.4

   if (associated(cv%r)) &
      deallocate(cv%r)

   allocate(cv%r(size(cv%i)), stat = error)
   if (error.ne.0) &
      NCSU_OUT_OF_MEMORY

   cv%r = ZERO
   do error = 1, size(cv%i)
      if (cv%i(error).gt.0) &
         cv%r(error) = amass(cv%i(error))
   end do

   call com_init_weights(cv%r)

end subroutine colvar_bootstrap

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

subroutine print_details(cv, lun)

#ifndef NCSU_DISABLE_ASSERT
   use ncsu_utils
   use ncsu_sander_proxy
#endif /* NCSU_DISABLE_ASSERT */

   use ncsu_colvar_type
   use ncsu_colvar_utils

   implicit none

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

   ncsu_assert(is_master())
   ncsu_assert(cv%type == COLVAR_COM_TORSION)
   ncsu_assert(associated(cv%i))

   call com_print_i(cv%i, lun)

end subroutine print_details

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

end module ncsu_cv_COM_TORSION