mdtools/libmdCode/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.18 2003-03-13 15:28:43 gezelter Exp $, $Date: 2003-03-13 15:28:43 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $


module do_Forces
  use simulation
  use definitions
  use atype_module
  use neighborLists  
  use lj
  use sticky_pair
  use dipole_dipole
  use reaction_field

#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

  logical, save :: do_forces_initialized = .false.
  logical, save :: FF_uses_LJ
  logical, save :: FF_uses_sticky
  logical, save :: FF_uses_dipoles
  logical, save :: FF_uses_RF
  logical, save :: FF_uses_GB
  logical, save :: FF_uses_EAM

  public :: init_FF
  public :: do_force_loop

contains

  subroutine init_FF(LJ_mix_policy, use_RF_c, thisStat)
    logical(kind = 2), intent(in) :: use_RF_c
    logical :: use_RF_f
    integer, intent(out) :: thisStat   
    integer :: my_status, nMatches
    character(len = 100) :: LJ_mix_Policy
    integer, pointer :: MatchList(:)
    
    !! Fortran's version of a cast:
    use_RF_f = use_RF_c

    !! assume things are copacetic, unless they aren't
    thisStat = 0
    
    !! 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_dipoles = .false.
    FF_uses_GB = .false.
    FF_uses_EAM = .false.
    
    call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
    deallocate(MatchList)
    if (nMatches .gt. 0) FF_uses_LJ = .true.

    call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
    deallocate(MatchList)
    if (nMatches .gt. 0) FF_uses_dipoles = .true.

    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
         MatchList) 
    deallocate(MatchList)
    if (nMatches .gt. 0) FF_uses_Sticky = .true.
    
    call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
    deallocate(MatchList)
    if (nMatches .gt. 0) FF_uses_GB = .true.

    call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
    deallocate(MatchList)
    if (nMatches .gt. 0) FF_uses_EAM = .true.

    !! check to make sure the use_RF setting makes sense
    if (use_RF_f) then
       if (FF_uses_dipoles) then
          FF_uses_RF = .true.
          call initialize_rf()
       else
          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
          thisStat = -1
          return
       endif
    endif
    
    call init_lj_FF(LJ_mix_Policy, my_status)
    if (my_status /= 0) then
       thisStat = -1
       return
    end if
    
    call check_sticky_FF(my_status)
    if (my_status /= 0) then
       thisStat = -1
       return
    end if
    
    do_forces_initialized = .true.    
    
  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,getNlocal()) :: q
    !! Rotation Matrix for each long range particle in simulation.
    real( kind = dp), dimension(9,getNlocal()) :: A    
    !! 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 ) :: 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
#endif
    integer :: nlocal
    integer :: natoms    
    logical :: update_nlist   
    integer :: i, j, jbeg, jend, jnab
    integer :: nlist
    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
    real(kind=dp),dimension(3) :: d
    real(kind=dp) :: rfpot, mu_i, virial
    integer :: me_i
    logical :: is_dp_i
    integer :: neighborListSize
    integer :: listerror, error
    integer :: localError

    !! initialize local variables  

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

    call getRcut(rcut,rc2=rcutsq)
    call getRlist(rlist,rlistsq)
    
    call check_initialization(localError)
    if ( localError .ne. 0 ) then
       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. SimUsesDirectionalAtoms()) 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
    
    if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
       !! See if we need to update neighbor lists
       call checkNeighborList(nlocal, q, rcut, rlist, 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
    else
       !! See if we need to update neighbor lists
       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)   
    endif
    
#ifdef IS_MPI
    
    if (update_nlist) then
       
       !! save current configuration, construct neighbor list, 
       !! and calculate forces
       call save_neighborList(q)
       
       neighborListSize = getNeighborListSize()
       nlist = 0       
       
       do i = 1, nrow
          point(i) = nlist + 1
          
          inner: do j = 1, ncol
             
             if (checkExcludes(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
                endif
                
                list(nlist) = j
                                
                if (rijsq <  rcutsq) then
                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
                endif
             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)

             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
          
          inner: do j = i+1, natoms
             
             if (checkExcludes(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 expandList(natoms, listerror)
                   if (listerror /= 0) then
                      error = -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(i, j, rijsq, d, do_pot, do_stress)
                endif
             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)

             enddo
          endif
       enddo
    endif
    
#endif 
    
    ! phew, done with main loop.
    
#ifdef IS_MPI
    !!distribute forces
    
    call scatter(f_Row,f,plan_row3d)
    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. SimUsesDirectionalAtoms()) then
       call scatter(t_Row,t,plan_row3d)
       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. SimRequiresPostpairCalc()) then
       
       if (FF_uses_RF .and. SimUsesRF()) 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, getNlocal()

             rfpot = 0.0_DP
#ifdef IS_MPI
             me_i = atid_row(i)
#else
             me_i = atid(i)
#endif
             call getElementProperty(atypes, me_i, "is_DP", is_DP_i)       
             if ( is_DP_i ) then
                call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
                !! 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_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, tau_Temp,9,mpi_double_precision,mpi_sum, &
            mpi_comm_world,mpi_err)
       call mpi_allreduce(virial, virial_Temp,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)

    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3,getNlocal()) :: u_l
    real (kind=dp), dimension(9,getNlocal()) :: A
    real (kind=dp), dimension(3,getNlocal()) :: f
    real (kind=dp), dimension(3,getNlocal()) :: 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_LJ_i, is_LJ_j
    logical :: is_DP_i, is_DP_j
    logical :: is_Sticky_i, is_Sticky_j
    integer :: me_i, me_j

    r = sqrt(rijsq)
    
#ifdef IS_MPI

    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. SimUsesLJ()) then
       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
       
       if ( is_LJ_i .and. is_LJ_j ) &
            call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
    endif
       

    if (FF_uses_dipoles .and. SimUsesDipoles()) then
       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
       
       if ( is_DP_i .and. is_DP_j ) then
          
          call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
          
          if (FF_uses_RF .and. SimUsesRF()) 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. SimUsesSticky()) then

       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
       
       if ( is_Sticky_i .and. is_Sticky_j ) then
          call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
               do_pot, do_stress)
       endif
    endif
      
  end subroutine do_pair


  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)

    d(1:3) = q_i(1:3) - q_j(1:3)
    
    ! Wrap back into periodic box if necessary
    if ( SimUsesPBC() ) then
       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,box(1:3)) * &
            int(abs(d(1:3)/box(1:3) + 0.5_dp))
    endif
    
    r_sq = dot_product(d,d)
        
  end subroutine get_interatomic_vector

  subroutine check_initialization(error)
    integer, intent(out) :: error
    
    error = 0
    ! Make sure we are properly initialized.
    if (.not. do_forces_initialized) then
       write(default_error,*) "ERROR: do_Forces has not been initialized!"
       error = -1
       return
    endif

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

  
  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

    rf = 0.0_dp
    tau_Temp = 0.0_dp
    virial_Temp = 0.0_dp
    
  end subroutine zero_work_arrays
  

  !! 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 checkExcludes(atom1,atom2) result(do_cycle)
    
    integer,intent(in) :: atom1
    integer,intent(in), optional :: atom2
    logical :: do_cycle
    integer :: unique_id_1, unique_id_2
    integer :: i, j

    do_cycle = .false. 
    
#ifdef IS_MPI
    unique_id_1 = tagRow(atom1)
#else
    unique_id_1 = tag(atom1)
#endif
    
    !! Check global excludes first
    if (.not. present(atom2)) then
       do i = 1, nExcludes_global
          if (excludesGlobal(i) == unique_id_1) then
             do_cycle = .true.
             return
          end if
       end do
       return !! return after checking globals
    end if

    !! we return if atom2 not present here.
    
#ifdef IS_MPI
    unique_id_2 = tagColumn(atom2)
#else
    unique_id_2 = tag(atom2)
#endif
    
    if (unique_id_1 == unique_id_2) then
       do_cycle = .true.
       return
    end if
    
    if (unique_id_1 < unique_id_2) then
       if (mod(unique_id_1 + unique_id_2,2) == 0) do_cycle = .true.
       return
    else                
       if (mod(unique_id_1 + unique_id_2,2) == 1) do_cycle = .true.
    endif
    
    do i = 1, nExcludes_local
       if ((unique_id_1 == excludesLocal(1,i)) .and.  &
            (excludesLocal(2,i) < 0)) then
          do_cycle = .true.
          return
       end if
    end do
    
  end function checkExcludes

  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
  
end module do_Forces
Revision:	332
Committed:	Thu Mar 13 15:28:43 2003 UTC (21 years, 5 months ago) by gezelter
File size:	18250 byte(s)
Log Message:	initialization routines, bug fixes, exclude lists
#	User	Rev	Content
1	chuckv	306	!! do_Forces.F90
2			!! module do_Forces
3	chuckv	292	!! Calculates Long Range forces.
4	chuckv	306
5	chuckv	292	!! @author Charles F. Vardeman II
6			!! @author Matthew Meineke
7	gezelter	332	!! @version $Id: do_Forces.F90,v 1.18 2003-03-13 15:28:43 gezelter Exp $, $Date: 2003-03-13 15:28:43 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
8	chuckv	292
9
10
11			module do_Forces
12			use simulation
13			use definitions
14	gezelter	328	use atype_module
15	gezelter	325	use neighborLists
16	gezelter	330	use lj
17			use sticky_pair
18	gezelter	309	use dipole_dipole
19	gezelter	332	use reaction_field
20	chuckv	292
21			#ifdef IS_MPI
22			use mpiSimulation
23			#endif
24			implicit none
25			PRIVATE
26
27	gezelter	329	logical, save :: do_forces_initialized = .false.
28			logical, save :: FF_uses_LJ
29			logical, save :: FF_uses_sticky
30			logical, save :: FF_uses_dipoles
31			logical, save :: FF_uses_RF
32			logical, save :: FF_uses_GB
33			logical, save :: FF_uses_EAM
34	chuckv	292
35	gezelter	329	public :: init_FF
36	gezelter	330	public :: do_force_loop
37	gezelter	329
38	chuckv	292	contains
39
40	gezelter	332	subroutine init_FF(LJ_mix_policy, use_RF_c, thisStat)
41			logical(kind = 2), intent(in) :: use_RF_c
42			logical :: use_RF_f
43	chuckv	331	integer, intent(out) :: thisStat
44	gezelter	332	integer :: my_status, nMatches
45			character(len = 100) :: LJ_mix_Policy
46			integer, pointer :: MatchList(:)
47
48			!! Fortran's version of a cast:
49			use_RF_f = use_RF_c
50	gezelter	329
51	gezelter	332	!! assume things are copacetic, unless they aren't
52	chuckv	331	thisStat = 0
53	gezelter	332
54			!! init_FF is called after all of the atom types have been
55			!! defined in atype_module using the new_atype subroutine.
56			!!
57			!! this will scan through the known atypes and figure out what
58			!! interactions are used by the force field.
59
60			FF_uses_LJ = .false.
61			FF_uses_sticky = .false.
62			FF_uses_dipoles = .false.
63			FF_uses_GB = .false.
64			FF_uses_EAM = .false.
65
66			call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
67			deallocate(MatchList)
68			if (nMatches .gt. 0) FF_uses_LJ = .true.
69
70			call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
71			deallocate(MatchList)
72			if (nMatches .gt. 0) FF_uses_dipoles = .true.
73
74			call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
75			MatchList)
76			deallocate(MatchList)
77			if (nMatches .gt. 0) FF_uses_Sticky = .true.
78
79			call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
80			deallocate(MatchList)
81			if (nMatches .gt. 0) FF_uses_GB = .true.
82
83			call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
84			deallocate(MatchList)
85			if (nMatches .gt. 0) FF_uses_EAM = .true.
86
87			!! check to make sure the use_RF setting makes sense
88			if (use_RF_f) then
89			if (FF_uses_dipoles) then
90			FF_uses_RF = .true.
91			call initialize_rf()
92			else
93			write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
94			thisStat = -1
95			return
96			endif
97			endif
98
99			call init_lj_FF(LJ_mix_Policy, my_status)
100	gezelter	329	if (my_status /= 0) then
101			thisStat = -1
102			return
103			end if
104
105			call check_sticky_FF(my_status)
106			if (my_status /= 0) then
107			thisStat = -1
108			return
109			end if
110
111			do_forces_initialized = .true.
112
113			end subroutine init_FF
114
115
116
117	gezelter	317	!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
118			!------------------------------------------------------------->
119	gezelter	325	subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
120	gezelter	329	error)
121	gezelter	317	!! Position array provided by C, dimensioned by getNlocal
122	chuckv	292	real ( kind = dp ), dimension(3,getNlocal()) :: q
123	gezelter	317	!! Rotation Matrix for each long range particle in simulation.
124	gezelter	325	real( kind = dp), dimension(9,getNlocal()) :: A
125	gezelter	317	!! Unit vectors for dipoles (lab frame)
126	chuckv	292	real( kind = dp ), dimension(3,getNlocal()) :: u_l
127	gezelter	317	!! Force array provided by C, dimensioned by getNlocal
128	chuckv	292	real ( kind = dp ), dimension(3,getNlocal()) :: f
129	gezelter	317	!! Torsion array provided by C, dimensioned by getNlocal
130	gezelter	325	real( kind = dp ), dimension(3,getNlocal()) :: t
131	gezelter	317	!! Stress Tensor
132	gezelter	325	real( kind = dp), dimension(9) :: tau
133			real ( kind = dp ) :: pot
134			logical ( kind = 2) :: do_pot_c, do_stress_c
135	gezelter	329	logical :: do_pot
136			logical :: do_stress
137	gezelter	317	#ifdef IS_MPI
138			real( kind = DP ) :: pot_local
139	gezelter	329	integer :: nrow
140			integer :: ncol
141	gezelter	317	#endif
142	gezelter	330	integer :: nlocal
143	gezelter	329	integer :: natoms
144	gezelter	325	logical :: update_nlist
145			integer :: i, j, jbeg, jend, jnab
146	gezelter	317	integer :: nlist
147	gezelter	325	real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut
148	gezelter	330	real(kind=dp),dimension(3) :: d
149			real(kind=dp) :: rfpot, mu_i, virial
150			integer :: me_i
151			logical :: is_dp_i
152	gezelter	317	integer :: neighborListSize
153	gezelter	330	integer :: listerror, error
154	chuckv	331	integer :: localError
155	chuckv	292
156	gezelter	329	!! initialize local variables
157	gezelter	325
158	chuckv	292	#ifdef IS_MPI
159	gezelter	329	nlocal = getNlocal()
160			nrow = getNrow(plan_row)
161			ncol = getNcol(plan_col)
162			#else
163			nlocal = getNlocal()
164			natoms = nlocal
165	chuckv	292	#endif
166	gezelter	329
167	chuckv	331	call getRcut(rcut,rc2=rcutsq)
168	gezelter	317	call getRlist(rlist,rlistsq)
169
170	chuckv	331	call check_initialization(localError)
171			if ( localError .ne. 0 ) then
172			error = -1
173			return
174			end if
175	gezelter	329	call zero_work_arrays()
176
177			do_pot = do_pot_c
178			do_stress = do_stress_c
179
180			! Gather all information needed by all force loops:
181	gezelter	317
182	gezelter	329	#ifdef IS_MPI
183	chuckv	292
184	gezelter	309	call gather(q,q_Row,plan_row3d)
185			call gather(q,q_Col,plan_col3d)
186	gezelter	329
187			if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
188			call gather(u_l,u_l_Row,plan_row3d)
189			call gather(u_l,u_l_Col,plan_col3d)
190
191			call gather(A,A_Row,plan_row_rotation)
192			call gather(A,A_Col,plan_col_rotation)
193			endif
194	gezelter	317
195	gezelter	329	#endif
196	gezelter	317
197	gezelter	329	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
198			!! See if we need to update neighbor lists
199			call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
200			!! if_mpi_gather_stuff_for_prepair
201			!! do_prepair_loop_if_needed
202			!! if_mpi_scatter_stuff_from_prepair
203			!! if_mpi_gather_stuff_from_prepair_to_main_loop
204			else
205			!! See if we need to update neighbor lists
206			call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
207			endif
208
209	gezelter	297	#ifdef IS_MPI
210	chuckv	292
211	gezelter	297	if (update_nlist) then
212
213	gezelter	329	!! save current configuration, construct neighbor list,
214			!! and calculate forces
215	gezelter	297	call save_neighborList(q)
216
217			neighborListSize = getNeighborListSize()
218	gezelter	329	nlist = 0
219	gezelter	297
220			do i = 1, nrow
221			point(i) = nlist + 1
222
223			inner: do j = 1, ncol
224
225	chuckv	331	if (checkExcludes(i,j)) cycle inner
226	gezelter	329
227	gezelter	317	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
228	gezelter	330
229	gezelter	297	if (rijsq < rlistsq) then
230
231			nlist = nlist + 1
232
233			if (nlist > neighborListSize) then
234	gezelter	325	call expandNeighborList(nlocal, listerror)
235	gezelter	297	if (listerror /= 0) then
236	gezelter	329	error = -1
237	gezelter	297	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
238			return
239			end if
240			endif
241
242			list(nlist) = j
243	gezelter	317
244	gezelter	297	if (rijsq < rcutsq) then
245	gezelter	329	call do_pair(i, j, rijsq, d, do_pot, do_stress)
246	gezelter	297	endif
247			endif
248			enddo inner
249			enddo
250	chuckv	292
251	gezelter	297	point(nrow + 1) = nlist + 1
252
253	gezelter	329	else !! (of update_check)
254	chuckv	292
255	gezelter	297	! use the list to find the neighbors
256			do i = 1, nrow
257			JBEG = POINT(i)
258			JEND = POINT(i+1) - 1
259			! check thiat molecule i has neighbors
260			if (jbeg .le. jend) then
261	gezelter	317
262	gezelter	297	do jnab = jbeg, jend
263			j = list(jnab)
264	gezelter	317
265			call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
266	gezelter	329	call do_pair(i, j, rijsq, d, do_pot, do_stress)
267	gezelter	317
268	gezelter	297	enddo
269			endif
270			enddo
271			endif
272	gezelter	329
273	gezelter	297	#else
274
275			if (update_nlist) then
276	chuckv	292
277	gezelter	297	! save current configuration, contruct neighbor list,
278			! and calculate forces
279			call save_neighborList(q)
280
281			neighborListSize = getNeighborListSize()
282			nlist = 0
283	gezelter	329
284	gezelter	297	do i = 1, natoms-1
285			point(i) = nlist + 1
286	gezelter	329
287	gezelter	297	inner: do j = i+1, natoms
288	gezelter	329
289	chuckv	331	if (checkExcludes(i,j)) cycle inner
290	gezelter	329
291	gezelter	317	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
292	chuckv	295
293	gezelter	297	if (rijsq < rlistsq) then
294	gezelter	317
295	gezelter	297	nlist = nlist + 1
296
297			if (nlist > neighborListSize) then
298	gezelter	325	call expandList(natoms, listerror)
299	gezelter	297	if (listerror /= 0) then
300	gezelter	329	error = -1
301	gezelter	297	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
302			return
303			end if
304			endif
305
306			list(nlist) = j
307	gezelter	329
308	gezelter	297	if (rijsq < rcutsq) then
309	gezelter	329	call do_pair(i, j, rijsq, d, do_pot, do_stress)
310	gezelter	297	endif
311			endif
312	chuckv	292	enddo inner
313	gezelter	297	enddo
314
315			point(natoms) = nlist + 1
316
317			else !! (update)
318
319			! use the list to find the neighbors
320	gezelter	330	do i = 1, natoms-1
321	gezelter	297	JBEG = POINT(i)
322			JEND = POINT(i+1) - 1
323			! check thiat molecule i has neighbors
324			if (jbeg .le. jend) then
325
326			do jnab = jbeg, jend
327			j = list(jnab)
328	gezelter	317
329			call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
330	gezelter	329	call do_pair(i, j, rijsq, d, do_pot, do_stress)
331	gezelter	317
332	gezelter	297	enddo
333			endif
334			enddo
335			endif
336	gezelter	317
337	gezelter	297	#endif
338	gezelter	317
339	gezelter	329	! phew, done with main loop.
340	gezelter	317
341	chuckv	292	#ifdef IS_MPI
342	gezelter	329	!!distribute forces
343	gezelter	317
344	gezelter	309	call scatter(f_Row,f,plan_row3d)
345			call scatter(f_Col,f_temp,plan_col3d)
346	chuckv	295	do i = 1,nlocal
347	gezelter	309	f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
348	chuckv	295	end do
349	gezelter	329
350			if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
351	gezelter	309	call scatter(t_Row,t,plan_row3d)
352			call scatter(t_Col,t_temp,plan_col3d)
353	gezelter	329
354	chuckv	295	do i = 1,nlocal
355	gezelter	309	t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
356	chuckv	295	end do
357			endif
358	gezelter	309
359	chuckv	295	if (do_pot) then
360			! scatter/gather pot_row into the members of my column
361	gezelter	309	call scatter(pot_Row, pot_Temp, plan_row)
362	chuckv	295
363			! scatter/gather pot_local into all other procs
364			! add resultant to get total pot
365			do i = 1, nlocal
366	gezelter	309	pot_local = pot_local + pot_Temp(i)
367	chuckv	295	enddo
368
369	gezelter	309	pot_Temp = 0.0_DP
370
371			call scatter(pot_Col, pot_Temp, plan_col)
372	chuckv	295	do i = 1, nlocal
373	gezelter	309	pot_local = pot_local + pot_Temp(i)
374	chuckv	295	enddo
375
376	gezelter	330	endif
377			#endif
378
379	gezelter	329	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
380
381			if (FF_uses_RF .and. SimUsesRF()) then
382
383			#ifdef IS_MPI
384			call scatter(rf_Row,rf,plan_row3d)
385			call scatter(rf_Col,rf_Temp,plan_col3d)
386			do i = 1,nlocal
387			rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
388			end do
389			#endif
390
391			do i = 1, getNlocal()
392
393	gezelter	330	rfpot = 0.0_DP
394	gezelter	329	#ifdef IS_MPI
395			me_i = atid_row(i)
396			#else
397			me_i = atid(i)
398			#endif
399			call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
400			if ( is_DP_i ) then
401			call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
402			!! The reaction field needs to include a self contribution
403			!! to the field:
404			call accumulate_self_rf(i, mu_i, u_l)
405			!! Get the reaction field contribution to the
406			!! potential and torques:
407	gezelter	330	call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
408	gezelter	329	#ifdef IS_MPI
409			pot_local = pot_local + rfpot
410			#else
411			pot = pot + rfpot
412			#endif
413			endif
414			enddo
415			endif
416			endif
417
418
419			#ifdef IS_MPI
420
421			if (do_pot) then
422	gezelter	309	pot = pot_local
423	gezelter	329	!! we assume the c code will do the allreduce to get the total potential
424			!! we could do it right here if we needed to...
425	chuckv	295	endif
426
427	gezelter	329	if (do_stress) then
428	gezelter	330	call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
429	gezelter	309	mpi_comm_world,mpi_err)
430	gezelter	330	call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
431	gezelter	309	mpi_comm_world,mpi_err)
432			endif
433	chuckv	295
434	gezelter	329	#else
435	chuckv	295
436	gezelter	329	if (do_stress) then
437	gezelter	309	tau = tau_Temp
438			virial = virial_Temp
439	chuckv	295	endif
440
441	gezelter	329	#endif
442
443	chuckv	295	end subroutine do_force_loop
444
445	gezelter	329	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress)
446	chuckv	295
447	gezelter	330	real( kind = dp ) :: pot
448			real( kind = dp ), dimension(3,getNlocal()) :: u_l
449			real (kind=dp), dimension(9,getNlocal()) :: A
450			real (kind=dp), dimension(3,getNlocal()) :: f
451			real (kind=dp), dimension(3,getNlocal()) :: t
452
453	gezelter	329	logical, intent(inout) :: do_pot, do_stress
454	gezelter	317	integer, intent(in) :: i, j
455	chuckv	331	real ( kind = dp ), intent(inout) :: rijsq
456	gezelter	317	real ( kind = dp ) :: r
457			real ( kind = dp ), intent(inout) :: d(3)
458			logical :: is_LJ_i, is_LJ_j
459			logical :: is_DP_i, is_DP_j
460			logical :: is_Sticky_i, is_Sticky_j
461			integer :: me_i, me_j
462	chuckv	295
463	gezelter	330	r = sqrt(rijsq)
464
465	gezelter	302	#ifdef IS_MPI
466
467	gezelter	317	me_i = atid_row(i)
468			me_j = atid_col(j)
469	chuckv	295
470	gezelter	317	#else
471	chuckv	295
472	gezelter	317	me_i = atid(i)
473			me_j = atid(j)
474	gezelter	302
475	gezelter	317	#endif
476	gezelter	302
477
478	gezelter	329	if (FF_uses_LJ .and. SimUsesLJ()) then
479			call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
480			call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
481
482			if ( is_LJ_i .and. is_LJ_j ) &
483			call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
484			endif
485
486	gezelter	302
487	gezelter	330	if (FF_uses_dipoles .and. SimUsesDipoles()) then
488	gezelter	329	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
489			call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
490
491			if ( is_DP_i .and. is_DP_j ) then
492
493			call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
494
495			if (FF_uses_RF .and. SimUsesRF()) then
496
497			call accumulate_rf(i, j, r, u_l)
498			call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
499
500			endif
501
502	gezelter	297	endif
503			endif
504
505	gezelter	329	if (FF_uses_Sticky .and. SimUsesSticky()) then
506	gezelter	317
507	gezelter	329	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
508			call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
509
510			if ( is_Sticky_i .and. is_Sticky_j ) then
511			call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
512			do_pot, do_stress)
513			endif
514	gezelter	297	endif
515	gezelter	309
516	chuckv	295	end subroutine do_pair
517	chuckv	292
518
519	gezelter	317	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
520
521	chuckv	295	real (kind = dp), dimension(3) :: q_i
522			real (kind = dp), dimension(3) :: q_j
523			real ( kind = dp ), intent(out) :: r_sq
524			real( kind = dp ) :: d(3)
525
526			d(1:3) = q_i(1:3) - q_j(1:3)
527	gezelter	317
528			! Wrap back into periodic box if necessary
529	gezelter	330	if ( SimUsesPBC() ) then
530			d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,box(1:3)) * &
531			int(abs(d(1:3)/box(1:3) + 0.5_dp))
532	gezelter	317	endif
533	chuckv	292
534	chuckv	295	r_sq = dot_product(d,d)
535	gezelter	317
536	chuckv	295	end subroutine get_interatomic_vector
537	gezelter	329
538			subroutine check_initialization(error)
539			integer, intent(out) :: error
540	gezelter	332
541	gezelter	329	error = 0
542			! Make sure we are properly initialized.
543	gezelter	332	if (.not. do_forces_initialized) then
544	gezelter	329	write(default_error,*) "ERROR: do_Forces has not been initialized!"
545			error = -1
546			return
547			endif
548	gezelter	332
549	gezelter	329	#ifdef IS_MPI
550			if (.not. isMPISimSet()) then
551			write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
552			error = -1
553			return
554			endif
555			#endif
556	gezelter	332
557	gezelter	329	return
558			end subroutine check_initialization
559
560	gezelter	317
561	gezelter	325	subroutine zero_work_arrays()
562
563			#ifdef IS_MPI
564	chuckv	292
565	gezelter	325	q_Row = 0.0_dp
566			q_Col = 0.0_dp
567	chuckv	295
568	gezelter	325	u_l_Row = 0.0_dp
569			u_l_Col = 0.0_dp
570	chuckv	295
571	gezelter	325	A_Row = 0.0_dp
572			A_Col = 0.0_dp
573	chuckv	295
574	gezelter	325	f_Row = 0.0_dp
575			f_Col = 0.0_dp
576			f_Temp = 0.0_dp
577
578			t_Row = 0.0_dp
579			t_Col = 0.0_dp
580			t_Temp = 0.0_dp
581	chuckv	295
582	gezelter	325	pot_Row = 0.0_dp
583			pot_Col = 0.0_dp
584			pot_Temp = 0.0_dp
585	chuckv	292
586	gezelter	332	rf_Row = 0.0_dp
587			rf_Col = 0.0_dp
588			rf_Temp = 0.0_dp
589
590	chuckv	295	#endif
591	chuckv	292
592	gezelter	332	rf = 0.0_dp
593	gezelter	325	tau_Temp = 0.0_dp
594			virial_Temp = 0.0_dp
595
596			end subroutine zero_work_arrays
597
598	chuckv	292
599	gezelter	330	!! Function to properly build neighbor lists in MPI using newtons 3rd law.
600			!! We don't want 2 processors doing the same i j pair twice.
601			!! Also checks to see if i and j are the same particle.
602
603	chuckv	306	function checkExcludes(atom1,atom2) result(do_cycle)
604	gezelter	332
605	chuckv	306	integer,intent(in) :: atom1
606			integer,intent(in), optional :: atom2
607	chuckv	295	logical :: do_cycle
608	gezelter	332	integer :: unique_id_1, unique_id_2
609	gezelter	330	integer :: i, j
610	gezelter	332
611	chuckv	306	do_cycle = .false.
612	gezelter	330
613	chuckv	306	#ifdef IS_MPI
614	gezelter	332	unique_id_1 = tagRow(atom1)
615	chuckv	306	#else
616	gezelter	332	unique_id_1 = tag(atom1)
617	chuckv	306	#endif
618	gezelter	330
619			!! Check global excludes first
620	chuckv	306	if (.not. present(atom2)) then
621	gezelter	330	do i = 1, nExcludes_global
622	gezelter	332	if (excludesGlobal(i) == unique_id_1) then
623	chuckv	306	do_cycle = .true.
624			return
625			end if
626			end do
627			return !! return after checking globals
628			end if
629
630	gezelter	330	!! we return if atom2 not present here.
631
632	gezelter	332	#ifdef IS_MPI
633			unique_id_2 = tagColumn(atom2)
634			#else
635			unique_id_2 = tag(atom2)
636			#endif
637
638			if (unique_id_1 == unique_id_2) then
639	chuckv	295	do_cycle = .true.
640			return
641			end if
642	gezelter	330
643	gezelter	332	if (unique_id_1 < unique_id_2) then
644			if (mod(unique_id_1 + unique_id_2,2) == 0) do_cycle = .true.
645	chuckv	295	return
646			else
647	gezelter	332	if (mod(unique_id_1 + unique_id_2,2) == 1) do_cycle = .true.
648	chuckv	295	endif
649	gezelter	332
650	gezelter	330	do i = 1, nExcludes_local
651	gezelter	332	if ((unique_id_1 == excludesLocal(1,i)) .and. &
652			(excludesLocal(2,i) < 0)) then
653	chuckv	306	do_cycle = .true.
654			return
655			end if
656			end do
657	gezelter	330
658	chuckv	306	end function checkExcludes
659
660	gezelter	329	function FF_UsesDirectionalAtoms() result(doesit)
661			logical :: doesit
662			doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
663			FF_uses_GB .or. FF_uses_RF
664			end function FF_UsesDirectionalAtoms
665
666			function FF_RequiresPrepairCalc() result(doesit)
667			logical :: doesit
668			doesit = FF_uses_EAM
669			end function FF_RequiresPrepairCalc
670
671			function FF_RequiresPostpairCalc() result(doesit)
672			logical :: doesit
673			doesit = FF_uses_RF
674			end function FF_RequiresPostpairCalc
675
676	chuckv	292	end module do_Forces