mdtools/md_code/lj_FF.F90

!! Calculates Long Range forces Lennard-Jones interactions.
!! Corresponds to the force field defined in lj_FF.cpp 
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: lj_FF.F90,v 1.10 2003-01-24 21:36:52 chuckv Exp $, $Date: 2003-01-24 21:36:52 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $


module lj_ff
  use simulation
  use definitions
#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

!! Number of lj_atypes in lj_atype_list
  integer, save :: n_lj_atypes = 0

!! Starting Size for ljMixed Array
  integer, parameter :: ljMixed_blocksize = 10

!! Basic atom type for a Lennard-Jones Atom.
  type, public :: lj_atype
     private
     sequence
!! Unique number for place in linked list
     integer :: atype_number = 0
!! Unique indentifier number (ie atomic no, etc)
     integer :: atype_ident = 0
!! Mass of Particle
     real ( kind = dp )  :: mass = 0.0_dp
!! Lennard-Jones epslon
     real ( kind = dp )  :: epsilon = 0.0_dp
!! Lennard-Jones Sigma
     real ( kind = dp )  :: sigma = 0.0_dp
!! Lennard-Jones Sigma Squared
     real ( kind = dp )  :: sigma2 = 0.0_dp
!! Lennard-Jones Sigma to sixth
     real ( kind = dp )  :: sigma6 = 0.0_dp
!! Pointer for linked list creation
     type (lj_atype), pointer :: next => null()
  end type lj_atype

!! Pointer type for atype ident array
  type, public :: lj_atypePtr
     type (lj_atype), pointer :: this => null()
  end type lj_atypePtr

!! Global list of lj atypes in simulation
  type (lj_atype), pointer :: lj_atype_list => null()
!! LJ mixing array  
  type (lj_atype), dimension(:,:), pointer :: ljMixed => null() 
!! identity pointer list for force loop.
  type (lj_atypePtr), dimension(:), pointer :: identPtrList => null()


!! Neighbor list and commom arrays
!! This should eventally get moved to a neighbor list type
  integer, allocatable, dimension(:) :: point
  integer, allocatable, dimension(:) :: list
  integer, parameter :: listMultiplier = 80
  integer :: nListAllocs = 0
  integer, parameter :: maxListAllocs = 5

#ifdef IS_MPI
! Universal routines: All types of force calculations will need these arrays
! Arrays specific to a type of force calculation should be declared in that module.
!! Row position array.
  real( kind = dp ), allocatable, dimension(:,:) :: qRow
!! Column position array.
  real( kind = dp ), allocatable, dimension(:,:) :: qColumn

!! Row Force array.
  real( kind = dp ), allocatable, dimension(:,:) :: fRow
!! Column Force Array.
  real( kind = dp ), allocatable, dimension(:,:) :: fColumn

!! Row potential energy array
  real( kind = dp ), allocatable, dimension(:) :: eRow
!! Column potential energy array
  real( kind = dp ), allocatable, dimension(:) :: eColumn


!! Row lj_atype pointer list
  type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null()
!! Column lj_atype pointer list
  type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null()
#endif


!! Logical has lj force field module been initialized?
  logical, save :: isljFFinit = .false.


!! Public methods and data
  public :: new_lj_atype 
  public :: do_lj_ff
  public :: getLjPot
  public :: init_ljFF

  
contains

  subroutine new_lj_atype(ident,mass,epsilon,sigma,status)
    real( kind = dp ), intent(in) :: mass
    real( kind = dp ), intent(in) :: epsilon
    real( kind = dp ), intent(in) :: sigma
    integer, intent(in) :: ident
    integer, intent(out) :: status

    type (lj_atype), pointer :: this_lj_atype
    type (lj_atype), pointer :: lj_atype_ptr

    type (lj_atype), dimension(:,:), pointer :: thisMix
    type (lj_atype), dimension(:,:), pointer :: oldMix
    integer :: alloc_error
    integer :: atype_counter = 0
    integer :: alloc_size

    status = 0


! allocate a new atype    
    allocate(this_lj_atype,stat=alloc_error)
    if (alloc_error /= 0 ) then
       status = -1
       return
    end if

! assign our new lj_atype information
    this_lj_atype%mass        = mass
    this_lj_atype%epsilon     = epsilon
    this_lj_atype%sigma       = sigma
    this_lj_atype%sigma2      = sigma * sigma
    this_lj_atype%sigma6      = this_lj_atype%sigma2 * this_lj_atype%sigma2 &
         * this_lj_atype%sigma2
! assume that this atype will be successfully added
    this_lj_atype%atype_ident = ident
    this_lj_atype%atype_number = n_lj_atypes + 1


! First time through allocate a array of size ljMixed_blocksize
    if(.not. associated(ljMixed)) then
       allocate(thisMix(ljMixed_blocksize,ljMixed_blocksize))
       if (alloc_error /= 0 ) then
          status = -1
          return
       end if
       ljMixed => thisMix
! If we have outgrown ljMixed_blocksize, allocate a new matrix twice the size and
! point ljMix at the new matrix.
    else if( (n_lj_atypes + 1) > size(ljMixed)) then
       alloc_size = 2*size(ljMixed)
       allocate(thisMix(alloc_size,alloc_size))
       if (alloc_error /= 0 ) then
          status = -1
          return
       end if
! point oldMix at old ljMixed array
       oldMix => ljMixed
! Copy oldMix into new Mixed array      
       thisMix = oldMix
! Point ljMixed at new array
       ljMixed => thisMix
! Free old array so we don't have a memory leak
       deallocate(oldMix)
    endif


! Find bottom of atype master list
! if lj_atype_list is null then we are at the top of the list.
    if (.not. associated(lj_atype_list)) then

       write(*,*) "Associating lists for first time"
       

       this_lj_atype%atype_number = 1
       
       ljMixed(n_lj_atypes,n_lj_atypes)%sigma  = this_lj_atype%sigma
       
       ljMixed(n_lj_atypes,n_lj_atypes)%sigma2  = ljMixed(n_lj_atypes,n_lj_atypes)%sigma & 
            * ljMixed(n_lj_atypes,n_lj_atypes)%sigma
       
       ljMixed(n_lj_atypes,n_lj_atypes)%sigma6 = ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 &
            * ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 &
            * ljMixed(n_lj_atypes,n_lj_atypes)%sigma2
       
       ljMixed(n_lj_atypes,n_lj_atypes)%epsilon = this_lj_atype%epsilon
       
       lj_atype_list => this_lj_atype
       
       write(*,*) "lj_atype_list first time through -- ident", lj_atype_list%atype_ident
    else ! we need to find the bottom of the list to insert new atype
       write(*,*) "Second time through"
       lj_atype_ptr => lj_atype_list
       write(*,*) "lj_atype_ptr to next"
       atype_counter = 1
       find_end: do
          if (.not. associated(lj_atype_ptr)) exit find_end

          write(*,*) "Inside find_end-this ident", lj_atype_ptr%atype_ident
! Set up mixing for new atype and current atype in list
       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma  =  &
            calcLJMix("sigma",this_lj_atype%sigma, &
            lj_atype_ptr%sigma)

       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2  = &
            ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma & 
            * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma

       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma6 = &
            ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 &
            * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 &
            * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2

       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%epsilon = &
            calcLJMix("epsilon",this_lj_atype%epsilon, &
            lj_atype_ptr%epsilon)

! Advance to next pointer
       lj_atype_ptr => lj_atype_ptr%next           
       atype_counter = atype_counter + 1          
       end do find_end

       lj_atype_ptr => this_lj_atype
       write(*,*) "just added: ", lj_atype_ptr%atype_ident
    end if

    write(*,*) "lj_atype_list now contains.."
    
    lj_atype_ptr => lj_atype_list
    do while (associated(lj_atype_ptr))
       write(*,*) "lj_atype_list contains ident", lj_atype_ptr%atype_ident
       lj_atype_ptr => lj_atype_ptr%next
    end do
    
! Insert new atype at end of list

! Increment number of atypes

    n_lj_atypes = n_lj_atypes + 1

! Set up self mixing

  end subroutine new_lj_atype


  subroutine init_ljFF(nComponents,ident, status)
 !! Number of components in ident array
    integer, intent(inout) :: nComponents
!! Array of identities nComponents long corresponding to
!! ljatype ident.
    integer, dimension(nComponents),intent(inout) :: ident
!!  Result status, success = 0, error = -1
    integer, intent(out) :: Status

    integer :: alloc_stat

    integer :: thisStat
    integer :: i
#ifdef IS_MPI
    integer, allocatable, dimension(:) :: identRow
    integer, allocatable, dimension(:) :: identCol
    integer :: nrow
    integer :: ncol
    integer :: alloc_stat
#endif
    status = 0

    write(*,*) "Initializing ljFF"
!! if were're not in MPI, we just update ljatypePtrList
#ifndef IS_MPI
    call new_ljatypePtrList(nComponents,ident,identPtrList,thisStat)
    if ( thisStat /= 0 ) then
       status = -1
       return
    endif
    write(*,*) "component pointer list initialized"
!! Allocate pointer lists
    allocate(point(nComponents),stat=alloc_stat)
    if (alloc_stat /=0) then
       status = -1
       return
    endif

    allocate(list(nComponents*listMultiplier),stat=alloc_stat)
    if (alloc_stat /=0) then
       status = -1
       return
    endif
    
! if were're in MPI, we also have to worry about row and col lists    
#else
! We can only set up forces if mpiSimulation has been setup.
    if (.not. isMPISimSet()) then
       status = -1
       return
    endif
    nrow = getNrow()
    ncol = getNcol()
!! Allocate temperary arrays to hold gather information
    allocate(identRow(nrow),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       status = -1
       return
    endif

    allocate(identCol(ncol),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       status = -1
       return
    endif
!! Gather idents into row and column idents
    call gather(ident,identRow,plan_row)
    call gather(ident,identCol,plan_col)

!! Create row and col pointer lists
    call new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat)
    if (thisStat /= 0 ) then
       status = -1
       return
    endif

    call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat)
    if (thisStat /= 0 ) then
       status = -1
       return
    endif

!! free temporary ident arrays
    deallocate(identCol)
    deallocate(identRow)

!! Allocate Simulation arrays
!! NOTE: This bit of code should be fixed, it can cause large
!! memory fragmentation if call repeatedly

    if (.not.allocated(qRow)) then
       allocate(qRow(3,nrow),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(qrow)
       allocate(qRow(3,nrow),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif

    if (.not.allocated(3,qCol)) then
       allocate(qCol(ncol),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(qCol)
       allocate(qCol(3,ncol),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif

    if (.not.allocated(fRow)) then
       allocate(fRow(3,nrow),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(fRow)
       allocate(fRow(3,nrow),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif

    if (.not.allocated(fCol)) then
       allocate(fCol(3,ncol),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(fCol)
       allocate(fCol(3,ncol),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif

    if (.not.allocated(eRow)) then
       allocate(eRow(nrow),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(eRow)
       allocate(eRow(nrow),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif

    if (.not.allocated(eCol)) then
       allocate(eCol(ncol),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(eCol)
       allocate(eCol(ncol),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif

    if (.not.allocated(eTemp)) then
       allocate(eTemp(nComponents),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    else
       deallocate(eTemp)
       allocate(eTemp(nComponets),stat=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1
          return
       endif
    endif


!! Allocate neighbor lists for mpi simulations.
    if (.not. allocated(point)) then
       allocate(point(nrow),stat=alloc_stat)
       if (alloc_stat /=0) then
          status = -1
          return
       endif

       allocate(list(ncol*listMultiplier),stat=alloc_stat)
       if (alloc_stat /=0) then
          status = -1
          return
       endif
    else
       deallocate(list)
       deallocate(point)
       allocate(point(nrow),stat=alloc_stat)
       if (alloc_stat /=0) then
          status = -1
          return
       endif

       allocate(list(ncol*listMultiplier),stat=alloc_stat)
       if (alloc_stat /=0) then
          status = -1
          return
       endif
    endif

#endif
    
    
    do i = 1, nComponents
       write(*,*) "Ident pointer list ident: ", identPtrList(i)%this%atype_ident
    end do

    isljFFinit = .true.

  end subroutine init_ljFF


!! Takes an ident array and creates an atype pointer list
!! based on those identities
  subroutine new_ljatypePtrList(mysize,ident,PtrList,status)
    integer, intent(in) :: mysize
    integer,dimension(:), intent(in) :: ident
    integer, optional :: status
    type(lj_atypePtr), dimension(:), pointer :: PtrList

    integer :: thisIdent
    integer :: i
    integer :: alloc_error
    type (lj_atype), pointer :: tmpPtr

    if (present(status)) status = 0

    write(*,*) "Creating new ljatypePtrList...."
! First time through, allocate list
    if (.not.associated(PtrList)) then
       write(*,*) "allocating pointer list...."
       allocate(PtrList(mysize))
    else
! We want to creat a new ident list so free old list
       deallocate(PtrList)
       allocate(PtrList(mysize))
    endif

! Match pointer list
    write(*,*) "Doing ident search"
    do i = 1, mysize
       thisIdent = ident(i)
       write(*,*) "Calling getLJatype for ident ", thisIdent
       call getLJatype(thisIdent,tmpPtr)
 
      if (.not. associated(tmpPtr)) then
         write(*,*) "tmpptr was not associated"
          status = -1
          return
       endif
       
       PtrList(i)%this => tmpPtr
    end do

  end subroutine new_ljatypePtrList

!! Finds a lj_atype based upon numerical ident 
!! returns a null pointer if error
  subroutine getLJatype(ident,ljAtypePtr)
    integer, intent(in) :: ident
    type (lj_atype), pointer :: ljAtypePtr    
    type (lj_atype), pointer :: tmplj_atype_ptr => null()

    write(*,*) "Finding ljAtype for ident ", ident
    ljAtypePtr => null()

    if(.not. associated(lj_atype_list)) return

! Point at head of list.
    write(*,*) "Searching at head of list"
    tmplj_atype_ptr => lj_atype_list
    find_ident: do 
       if (.not.associated(tmplj_atype_ptr)) then
          write(*,*) "Find_ident, pointer not associated"
          exit find_ident
       else if( tmplj_atype_ptr%atype_ident == ident) then
          ljAtypePtr => tmplj_atype_ptr
          exit find_ident
       endif
       tmplj_atype_ptr => tmplj_atype_ptr%next
    end do find_ident


  end subroutine getLJatype


!! FORCE routine Calculates Lennard Jones forces.
!------------------------------------------------------------->
  subroutine do_lj_ff(q,f,potE,do_pot)
!! Position array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,getNlocal()) :: q
!! Force array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,getNlocal()) :: f
    real ( kind = dp ) :: potE
    logical ( kind = 2) :: do_pot
    
    type(lj_atype), pointer :: ljAtype_i
    type(lj_atype), pointer :: ljAtype_j

#ifdef MPI
  real( kind = DP ), dimension(3,getNcol()) :: efr
  real( kind = DP ) :: pot_local
#else
  real( kind = DP ), dimension(3,getNlocal()) :: efr
#endif
  
  real( kind = DP )   :: pe
  logical             :: update_nlist


  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
  integer :: nlist
  integer :: j_start
  integer :: tag_i,tag_j
  real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
  real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut

! a rig that need to be fixed.
#ifdef IS_MPI
  logical :: newtons_thrd = .true.
  real( kind = dp ) :: pe_local
  integer :: nlocal
#endif
  integer :: nrow
  integer :: ncol
  integer :: natoms

! Make sure we are properly initialized.
  if (.not. isljFFInit) then
     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
     return
  endif
#ifdef IS_MPI
    if (.not. isMPISimSet()) then
     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
     return
  endif
#endif

!! initialize local variables  
  natoms = getNlocal()
  call getRcut(rcut,rcut2=rcutsq)
  call getRlist(rlist,rlistsq)
#ifndef IS_MPI
  nrow = natoms - 1
  ncol = natoms
#else
  nrow = getNrow(plan_row)
  ncol = getNcol(plan_col)
  nlocal = natoms
  j_start = 1
#endif

  
!! See if we need to update neighbor lists
  call check(q,update_nlist)

!--------------WARNING...........................
! Zero variables, NOTE:::: Forces are zeroed in C
! Zeroing them here could delete previously computed
! Forces.
!------------------------------------------------
#ifndef IS_MPI
!  nloops = nloops + 1
  pe = 0.0E0_DP
 
#else
    f_row = 0.0E0_DP
    f_col = 0.0E0_DP

    pe_local = 0.0E0_DP

    eRow = 0.0E0_DP
    eCol = 0.0E0_DP
    eTemp = 0.0E0_DP
#endif
    efr = 0.0E0_DP

! communicate MPI positions
#ifdef MPI    
    call gather(q,qRow,plan_row3)
    call gather(q,qCol,plan_col3)
#endif


  if (update_nlist) then

     ! save current configuration, contruct neighbor list, 
     ! and calculate forces
     call save_nlist(q)
     
     nlist = 0
     
     
     do i = 1, nrow
        point(i) = nlist + 1
#ifdef MPI
        ljAtype_i => identPtrListRow(i)%this
        tag_i = tagRow(i)
        rxi = qRow(1,i)
        ryi = qRow(2,i)
        rzi = qRow(3,i)
#else
        ljAtype_i   => identPtrList(i)%this
        j_start = i + 1
        rxi = q(1,i)
        ryi = q(2,i)
        rzi = q(3,i)
#endif

        inner: do j = j_start, ncol
#ifdef MPI
! Assign identity pointers and tags
           ljAtype_j => identPtrListColumn(j)%this
           tag_j = tagCol(j)
           if (newtons_thrd) then
              if (tag_i <= tag_j) then
                 if (mod(tag_i + tag_j,2) == 0) cycle inner
              else                
                 if (mod(tag_i + tag_j,2) == 1) cycle inner
              endif
           endif

           rxij = wrap(rxi - qCol(1,j), 1)
           ryij = wrap(ryi - qCol(2,j), 2)
           rzij = wrap(rzi - qCol(3,j), 3)
#else           
           ljAtype_j   => identPtrList(j)%this
           rxij = wrap(rxi - q(1,j), 1)
           ryij = wrap(ryi - q(2,j), 2)
           rzij = wrap(rzi - q(3,j), 3)
       
#endif           
           rijsq = rxij*rxij + ryij*ryij + rzij*rzij

#ifdef MPI
             if (rijsq <=  rlistsq .AND. &
                  tag_j /= tag_i) then
#else
             if (rijsq <  rlistsq) then
#endif
             
              nlist = nlist + 1
              if (nlist > size(list)) then
!!  "Change how nlist size is done"
                 write(DEFAULT_ERROR,*) "ERROR: nlist > list size"
              endif
              list(nlist) = j

              
              if (rijsq <  rcutsq) then
                
                 r = dsqrt(rijsq)
       
                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
       
#ifdef MPI
                eRow(i) = eRow(i) + pot*0.5
                eCol(i) = eCol(i) + pot*0.5
#else
                pe = pe + pot 
#endif                 
             
                 efr(1,j) = -rxij
                 efr(2,j) = -ryij
                 efr(3,j) = -rzij

                 do dim = 1, 3   

             
                    drdx1 = efr(dim,j) / r
                    ftmp = dudr * drdx1


#ifdef MPI
                    fCol(dim,j) = fCol(dim,j) - ftmp
                    fRow(dim,i) = fRow(dim,i) + ftmp
#else                    
             
                    f(dim,j) = f(dim,j) - ftmp
                    f(dim,i) = f(dim,i) + ftmp

#endif                    
                 enddo
              endif
           endif
        enddo inner
     enddo 

#ifdef MPI
     point(nrow + 1) = nlist + 1
#else
     point(natoms) = nlist + 1
#endif

  else

     ! use the list to find the neighbors
     do i = 1, nrow
        JBEG = POINT(i)
        JEND = POINT(i+1) - 1
        ! check thiat molecule i has neighbors
        if (jbeg .le. jend) then
#ifdef MPI
           ljAtype_i => identPtrListRow(i)%this
           rxi = qRow(1,i)
           ryi = qRow(2,i)
           rzi = qRow(3,i)
#else
           ljAtype_i   => identPtrList(i)%this
           rxi = q(1,i)
           ryi = q(2,i)
           rzi = q(3,i)
#endif
           do jnab = jbeg, jend
              j = list(jnab)
#ifdef MPI
              ljAtype_j = identPtrListColumn(j)%this
              rxij = wrap(rxi - q_col(1,j), 1)
              ryij = wrap(ryi - q_col(2,j), 2)
              rzij = wrap(rzi - q_col(3,j), 3)
#else
              ljAtype_j = identPtrList(j)%this
              rxij = wrap(rxi - q(1,j), 1)
              ryij = wrap(ryi - q(2,j), 2)
              rzij = wrap(rzi - q(3,j), 3)
#endif
              rijsq = rxij*rxij + ryij*ryij + rzij*rzij
              
              if (rijsq <  rcutsq) then

                 r = dsqrt(rijsq)
                 
                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
#ifdef MPI
                eRow(i) = eRow(i) + pot*0.5
                eCol(i) = eCol(i) + pot*0.5
#else
                pe = pe + pot 
#endif                 

                 
                 efr(1,j) = -rxij
                 efr(2,j) = -ryij
                 efr(3,j) = -rzij

                 do dim = 1, 3                        
                    
                    drdx1 = efr(dim,j) / r
                    ftmp = dudr * drdx1
#ifdef MPI
                    fCol(dim,j) = fCol(dim,j) - ftmp
                    fRow(dim,i) = fRow(dim,i) + ftmp
#else                    
                    f(dim,j) = f(dim,j) - ftmp
                    f(dim,i) = f(dim,i) + ftmp
#endif                    
                 enddo
              endif
           enddo
        endif
     enddo
  endif


#ifdef MPI
    !!distribute forces
    call scatter(fRow,f,plan_row3)

    call scatter(fCol,f_tmp,plan_col3)

    do i = 1,nlocal
       f(1:3,i) = f(1:3,i) + f_tmp(1:3,i)
    end do


    if (do_pot) then
! scatter/gather pot_row into the members of my column
       call scatter(e_row,e_tmp,plan_row)
       
       ! scatter/gather pot_local into all other procs
       ! add resultant to get total pot
       do i = 1, nlocal
          pe_local = pe_local + e_tmp(i)
       enddo
       if (newtons_thrd) then
          e_tmp = 0.0E0_DP
          call scatter(e_col,e_tmp,plan_col)
          do i = 1, nlocal
             pe_local = pe_local + e_tmp(i)
          enddo
       endif
    endif
    potE = pe_local
#endif

    potE = pe


  end subroutine do_lj_ff

!! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
!! derivatives.
  subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status)
! arguments
!! Length of vector between particles
    real( kind = dp ), intent(in)  :: r
!! Potential Energy
    real( kind = dp ), intent(out) :: pot
!! Derivatve wrt postion
    real( kind = dp ), intent(out) :: dudr
!! Second Derivative, optional, used mainly for normal mode calculations.
    real( kind = dp ), intent(out), optional :: d2
    
    type (lj_atype), pointer :: atype1
    type (lj_atype), pointer :: atype2

    integer, intent(out), optional :: status

! local Variables
    real( kind = dp ) :: sigma
    real( kind = dp ) :: sigma2
    real( kind = dp ) :: sigma6
    real( kind = dp ) :: epsilon

    real( kind = dp ) :: rcut
    real( kind = dp ) :: rcut2
    real( kind = dp ) :: rcut6
    real( kind = dp ) :: r2
    real( kind = dp ) :: r6

    real( kind = dp ) :: t6
    real( kind = dp ) :: t12
    real( kind = dp ) :: tp6
    real( kind = dp ) :: tp12
    real( kind = dp ) :: delta

    logical :: doSec = .false.

    integer :: errorStat

!! Optional Argument Checking
! Check to see if we need to do second derivatives
    
    if (present(d2))     doSec = .true.
    if (present(status)) status = 0

! Look up the correct parameters in the mixing matrix
    sigma    = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
    sigma2   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
    sigma6   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
    epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
 

    call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
    
    r2 = r * r
    r6 = r2 * r2 * r2

    t6  = sigma6/ r6
    t12 = t6 * t6


    tp6 = sigma6 / rcut6
    tp12 = tp6*tp6
                                                                               
    delta = -4.0E0_DP*epsilon * (tp12 - tp6)
                                                                               
    if (r.le.rcut) then
       pot = 4.0E0_DP * epsilon * (t12 - t6) + delta
       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / r
       if(doSec)  d2 = 24.0E0_DP * epsilon * (26.0E0_DP*t12 - 7.0E0_DP*t6)/r/r
    else
       pot = 0.0E0_DP
       dudr = 0.0E0_DP
       if(doSec) d2 = 0.0E0_DP
    endif
                                                                               
  return


  end subroutine getLjPot


!! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array.
  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
    character(len=*) :: thisParam
    real(kind = dp)  :: param1
    real(kind = dp)  :: param2
    real(kind = dp ) :: myMixParam
    integer, optional :: status


    myMixParam = 0.0_dp

    if (present(status)) status = 0

    select case (thisParam)

    case ("sigma")
       myMixParam = 0.5_dp * (param1 + param2)
    case ("epsilon")
       myMixParam = sqrt(param1 * param2)
    case default
       status = -1
    end select

  end function calcLJMix


end module lj_ff
Revision:	246
Committed:	Fri Jan 24 21:36:52 2003 UTC (21 years, 5 months ago) by chuckv
File size:	27267 byte(s)
Log Message:	lj_FF.F90 and simulation_module now compile. Logical bug still exists in lj_FF.F90 lj_atypes_list.