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.12 2003-03-12 15:39:01 gezelter Exp $, $Date: 2003-03-12 15:39:01 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $


module do_Forces
  use simulation
  use definitions
  use forceGlobals
  use atype_typedefs
  use neighborLists

  
  use lj
  use sticky_FF
  use dipole_dipole
  use gb_FF

#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

public :: do_force_loop

contains

  !! Does force loop over i,j pairs. Calls do_pair to calculates 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

#ifdef IS_MPI
    real( kind = DP ) :: pot_local
#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

    ! 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)
    
    !! 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,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 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(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(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)

    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)

       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)
    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_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


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