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.16 2003-03-12 23:15:46 gezelter Exp $, $Date: 2003-03-12 23:15:46 $, $Name: not supported by cvs2svn $, $Revision: 1.16 $


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

#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE

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


  public :: init_FF
  public :: do_force_loop

contains

  subroutine init_FF(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 :: nrow
    integer :: ncol
#endif
    integer :: nlocal
    integer :: natoms    
    logical :: update_nlist   
    integer :: i, j, jbeg, jend, jnab
    integer :: nlist
    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
    real(kind=dp),dimension(3) :: d
    real(kind=dp) :: rfpot, mu_i, virial
    integer :: me_i
    logical :: is_dp_i
    integer :: neighborListSize
    integer :: listerror, error

    !! 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 (checkExcludes(i,j)) cycle inner:
             
             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
             
             if (rijsq <  rlistsq) then             
                
                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandNeighborList(nlocal, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                endif
                
                list(nlist) = j
                                
                if (rijsq <  rcutsq) then
                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
                endif
             endif
          enddo inner
       enddo

       point(nrow + 1) = nlist + 1
       
    else  !! (of update_check)

       ! use the list to find the neighbors
       do i = 1, nrow
          JBEG = POINT(i)
          JEND = POINT(i+1) - 1
          ! check thiat molecule i has neighbors
          if (jbeg .le. jend) then
             
             do jnab = jbeg, jend
                j = list(jnab)

                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
                call do_pair(i, j, rijsq, d, do_pot, do_stress)

             enddo
          endif
       enddo
    endif
    
#else
    
    if (update_nlist) then
       
       ! save current configuration, contruct neighbor list, 
       ! and calculate forces
       call save_neighborList(q)
       
       neighborListSize = getNeighborListSize()
       nlist = 0
       
       do i = 1, natoms-1
          point(i) = nlist + 1
          
          inner: do j = i+1, natoms
             
             if (checkExcludes(i,j)) cycle inner:
             
             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
           
             if (rijsq <  rlistsq) then
                
                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandList(natoms, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                endif
                
                list(nlist) = j
                
                if (rijsq <  rcutsq) then
                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
                endif
             endif
          enddo inner
       enddo
       
       point(natoms) = nlist + 1
       
    else !! (update)
       
       ! use the list to find the neighbors
       do i = 1, natoms-1
          JBEG = POINT(i)
          JEND = POINT(i+1) - 1
          ! check thiat molecule i has neighbors
          if (jbeg .le. jend) then
             
             do jnab = jbeg, jend
                j = list(jnab)

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

             enddo
          endif
       enddo
    endif
    
#endif 
    
    ! phew, done with main loop.
    
#ifdef IS_MPI
    !!distribute forces
    
    call scatter(f_Row,f,plan_row3d)
    call scatter(f_Col,f_temp,plan_col3d)
    do i = 1,nlocal
       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
    end do
    
    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
       call scatter(t_Row,t,plan_row3d)
       call scatter(t_Col,t_temp,plan_col3d)
       
       do i = 1,nlocal
          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
       end do
    endif
    
    if (do_pot) then
       ! scatter/gather pot_row into the members of my column
       call scatter(pot_Row, pot_Temp, plan_row)
       
       ! scatter/gather pot_local into all other procs
       ! add resultant to get total pot
       do i = 1, nlocal
          pot_local = pot_local + pot_Temp(i)
       enddo

       pot_Temp = 0.0_DP

       call scatter(pot_Col, pot_Temp, plan_col)
       do i = 1, nlocal
          pot_local = pot_local + pot_Temp(i)
       enddo
       
    endif    
#endif

    if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
       
       if (FF_uses_RF .and. SimUsesRF()) then
          
#ifdef IS_MPI
          call scatter(rf_Row,rf,plan_row3d)
          call scatter(rf_Col,rf_Temp,plan_col3d)
          do i = 1,nlocal
             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
          end do
#endif
          
          do i = 1, getNlocal()

             rfpot = 0.0_DP
#ifdef IS_MPI
             me_i = atid_row(i)
#else
             me_i = atid(i)
#endif
             call getElementProperty(atypes, me_i, "is_DP", is_DP_i)       
             if ( is_DP_i ) then
                call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
                !! The reaction field needs to include a self contribution 
                !! to the field:
                call accumulate_self_rf(i, mu_i, u_l)             
                !! Get the reaction field contribution to the 
                !! potential and torques:
                call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#ifdef IS_MPI
                pot_local = pot_local + rfpot
#else
                pot = pot + rfpot
#endif
             endif             
          enddo
       endif
    endif


#ifdef IS_MPI

    if (do_pot) then
       pot = pot_local
       !! we assume the c code will do the allreduce to get the total potential
       !! we could do it right here if we needed to...
    endif

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

#else

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

#endif
    
  end subroutine do_force_loop


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

    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3,getNlocal()) :: u_l
    real (kind=dp), dimension(9,getNlocal()) :: A
    real (kind=dp), dimension(3,getNlocal()) :: f
    real (kind=dp), dimension(3,getNlocal()) :: t

    logical, intent(inout) :: do_pot, do_stress
    integer, intent(in) :: i, j
    real ( kind = dp ), intent(in)    :: rijsq
    real ( kind = dp )                :: r
    real ( kind = dp ), intent(inout) :: d(3)
    logical :: is_LJ_i, is_LJ_j
    logical :: is_DP_i, is_DP_j
    logical :: is_Sticky_i, is_Sticky_j
    integer :: me_i, me_j

    r = sqrt(rijsq)
    
#ifdef IS_MPI

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

#else

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

#endif


    if (FF_uses_LJ .and. SimUsesLJ()) then
       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
       
       if ( is_LJ_i .and. is_LJ_j ) &
            call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
    endif
       

    if (FF_uses_dipoles .and. SimUsesDipoles()) then
       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
       
       if ( is_DP_i .and. is_DP_j ) then
          
          call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
          
          if (FF_uses_RF .and. SimUsesRF()) then
             
             call accumulate_rf(i, j, r, u_l)
             call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
             
          endif
          
       endif
    endif

    if (FF_uses_Sticky .and. SimUsesSticky()) then

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


  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
    
    real (kind = dp), dimension(3) :: q_i
    real (kind = dp), dimension(3) :: q_j
    real ( kind = dp ), intent(out) :: r_sq
    real( kind = dp ) :: d(3)

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

  subroutine check_initialization(error)
    integer, intent(out) :: error

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

    return
  end subroutine check_initialization

  
  subroutine zero_work_arrays()
    
#ifdef IS_MPI

    q_Row = 0.0_dp
    q_Col = 0.0_dp   
    
    u_l_Row = 0.0_dp
    u_l_Col = 0.0_dp
    
    A_Row = 0.0_dp
    A_Col = 0.0_dp
    
    f_Row = 0.0_dp
    f_Col = 0.0_dp
    f_Temp = 0.0_dp
      
    t_Row = 0.0_dp
    t_Col = 0.0_dp
    t_Temp = 0.0_dp
 
    pot_Row = 0.0_dp
    pot_Col = 0.0_dp
    pot_Temp = 0.0_dp

#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, j
    !--------------- 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, nExcludes_global
          if (excludeGlobal(i) == tag_i) then
             do_cycle = .true.
             return
          end if
       end do
       return !! return after checking globals
    end if

    !! we return if atom2 not present here.
    tag_j = tagColumn(atom2)
    
    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, nExcludes_local
       if ((tag_i == excludesLocal(1,i)) .and. (excludesLocal(2,i) < 0)) then
          do_cycle = .true.
          return
       end if
    end do
    
    
  end function checkExcludes

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