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.13 2003-03-12 19:10:54 gezelter Exp $, $Date: 2003-03-12 19:10:54 $, $Name: not supported by cvs2svn $, $Revision: 1.13 $


module do_Forces
  use simulation
  use definitions
  use forceGlobals
  use atype_typedefs
  use neighborLists  
  use lj_FF
  use sticky_FF
  use dipole_dipole
  use gb_FF

#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

  public :: do_force_loop

  logical :: do_pot
  logical :: do_stress

contains

  !! 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, &
       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    
    !! 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
    integer :: FFerror

#ifdef IS_MPI
    real( kind = DP ) :: pot_local
#endif
    
    logical :: update_nlist   
    integer :: i, j, jbeg, jend, jnab
    integer :: nlist
    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut

#ifdef IS_MPI
    integer :: nlocal
#endif
    integer :: nrow
    integer :: ncol
    integer :: natoms
    integer :: neighborListSize
    integer :: listerror
    FFerror = 0

    do_pot = do_pot_c
    do_stress = do_stress_c

    ! 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)
    
    !! See if we need to update neighbor lists
    call checkNeighborList(natoms, q, rcut, rlist, 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,q_Row,plan_row3d)
    call gather(q,q_Col,plan_col3d)
    
    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


#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
          
          inner: do j = 1, ncol
             
             if (check_exclude(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
                      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(i, j, rijsq, d)
                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
             
             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)

             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 (check_exclude(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
                      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(i, j, rijsq, d)
                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
             
             do jnab = jbeg, jend
                j = list(jnab)

                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
                call do_pair(i, j, rijsq, d)

             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(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 (doTorque()) 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
       
       pot = pot_local
    endif

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

#endif

    if (doStress()) then
       tau = tau_Temp
       virial = virial_Temp
    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(i, j, rijsq, d)

    integer, intent(in) :: i, j
    real ( kind = dp ), intent(in)    :: rijsq
    real ( kind = dp )                :: r
    real ( kind = dp ), intent(inout) :: d(3)

    r = sqrt(rijsq)
    
    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

#ifdef IS_MPI

    me_i = atid_row(i)
    me_j = atid_col(j)

#else

    me_i = atid(i)
    me_j = atid(j)

#endif

    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)

    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 (do_reaction_field) then
          call accumulate_rf(i, j, r)
       endif

    endif

    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

      
  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 ( isPBC() ) 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)
    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

#endif

    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)
!--------------- Arguments--------------------------
! Index i
    integer,intent(in) :: atom1
! Index j
    integer,intent(in), optional :: atom2
! Result do_cycle
    logical :: do_cycle
!--------------- Local variables--------------------
    integer :: tag_i
    integer :: tag_j
    integer :: i
!--------------- END DECLARATIONS------------------   
    do_cycle = .false. 

#ifdef IS_MPI
    tag_i = tagRow(atom1)
#else
    tag_i = tag(atom1)
#endif

!! Check global excludes first
    if (.not. present(atom2)) then
       do i = 1,nGlobalExcludes
          if (excludeGlobal(i) == tag_i) then
             do_cycle = .true.
             return
          end if
       end do
       return !! return after checking globals
    end if

!! we return if j not present here.
    tag_j = tagColumn(j)

 
    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


    do i = 1, nLocalExcludes
       if (tag_i = excludes(1,i) .and. excludes(2,i) < 0) then
          do_cycle = .true.
          return
       end if
    end do
       

  end function checkExcludes


end module do_Forces
Revision:	325
Committed:	Wed Mar 12 19:10:54 2003 UTC (21 years, 5 months ago) by gezelter
File size:	13130 byte(s)
Log Message:	Massive rewrite
#	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	325	!! @version $Id: do_Forces.F90,v 1.13 2003-03-12 19:10:54 gezelter Exp $, $Date: 2003-03-12 19:10:54 $, $Name: not supported by cvs2svn $, $Revision: 1.13 $
8	chuckv	292
9
10
11			module do_Forces
12			use simulation
13			use definitions
14	chuckv	298	use forceGlobals
15			use atype_typedefs
16	gezelter	325	use neighborLists
17			use lj_FF
18	chuckv	292	use sticky_FF
19	gezelter	309	use dipole_dipole
20	chuckv	292	use gb_FF
21
22			#ifdef IS_MPI
23			use mpiSimulation
24			#endif
25			implicit none
26			PRIVATE
27
28	gezelter	325	public :: do_force_loop
29	chuckv	292
30	gezelter	325	logical :: do_pot
31			logical :: do_stress
32
33	chuckv	292	contains
34
35	gezelter	317	!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
36			!------------------------------------------------------------->
37	gezelter	325	subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
38			FFerror)
39	gezelter	317	!! Position array provided by C, dimensioned by getNlocal
40	chuckv	292	real ( kind = dp ), dimension(3,getNlocal()) :: q
41	gezelter	317	!! Rotation Matrix for each long range particle in simulation.
42	gezelter	325	real( kind = dp), dimension(9,getNlocal()) :: A
43	gezelter	317	!! Unit vectors for dipoles (lab frame)
44	chuckv	292	real( kind = dp ), dimension(3,getNlocal()) :: u_l
45	gezelter	317	!! Force array provided by C, dimensioned by getNlocal
46	chuckv	292	real ( kind = dp ), dimension(3,getNlocal()) :: f
47	gezelter	317	!! Torsion array provided by C, dimensioned by getNlocal
48	gezelter	325	real( kind = dp ), dimension(3,getNlocal()) :: t
49	gezelter	317	!! Stress Tensor
50	gezelter	325	real( kind = dp), dimension(9) :: tau
51			real ( kind = dp ) :: pot
52			logical ( kind = 2) :: do_pot_c, do_stress_c
53	chuckv	292	integer :: FFerror
54
55	gezelter	317	#ifdef IS_MPI
56			real( kind = DP ) :: pot_local
57			#endif
58	chuckv	292
59	gezelter	325	logical :: update_nlist
60			integer :: i, j, jbeg, jend, jnab
61	gezelter	317	integer :: nlist
62	gezelter	325	real( kind = DP ) :: rijsq, rlistsq, rcutsq, rlist, rcut
63	chuckv	292
64			#ifdef IS_MPI
65	gezelter	317	integer :: nlocal
66	chuckv	292	#endif
67	gezelter	317	integer :: nrow
68			integer :: ncol
69			integer :: natoms
70			integer :: neighborListSize
71			integer :: listerror
72			FFerror = 0
73	chuckv	292
74	gezelter	325	do_pot = do_pot_c
75			do_stress = do_stress_c
76
77	gezelter	317	! Make sure we are properly initialized.
78			if (.not. isFFInit) then
79			write(default_error,*) "ERROR: lj_FF has not been properly initialized"
80			FFerror = -1
81			return
82			endif
83	chuckv	292	#ifdef IS_MPI
84			if (.not. isMPISimSet()) then
85	gezelter	317	write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
86			FFerror = -1
87			return
88			endif
89	chuckv	292	#endif
90	gezelter	317
91			!! initialize local variables
92			natoms = getNlocal()
93			call getRcut(rcut,rcut2=rcutsq)
94			call getRlist(rlist,rlistsq)
95
96			!! See if we need to update neighbor lists
97	gezelter	325	call checkNeighborList(natoms, q, rcut, rlist, update_nlist)
98	gezelter	317
99			!--------------WARNING...........................
100			! Zero variables, NOTE:::: Forces are zeroed in C
101			! Zeroing them here could delete previously computed
102			! Forces.
103			!------------------------------------------------
104			call zero_module_variables()
105	chuckv	292
106	gezelter	317	! communicate MPI positions
107	chuckv	292	#ifdef IS_MPI
108	gezelter	309	call gather(q,q_Row,plan_row3d)
109			call gather(q,q_Col,plan_col3d)
110	gezelter	317
111	gezelter	309	call gather(u_l,u_l_Row,plan_row3d)
112			call gather(u_l,u_l_Col,plan_col3d)
113	gezelter	317
114	gezelter	309	call gather(A,A_Row,plan_row_rotation)
115			call gather(A,A_Col,plan_col_rotation)
116	chuckv	292	#endif
117
118
119	gezelter	297	#ifdef IS_MPI
120	chuckv	292
121	gezelter	297	if (update_nlist) then
122
123			! save current configuration, contruct neighbor list,
124			! and calculate forces
125			call save_neighborList(q)
126
127			neighborListSize = getNeighborListSize()
128			nlist = 0
129
130			nrow = getNrow(plan_row)
131			ncol = getNcol(plan_col)
132			nlocal = getNlocal()
133
134			do i = 1, nrow
135			point(i) = nlist + 1
136
137			inner: do j = 1, ncol
138
139	gezelter	317	if (check_exclude(i,j)) cycle inner:
140
141			call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
142
143	gezelter	297	if (rijsq < rlistsq) then
144
145			nlist = nlist + 1
146
147			if (nlist > neighborListSize) then
148	gezelter	325	call expandNeighborList(nlocal, listerror)
149	gezelter	297	if (listerror /= 0) then
150			FFerror = -1
151			write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
152			return
153			end if
154			endif
155
156			list(nlist) = j
157	gezelter	317
158	gezelter	297	if (rijsq < rcutsq) then
159	gezelter	317	call do_pair(i, j, rijsq, d)
160	gezelter	297	endif
161			endif
162			enddo inner
163			enddo
164	chuckv	292
165	gezelter	297	point(nrow + 1) = nlist + 1
166
167			else !! (update)
168	chuckv	292
169	gezelter	297	! use the list to find the neighbors
170			do i = 1, nrow
171			JBEG = POINT(i)
172			JEND = POINT(i+1) - 1
173			! check thiat molecule i has neighbors
174			if (jbeg .le. jend) then
175	gezelter	317
176	gezelter	297	do jnab = jbeg, jend
177			j = list(jnab)
178	gezelter	317
179			call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
180			call do_pair(i, j, rijsq, d)
181
182	gezelter	297	enddo
183			endif
184			enddo
185			endif
186
187			#else
188
189			if (update_nlist) then
190	chuckv	292
191	gezelter	297	! save current configuration, contruct neighbor list,
192			! and calculate forces
193			call save_neighborList(q)
194
195			neighborListSize = getNeighborListSize()
196			nlist = 0
197	gezelter	317
198	gezelter	297	do i = 1, natoms-1
199			point(i) = nlist + 1
200
201			inner: do j = i+1, natoms
202	gezelter	317
203			if (check_exclude(i,j)) cycle inner:
204
205			call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
206	chuckv	295
207	gezelter	297	if (rijsq < rlistsq) then
208	gezelter	317
209	gezelter	297	nlist = nlist + 1
210
211			if (nlist > neighborListSize) then
212	gezelter	325	call expandList(natoms, listerror)
213	gezelter	297	if (listerror /= 0) then
214			FFerror = -1
215			write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
216			return
217			end if
218			endif
219
220			list(nlist) = j
221	gezelter	317
222	gezelter	297	if (rijsq < rcutsq) then
223	gezelter	317	call do_pair(i, j, rijsq, d)
224	gezelter	297	endif
225			endif
226	chuckv	292	enddo inner
227	gezelter	297	enddo
228
229			point(natoms) = nlist + 1
230
231			else !! (update)
232
233			! use the list to find the neighbors
234			do i = 1, nrow
235			JBEG = POINT(i)
236			JEND = POINT(i+1) - 1
237			! check thiat molecule i has neighbors
238			if (jbeg .le. jend) then
239
240			do jnab = jbeg, jend
241			j = list(jnab)
242	gezelter	317
243			call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
244			call do_pair(i, j, rijsq, d)
245
246	gezelter	297	enddo
247			endif
248			enddo
249			endif
250	gezelter	317
251	gezelter	297	#endif
252	gezelter	317
253
254	chuckv	292	#ifdef IS_MPI
255	gezelter	297	!! distribute all reaction field stuff (or anything for post-pair):
256			call scatter(rflRow,rflTemp1,plan_row3d)
257			call scatter(rflCol,rflTemp2,plan_col3d)
258			do i = 1,nlocal
259			rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
260			end do
261	chuckv	292	#endif
262	gezelter	317
263	gezelter	297	! This is the post-pair loop:
264			#ifdef IS_MPI
265	gezelter	317
266	gezelter	297	if (system_has_postpair_atoms) then
267			do i = 1, nlocal
268			Atype_i => identPtrListRow(i)%this
269			call do_postpair(i, Atype_i)
270			enddo
271			endif
272	gezelter	317
273	chuckv	292	#else
274	gezelter	317
275	gezelter	297	if (system_has_postpair_atoms) then
276			do i = 1, natoms
277			Atype_i => identPtr(i)%this
278			call do_postpair(i, Atype_i)
279			enddo
280			endif
281	gezelter	317
282	chuckv	292	#endif
283	gezelter	317
284	chuckv	292
285			#ifdef IS_MPI
286	chuckv	295	!!distribute forces
287
288	gezelter	309	call scatter(f_Row,f,plan_row3d)
289			call scatter(f_Col,f_temp,plan_col3d)
290	chuckv	295	do i = 1,nlocal
291	gezelter	309	f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
292	chuckv	295	end do
293
294	gezelter	309	if (doTorque()) then
295			call scatter(t_Row,t,plan_row3d)
296			call scatter(t_Col,t_temp,plan_col3d)
297	chuckv	295
298			do i = 1,nlocal
299	gezelter	309	t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
300	chuckv	295	end do
301			endif
302	gezelter	309
303	chuckv	295	if (do_pot) then
304			! scatter/gather pot_row into the members of my column
305	gezelter	309	call scatter(pot_Row, pot_Temp, plan_row)
306	chuckv	295
307			! scatter/gather pot_local into all other procs
308			! add resultant to get total pot
309			do i = 1, nlocal
310	gezelter	309	pot_local = pot_local + pot_Temp(i)
311	chuckv	295	enddo
312
313	gezelter	309	pot_Temp = 0.0_DP
314
315			call scatter(pot_Col, pot_Temp, plan_col)
316	chuckv	295	do i = 1, nlocal
317	gezelter	309	pot_local = pot_local + pot_Temp(i)
318	chuckv	295	enddo
319
320	gezelter	309	pot = pot_local
321	chuckv	295	endif
322
323	gezelter	309	if (doStress()) then
324			mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
325			mpi_comm_world,mpi_err)
326			mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
327			mpi_comm_world,mpi_err)
328			endif
329	chuckv	295
330	gezelter	309	#endif
331	chuckv	295
332	gezelter	309	if (doStress()) then
333			tau = tau_Temp
334			virial = virial_Temp
335	chuckv	295	endif
336
337			end subroutine do_force_loop
338
339
340			!! Calculate any pre-force loop components and update nlist if necessary.
341			subroutine do_preForce(updateNlist)
342			logical, intent(inout) :: updateNlist
343
344
345
346			end subroutine do_preForce
347
348			!! Calculate any post force loop components, i.e. reaction field, etc.
349			subroutine do_postForce()
350
351
352
353			end subroutine do_postForce
354
355	gezelter	317	subroutine do_pair(i, j, rijsq, d)
356	chuckv	295
357	gezelter	317	integer, intent(in) :: i, j
358			real ( kind = dp ), intent(in) :: rijsq
359			real ( kind = dp ) :: r
360			real ( kind = dp ), intent(inout) :: d(3)
361	chuckv	295
362	gezelter	317	r = sqrt(rijsq)
363	chuckv	295
364	gezelter	317	logical :: is_LJ_i, is_LJ_j
365			logical :: is_DP_i, is_DP_j
366			logical :: is_Sticky_i, is_Sticky_j
367			integer :: me_i, me_j
368	chuckv	295
369	gezelter	302	#ifdef IS_MPI
370
371	gezelter	317	me_i = atid_row(i)
372			me_j = atid_col(j)
373	chuckv	295
374	gezelter	317	#else
375	chuckv	295
376	gezelter	317	me_i = atid(i)
377			me_j = atid(j)
378	gezelter	302
379	gezelter	317	#endif
380	gezelter	302
381	gezelter	317	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
382			call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
383	gezelter	302
384	gezelter	325	if ( is_LJ_i .and. is_LJ_j ) &
385			call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
386	gezelter	302
387	gezelter	317	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
388			call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
389	gezelter	302
390	gezelter	317	if ( is_DP_i .and. is_DP_j ) then
391	gezelter	302
392	gezelter	325	call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
393	chuckv	295
394	gezelter	297	if (do_reaction_field) then
395	gezelter	323	call accumulate_rf(i, j, r)
396	gezelter	297	endif
397	chuckv	295
398	gezelter	297	endif
399
400	gezelter	317	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
401			call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
402
403			if ( is_Sticky_i .and. is_Sticky_j ) then
404	gezelter	325	call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, do_pot, do_stress)
405	gezelter	297	endif
406
407	gezelter	309
408	chuckv	295	end subroutine do_pair
409	chuckv	292
410
411	gezelter	317	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
412
413	chuckv	295	real (kind = dp), dimension(3) :: q_i
414			real (kind = dp), dimension(3) :: q_j
415			real ( kind = dp ), intent(out) :: r_sq
416			real( kind = dp ) :: d(3)
417
418			d(1:3) = q_i(1:3) - q_j(1:3)
419	gezelter	317
420			! Wrap back into periodic box if necessary
421			if ( isPBC() ) then
422	chuckv	295	d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
423			int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
424	gezelter	317	endif
425	chuckv	292
426	chuckv	295	r_sq = dot_product(d,d)
427	gezelter	317
428	chuckv	295	end subroutine get_interatomic_vector
429	gezelter	317
430	gezelter	325	subroutine zero_work_arrays()
431
432			#ifdef IS_MPI
433	chuckv	292
434	gezelter	325	q_Row = 0.0_dp
435			q_Col = 0.0_dp
436	chuckv	295
437	gezelter	325	u_l_Row = 0.0_dp
438			u_l_Col = 0.0_dp
439	chuckv	295
440	gezelter	325	A_Row = 0.0_dp
441			A_Col = 0.0_dp
442	chuckv	295
443	gezelter	325	f_Row = 0.0_dp
444			f_Col = 0.0_dp
445			f_Temp = 0.0_dp
446
447			t_Row = 0.0_dp
448			t_Col = 0.0_dp
449			t_Temp = 0.0_dp
450	chuckv	295
451	gezelter	325	pot_Row = 0.0_dp
452			pot_Col = 0.0_dp
453			pot_Temp = 0.0_dp
454	chuckv	292
455	chuckv	295	#endif
456	chuckv	292
457	gezelter	325	tau_Temp = 0.0_dp
458			virial_Temp = 0.0_dp
459
460			end subroutine zero_work_arrays
461
462	chuckv	292
463	chuckv	295	!! Function to properly build neighbor lists in MPI using newtons 3rd law.
464			!! We don't want 2 processors doing the same i j pair twice.
465			!! Also checks to see if i and j are the same particle.
466	chuckv	306	function checkExcludes(atom1,atom2) result(do_cycle)
467	chuckv	295	!--------------- Arguments--------------------------
468			! Index i
469	chuckv	306	integer,intent(in) :: atom1
470	chuckv	295	! Index j
471	chuckv	306	integer,intent(in), optional :: atom2
472	chuckv	295	! Result do_cycle
473			logical :: do_cycle
474			!--------------- Local variables--------------------
475			integer :: tag_i
476			integer :: tag_j
477	chuckv	306	integer :: i
478			!--------------- END DECLARATIONS------------------
479			do_cycle = .false.
480
481			#ifdef IS_MPI
482			tag_i = tagRow(atom1)
483			#else
484			tag_i = tag(atom1)
485			#endif
486
487			!! Check global excludes first
488			if (.not. present(atom2)) then
489			do i = 1,nGlobalExcludes
490			if (excludeGlobal(i) == tag_i) then
491			do_cycle = .true.
492			return
493			end if
494			end do
495			return !! return after checking globals
496			end if
497
498			!! we return if j not present here.
499	chuckv	295	tag_j = tagColumn(j)
500	chuckv	292
501	chuckv	306
502	chuckv	292
503	chuckv	295	if (tag_i == tag_j) then
504			do_cycle = .true.
505			return
506			end if
507
508			if (tag_i < tag_j) then
509			if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
510			return
511			else
512			if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
513			endif
514
515	chuckv	306
516
517			do i = 1, nLocalExcludes
518			if (tag_i = excludes(1,i) .and. excludes(2,i) < 0) then
519			do_cycle = .true.
520			return
521			end if
522			end do
523
524
525			end function checkExcludes
526
527
528	chuckv	292	end module do_Forces