!--------------------------------------------------------------------
! This module contains the data structure for Lennard-Jones type 
! interactions
!--------------------------------------------------------------------
Module lj
  Use defaults, Only: Rdbl, strLen, kcalmole_kb
  Use vector3D, Only: VecType, Assignment(=), Operator(+), &
	Operator(*), vector3D_getnormsq, vector3D_mult_scalar

  Implicit None
  Save

  Private
  Public :: LJ_Params, lj_init, lj_getinteraction, lj_display, &
	lj_cleanup

  Type LJ_Params
    Character(len=strLen) :: atom_name1, atom_name2
    Real(kind=RDbl)     :: sig, eps
    Real(kind=RDbl)     :: A, B      ! A=4*eps*sig^12, B=4*eps*sig^6
    Real(kind=RDbl)     :: cutrad,cutrad2
    Real(kind=RDbl)     :: lowcut    ! lower bound to prevent overflow
  End Type LJ_Params

Contains
  !------------------------------------------------------------------
  ! Initialize the interaction information from the file line.  This
  ! module "knows" how the information is arranged in the 
  ! paramsFilename and is able to read.  This way the file format 
  ! can be tailored for each specific forcefield initialization.
  !------------------------------------------------------------------
  Subroutine lj_init(params, paramsFilename)
    Type(LJ_Params), Intent(INOUT)            :: params
    Character(*), Intent(IN)                  :: paramsFilename

    Integer                                   :: ios, paramunit
    Character(len=strLen)                     :: paramFilename
    Character(len=strLen), Dimension(10)      :: fields,chunks

    !** Open the file
    paramunit = 13
    Open(unit=paramunit, file=paramFilename, status='old', IOSTAT=ios)
    If (ios /= 0) Then
      Write(0,'(2a,i4,2a)') __FILE__,": ",__LINE__, &
          '  Could not open file ',Trim(paramFilename)
      Stop
    End If

    !** Read the various parameters
    Read(paramunit, *) params%atom_name1, params%atom_name2
    Read(paramunit, *) params%sig
    Read(paramunit, *) params%eps
    Read(paramunit, *) params%lowcut
    Read(paramunit, *) params%cutrad
    params%cutrad2 = params%cutrad**2

    Close(paramunit)

    !** Generate the AB parameters
    Call lj_calc_AB(params)
  End Subroutine lj_init

  !------------------------------------------------------------------
  ! Calculate the interaction between two atoms
  !------------------------------------------------------------------
  Subroutine lj_getinteraction(params, sepvec, pot, fvec, ierr)
    Type(LJ_Params), Intent(IN)      :: params
    Type(VecType), Intent(IN)        :: sepvec
    Real(kind = RDbl), Intent(OUT)   :: pot
    Type(VecType), Intent(OUT)       :: fvec
    Integer, Intent(out)             :: ierr

    Real(kind = RDbl)   :: A,B,r2i,r6i,r12i,r2

    pot = 0.0_RDbl
    fvec = 0.0_RDbl       ! Using the overloaded '=' operator
    ierr = 0

    !** Get the square of the distance
    r2 = vector3D_getnormsq(sepvec)

    !**Check if it is within the cut-off radius
    If (r2 > params%cutrad2) Return

    !** Make sure the atoms are not too close
    If (r2 < params%lowcut) Then
      ierr = 1
      Return
    End If

    A = params%A
    B = params%B

    r2i = 1.0_RDbl/r2
    r6i = r2i*r2i*r2i
    r12i= r6i*r6i

    fvec = sepvec  * Real( (12.0_RDbl*A*r12i - 6.0_RDbl*B*r6i)*r2i )  
    pot = A*r12i - B*r6i
  End Subroutine lj_getinteraction

  !------------------------------------------------------------------
  ! This routine calculates A, B parameters from sigma
  ! and epsilon.  A = 4*eps*(sig**12), B = 4*eps*(sig**6)
  !------------------------------------------------------------------
  Subroutine lj_calc_AB(obj)
    Implicit None
    Type(LJ_Params), Intent(INOUT)         :: obj

    obj%B = 4.0_Rdbl*obj%eps*(obj%sig**6)*kcalmole_kb
    obj%A = obj%b*obj%sig**6
  End Subroutine lj_calc_AB

  !------------------------------------------------------------------
  ! This routine calculates the sigma and epsilon values from A and B
  !------------------------------------------------------------------
  Subroutine lj_calc_sigeps(obj)
    Implicit None
    Type(LJ_Params), Intent(INOUT)   :: obj

    Real(kind=RDbl)    :: sig6

    sig6    = (obj%A/obj%B)
    obj%eps = (obj%B/(4.0_Rdbl*sig6*kcalmole_kb))
    obj%sig = sig6**(1/6.0)
  End Subroutine lj_calc_sigeps

  !------------------------------------------------------------------
  ! Display the fields of the LJ structure
  !------------------------------------------------------------------
  Subroutine lj_display(params, optunitno)
    Type(LJ_Params), Intent(in)    :: params
    Integer, Optional, Intent(in)  :: optunitno
    
    Integer    :: unitno

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

    Write(unitno,'(4x,3a)')"Atoms   : ", Trim(params%atom_name1), &
        Trim(params%atom_name2)
    Write(unitno,'(4x,a,f8.3)') "Sigma   : ", params%sig
    Write(unitno,'(4x,a,f8.3)') "Epsilon : ", params%eps
    Write(unitno,'(4x,a,e12.4)')"Sigma   : ", params%A
    Write(unitno,'(4x,a,e12.4)')"Epsilon : ", params%B
  End Subroutine lj_display

  !----------------------------------------------------------
  ! Nothing to do here
  !----------------------------------------------------------
  Subroutine lj_cleanup(params)
    Type(LJ_Params), Intent(inout) :: params
  End Subroutine lj_cleanup

End Module lj