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.15 2003-03-12 22:27:59 gezelter Exp $, $Date: 2003-03-12 22:27:59 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $


module do_Forces
  use simulation
  use definitions
  use atype_module
  use neighborLists  
  use lj_FF
  use sticky_FF
  use dipole_dipole
  use gb_FF

#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_forces

contains

  subroutine init_FF(thisStat)
   
    integer, intent(out) :: my_status 
    integer :: thisStat = 0

    ! be a smarter subroutine.

    
    call init_lj_FF(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 :: nlocal
    integer :: nrow
    integer :: ncol
#endif
    integer :: natoms    
    logical :: update_nlist   
    integer :: i, j, jbeg, jend, jnab
    integer :: nlist
    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
    integer :: neighborListSize
    integer :: listerror

    !! 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,rcut2=rcutsq)
    call getRlist(rlist,rlistsq)
    
    call check_initialization()
    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 (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
                      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 (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
                      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, 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, 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
    
    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()

#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(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
       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

#else

    if (do_stress) then
       tau = tau_Temp
       virial = virial_Temp
    endif

#endif
    
  end subroutine do_force_loop


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


  end subroutine do_preForce

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


  end subroutine do_postForce

  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress)

    logical, intent(inout) :: do_pot, do_stress
    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


    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_DP .and. SimUsesDP()) 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 ( 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 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

#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

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