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.17 2003-03-13 00:33:18 chuckv Exp $, $Date: 2003-03-13 00:33:18 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $


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) :: thisStat   
    integer :: my_status 
    character(len = 100) :: mix_Policy

    ! be a smarter subroutine.
    mix_Policy = "FIXME"
    thisStat = 0
    call init_lj_FF(mix_Policy,my_status)
    if (my_status /= 0) then
       thisStat = -1
       return
    end if
    
    call check_sticky_FF(my_status)
    if (my_status /= 0) then
       thisStat = -1
       return
    end if
    
    do_forces_initialized = .true.    
    
  end subroutine init_FF


  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
  !------------------------------------------------------------->
  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
       error)
    !! Position array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,getNlocal()) :: q
    !! Rotation Matrix for each long range particle in simulation.
    real( kind = dp), dimension(9,getNlocal()) :: A    
    !! Unit vectors for dipoles (lab frame)
    real( kind = dp ), dimension(3,getNlocal()) :: u_l
    !! Force array provided by C, dimensioned by getNlocal
    real ( kind = dp ), dimension(3,getNlocal()) :: f
    !! Torsion array provided by C, dimensioned by getNlocal
    real( kind = dp ), dimension(3,getNlocal()) :: t    
    !! Stress Tensor
    real( kind = dp), dimension(9) :: tau   
    real ( kind = dp ) :: pot
    logical ( kind = 2) :: do_pot_c, do_stress_c
    logical :: do_pot
    logical :: do_stress
#ifdef IS_MPI
    real( kind = DP ) :: pot_local
    integer :: nrow
    integer :: ncol
#endif
    integer :: nlocal
    integer :: natoms    
    logical :: update_nlist   
    integer :: i, j, jbeg, jend, jnab
    integer :: nlist
    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
    real(kind=dp),dimension(3) :: d
    real(kind=dp) :: rfpot, mu_i, virial
    integer :: me_i
    logical :: is_dp_i
    integer :: neighborListSize
    integer :: listerror, error
    integer :: localError

    !! initialize local variables  

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

    call getRcut(rcut,rc2=rcutsq)
    call getRlist(rlist,rlistsq)
    
    call check_initialization(localError)
    if ( localError .ne. 0 ) then
       error = -1
       return
    end if
    call zero_work_arrays()

    do_pot = do_pot_c
    do_stress = do_stress_c

    ! Gather all information needed by all force loops:
    
#ifdef IS_MPI    

    call gather(q,q_Row,plan_row3d)
    call gather(q,q_Col,plan_col3d)
        
    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
       call gather(u_l,u_l_Row,plan_row3d)
       call gather(u_l,u_l_Col,plan_col3d)
       
       call gather(A,A_Row,plan_row_rotation)
       call gather(A,A_Col,plan_col_rotation)
    endif
    
#endif
    
    if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
       !! See if we need to update neighbor lists
       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)   
       !! if_mpi_gather_stuff_for_prepair
       !! do_prepair_loop_if_needed
       !! if_mpi_scatter_stuff_from_prepair
       !! if_mpi_gather_stuff_from_prepair_to_main_loop
    else
       !! See if we need to update neighbor lists
       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)   
    endif
    
#ifdef IS_MPI
    
    if (update_nlist) then
       
       !! save current configuration, construct neighbor list, 
       !! and calculate forces
       call save_neighborList(q)
       
       neighborListSize = getNeighborListSize()
       nlist = 0       
       
       do i = 1, nrow
          point(i) = nlist + 1
          
          inner: do j = 1, ncol
             
             if (checkExcludes(i,j)) cycle inner
             
             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
             
             if (rijsq <  rlistsq) then             
                
                nlist = nlist + 1
                
                if (nlist > neighborListSize) then
                   call expandNeighborList(nlocal, listerror)
                   if (listerror /= 0) then
                      error = -1
                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
                      return
                   end if
                endif
                
                list(nlist) = j
                                
                if (rijsq <  rcutsq) then
                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
                endif
             endif
          enddo inner
       enddo

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

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

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

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

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

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

       pot_Temp = 0.0_DP

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

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

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


#ifdef IS_MPI

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

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

#else

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

#endif
    
  end subroutine do_force_loop


!! 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(inout)    :: rijsq
    real ( kind = dp )                :: r
    real ( kind = dp ), intent(inout) :: d(3)
    logical :: is_LJ_i, is_LJ_j
    logical :: is_DP_i, is_DP_j
    logical :: is_Sticky_i, is_Sticky_j
    integer :: me_i, me_j

    r = sqrt(rijsq)
    
#ifdef IS_MPI

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

#else

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

#endif


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

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

    if (FF_uses_Sticky .and. SimUsesSticky()) then

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


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

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

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

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

    return
  end subroutine check_initialization

  
  subroutine zero_work_arrays()
    
#ifdef IS_MPI

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

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