mdtools/libmdCode/do_Forces.F90

!! Calculates Long Range forces.
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: do_Forces.F90,v 1.4 2003-03-06 22:08:29 gezelter Exp $, $Date: 2003-03-06 22:08:29 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $


module do_Forces
  use simulation
  use definitions
  use generic_atypes
  use neighborLists
  
  use lj_FF
  use sticky_FF
  use dp_FF
  use gb_FF

#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

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

!! Global list of lj atypes in simulation
  type (atype), pointer :: ListHead => null()
  type (atype), pointer :: ListTail => null()


  logical, save :: firstTime = .True.

!! Atype identity pointer lists
#ifdef IS_MPI
!! Row lj_atype pointer list
  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
!! Column lj_atype pointer list
  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
#else
  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
#endif


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

!! Use periodic boundry conditions
  logical :: wrap = .false.

!! Potential energy global module variables
#ifdef IS_MPI
  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp

  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp

  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp

  real(kind = dp), dimension(9,getNrow(plan_row)) :: ARow = 0.0_dp
  real(kind = dp), dimension(9,getNcol(plan_col)) :: ACol = 0.0_dp  

  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: fTemp1 = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: tTemp1 = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: fTemp2 = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: tTemp2 = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp

  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp

  real(kind = dp), dimension(3,getNrow(plan_row)) :: rflRow = 0.0_dp
  real(kind = dp), dimension(3,getNcol(plan_col)) :: rflCol = 0.0_dp
  real(kind = dp), dimension(3,getNlocal()) :: rflTemp = 0.0_dp

  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp

  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
#endif
  real(kind = dp) :: pe = 0.0_dp
  real(kind = dp), dimension(3,natoms) :: fTemp = 0.0_dp
  real(kind = dp), dimension(3,natoms) :: tTemp = 0.0_dp
  real(kind = dp), dimension(3,natoms) :: rflTemp = 0.0_dp
  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp

  logical :: do_preForce  = .false.
  logical :: do_postForce = .false.


!! Public methods and data
  public :: new_atype 
  public :: do_forceLoop
  public :: init_FF

  
contains

!! Adds a new lj_atype to the list. 
  subroutine new_atype(ident,mass,epsilon,sigma, & 
       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
    real( kind = dp ), intent(in) :: mass
    real( kind = dp ), intent(in) :: epsilon
    real( kind = dp ), intent(in) :: sigma
    real( kind = dp ), intent(in) :: w0
    real( kind = dp ), intent(in) :: v0
    real( kind = dp ), intent(in) :: dipoleMoment

    integer, intent(in) :: ident
    integer, intent(out) :: status
    integer, intent(in) :: is_Sticky
    integer, intent(in) :: is_DP
    integer, intent(in) :: is_GB
    integer, intent(in) :: is_LJ


    type (atype), pointer :: the_new_atype
    integer :: alloc_error
    integer :: atype_counter = 0
    integer :: alloc_size
    integer :: err_stat
    status = 0


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

! assign our new atype information
    the_new_atype%mass        = mass
    the_new_atype%epsilon     = epsilon
    the_new_atype%sigma       = sigma
    the_new_atype%sigma2      = sigma * sigma
    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
         * the_new_atype%sigma2
    the_new_atype%w0       = w0
    the_new_atype%v0       = v0
    the_new_atype%dipoleMoment       = dipoleMoment
 
    
! assume that this atype will be successfully added
    the_new_atype%atype_ident = ident
    the_new_atype%atype_number = n_lj_atypes + 1
    
    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.

    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
    if (err_stat /= 0 ) then 
       status = -1
       return
    endif

    n_atypes = n_atypes + 1


  end subroutine new_atype


  subroutine init_FF(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

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

!! if were're not in MPI, we just update ljatypePtrList
#ifndef IS_MPI
    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
    if ( thisStat /= 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
       write(default_error,*) "MPI is not set"
       status = -1
       return
    endif
    nrow = getNrow(plan_row)
    ncol = getNcol(plan_col)
    mynode = getMyNode()
!! 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 create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
    if (thisStat /= 0 ) then
       status = -1
       return
    endif
  
    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
    if (thisStat /= 0 ) then
       status = -1
       return
    endif

!! free temporary ident arrays
    if (allocated(identCol)) then
       deallocate(identCol)
    end if
    if (allocated(identCol)) then
       deallocate(identRow)
    endif

#endif
    
    call initForce_Modules(thisStat)
    if (thisStat /= 0) then
       status = -1
       return
    endif

!! Create neighbor lists
    call expandList(thisStat)
    if (thisStat /= 0) then
       status = -1
       return
    endif

    isFFinit = .true.


  end subroutine init_FF


  subroutine initForce_Modules(thisStat)
    integer, intent(out) :: thisStat
    integer :: my_status
    
    thisStat = 0
    call init_lj_FF(ListHead,my_status)
    if (my_status /= 0) then
       thisStat = -1
       return
    end if

  end subroutine initForce_Modules


!! FORCE routine Calculates Lennard Jones forces.
!------------------------------------------------------------->
  subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
!! Position array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,getNlocal()) :: q
  !! Rotation Matrix for each long range particle in simulation.
    real( kind = dp), dimension(9,getNlocal()) :: A

  !! Magnitude dipole moment
    real( kind = dp ), dimension(3,getNlocal()) :: mu
  !! Unit vectors for dipoles (lab frame)
    real( kind = dp ), dimension(3,getNlocal()) :: u_l
!! Force array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,getNlocal()) :: f
!! Torsion array provided by C, dimensioned by getNlocal
    real( kind = dp ), dimension(3,getNlocal()) :: t

!! Stress Tensor
    real( kind = dp), dimension(9) :: tau

    real ( kind = dp ) :: potE
    logical ( kind = 2) :: do_pot
    integer :: FFerror

    
    type(atype), pointer :: Atype_i
    type(atype), pointer :: Atype_j


#ifdef IS_MPI
  real( kind = DP ) :: pot_local

!! Local arrays needed for MPI

 #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
 
  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp

  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut

  real( kind = DP ) :: dielectric = 0.0_dp

! a rig that need to be fixed.
#ifdef IS_MPI
  real( kind = dp ) :: pe_local
  integer :: nlocal
#endif
  integer :: nrow
  integer :: ncol
  integer :: natoms
  integer :: neighborListSize
  integer :: listerror
!! should we calculate the stress tensor
  logical  :: do_stress = .false.


  FFerror = 0

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

!! initialize local variables  
  natoms = getNlocal()
  call getRcut(rcut,rcut2=rcutsq)
  call getRlist(rlist,rlistsq)

!! Find ensemble
  if (isEnsemble("NPT")) do_stress = .true.
!! set to wrap
  if (isPBC()) wrap = .true.


!! 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.
!------------------------------------------------
  call zero_module_variables()


! communicate MPI positions
#ifdef IS_MPI    
    call gather(q,qRow,plan_row3d)
    call gather(q,qCol,plan_col3d)

    call gather(mu,muRow,plan_row3d)
    call gather(mu,muCol,plan_col3d)

    call gather(u_l,u_lRow,plan_row3d)
    call gather(u_l,u_lCol,plan_col3d)

    call gather(A,ARow,plan_row_rotation)
    call gather(A,ACol,plan_col_rotation)
#endif


#ifdef IS_MPI
    
    if (update_nlist) then
       
       ! save current configuration, contruct neighbor list, 
       ! and calculate forces
       call save_neighborList(q)
       
       neighborListSize = getNeighborListSize()
       nlist = 0
       
       nrow = getNrow(plan_row)
       ncol = getNcol(plan_col)
       nlocal = getNlocal()
       
       do i = 1, nrow
          point(i) = nlist + 1
          Atype_i => identPtrListRow(i)%this
          
          inner: do j = 1, ncol
             Atype_j => identPtrListColumn(j)%this
             
             call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
                  rxij,ryij,rzij,rijsq,r)
             
             ! skip the loop if the atoms are identical
             if (mpi_cycle_jLoop(i,j)) cycle inner:
             
             if (rijsq <  rlistsq) then             
                
                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandList(listerror)
                   if (listerror /= 0) then
                      FFerror = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                endif
                
                list(nlist) = j
                
                
                if (rijsq <  rcutsq) then
                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
                endif
             endif
          enddo inner
       enddo

       point(nrow + 1) = nlist + 1
       
    else !! (update)

       ! 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

             Atype_i => identPtrListRow(i)%this
             do jnab = jbeg, jend
                j = list(jnab)
                Atype_j = identPtrListColumn(j)%this
                call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
                     rxij,ryij,rzij,rijsq,r)
                
                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
             enddo
          endif
       enddo
    endif

#else
    
    if (update_nlist) then
       
       ! save current configuration, contruct neighbor list, 
       ! and calculate forces
       call save_neighborList(q)
       
       neighborListSize = getNeighborListSize()
       nlist = 0
       
    
       do i = 1, natoms-1
          point(i) = nlist + 1
          Atype_i   => identPtrList(i)%this

          inner: do j = i+1, natoms
             Atype_j   => identPtrList(j)%this
             call get_interatomic_vector(i,j,q(:,i),q(:,j),&
                  rxij,ryij,rzij,rijsq,r)
           
             if (rijsq <  rlistsq) then

                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandList(listerror)
                   if (listerror /= 0) then
                      FFerror = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                endif
                
                list(nlist) = j

    
                if (rijsq <  rcutsq) then
                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
                endif
             endif
          enddo inner
       enddo
       
       point(natoms) = nlist + 1
       
    else !! (update)
       
       ! 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
             
             Atype_i => identPtrList(i)%this
             do jnab = jbeg, jend
                j = list(jnab)
                Atype_j = identPtrList(j)%this
                call get_interatomic_vector(i,j,q(:,i),q(:,j),&
                     rxij,ryij,rzij,rijsq,r)
                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
             enddo
          endif
       enddo
    endif

#endif 


#ifdef IS_MPI
    !! distribute all reaction field stuff (or anything for post-pair):
    call scatter(rflRow,rflTemp1,plan_row3d)
    call scatter(rflCol,rflTemp2,plan_col3d)
    do i = 1,nlocal
       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
    end do
#endif

! This is the post-pair loop:
#ifdef IS_MPI

    if (system_has_postpair_atoms) then
       do i = 1, nlocal
          Atype_i => identPtrListRow(i)%this
          call do_postpair(i, Atype_i)
       enddo
    endif

#else

    if (system_has_postpair_atoms) then
       do i = 1, natoms
          Atype_i => identPtr(i)%this
          call do_postpair(i, Atype_i)
       enddo
    endif

#endif


#ifdef IS_MPI
    !!distribute forces

    call scatter(fRow,fTemp1,plan_row3d)
    call scatter(fCol,fTemp2,plan_col3d)


    do i = 1,nlocal
       fTemp(1:3,i) = fTemp1(1:3,i) + fTemp2(1:3,i)
    end do

    if (do_torque) then
       call scatter(tRow,tTemp1,plan_row3d)
       call scatter(tCol,tTemp2,plan_col3d)
    
       do i = 1,nlocal
          tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(1:3,i)
       end do
    endif

    if (do_pot) then
       ! scatter/gather pot_row into the members of my column
       call scatter(eRow,eTemp,plan_row)
       
       ! scatter/gather pot_local into all other procs
       ! add resultant to get total pot
       do i = 1, nlocal
          pe_local = pe_local + eTemp(i)
       enddo

       eTemp = 0.0E0_DP
       call scatter(eCol,eTemp,plan_col)
       do i = 1, nlocal
          pe_local = pe_local + eTemp(i)
       enddo
       
       pe = pe_local
    endif
#else
! Copy local array into return array for c
    f = f+fTemp
    t = t+tTemp
#endif

    potE = pe


    if (do_stress) then
#ifdef IS_MPI
       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
            mpi_comm_world,mpi_err)
#else
       tau = tauTemp
#endif       
    endif

  end subroutine do_force_loop


!! Calculate any pre-force loop components and update nlist if necessary.
  subroutine do_preForce(updateNlist)
    logical, intent(inout) :: updateNlist


  end subroutine do_preForce


!! Calculate any post force loop components, i.e. reaction field, etc.
  subroutine do_postForce()


  end subroutine do_postForce


  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
    type (atype ), pointer, intent(inout) :: atype_i
    type (atype ), pointer, intent(inout) :: atype_j
    integer :: i
    integer :: j
    real ( kind = dp ), intent(inout) :: rx_ij
    real ( kind = dp ), intent(inout) :: ry_ij
    real ( kind = dp ), intent(inout) :: rz_ij


    real( kind = dp ) :: fx = 0.0_dp 
    real( kind = dp ) :: fy = 0.0_dp
    real( kind = dp ) :: fz = 0.0_dp  
   
    real( kind = dp ) ::  drdx = 0.0_dp
    real( kind = dp ) ::  drdy = 0.0_dp
    real( kind = dp ) ::  drdz = 0.0_dp
    

    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
    endif

    if (Atype_i%is_dp .and. Atype_j%is_dp) then

#ifdef IS_MPI
       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
            ulRow(:,i), ulCol(:,j), rt, rrf, pot)
#else
       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
            ul(:,i), ul(:,j), rt, rrf, pot)
#endif

       if (do_reaction_field) then
#ifdef IS_MPI
          call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
               ulRow(:i), ulCol(:,j), rt, rrf)
#else
          call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
               ul(:,i), ul(:,j), rt, rrf)
#endif
       endif


    endif

    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
       call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
    endif

       
#ifdef IS_MPI
                eRow(i) = eRow(i) + pot*0.5
                eCol(i) = eCol(i) + pot*0.5
#else
                    pe = pe + pot 
#endif                 
             
                drdx = -rxij / r
                drdy = -ryij / r
                drdz = -rzij / r
                
                fx = dudr * drdx
                fy = dudr * drdy
                fz = dudr * drdz


#ifdef IS_MPI
                fCol(1,j) = fCol(1,j) - fx 
                fCol(2,j) = fCol(2,j) - fy
                fCol(3,j) = fCol(3,j) - fz
                
                fRow(1,j) = fRow(1,j) + fx 
                fRow(2,j) = fRow(2,j) + fy
                fRow(3,j) = fRow(3,j) + fz
#else
                fTemp(1,j) = fTemp(1,j) - fx
                fTemp(2,j) = fTemp(2,j) - fy
                fTemp(3,j) = fTemp(3,j) - fz
                fTemp(1,i) = fTemp(1,i) + fx
                fTemp(2,i) = fTemp(2,i) + fy
                fTemp(3,i) = fTemp(3,i) + fz
#endif
                
                if (do_stress) then
                   tauTemp(1) = tauTemp(1) + fx * rxij
                   tauTemp(2) = tauTemp(2) + fx * ryij
                   tauTemp(3) = tauTemp(3) + fx * rzij
                   tauTemp(4) = tauTemp(4) + fy * rxij
                   tauTemp(5) = tauTemp(5) + fy * ryij
                   tauTemp(6) = tauTemp(6) + fy * rzij
                   tauTemp(7) = tauTemp(7) + fz * rxij
                   tauTemp(8) = tauTemp(8) + fz * ryij
                   tauTemp(9) = tauTemp(9) + fz * rzij
                endif


  end subroutine do_pair


  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
!---------------- Arguments------------------------------- 
   !! index i

    !! Position array
    real (kind = dp), dimension(3) :: q_i
    real (kind = dp), dimension(3) :: q_j
    !! x component of vector between i and j
    real ( kind = dp ), intent(out)  :: rx_ij 
    !! y component of vector between i and j
    real ( kind = dp ), intent(out)  :: ry_ij 
    !! z component of vector between i and j
    real ( kind = dp ), intent(out)  :: rz_ij 
    !! magnitude of r squared
    real ( kind = dp ), intent(out) :: r_sq
    !! magnitude of vector r between atoms i and j.
    real ( kind = dp ), intent(out) :: r_ij
    !! wrap into periodic box.
    logical, intent(in) :: wrap

!--------------- Local Variables---------------------------
    !! Distance between i and j
    real( kind = dp ) :: d(3)
!---------------- END DECLARATIONS-------------------------


! Find distance between i and j
    d(1:3) = q_i(1:3) - q_j(1:3)

! Wrap back into periodic box if necessary
    if ( wrap ) then
       d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
            int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
    end if
    
!   Find Magnitude of the vector
    r_sq = dot_product(d,d)
    r_ij = sqrt(r_sq)

!   Set each component for force calculation
    rx_ij = d(1)
    ry_ij = d(2)
    rz_ij = d(3)


  end subroutine get_interatomic_vector

  subroutine zero_module_variables()

#ifndef IS_MPI

    pe = 0.0E0_DP
    tauTemp = 0.0_dp
    fTemp = 0.0_dp
    tTemp = 0.0_dp
#else
    qRow = 0.0_dp
    qCol = 0.0_dp
    
    muRow = 0.0_dp
    muCol = 0.0_dp
    
    u_lRow = 0.0_dp
    u_lCol = 0.0_dp
    
    ARow = 0.0_dp
    ACol = 0.0_dp
    
    fRow = 0.0_dp
    fCol = 0.0_dp
    
  
    tRow = 0.0_dp
    tCol = 0.0_dp
 
  
    eRow = 0.0_dp
    eCol = 0.0_dp
    eTemp = 0.0_dp
#endif

  end subroutine zero_module_variables

#ifdef IS_MPI
!! Function to properly build neighbor lists in MPI using newtons 3rd law.
!! We don't want 2 processors doing the same i j pair twice.
!! Also checks to see if i and j are the same particle.
  function mpi_cycle_jLoop(i,j) result(do_cycle)
!--------------- Arguments--------------------------
! Index i
    integer,intent(in) :: i
! Index j
    integer,intent(in) :: j
! Result do_cycle
    logical :: do_cycle
!--------------- Local variables--------------------
    integer :: tag_i
    integer :: tag_j
!--------------- END DECLARATIONS------------------    
    tag_i = tagRow(i)
    tag_j = tagColumn(j)

    do_cycle = .false.

    if (tag_i == tag_j) then
       do_cycle = .true.
       return
    end if

    if (tag_i < tag_j) then
       if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
       return
    else                
       if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
    endif
  end function mpi_cycle_jLoop
#endif

end module do_Forces
Revision:	297
Committed:	Thu Mar 6 22:08:29 2003 UTC (21 years, 6 months ago) by gezelter
File size:	23620 byte(s)
Log Message:	Dan Goes Crazy