ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90 (file contents):
Revision 292 by chuckv, Thu Mar 6 14:52:44 2003 UTC vs.
Revision 334 by gezelter, Thu Mar 13 17:45:54 2003 UTC

# Line 1 | Line 1
1 + !! do_Forces.F90
2 + !! module do_Forces
3   !! Calculates Long Range forces.
4 +
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.1 2003-03-06 14:52:44 chuckv Exp $, $Date: 2003-03-06 14:52:44 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
7 > !! @version $Id: do_Forces.F90,v 1.19 2003-03-13 17:45:54 gezelter Exp $, $Date: 2003-03-13 17:45:54 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
8  
9  
10  
11   module do_Forces
12    use simulation
13    use definitions
14 <  use generic_atypes
15 <  use neighborLists
16 <  
17 <  use lj_FF
18 <  use sticky_FF
19 <  use dp_FF
17 <  use gb_FF
14 >  use atype_module
15 >  use neighborLists  
16 >  use lj
17 >  use sticky_pair
18 >  use dipole_dipole
19 >  use reaction_field
20  
21   #ifdef IS_MPI
22    use mpiSimulation
# Line 22 | Line 24 | module do_Forces
24    implicit none
25    PRIVATE
26  
27 < !! Number of lj_atypes in lj_atype_list
28 <  integer, save :: n_atypes = 0
27 >  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  
28 – !! Global list of lj atypes in simulation
29 –  type (atype), pointer :: ListHead => null()
30 –  type (atype), pointer :: ListTail => null()
31 –
32 –
33 – !! LJ mixing array  
34 – !  type (lj_atype), dimension(:,:), pointer :: ljMixed => null()
35 –
36 –
37 –  logical, save :: firstTime = .True.
38 –
39 – !! Atype identity pointer lists
40 – #ifdef IS_MPI
41 – !! Row lj_atype pointer list
42 –  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
43 – !! Column lj_atype pointer list
44 –  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
45 – #else
46 –  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
47 – #endif
48 –
49 –
50 – !! Logical has lj force field module been initialized?
51 –  logical, save :: isFFinit = .false.
52 –
53 –
54 – !! Public methods and data
55 –  public :: new_atype
56 –  public :: do_lj_ff
57 –  public :: getLjPot
35    public :: init_FF
36 +  public :: do_force_loop
37  
60 –  
61 –
62 –
38   contains
39  
40 < !! Adds a new lj_atype to the list.
41 <  subroutine new_atype(ident,mass,epsilon,sigma, &
42 <       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
43 <    real( kind = dp ), intent(in) :: mass
44 <    real( kind = dp ), intent(in) :: epsilon
45 <    real( kind = dp ), intent(in) :: sigma
46 <    real( kind = dp ), intent(in) :: w0
47 <    real( kind = dp ), intent(in) :: v0
48 <    real( kind = dp ), intent(in) :: dipoleMoment
40 >  subroutine init_FF(LJ_mix_policy, use_RF_c, thisStat)
41 >    logical(kind = 2), intent(in) :: use_RF_c
42 >    logical :: use_RF_f
43 >    integer, intent(out) :: thisStat  
44 >    integer :: my_status, nMatches
45 >    character(len = 100) :: LJ_mix_Policy
46 >    integer, pointer :: MatchList(:)
47 >    
48 >    !! Fortran's version of a cast:
49 >    use_RF_f = use_RF_c
50  
51 <    integer, intent(in) :: ident
52 <    integer, intent(out) :: status
53 <    integer, intent(in) :: is_Sticky
54 <    integer, intent(in) :: is_DP
55 <    integer, intent(in) :: is_GB
56 <    integer, intent(in) :: is_LJ
51 >    !! assume things are copacetic, unless they aren't
52 >    thisStat = 0
53 >    
54 >    !! init_FF is called *after* all of the atom types have been
55 >    !! defined in atype_module using the new_atype subroutine.
56 >    !!
57 >    !! this will scan through the known atypes and figure out what
58 >    !! interactions are used by the force field.    
59  
60 +    FF_uses_LJ = .false.
61 +    FF_uses_sticky = .false.
62 +    FF_uses_dipoles = .false.
63 +    FF_uses_GB = .false.
64 +    FF_uses_EAM = .false.
65 +    
66 +    call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
67 +    deallocate(MatchList)
68 +    if (nMatches .gt. 0) FF_uses_LJ = .true.
69  
70 <    type (atype), pointer :: the_new_atype
71 <    integer :: alloc_error
72 <    integer :: atype_counter = 0
86 <    integer :: alloc_size
87 <    integer :: err_stat
88 <    status = 0
70 >    call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
71 >    deallocate(MatchList)
72 >    if (nMatches .gt. 0) FF_uses_dipoles = .true.
73  
74 +    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
75 +         MatchList)
76 +    deallocate(MatchList)
77 +    if (nMatches .gt. 0) FF_uses_Sticky = .true.
78 +    
79 +    call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
80 +    deallocate(MatchList)
81 +    if (nMatches .gt. 0) FF_uses_GB = .true.
82  
83 +    call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
84 +    deallocate(MatchList)
85 +    if (nMatches .gt. 0) FF_uses_EAM = .true.
86  
87 < ! allocate a new atype    
88 <    allocate(the_new_atype,stat=alloc_error)
89 <    if (alloc_error /= 0 ) then
90 <       status = -1
87 >    !! check to make sure the use_RF setting makes sense
88 >    if (use_RF_f) then
89 >       if (FF_uses_dipoles) then
90 >          FF_uses_RF = .true.
91 >          call initialize_rf()
92 >       else
93 >          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
94 >          thisStat = -1
95 >          return
96 >       endif
97 >    endif
98 >    
99 >    call init_lj_FF(LJ_mix_Policy, my_status)
100 >    if (my_status /= 0) then
101 >       thisStat = -1
102         return
103      end if
98 –
99 – ! assign our new atype information
100 –    the_new_atype%mass        = mass
101 –    the_new_atype%epsilon     = epsilon
102 –    the_new_atype%sigma       = sigma
103 –    the_new_atype%sigma2      = sigma * sigma
104 –    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
105 –         * the_new_atype%sigma2
106 –    the_new_atype%w0       = w0
107 –    the_new_atype%v0       = v0
108 –    the_new_atype%dipoleMoment       = dipoleMoment
109 –
104      
105 < ! assume that this atype will be successfully added
106 <    the_new_atype%atype_ident = ident
107 <    the_new_atype%atype_number = n_lj_atypes + 1
105 >    call check_sticky_FF(my_status)
106 >    if (my_status /= 0) then
107 >       thisStat = -1
108 >       return
109 >    end if
110      
111 <    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
112 <    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
113 <    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
118 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
111 >    do_forces_initialized = .true.    
112 >    
113 >  end subroutine init_FF
114  
120 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
121 –    if (err_stat /= 0 ) then
122 –       status = -1
123 –       return
124 –    endif
115  
126 –    n_atypes = n_atypes + 1
116  
117 <
118 <  end subroutine new_atype
119 <
120 <
121 <  subroutine init_FF(nComponents,ident, status)
122 < !! Number of components in ident array
123 <    integer, intent(inout) :: nComponents
124 < !! Array of identities nComponents long corresponding to
125 < !! ljatype ident.
126 <    integer, dimension(nComponents),intent(inout) :: ident
127 < !!  Result status, success = 0, error = -1
128 <    integer, intent(out) :: Status
129 <
130 <    integer :: alloc_stat
131 <
132 <    integer :: thisStat
133 <    integer :: i
134 <
135 <    integer :: myNode
117 >  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
118 >  !------------------------------------------------------------->
119 >  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
120 >       error)
121 >    !! Position array provided by C, dimensioned by getNlocal
122 >    real ( kind = dp ), dimension(3,getNlocal()) :: q
123 >    !! Rotation Matrix for each long range particle in simulation.
124 >    real( kind = dp), dimension(9,getNlocal()) :: A    
125 >    !! Unit vectors for dipoles (lab frame)
126 >    real( kind = dp ), dimension(3,getNlocal()) :: u_l
127 >    !! Force array provided by C, dimensioned by getNlocal
128 >    real ( kind = dp ), dimension(3,getNlocal()) :: f
129 >    !! Torsion array provided by C, dimensioned by getNlocal
130 >    real( kind = dp ), dimension(3,getNlocal()) :: t    
131 >    !! Stress Tensor
132 >    real( kind = dp), dimension(9) :: tau  
133 >    real ( kind = dp ) :: pot
134 >    logical ( kind = 2) :: do_pot_c, do_stress_c
135 >    logical :: do_pot
136 >    logical :: do_stress
137   #ifdef IS_MPI
138 <    integer, allocatable, dimension(:) :: identRow
149 <    integer, allocatable, dimension(:) :: identCol
138 >    real( kind = DP ) :: pot_local
139      integer :: nrow
140      integer :: ncol
141   #endif
142 <    status = 0
143 <  
142 >    integer :: nlocal
143 >    integer :: natoms    
144 >    logical :: update_nlist  
145 >    integer :: i, j, jbeg, jend, jnab
146 >    integer :: nlist
147 >    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
148 >    real(kind=dp),dimension(3) :: d
149 >    real(kind=dp) :: rfpot, mu_i, virial
150 >    integer :: me_i
151 >    logical :: is_dp_i
152 >    integer :: neighborListSize
153 >    integer :: listerror, error
154 >    integer :: localError
155  
156 <    
156 >    !! initialize local variables  
157  
158 < !! if were're not in MPI, we just update ljatypePtrList
159 < #ifndef IS_MPI
160 <    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
161 <    if ( thisStat /= 0 ) then
162 <       status = -1
163 <       return
164 <    endif
165 <
166 <    
167 < ! if were're in MPI, we also have to worry about row and col lists    
158 > #ifdef IS_MPI
159 >    nlocal = getNlocal()
160 >    nrow   = getNrow(plan_row)
161 >    ncol   = getNcol(plan_col)
162   #else
163 <  
164 < ! We can only set up forces if mpiSimulation has been setup.
165 <    if (.not. isMPISimSet()) then
172 <       write(default_error,*) "MPI is not set"
173 <       status = -1
174 <       return
175 <    endif
176 <    nrow = getNrow(plan_row)
177 <    ncol = getNcol(plan_col)
178 <    mynode = getMyNode()
179 < !! Allocate temperary arrays to hold gather information
180 <    allocate(identRow(nrow),stat=alloc_stat)
181 <    if (alloc_stat /= 0 ) then
182 <       status = -1
183 <       return
184 <    endif
163 >    nlocal = getNlocal()
164 >    natoms = nlocal
165 > #endif
166  
167 <    allocate(identCol(ncol),stat=alloc_stat)
168 <    if (alloc_stat /= 0 ) then
188 <       status = -1
189 <       return
190 <    endif
191 <
192 < !! Gather idents into row and column idents
193 <
194 <    call gather(ident,identRow,plan_row)
195 <    call gather(ident,identCol,plan_col)
167 >    call getRcut(rcut,rc2=rcutsq)
168 >    call getRlist(rlist,rlistsq)
169      
170 <  
171 < !! Create row and col pointer lists
172 <  
200 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
201 <    if (thisStat /= 0 ) then
202 <       status = -1
170 >    call check_initialization(localError)
171 >    if ( localError .ne. 0 ) then
172 >       error = -1
173         return
204 –    endif
205 –  
206 –    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
207 –    if (thisStat /= 0 ) then
208 –       status = -1
209 –       return
210 –    endif
211 –
212 – !! free temporary ident arrays
213 –    if (allocated(identCol)) then
214 –       deallocate(identCol)
174      end if
175 <    if (allocated(identCol)) then
217 <       deallocate(identRow)
218 <    endif
175 >    call zero_work_arrays()
176  
177 < #endif
178 <    
222 <    call initForce_Modules(thisStat)
223 <    if (thisStat /= 0) then
224 <       status = -1
225 <       return
226 <    endif
227 <
228 < !! Create neighbor lists
229 <    call expandList(thisStat)
230 <    if (thisStat /= 0) then
231 <       status = -1
232 <       return
233 <    endif
177 >    do_pot = do_pot_c
178 >    do_stress = do_stress_c
179  
180 <    isFFinit = .true.
236 <
237 <
238 <  end subroutine init_FF
239 <
240 <
241 <
242 <
243 <  subroutine initForce_Modules(thisStat)
244 <    integer, intent(out) :: thisStat
245 <    integer :: my_status
180 >    ! Gather all information needed by all force loops:
181      
182 <    thisStat = 0
248 <    call init_lj_FF(ListHead,my_status)
249 <    if (my_status /= 0) then
250 <       thisStat = -1
251 <       return
252 <    end if
182 > #ifdef IS_MPI    
183  
184 <  end subroutine initForce_Modules
185 <
186 <
187 <
188 <
189 < !! FORCE routine Calculates Lennard Jones forces.
190 < !------------------------------------------------------------->
191 <  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
192 < !! Position array provided by C, dimensioned by getNlocal
193 <    real ( kind = dp ), dimension(3,getNlocal()) :: q
264 <  !! Rotation Matrix for each long range particle in simulation.
265 <    real( kind = dp), dimension(9,getNlocal()) :: A
266 <
267 <  !! Magnitude dipole moment
268 <    real( kind = dp ), dimension(3,getNlocal()) :: mu
269 <  !! Unit vectors for dipoles (lab frame)
270 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
271 < !! Force array provided by C, dimensioned by getNlocal
272 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
273 < !! Torsion array provided by C, dimensioned by getNlocal
274 <    real( kind = dp ), dimension(3,getNlocal()) :: t
275 <
276 < !! Stress Tensor
277 <    real( kind = dp), dimension(9) :: tau
278 <    real( kind = dp), dimension(9) :: tauTemp
279 <    real ( kind = dp ) :: potE
280 <    logical ( kind = 2) :: do_pot
281 <    integer :: FFerror
282 <
184 >    call gather(q,q_Row,plan_row3d)
185 >    call gather(q,q_Col,plan_col3d)
186 >        
187 >    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
188 >       call gather(u_l,u_l_Row,plan_row3d)
189 >       call gather(u_l,u_l_Col,plan_col3d)
190 >      
191 >       call gather(A,A_Row,plan_row_rotation)
192 >       call gather(A,A_Col,plan_col_rotation)
193 >    endif
194      
284 –    type(atype), pointer :: Atype_i
285 –    type(atype), pointer :: Atype_j
286 –
287 –
288 –
289 –
290 –  
291 –
292 – #ifdef IS_MPI
293 –  real( kind = DP ) :: pot_local
294 –
295 – !! Local arrays needed for MPI
296 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
297 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
298 –
299 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
300 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
301 –
302 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
303 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
304 –
305 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: ARow = 0.0_dp
306 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp
307 –
308 –  
309 –
310 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
311 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
312 –  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp = 0.0_dp
313 –
314 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
315 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
316 –  real(kind = dp), dimension(3,getNlocal()) :: tMPITemp = 0.0_dp
317 –
318 –
319 –  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
320 –  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
321 –
322 –  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
323 –
324 – #endif
325 –
326 –
327 –
328 –  real( kind = DP )   :: pe
329 –  logical             :: update_nlist
330 –
331 –
332 –  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
333 –  integer :: nlist
334 –  integer :: j_start
335 –  integer :: tag_i,tag_j
336 –  real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
337 –  real( kind = dp ) :: fx,fy,fz
338 –  real( kind = DP ) ::  drdx, drdy, drdz
339 –  real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
340 –  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
341 –
342 –  real( kind = DP ) :: dielectric = 0.0_dp
343 –
344 – ! a rig that need to be fixed.
345 – #ifdef IS_MPI
346 –  logical :: newtons_thrd = .true.
347 –  real( kind = dp ) :: pe_local
348 –  integer :: nlocal
195   #endif
350 –  integer :: nrow
351 –  integer :: ncol
352 –  integer :: natoms
353 –  integer :: neighborListSize
354 –  integer :: listerror
355 – !! should we calculate the stress tensor
356 –  logical  :: do_stress = .false.
357 –
358 –
359 –  FFerror = 0
360 –
361 – ! Make sure we are properly initialized.
362 –  if (.not. isFFInit) then
363 –     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
364 –     FFerror = -1
365 –     return
366 –  endif
367 – #ifdef IS_MPI
368 –    if (.not. isMPISimSet()) then
369 –     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
370 –     FFerror = -1
371 –     return
372 –  endif
373 – #endif
374 –
375 – !! initialize local variables  
376 –  natoms = getNlocal()
377 –  call getRcut(rcut,rcut2=rcutsq)
378 –  call getRlist(rlist,rlistsq)
379 – !! Find ensemble
380 –  if (isEnsemble("NPT")) do_stress = .true.
381 –
382 – #ifndef IS_MPI
383 –  nrow = natoms - 1
384 –  ncol = natoms
385 – #else
386 –  nrow = getNrow(plan_row)
387 –  ncol = getNcol(plan_col)
388 –  nlocal = natoms
389 –  j_start = 1
390 – #endif
391 –
392 –  
393 – !! See if we need to update neighbor lists
394 –  call check(q,update_nlist)
395 – !  if (firstTime) then
396 – !     update_nlist = .true.
397 – !     firstTime = .false.
398 – !  endif
399 –
400 – !--------------WARNING...........................
401 – ! Zero variables, NOTE:::: Forces are zeroed in C
402 – ! Zeroing them here could delete previously computed
403 – ! Forces.
404 – !------------------------------------------------
405 – #ifndef IS_MPI
406 – !  nloops = nloops + 1
407 –  pe = 0.0E0_DP
408 –
409 – #else
410 –    fRow = 0.0E0_DP
411 –    fCol = 0.0E0_DP
412 –
413 –    pe_local = 0.0E0_DP
414 –
415 –    eRow = 0.0E0_DP
416 –    eCol = 0.0E0_DP
417 –    eTemp = 0.0E0_DP
418 – #endif
419 –
420 – ! communicate MPI positions
421 – #ifdef IS_MPI    
422 –    call gather(q,qRow,plan_row3d)
423 –    call gather(q,qCol,plan_col3d)
424 –
425 –    call gather(mu,muRow,plan_row3d)
426 –    call gather(mu,muCol,plan_col3d)
427 –
428 –    call gather(u_l,u_lRow,plan_row3d)
429 –    call gather(u_l,u_lCol,plan_col3d)
430 –
431 –    call gather(A,ARow,plan_row_rotation)
432 –    call gather(A,ACol,plan_col_rotation)
433 –
434 – #endif
435 –
436 –
437 –  if (update_nlist) then
438 –
439 –     ! save current configuration, contruct neighbor list,
440 –     ! and calculate forces
441 –     call save_neighborList(q)
442 –    
443 –     neighborListSize = getNeighborListSize()
444 –     nlist = 0
445 –    
196      
197 <
198 <     do i = 1, nrow
199 <        point(i) = nlist + 1
197 >    if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
198 >       !! See if we need to update neighbor lists
199 >       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
200 >       !! if_mpi_gather_stuff_for_prepair
201 >       !! do_prepair_loop_if_needed
202 >       !! if_mpi_scatter_stuff_from_prepair
203 >       !! if_mpi_gather_stuff_from_prepair_to_main_loop
204 >    else
205 >       !! See if we need to update neighbor lists
206 >       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
207 >    endif
208 >    
209   #ifdef IS_MPI
451 –        ljAtype_i => identPtrListRow(i)%this
452 –        tag_i = tagRow(i)
453 –        rxi = qRow(1,i)
454 –        ryi = qRow(2,i)
455 –        rzi = qRow(3,i)
456 – #else
457 –        ljAtype_i   => identPtrList(i)%this
458 –        j_start = i + 1
459 –        rxi = q(1,i)
460 –        ryi = q(2,i)
461 –        rzi = q(3,i)
462 – #endif
463 –
464 –        inner: do j = j_start, ncol
465 – #ifdef IS_MPI
466 – ! Assign identity pointers and tags
467 –           ljAtype_j => identPtrListColumn(j)%this
468 –           tag_j = tagColumn(j)
469 –           if (newtons_thrd) then
470 –              if (tag_i <= tag_j) then
471 –                 if (mod(tag_i + tag_j,2) == 0) cycle inner
472 –              else                
473 –                 if (mod(tag_i + tag_j,2) == 1) cycle inner
474 –              endif
475 –           endif
476 –
477 –           rxij = wrap(rxi - qCol(1,j), 1)
478 –           ryij = wrap(ryi - qCol(2,j), 2)
479 –           rzij = wrap(rzi - qCol(3,j), 3)
480 – #else          
481 –           ljAtype_j   => identPtrList(j)%this
482 –           rxij = wrap(rxi - q(1,j), 1)
483 –           ryij = wrap(ryi - q(2,j), 2)
484 –           rzij = wrap(rzi - q(3,j), 3)
485 –      
486 – #endif          
487 –           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
488 –
489 – #ifdef IS_MPI
490 –             if (rijsq <=  rlistsq .AND. &
491 –                  tag_j /= tag_i) then
492 – #else
493 –          
494 –             if (rijsq <  rlistsq) then
495 – #endif
496 –            
497 –              nlist = nlist + 1
498 –              if (nlist > neighborListSize) then
499 –                 call expandList(listerror)
500 –                 if (listerror /= 0) then
501 –                    FFerror = -1
502 –                    write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
503 –                    return
504 –                 end if
505 –              endif
506 –              list(nlist) = j
507 –
210      
211 <              if (rijsq <  rcutsq) then
510 <                
511 <                 r = dsqrt(rijsq)
211 >    if (update_nlist) then
212        
213 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
213 >       !! save current configuration, construct neighbor list,
214 >       !! and calculate forces
215 >       call saveNeighborList(q)
216        
217 < #ifdef IS_MPI
218 <                eRow(i) = eRow(i) + pot*0.5
219 <                eCol(i) = eCol(i) + pot*0.5
220 < #else
221 <                    pe = pe + pot
222 < #endif                
223 <            
224 <                drdx = -rxij / r
225 <                drdy = -ryij / r
226 <                drdz = -rzij / r
217 >       neighborListSize = getNeighborListSize()
218 >       nlist = 0      
219 >      
220 >       do i = 1, nrow
221 >          point(i) = nlist + 1
222 >          
223 >          inner: do j = 1, ncol
224 >            
225 >             if (skipThisPair(i,j)) cycle inner
226 >            
227 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
228 >            
229 >             if (rijsq <  rlistsq) then            
230                  
231 <                fx = dudr * drdx
527 <                fy = dudr * drdy
528 <                fz = dudr * drdz
231 >                nlist = nlist + 1
232                  
233 < #ifdef IS_MPI
234 <                fCol(1,j) = fCol(1,j) - fx
235 <                fCol(2,j) = fCol(2,j) - fy
236 <                fCol(3,j) = fCol(3,j) - fz
233 >                if (nlist > neighborListSize) then
234 >                   call expandNeighborList(nlocal, listerror)
235 >                   if (listerror /= 0) then
236 >                      error = -1
237 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
238 >                      return
239 >                   end if
240 >                endif
241                  
242 <                fRow(1,j) = fRow(1,j) + fx
243 <                fRow(2,j) = fRow(2,j) + fy
244 <                fRow(3,j) = fRow(3,j) + fz
245 < #else
539 <                f(1,j) = f(1,j) - fx
540 <                f(2,j) = f(2,j) - fy
541 <                f(3,j) = f(3,j) - fz
542 <                f(1,i) = f(1,i) + fx
543 <                f(2,i) = f(2,i) + fy
544 <                f(3,i) = f(3,i) + fz
545 < #endif
546 <                
547 <                if (do_stress) then
548 <                   tauTemp(1) = tauTemp(1) + fx * rxij
549 <                   tauTemp(2) = tauTemp(2) + fx * ryij
550 <                   tauTemp(3) = tauTemp(3) + fx * rzij
551 <                   tauTemp(4) = tauTemp(4) + fy * rxij
552 <                   tauTemp(5) = tauTemp(5) + fy * ryij
553 <                   tauTemp(6) = tauTemp(6) + fy * rzij
554 <                   tauTemp(7) = tauTemp(7) + fz * rxij
555 <                   tauTemp(8) = tauTemp(8) + fz * ryij
556 <                   tauTemp(9) = tauTemp(9) + fz * rzij
242 >                list(nlist) = j
243 >                                
244 >                if (rijsq <  rcutsq) then
245 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
246                  endif
247               endif
248            enddo inner
249 <     enddo
249 >       enddo
250  
251 < #ifdef IS_MPI
252 <     point(nrow + 1) = nlist + 1
253 < #else
565 <     point(natoms) = nlist + 1
566 < #endif
251 >       point(nrow + 1) = nlist + 1
252 >      
253 >    else  !! (of update_check)
254  
255 <  else
255 >       ! use the list to find the neighbors
256 >       do i = 1, nrow
257 >          JBEG = POINT(i)
258 >          JEND = POINT(i+1) - 1
259 >          ! check thiat molecule i has neighbors
260 >          if (jbeg .le. jend) then
261 >            
262 >             do jnab = jbeg, jend
263 >                j = list(jnab)
264  
265 <     ! use the list to find the neighbors
266 <     do i = 1, nrow
572 <        JBEG = POINT(i)
573 <        JEND = POINT(i+1) - 1
574 <        ! check thiat molecule i has neighbors
575 <        if (jbeg .le. jend) then
576 < #ifdef IS_MPI
577 <           ljAtype_i => identPtrListRow(i)%this
578 <           rxi = qRow(1,i)
579 <           ryi = qRow(2,i)
580 <           rzi = qRow(3,i)
581 < #else
582 <           ljAtype_i   => identPtrList(i)%this
583 <           rxi = q(1,i)
584 <           ryi = q(2,i)
585 <           rzi = q(3,i)
586 < #endif
587 <           do jnab = jbeg, jend
588 <              j = list(jnab)
589 < #ifdef IS_MPI
590 <              ljAtype_j = identPtrListColumn(j)%this
591 <              rxij = wrap(rxi - qCol(1,j), 1)
592 <              ryij = wrap(ryi - qCol(2,j), 2)
593 <              rzij = wrap(rzi - qCol(3,j), 3)
594 < #else
595 <              ljAtype_j = identPtrList(j)%this
596 <              rxij = wrap(rxi - q(1,j), 1)
597 <              ryij = wrap(ryi - q(2,j), 2)
598 <              rzij = wrap(rzi - q(3,j), 3)
599 < #endif
600 <              rijsq = rxij*rxij + ryij*ryij + rzij*rzij
601 <              
602 <              if (rijsq <  rcutsq) then
265 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
266 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
267  
268 <                 r = dsqrt(rijsq)
269 <                
270 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
271 < #ifdef IS_MPI
272 <                eRow(i) = eRow(i) + pot*0.5
609 <                eCol(i) = eCol(i) + pot*0.5
268 >             enddo
269 >          endif
270 >       enddo
271 >    endif
272 >    
273   #else
274 <                pe = pe + pot
275 < #endif                
276 <  
277 <                drdx = -rxij / r
278 <                drdy = -ryij / r
279 <                drdz = -rzij / r
274 >    
275 >    if (update_nlist) then
276 >      
277 >       ! save current configuration, contruct neighbor list,
278 >       ! and calculate forces
279 >       call saveNeighborList(q)
280 >      
281 >       neighborListSize = getNeighborListSize()
282 >       nlist = 0
283 >      
284 >       do i = 1, natoms-1
285 >          point(i) = nlist + 1
286 >          
287 >          inner: do j = i+1, natoms
288 >            
289 >             if (skipThisPair(i,j)) cycle inner
290 >            
291 >             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
292 >          
293 >             if (rijsq <  rlistsq) then
294                  
295 <                fx = dudr * drdx
619 <                fy = dudr * drdy
620 <                fz = dudr * drdz
295 >                nlist = nlist + 1
296                  
297 < #ifdef IS_MPI
298 <                fCol(1,j) = fCol(1,j) - fx
299 <                fCol(2,j) = fCol(2,j) - fy
300 <                fCol(3,j) = fCol(3,j) - fz
301 <                
302 <                fRow(1,j) = fRow(1,j) + fx
303 <                fRow(2,j) = fRow(2,j) + fy
629 <                fRow(3,j) = fRow(3,j) + fz
630 < #else
631 <                f(1,j) = f(1,j) - fx
632 <                f(2,j) = f(2,j) - fy
633 <                f(3,j) = f(3,j) - fz
634 <                f(1,i) = f(1,i) + fx
635 <                f(2,i) = f(2,i) + fy
636 <                f(3,i) = f(3,i) + fz
637 < #endif
638 <                
639 <                if (do_stress) then
640 <                   tauTemp(1) = tauTemp(1) + fx * rxij
641 <                   tauTemp(2) = tauTemp(2) + fx * ryij
642 <                   tauTemp(3) = tauTemp(3) + fx * rzij
643 <                   tauTemp(4) = tauTemp(4) + fy * rxij
644 <                   tauTemp(5) = tauTemp(5) + fy * ryij
645 <                   tauTemp(6) = tauTemp(6) + fy * rzij
646 <                   tauTemp(7) = tauTemp(7) + fz * rxij
647 <                   tauTemp(8) = tauTemp(8) + fz * ryij
648 <                   tauTemp(9) = tauTemp(9) + fz * rzij
297 >                if (nlist > neighborListSize) then
298 >                   call expandList(natoms, listerror)
299 >                   if (listerror /= 0) then
300 >                      error = -1
301 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
302 >                      return
303 >                   end if
304                  endif
305                  
306 +                list(nlist) = j
307                  
308 +                if (rijsq <  rcutsq) then
309 +                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
310 +                endif
311               endif
312 <          enddo
313 <       endif
314 <    enddo
315 < endif
316 <
312 >          enddo inner
313 >       enddo
314 >      
315 >       point(natoms) = nlist + 1
316 >      
317 >    else !! (update)
318 >      
319 >       ! use the list to find the neighbors
320 >       do i = 1, natoms-1
321 >          JBEG = POINT(i)
322 >          JEND = POINT(i+1) - 1
323 >          ! check thiat molecule i has neighbors
324 >          if (jbeg .le. jend) then
325 >            
326 >             do jnab = jbeg, jend
327 >                j = list(jnab)
328  
329 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
330 +                call do_pair(i, j, rijsq, d, do_pot, do_stress)
331  
332 +             enddo
333 +          endif
334 +       enddo
335 +    endif
336 +    
337 + #endif
338 +    
339 +    ! phew, done with main loop.
340 +    
341   #ifdef IS_MPI
342      !!distribute forces
662 –
663 –    call scatter(fRow,f,plan_row3d)
664 –    call scatter(fCol,fMPITemp,plan_col3d)
665 –
666 –    do i = 1,nlocal
667 –       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
668 –    end do
669 –
670 –
671 –    call scatter(tRow,t,plan_row3d)
672 –    call scatter(tCol,tMPITemp,plan_col3d)
343      
344 +    call scatter(f_Row,f,plan_row3d)
345 +    call scatter(f_Col,f_temp,plan_col3d)
346      do i = 1,nlocal
347 <       t(1:3,i) = t(1:3,i) + tMPITemp(1:3,i)
347 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
348      end do
349 <
350 <
349 >    
350 >    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
351 >       call scatter(t_Row,t,plan_row3d)
352 >       call scatter(t_Col,t_temp,plan_col3d)
353 >      
354 >       do i = 1,nlocal
355 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
356 >       end do
357 >    endif
358 >    
359      if (do_pot) then
360         ! scatter/gather pot_row into the members of my column
361 <       call scatter(eRow,eTemp,plan_row)
361 >       call scatter(pot_Row, pot_Temp, plan_row)
362        
363         ! scatter/gather pot_local into all other procs
364         ! add resultant to get total pot
365         do i = 1, nlocal
366 <          pe_local = pe_local + eTemp(i)
366 >          pot_local = pot_local + pot_Temp(i)
367         enddo
688 –       if (newtons_thrd) then
689 –          eTemp = 0.0E0_DP
690 –          call scatter(eCol,eTemp,plan_col)
691 –          do i = 1, nlocal
692 –             pe_local = pe_local + eTemp(i)
693 –          enddo
694 –       endif
695 –       pe = pe_local
696 –    endif
368  
369 +       pot_Temp = 0.0_DP
370 +
371 +       call scatter(pot_Col, pot_Temp, plan_col)
372 +       do i = 1, nlocal
373 +          pot_local = pot_local + pot_Temp(i)
374 +       enddo
375 +      
376 +    endif    
377   #endif
378  
379 <    potE = pe
379 >    if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
380 >      
381 >       if (FF_uses_RF .and. SimUsesRF()) then
382 >          
383 > #ifdef IS_MPI
384 >          call scatter(rf_Row,rf,plan_row3d)
385 >          call scatter(rf_Col,rf_Temp,plan_col3d)
386 >          do i = 1,nlocal
387 >             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
388 >          end do
389 > #endif
390 >          
391 >          do i = 1, getNlocal()
392 >
393 >             rfpot = 0.0_DP
394 > #ifdef IS_MPI
395 >             me_i = atid_row(i)
396 > #else
397 >             me_i = atid(i)
398 > #endif
399 >             call getElementProperty(atypes, me_i, "is_DP", is_DP_i)      
400 >             if ( is_DP_i ) then
401 >                call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
402 >                !! The reaction field needs to include a self contribution
403 >                !! to the field:
404 >                call accumulate_self_rf(i, mu_i, u_l)            
405 >                !! Get the reaction field contribution to the
406 >                !! potential and torques:
407 >                call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
408 > #ifdef IS_MPI
409 >                pot_local = pot_local + rfpot
410 > #else
411 >                pot = pot + rfpot
412 > #endif
413 >             endif            
414 >          enddo
415 >       endif
416 >    endif
417  
418  
703 –    if (do_stress) then
419   #ifdef IS_MPI
420 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
420 >
421 >    if (do_pot) then
422 >       pot = pot_local
423 >       !! we assume the c code will do the allreduce to get the total potential
424 >       !! we could do it right here if we needed to...
425 >    endif
426 >
427 >    if (do_stress) then
428 >       call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
429              mpi_comm_world,mpi_err)
430 < #else
431 <       tau = tauTemp
709 < #endif      
430 >       call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
431 >            mpi_comm_world,mpi_err)
432      endif
433  
434 + #else
435  
436 +    if (do_stress) then
437 +       tau = tau_Temp
438 +       virial = virial_Temp
439 +    endif
440 +
441 + #endif
442 +    
443    end subroutine do_force_loop
444  
445 +  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress)
446  
447 +    real( kind = dp ) :: pot
448 +    real( kind = dp ), dimension(3,getNlocal()) :: u_l
449 +    real (kind=dp), dimension(9,getNlocal()) :: A
450 +    real (kind=dp), dimension(3,getNlocal()) :: f
451 +    real (kind=dp), dimension(3,getNlocal()) :: t
452  
453 +    logical, intent(inout) :: do_pot, do_stress
454 +    integer, intent(in) :: i, j
455 +    real ( kind = dp ), intent(inout)    :: rijsq
456 +    real ( kind = dp )                :: r
457 +    real ( kind = dp ), intent(inout) :: d(3)
458 +    logical :: is_LJ_i, is_LJ_j
459 +    logical :: is_DP_i, is_DP_j
460 +    logical :: is_Sticky_i, is_Sticky_j
461 +    integer :: me_i, me_j
462 +
463 +    r = sqrt(rijsq)
464 +    
465 + #ifdef IS_MPI
466 +
467 +    me_i = atid_row(i)
468 +    me_j = atid_col(j)
469 +
470 + #else
471 +
472 +    me_i = atid(i)
473 +    me_j = atid(j)
474 +
475 + #endif
476 +
477 +
478 +    if (FF_uses_LJ .and. SimUsesLJ()) then
479 +       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
480 +       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
481 +      
482 +       if ( is_LJ_i .and. is_LJ_j ) &
483 +            call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
484 +    endif
485 +      
486 +
487 +    if (FF_uses_dipoles .and. SimUsesDipoles()) then
488 +       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
489 +       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
490 +      
491 +       if ( is_DP_i .and. is_DP_j ) then
492 +          
493 +          call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
494 +          
495 +          if (FF_uses_RF .and. SimUsesRF()) then
496 +            
497 +             call accumulate_rf(i, j, r, u_l)
498 +             call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
499 +            
500 +          endif
501 +          
502 +       endif
503 +    endif
504 +
505 +    if (FF_uses_Sticky .and. SimUsesSticky()) then
506 +
507 +       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
508 +       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
509 +      
510 +       if ( is_Sticky_i .and. is_Sticky_j ) then
511 +          call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
512 +               do_pot, do_stress)
513 +       endif
514 +    endif
515 +      
516 +  end subroutine do_pair
517 +
518 +
519 +  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
520 +    
521 +    real (kind = dp), dimension(3) :: q_i
522 +    real (kind = dp), dimension(3) :: q_j
523 +    real ( kind = dp ), intent(out) :: r_sq
524 +    real( kind = dp ) :: d(3)
525 +
526 +    d(1:3) = q_i(1:3) - q_j(1:3)
527 +    
528 +    ! Wrap back into periodic box if necessary
529 +    if ( SimUsesPBC() ) then
530 +       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,box(1:3)) * &
531 +            int(abs(d(1:3)/box(1:3) + 0.5_dp))
532 +    endif
533 +    
534 +    r_sq = dot_product(d,d)
535 +        
536 +  end subroutine get_interatomic_vector
537 +
538 +  subroutine check_initialization(error)
539 +    integer, intent(out) :: error
540 +    
541 +    error = 0
542 +    ! Make sure we are properly initialized.
543 +    if (.not. do_forces_initialized) then
544 +       write(default_error,*) "ERROR: do_Forces has not been initialized!"
545 +       error = -1
546 +       return
547 +    endif
548 +
549 + #ifdef IS_MPI
550 +    if (.not. isMPISimSet()) then
551 +       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
552 +       error = -1
553 +       return
554 +    endif
555 + #endif
556 +    
557 +    return
558 +  end subroutine check_initialization
559 +
560 +  
561 +  subroutine zero_work_arrays()
562 +    
563 + #ifdef IS_MPI
564 +
565 +    q_Row = 0.0_dp
566 +    q_Col = 0.0_dp  
567 +    
568 +    u_l_Row = 0.0_dp
569 +    u_l_Col = 0.0_dp
570 +    
571 +    A_Row = 0.0_dp
572 +    A_Col = 0.0_dp
573 +    
574 +    f_Row = 0.0_dp
575 +    f_Col = 0.0_dp
576 +    f_Temp = 0.0_dp
577 +      
578 +    t_Row = 0.0_dp
579 +    t_Col = 0.0_dp
580 +    t_Temp = 0.0_dp
581 +
582 +    pot_Row = 0.0_dp
583 +    pot_Col = 0.0_dp
584 +    pot_Temp = 0.0_dp
585 +
586 +    rf_Row = 0.0_dp
587 +    rf_Col = 0.0_dp
588 +    rf_Temp = 0.0_dp
589 +
590 + #endif
591 +
592 +    rf = 0.0_dp
593 +    tau_Temp = 0.0_dp
594 +    virial_Temp = 0.0_dp
595 +    
596 +  end subroutine zero_work_arrays
597 +  
598 +  function skipThisPair(atom1, atom2) result(skip_it)
599 +    
600 +    integer, intent(in) :: atom1
601 +    integer, intent(in), optional :: atom2
602 +    logical :: skip_it
603 +    integer :: unique_id_1, unique_id_2
604 +    integer :: i
605 +
606 +    skip_it = .false.
607 +
608 +    !! there are a number of reasons to skip a pair or a particle
609 +    !! mostly we do this to exclude atoms who are involved in short
610 +    !! range interactions (bonds, bends, torsions), but we also need
611 +    !! to exclude some overcounted interactions that result from
612 +    !! the parallel decomposition
613 +    
614 + #ifdef IS_MPI
615 +    !! in MPI, we have to look up the unique IDs for each atom
616 +    unique_id_1 = tagRow(atom1)
617 + #else
618 +    !! in the normal loop, the atom numbers are unique
619 +    unique_id_1 = atom1
620 + #endif
621 +    
622 +    !! We were called with only one atom, so just check the global exclude
623 +    !! list for this atom
624 +    if (.not. present(atom2)) then
625 +       do i = 1, nExcludes_global
626 +          if (excludesGlobal(i) == unique_id_1) then
627 +             skip_it = .true.
628 +             return
629 +          end if
630 +       end do
631 +       return
632 +    end if
633 +    
634 + #ifdef IS_MPI
635 +    unique_id_2 = tagColumn(atom2)
636 + #else
637 +    unique_id_2 = atom2
638 + #endif
639 +    
640 + #ifdef IS_MPI
641 +    !! this situation should only arise in MPI simulations
642 +    if (unique_id_1 == unique_id_2) then
643 +       skip_it = .true.
644 +       return
645 +    end if
646 +    
647 +    !! this prevents us from doing the pair on multiple processors
648 +    if (unique_id_1 < unique_id_2) then
649 +       if (mod(unique_id_1 + unique_id_2,2) == 0) skip_it = .true.
650 +       return
651 +    else                
652 +       if (mod(unique_id_1 + unique_id_2,2) == 1) skip_it = .true.
653 +    endif
654 + #endif
655 +
656 +    !! the rest of these situations can happen in all simulations:
657 +    do i = 1, nExcludes_global      
658 +       if ((excludesGlobal(i) == unique_id_1) .or. &
659 +            (excludesGlobal(i) == unique_id_2)) then
660 +          skip_it = .true.
661 +          return
662 +       endif
663 +    enddo
664 +    
665 +    do i = 1, nExcludes_local
666 +       if (excludesLocal(1,i) == unique_id_1) then
667 +          if (excludesLocal(2,i) == unique_id_2) then
668 +             skip_it = .true.
669 +             return
670 +          endif
671 +       else
672 +          if (excludesLocal(1,i) == unique_id_2) then
673 +             if (excludesLocal(2,i) == unique_id_1) then
674 +                skip_it = .true.
675 +                return
676 +             endif
677 +          endif
678 +       endif
679 +    end do
680 +    
681 +    return
682 +  end function skipThisPair
683 +
684 +  function FF_UsesDirectionalAtoms() result(doesit)
685 +    logical :: doesit
686 +    doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
687 +         FF_uses_GB .or. FF_uses_RF
688 +  end function FF_UsesDirectionalAtoms
689 +  
690 +  function FF_RequiresPrepairCalc() result(doesit)
691 +    logical :: doesit
692 +    doesit = FF_uses_EAM
693 +  end function FF_RequiresPrepairCalc
694 +  
695 +  function FF_RequiresPostpairCalc() result(doesit)
696 +    logical :: doesit
697 +    doesit = FF_uses_RF
698 +  end function FF_RequiresPostpairCalc
699 +  
700   end module do_Forces

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines