OOPSE/libmdtools/do_Forces.F90

!! do_Forces.F90
!! module do_Forces
!! Calculates Long Range forces.

!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: do_Forces.F90,v 1.48 2004-01-30 15:01:09 chrisfen Exp $, $Date: 2004-01-30 15:01:09 $, $Name: not supported by cvs2svn $, $Revision: 1.48 $

module do_Forces
  use force_globals
  use simulation
  use definitions
  use atype_module
  use neighborLists  
  use lj
  use sticky_pair
  use dipole_dipole
  use charge_charge
  use reaction_field
  use gb_pair
  use vector_class
  use eam
  use status
#ifdef IS_MPI
  use mpiSimulation
#endif

  implicit none
  PRIVATE

#define __FORTRAN90
#include "fForceField.h"

  logical, save :: haveRlist = .false.
  logical, save :: haveNeighborList = .false.
  logical, save :: havePolicies = .false.
  logical, save :: haveSIMvariables = .false.
  logical, save :: havePropertyMap = .false.
  logical, save :: haveSaneForceField = .false.
  logical, save :: FF_uses_LJ
  logical, save :: FF_uses_sticky
  logical, save :: FF_uses_charges
  logical, save :: FF_uses_dipoles
  logical, save :: FF_uses_RF
  logical, save :: FF_uses_GB
  logical, save :: FF_uses_EAM
  logical, save :: SIM_uses_LJ
  logical, save :: SIM_uses_sticky
  logical, save :: SIM_uses_charges
  logical, save :: SIM_uses_dipoles
  logical, save :: SIM_uses_RF
  logical, save :: SIM_uses_GB
  logical, save :: SIM_uses_EAM
  logical, save :: SIM_requires_postpair_calc
  logical, save :: SIM_requires_prepair_calc
  logical, save :: SIM_uses_directional_atoms
  logical, save :: SIM_uses_PBC

  real(kind=dp), save :: rlist, rlistsq

  public :: init_FF
  public :: do_force_loop
  public :: setRlistDF


#ifdef PROFILE
  public :: getforcetime
  real, save :: forceTime = 0
  real :: forceTimeInitial, forceTimeFinal
  integer :: nLoops
#endif

  type :: Properties
     logical :: is_lj     = .false.
     logical :: is_sticky = .false.
     logical :: is_dp     = .false.
     logical :: is_gb     = .false.
     logical :: is_eam    = .false.
     logical :: is_charge = .false.
     real(kind=DP) :: charge = 0.0_DP
     real(kind=DP) :: dipole_moment = 0.0_DP
  end type Properties

  type(Properties), dimension(:),allocatable :: PropertyMap

contains

  subroutine setRlistDF( this_rlist )
    
    real(kind=dp) :: this_rlist

    rlist = this_rlist
    rlistsq = rlist * rlist
    
    haveRlist = .true.

  end subroutine setRlistDF    

  subroutine createPropertyMap(status)
    integer :: nAtypes
    integer :: status
    integer :: i
    logical :: thisProperty
    real (kind=DP) :: thisDPproperty

    status = 0

    nAtypes = getSize(atypes)

    if (nAtypes == 0) then
       status = -1
       return
    end if
        
    if (.not. allocated(PropertyMap)) then
       allocate(PropertyMap(nAtypes))
    endif

    do i = 1, nAtypes
       call getElementProperty(atypes, i, "is_LJ", thisProperty)
       PropertyMap(i)%is_LJ = thisProperty

       call getElementProperty(atypes, i, "is_Charge", thisProperty)
       PropertyMap(i)%is_Charge = thisProperty
       
       if (thisProperty) then
          call getElementProperty(atypes, i, "charge", thisDPproperty)
          PropertyMap(i)%charge = thisDPproperty
       endif

       call getElementProperty(atypes, i, "is_DP", thisProperty)
       PropertyMap(i)%is_DP = thisProperty

       if (thisProperty) then
          call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
          PropertyMap(i)%dipole_moment = thisDPproperty
       endif

       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
       PropertyMap(i)%is_Sticky = thisProperty
       call getElementProperty(atypes, i, "is_GB", thisProperty)
       PropertyMap(i)%is_GB = thisProperty
       call getElementProperty(atypes, i, "is_EAM", thisProperty)
       PropertyMap(i)%is_EAM = thisProperty
    end do

    havePropertyMap = .true.

  end subroutine createPropertyMap

  subroutine setSimVariables()
    SIM_uses_LJ = SimUsesLJ()
    SIM_uses_sticky = SimUsesSticky()
    SIM_uses_charges = SimUsesCharges()
    SIM_uses_dipoles = SimUsesDipoles()
    SIM_uses_RF = SimUsesRF()
    SIM_uses_GB = SimUsesGB()
    SIM_uses_EAM = SimUsesEAM()
    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
    SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
    SIM_uses_PBC = SimUsesPBC()

    haveSIMvariables = .true.

    return
  end subroutine setSimVariables

  subroutine doReadyCheck(error)
    integer, intent(out) :: error

    integer :: myStatus

    error = 0
    
    if (.not. havePropertyMap) then 

       myStatus = 0

       call createPropertyMap(myStatus)

       if (myStatus .ne. 0) then
          write(default_error, *) 'createPropertyMap failed in do_Forces!'
          error = -1
          return
       endif
    endif

    if (.not. haveSIMvariables) then
       call setSimVariables()
    endif

    if (.not. haveRlist) then
       write(default_error, *) 'rList has not been set in do_Forces!'
       error = -1
       return
    endif

    if (SIM_uses_LJ .and. FF_uses_LJ) then
       if (.not. havePolicies) then
          write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
          error = -1
          return
       endif
    endif

    if (.not. haveNeighborList) then
       write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
       error = -1
       return
    end if

    if (.not. haveSaneForceField) then
       write(default_error, *) 'Force Field is not sane in do_Forces!'
       error = -1
       return
    end if

#ifdef IS_MPI
    if (.not. isMPISimSet()) then
       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
       error = -1
       return
    endif
#endif
    return
  end subroutine doReadyCheck
    

  subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)

    integer, intent(in) :: LJMIXPOLICY
    logical, intent(in) :: use_RF_c

    integer, intent(out) :: thisStat   
    integer :: my_status, nMatches
    integer, pointer :: MatchList(:) => null()
    real(kind=dp) :: rcut, rrf, rt, dielect

    !! assume things are copacetic, unless they aren't
    thisStat = 0

    !! Fortran's version of a cast:
    FF_uses_RF = use_RF_c
    
    !! init_FF is called *after* all of the atom types have been 
    !! defined in atype_module using the new_atype subroutine.
    !!
    !! this will scan through the known atypes and figure out what
    !! interactions are used by the force field.     
  
    FF_uses_LJ = .false.
    FF_uses_sticky = .false.
    FF_uses_charges = .false.
    FF_uses_dipoles = .false.
    FF_uses_GB = .false.
    FF_uses_EAM = .false.
    
    call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
    if (nMatches .gt. 0) FF_uses_LJ = .true.

    call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
    if (nMatches .gt. 0) FF_uses_charges = .true.
    
    call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
    if (nMatches .gt. 0) FF_uses_dipoles = .true.
    
    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
         MatchList) 
    if (nMatches .gt. 0) FF_uses_Sticky = .true.
    
    call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
    if (nMatches .gt. 0) FF_uses_GB = .true.
    
    call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
    if (nMatches .gt. 0) FF_uses_EAM = .true.
    
    !! Assume sanity (for the sake of argument)
    haveSaneForceField = .true.

    !! check to make sure the FF_uses_RF setting makes sense
    
    if (FF_uses_dipoles) then
       if (FF_uses_RF) then
          dielect = getDielect()
          call initialize_rf(dielect)
       endif
    else
       if (FF_uses_RF) then          
          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
          thisStat = -1
          haveSaneForceField = .false.
          return
       endif
    endif 

    if (FF_uses_LJ) then
       
       select case (LJMIXPOLICY)
       case (LB_MIXING_RULE)
          call init_lj_FF(LB_MIXING_RULE, my_status)             
       case (EXPLICIT_MIXING_RULE)
          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
       case default
          write(default_error,*) 'unknown LJ Mixing Policy!'
          thisStat = -1
          haveSaneForceField = .false.
          return             
       end select
       if (my_status /= 0) then
          thisStat = -1
          haveSaneForceField = .false.
          return
       end if
       havePolicies = .true.
    endif

    if (FF_uses_sticky) then
       call check_sticky_FF(my_status)
       if (my_status /= 0) then
          thisStat = -1
          haveSaneForceField = .false.
          return
       end if
    endif


    if (FF_uses_EAM) then
         call init_EAM_FF(my_status) 
       if (my_status /= 0) then
          write(default_error, *) "init_EAM_FF returned a bad status"
          thisStat = -1
          haveSaneForceField = .false.
          return
       end if
    endif

    if (FF_uses_GB) then
       call check_gb_pair_FF(my_status)
       if (my_status .ne. 0) then
          thisStat = -1
          haveSaneForceField = .false.
          return
       endif
    endif

    if (FF_uses_GB .and. FF_uses_LJ) then
    endif
    if (.not. haveNeighborList) then
       !! Create neighbor lists
       call expandNeighborList(nLocal, my_status)
       if (my_Status /= 0) then
          write(default_error,*) "SimSetup: ExpandNeighborList returned error."
          thisStat = -1
          return
       endif
       haveNeighborList = .true.
    endif
    
  end subroutine init_FF
  

  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
  !------------------------------------------------------------->
  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
       error)
    !! Position array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,nLocal) :: q
    !! Rotation Matrix for each long range particle in simulation.
    real( kind = dp), dimension(9,nLocal) :: A    
    !! Unit vectors for dipoles (lab frame)
    real( kind = dp ), dimension(3,nLocal) :: u_l
    !! Force array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,nLocal) :: f
    !! Torsion array provided by C, dimensioned by getNlocal
    real( kind = dp ), dimension(3,nLocal) :: t    

    !! Stress Tensor
    real( kind = dp), dimension(9) :: tau   
    real ( kind = dp ) :: pot
    logical ( kind = 2) :: do_pot_c, do_stress_c
    logical :: do_pot
    logical :: do_stress
#ifdef IS_MPI 
    real( kind = DP ) :: pot_local
    integer :: nrow
    integer :: ncol
    integer :: nprocs
#endif
    integer :: natoms    
    logical :: update_nlist   
    integer :: i, j, jbeg, jend, jnab
    integer :: nlist
    real( kind = DP ) ::  rijsq 
    real(kind=dp),dimension(3) :: d
    real(kind=dp) :: rfpot, mu_i, virial
    integer :: me_i, me_j
    logical :: is_dp_i
    integer :: neighborListSize
    integer :: listerror, error
    integer :: localError
    integer :: propPack_i, propPack_j

    real(kind=dp) :: listSkin = 1.0   

    !! initialize local variables  

#ifdef IS_MPI
    pot_local = 0.0_dp
    nrow   = getNrow(plan_row)
    ncol   = getNcol(plan_col)
#else
    natoms = nlocal
#endif

    call doReadyCheck(localError)
    if ( localError .ne. 0 ) then
       call handleError("do_force_loop", "Not Initialized")
       error = -1
       return
    end if
    call zero_work_arrays()

    do_pot = do_pot_c
    do_stress = do_stress_c

    ! Gather all information needed by all force loops:
    
#ifdef IS_MPI    

    call gather(q,q_Row,plan_row3d)
    call gather(q,q_Col,plan_col3d)
        
    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
       call gather(u_l,u_l_Row,plan_row3d)
       call gather(u_l,u_l_Col,plan_col3d)
       
       call gather(A,A_Row,plan_row_rotation)
       call gather(A,A_Col,plan_col_rotation)
    endif
    
#endif

!! Begin force loop timing:
#ifdef PROFILE
    call cpu_time(forceTimeInitial)
    nloops = nloops + 1
#endif
  
    if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
       !! See if we need to update neighbor lists
       call checkNeighborList(nlocal, q, listSkin, update_nlist)   
       !! if_mpi_gather_stuff_for_prepair
       !! do_prepair_loop_if_needed
       !! if_mpi_scatter_stuff_from_prepair
       !! if_mpi_gather_stuff_from_prepair_to_main_loop
     
!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
#ifdef IS_MPI
    
    if (update_nlist) then
       
       !! save current configuration, construct neighbor list, 
       !! and calculate forces
       call saveNeighborList(nlocal, q)
       
       neighborListSize = size(list)
       nlist = 0       
       
       do i = 1, nrow
          point(i) = nlist + 1
          
          prepair_inner: do j = 1, ncol
             
             if (skipThisPair(i,j)) cycle prepair_inner
             
             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
             
             if (rijsq < rlistsq) then             
                
                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandNeighborList(nlocal, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                   neighborListSize = size(list)
                endif
                
                list(nlist) = j
                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
             endif
          enddo prepair_inner
       enddo

       point(nrow + 1) = nlist + 1
       
    else  !! (of update_check)

       ! 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
             
             do jnab = jbeg, jend
                j = list(jnab)

                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
                     u_l, A, f, t, pot_local)

             enddo
          endif
       enddo
    endif
    
#else
    
    if (update_nlist) then
       
       ! save current configuration, contruct neighbor list, 
       ! and calculate forces
       call saveNeighborList(natoms, q)
       
       neighborListSize = size(list)
   
       nlist = 0

       do i = 1, natoms-1
          point(i) = nlist + 1
          
          prepair_inner: do j = i+1, natoms
             
             if (skipThisPair(i,j))  cycle prepair_inner
                          
             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
           

             if (rijsq < rlistsq) then

           
                nlist = nlist + 1
              
                if (nlist > neighborListSize) then
                   call expandNeighborList(natoms, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                   neighborListSize = size(list)
                endif
                
                list(nlist) = j
                
                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
                        u_l, A, f, t, pot)
                
             endif
          enddo prepair_inner
       enddo
       
       point(natoms) = nlist + 1
       
    else !! (update)
   
       ! use the list to find the neighbors
       do i = 1, natoms-1
          JBEG = POINT(i)
          JEND = POINT(i+1) - 1
          ! check thiat molecule i has neighbors
          if (jbeg .le. jend) then
             
             do jnab = jbeg, jend
                j = list(jnab)

                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
                     u_l, A, f, t, pot)

             enddo
          endif
       enddo
    endif    
#endif 
    !! Do rest of preforce calculations
    !! do necessary preforce calculations   
    call do_preforce(nlocal,pot)
   ! we have already updated the neighbor list set it to false...
   update_nlist = .false.
    else
       !! See if we need to update neighbor lists for non pre-pair
       call checkNeighborList(nlocal, q, listSkin, update_nlist)   
    endif
 

!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>


#ifdef IS_MPI
    
    if (update_nlist) then
       !! save current configuration, construct neighbor list, 
       !! and calculate forces
       call saveNeighborList(nlocal, q)
       
       neighborListSize = size(list)
       nlist = 0       
       
       do i = 1, nrow

          point(i) = nlist + 1
          
          inner: do j = 1, ncol
             
             if (skipThisPair(i,j)) cycle inner
             
             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
             
             if (rijsq < rlistsq) then             
                
                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandNeighborList(nlocal, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                   neighborListSize = size(list)
                endif
                
                list(nlist) = j
                                
                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
                     u_l, A, f, t, pot_local)
                
             endif
          enddo inner
       enddo

       point(nrow + 1) = nlist + 1
       
    else  !! (of update_check)

       ! 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
             
             do jnab = jbeg, jend
                j = list(jnab)

                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
                     u_l, A, f, t, pot_local)

             enddo
          endif
       enddo
    endif
    
#else
    
    if (update_nlist) then

       ! save current configuration, contruct neighbor list, 
       ! and calculate forces
       call saveNeighborList(natoms, q)
       
       neighborListSize = size(list)
   
       nlist = 0
       
       do i = 1, natoms-1
          point(i) = nlist + 1
          
          inner: do j = i+1, natoms
             
             if (skipThisPair(i,j))  cycle inner
                          
             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
           

             if (rijsq < rlistsq) then
                
                nlist = nlist + 1
              
                if (nlist > neighborListSize) then
                   call expandNeighborList(natoms, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                   neighborListSize = size(list)
                endif
                
                list(nlist) = j
                
                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
                        u_l, A, f, t, pot)
                
             endif
          enddo inner
       enddo
       
       point(natoms) = nlist + 1
       
    else !! (update)
       
       ! use the list to find the neighbors
       do i = 1, natoms-1
          JBEG = POINT(i)
          JEND = POINT(i+1) - 1
          ! check thiat molecule i has neighbors
          if (jbeg .le. jend) then
             
             do jnab = jbeg, jend
                j = list(jnab)

                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
                     u_l, A, f, t, pot)

             enddo
          endif
       enddo
    endif
    
#endif 
    
    ! phew, done with main loop.
 
!! Do timing
#ifdef PROFILE
    call cpu_time(forceTimeFinal)
    forceTime = forceTime + forceTimeFinal - forceTimeInitial
#endif


#ifdef IS_MPI
    !!distribute forces
   
    f_temp = 0.0_dp
    call scatter(f_Row,f_temp,plan_row3d)
    do i = 1,nlocal
       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
    end do

    f_temp = 0.0_dp
    call scatter(f_Col,f_temp,plan_col3d)
    do i = 1,nlocal
       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
    end do
    
    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
       t_temp = 0.0_dp
       call scatter(t_Row,t_temp,plan_row3d)
       do i = 1,nlocal
          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
       end do
       t_temp = 0.0_dp
       call scatter(t_Col,t_temp,plan_col3d)
       
       do i = 1,nlocal
          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
       end do
    endif
    
    if (do_pot) then
       ! scatter/gather pot_row into the members of my column
       call scatter(pot_Row, pot_Temp, plan_row)

       ! scatter/gather pot_local into all other procs
       ! add resultant to get total pot
       do i = 1, nlocal
          pot_local = pot_local + pot_Temp(i)
       enddo
       
       pot_Temp = 0.0_DP 

       call scatter(pot_Col, pot_Temp, plan_col)
       do i = 1, nlocal
          pot_local = pot_local + pot_Temp(i)
       enddo

    endif    
#endif

    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
       
       if (FF_uses_RF .and. SIM_uses_RF) then
 
#ifdef IS_MPI
          call scatter(rf_Row,rf,plan_row3d)
          call scatter(rf_Col,rf_Temp,plan_col3d)
          do i = 1,nlocal
             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
          end do
#endif
          
          do i = 1, nLocal

             rfpot = 0.0_DP
#ifdef IS_MPI
             me_i = atid_row(i)
#else
             me_i = atid(i)
#endif

             if (PropertyMap(me_i)%is_DP) then

                mu_i = PropertyMap(me_i)%dipole_moment

                !! The reaction field needs to include a self contribution 
                !! to the field:
                call accumulate_self_rf(i, mu_i, u_l)
                !! Get the reaction field contribution to the 
                !! potential and torques:
                call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#ifdef IS_MPI
                pot_local = pot_local + rfpot
#else
                pot = pot + rfpot
      
#endif
             endif             
          enddo
       endif
    endif


#ifdef IS_MPI

    if (do_pot) then
       pot = pot + pot_local
       !! we assume the c code will do the allreduce to get the total potential
       !! we could do it right here if we needed to...
    endif

    if (do_stress) then
      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
            mpi_comm_world,mpi_err)
       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
            mpi_comm_world,mpi_err)
    endif

#else

    if (do_stress) then
       tau = tau_Temp
       virial = virial_Temp
    endif
    
#endif
    
    
  end subroutine do_force_loop

  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)

    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3,nLocal) :: u_l
    real (kind=dp), dimension(9,nLocal) :: A
    real (kind=dp), dimension(3,nLocal) :: f
    real (kind=dp), dimension(3,nLocal) :: t

    logical, intent(inout) :: do_pot, do_stress
    integer, intent(in) :: i, j
    real ( kind = dp ), intent(inout)    :: rijsq
    real ( kind = dp )                :: r
    real ( kind = dp ), intent(inout) :: d(3)
    integer :: me_i, me_j

    r = sqrt(rijsq)

#ifdef IS_MPI
    if (tagRow(i) .eq. tagColumn(j)) then
       write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
    endif
    me_i = atid_row(i)
    me_j = atid_col(j)
#else
    me_i = atid(i)
    me_j = atid(j)
#endif
    
    if (FF_uses_LJ .and. SIM_uses_LJ) then
       
       if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
          call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
       endif
       
    endif
    
    if (FF_uses_charges .and. SIM_uses_charges) then
       
       if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
          call do_charge_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
       endif
       
    endif
    
    if (FF_uses_dipoles .and. SIM_uses_dipoles) then
       
       if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
          call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
               do_pot, do_stress)
          if (FF_uses_RF .and. SIM_uses_RF) then
             call accumulate_rf(i, j, r, u_l)
             call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
          endif          
       endif

    endif

    if (FF_uses_Sticky .and. SIM_uses_sticky) then

       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
          call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
               do_pot, do_stress)
       endif

    endif


    if (FF_uses_GB .and. SIM_uses_GB) then
       
       if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
          call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
               do_pot, do_stress)          
       endif

    endif
       
    if (FF_uses_EAM .and. SIM_uses_EAM) then
       
       if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
          call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
       endif
       
    endif

  end subroutine do_pair


  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
   real( kind = dp ) :: pot
   real( kind = dp ), dimension(3,nLocal) :: u_l
   real (kind=dp), dimension(9,nLocal) :: A
   real (kind=dp), dimension(3,nLocal) :: f
   real (kind=dp), dimension(3,nLocal) :: t
   
   logical, intent(inout) :: do_pot, do_stress
   integer, intent(in) :: i, j
   real ( kind = dp ), intent(inout)    :: rijsq
   real ( kind = dp )                :: r
   real ( kind = dp ), intent(inout) :: d(3)
   
   logical :: is_EAM_i, is_EAM_j
   
   integer :: me_i, me_j
   
   r = sqrt(rijsq)
   

#ifdef IS_MPI
   if (tagRow(i) .eq. tagColumn(j)) then
      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
   endif
   
   me_i = atid_row(i)
   me_j = atid_col(j)
   
#else
   
   me_i = atid(i)
   me_j = atid(j)
   
#endif
    
   if (FF_uses_EAM .and. SIM_uses_EAM) then

      if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
           call calc_EAM_prepair_rho(i, j, d, r, rijsq )

   endif
   
 end subroutine do_prepair


  subroutine do_preforce(nlocal,pot)
    integer :: nlocal
    real( kind = dp ) :: pot

    if (FF_uses_EAM .and. SIM_uses_EAM) then
       call calc_EAM_preforce_Frho(nlocal,pot)
    endif


  end subroutine do_preforce
  
  
  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
    
    real (kind = dp), dimension(3) :: q_i
    real (kind = dp), dimension(3) :: q_j
    real ( kind = dp ), intent(out) :: r_sq
    real( kind = dp ) :: d(3), scaled(3)
    integer i

    d(1:3) = q_j(1:3) - q_i(1:3)

    ! Wrap back into periodic box if necessary
    if ( SIM_uses_PBC ) then
       
       if( .not.boxIsOrthorhombic ) then
          ! calc the scaled coordinates.
          
          scaled = matmul(HmatInv, d)
          
          ! wrap the scaled coordinates

          scaled = scaled  - anint(scaled)
          

          ! calc the wrapped real coordinates from the wrapped scaled 
          ! coordinates

          d = matmul(Hmat,scaled)

       else
          ! calc the scaled coordinates.
          
          do i = 1, 3
             scaled(i) = d(i) * HmatInv(i,i)
             
             ! wrap the scaled coordinates
             
             scaled(i) = scaled(i) - anint(scaled(i))
             
             ! calc the wrapped real coordinates from the wrapped scaled 
             ! coordinates

             d(i) = scaled(i)*Hmat(i,i)
          enddo
       endif
       
    endif
    
    r_sq = dot_product(d,d)
    
  end subroutine get_interatomic_vector
  
  subroutine zero_work_arrays()
    
#ifdef IS_MPI

    q_Row = 0.0_dp
    q_Col = 0.0_dp   
    
    u_l_Row = 0.0_dp
    u_l_Col = 0.0_dp
    
    A_Row = 0.0_dp
    A_Col = 0.0_dp
    
    f_Row = 0.0_dp
    f_Col = 0.0_dp
    f_Temp = 0.0_dp
      
    t_Row = 0.0_dp
    t_Col = 0.0_dp
    t_Temp = 0.0_dp
 
    pot_Row = 0.0_dp
    pot_Col = 0.0_dp
    pot_Temp = 0.0_dp

    rf_Row = 0.0_dp
    rf_Col = 0.0_dp
    rf_Temp = 0.0_dp

#endif

 
    if (FF_uses_EAM .and. SIM_uses_EAM) then
       call clean_EAM()
    endif


    rf = 0.0_dp
    tau_Temp = 0.0_dp
    virial_Temp = 0.0_dp
  end subroutine zero_work_arrays
  
  function skipThisPair(atom1, atom2) result(skip_it)
    integer, intent(in) :: atom1
    integer, intent(in), optional :: atom2
    logical :: skip_it
    integer :: unique_id_1, unique_id_2
    integer :: me_i,me_j
    integer :: i

    skip_it = .false. 
    
    !! there are a number of reasons to skip a pair or a particle
    !! mostly we do this to exclude atoms who are involved in short
    !! range interactions (bonds, bends, torsions), but we also need 
    !! to exclude some overcounted interactions that result from
    !! the parallel decomposition
    
#ifdef IS_MPI
    !! in MPI, we have to look up the unique IDs for each atom
    unique_id_1 = tagRow(atom1)
#else
    !! in the normal loop, the atom numbers are unique
    unique_id_1 = atom1
#endif

    !! We were called with only one atom, so just check the global exclude
    !! list for this atom
    if (.not. present(atom2)) then
       do i = 1, nExcludes_global
          if (excludesGlobal(i) == unique_id_1) then
             skip_it = .true.
             return
          end if
       end do
       return 
    end if
    
#ifdef IS_MPI
    unique_id_2 = tagColumn(atom2)
#else
    unique_id_2 = atom2
#endif

#ifdef IS_MPI
    !! this situation should only arise in MPI simulations
    if (unique_id_1 == unique_id_2) then
       skip_it = .true.
       return
    end if
    
    !! this prevents us from doing the pair on multiple processors
    if (unique_id_1 < unique_id_2) then
       if (mod(unique_id_1 + unique_id_2,2) == 0) then
          skip_it = .true.
          return
       endif
    else                
       if (mod(unique_id_1 + unique_id_2,2) == 1) then 
          skip_it = .true.
          return
       endif
    endif
#endif
 
    !! the rest of these situations can happen in all simulations:
    do i = 1, nExcludes_global       
       if ((excludesGlobal(i) == unique_id_1) .or. &
            (excludesGlobal(i) == unique_id_2)) then
          skip_it = .true.
          return
       endif
    enddo

    do i = 1, nExcludes_local
       if (excludesLocal(1,i) == unique_id_1) then
          if (excludesLocal(2,i) == unique_id_2) then
             skip_it = .true.
             return
          endif
       else
          if (excludesLocal(1,i) == unique_id_2) then
             if (excludesLocal(2,i) == unique_id_1) then
                skip_it = .true.
                return
             endif
          endif
       endif
    end do
    
    return
  end function skipThisPair

  function FF_UsesDirectionalAtoms() result(doesit)
    logical :: doesit
    doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
         FF_uses_GB .or. FF_uses_RF
  end function FF_UsesDirectionalAtoms
  
  function FF_RequiresPrepairCalc() result(doesit)
    logical :: doesit
    doesit = FF_uses_EAM
  end function FF_RequiresPrepairCalc
  
  function FF_RequiresPostpairCalc() result(doesit)
    logical :: doesit
    doesit = FF_uses_RF
  end function FF_RequiresPostpairCalc
  
#ifdef PROFILE
  function getforcetime() result(totalforcetime)
    real(kind=dp) :: totalforcetime
    totalforcetime = forcetime
  end function getforcetime
#endif

!! This cleans componets of force arrays belonging only to fortran

end module do_Forces
Revision:	999
Committed:	Fri Jan 30 15:01:09 2004 UTC (20 years, 5 months ago) by chrisfen
File size:	32939 byte(s)
Log Message:	Substantial changes. OOPSE now has a working WATER.cpp forcefield and parser. This involved changes to WATER.cpp and ForceFields amoung other files. One important note: a hardwiring of LJ_rcut was made in calc_LJ_FF.F90. This will be removed on the next commit...