!------------------------------------------------------------------
! This is a "proxy" module to access different pairwise forcefields
! using manually implemented polymorphism.  It is used to provide
! a generic interface to initialize different forcefield models and 
! then get interactions from them.
! Lennard Jones potentials (lj) and Buckingham potentials (buck) are 
! accessible through this module. More potential models can be 
! added similarly.
!------------------------------------------------------------------
Module pairmodels
  Use vector3D, Only: VecType, vector3D_getnormsq
  Use lj, Only: LJ_Params, lj_init, lj_getinteraction, lj_display, &
	lj_cleanup
  Use buck, Only: BUCK_Params, buck_init, buck_getinteraction, &
        buck_display, buck_cleanup
  Use defaults, Only: RDbl
  Implicit None
  Save

  Private
  Public:: Pairmodels_Type, pairmodels_init, pairmodels_cleanup, &
	pairmodels_getinteraction, pairmodels_display
  
  Type Pairmodels_Type
    Type(LJ_Params), Pointer         :: lj
    Type(BUCK_Params), Pointer       :: buck
  End Type Pairmodels_Type

Contains
  !-------------------------------------------------------------------
  ! Initialize the parameters for atom types "atomtype1" and
  ! "atomtype2" to the "modeltype" of forcefield from file 
  ! "paramFilename". The parameters for all atoms are stored in 
  ! "pmodels".
  !-------------------------------------------------------------------
  Subroutine pairmodels_init(pmodel, modeltype, paramFilename)
    Type(Pairmodels_Type),Intent(inout) :: pmodel
    Character(*), Intent(in)            :: modeltype
    Character(*), Intent(in)            :: paramFilename

    !** Initialize the appropriate forcefield
    If (modeltype == 'LJ') Then
      Call lj_init(pmodel%lj, paramFilename)
    Else If (modeltype == 'BUCK') Then
      Call buck_init(pmodel%buck, paramFilename)
    End If
  End Subroutine pairmodels_init
  
  !-------------------------------------------------------------------
  ! Get the potential energy "pot" and force "f" between atom types 
  ! "atom1" and "atom2" separated by the vector "sepvec".  This is in 
  ! effect implementing dynamic polymorphism
  !-------------------------------------------------------------------
  Subroutine pairmodels_getinteraction(pmodel, sepvec, pot, f)
    Type(Pairmodels_Type),Intent(in) :: pmodel
    Type(VecType), Intent(in)        :: sepvec
    Type(VecType), Intent(out)        :: f 
    Real(kind=RDbl), Intent(out)     :: pot
    Integer :: ierr

    If (Associated(pmodel%lj)) Then
      Call lj_getinteraction(pmodel%lj, sepvec, pot, f, ierr)
    Else If (Associated(pmodel%buck)) Then
      Call buck_getinteraction(pmodel%buck, sepvec, pot, f, ierr)
    End If
  End Subroutine pairmodels_getinteraction
  
  !------------------------------------------------------------
  ! Displays the parameters 
  !------------------------------------------------------------
  Subroutine pairmodels_display(pmodel, optunitno)
    Type(Pairmodels_Type),Intent(in) :: pmodel
    Integer, Optional, Intent(in)  :: optunitno
    
    Integer    :: unitno

    If (Present(optunitno)) Then
      unitno = optunitno
    Else
      unitno = 6    !Standard output
    End If

    If (Associated(pmodel%lj)) Then
      Call lj_display(pmodel%lj, unitno)
    Else If (Associated(pmodel%buck)) Then
      Call buck_display(pmodel%buck, unitno)
    End If
  End Subroutine pairmodels_display

  !------------------------------------------------------------
  ! Call the cleanup routines for the different forcefields
  !------------------------------------------------------------
  Subroutine pairmodels_cleanup(pmodel)
    Type(Pairmodels_Type),Intent(inout) :: pmodel
        If (Associated(pmodel%lj)) Then
          Call lj_cleanup(pmodel%lj)
        Else If (Associated(pmodel%buck)) Then
          Call buck_cleanup(pmodel%buck)
        End If
  End Subroutine pairmodels_cleanup
End Module pairmodels