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 293 by mmeineke, Thu Mar 6 16:07:57 2003 UTC vs.
Revision 331 by chuckv, Thu Mar 13 00:33:18 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.2 2003-03-06 16:07:55 mmeineke Exp $, $Date: 2003-03-06 16:07:55 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
7 > !! @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  
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
16 <  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  
20   #ifdef IS_MPI
21    use mpiSimulation
# Line 22 | Line 23 | module do_Forces
23    implicit none
24    PRIVATE
25  
26 < !! Number of lj_atypes in lj_atype_list
27 <  integer, save :: n_atypes = 0
26 >  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  
28 – !! Global list of lj atypes in simulation
29 –  type (atype), pointer :: ListHead => null()
30 –  type (atype), pointer :: ListTail => null()
34  
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_forceLoop
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
69 <    real( kind = dp ), intent(in) :: epsilon
70 <    real( kind = dp ), intent(in) :: sigma
71 <    real( kind = dp ), intent(in) :: w0
72 <    real( kind = dp ), intent(in) :: v0
73 <    real( kind = dp ), intent(in) :: dipoleMoment
40 >  subroutine init_FF(thisStat)
41 >    integer, intent(out) :: thisStat  
42 >    integer :: my_status
43 >    character(len = 100) :: mix_Policy
44  
45 <    integer, intent(in) :: ident
46 <    integer, intent(out) :: status
47 <    integer, intent(in) :: is_Sticky
48 <    integer, intent(in) :: is_DP
49 <    integer, intent(in) :: is_GB
50 <    integer, intent(in) :: is_LJ
81 <
82 <
83 <    type (atype), pointer :: the_new_atype
84 <    integer :: alloc_error
85 <    integer :: atype_counter = 0
86 <    integer :: alloc_size
87 <    integer :: err_stat
88 <    status = 0
89 <
90 <
91 <
92 < ! allocate a new atype    
93 <    allocate(the_new_atype,stat=alloc_error)
94 <    if (alloc_error /= 0 ) then
95 <       status = -1
45 >    ! be a smarter subroutine.
46 >    mix_Policy = "FIXME"
47 >    thisStat = 0
48 >    call init_lj_FF(mix_Policy,my_status)
49 >    if (my_status /= 0) then
50 >       thisStat = -1
51         return
52      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 –
53      
54 < ! assume that this atype will be successfully added
55 <    the_new_atype%atype_ident = ident
56 <    the_new_atype%atype_number = n_lj_atypes + 1
54 >    call check_sticky_FF(my_status)
55 >    if (my_status /= 0) then
56 >       thisStat = -1
57 >       return
58 >    end if
59      
60 <    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
61 <    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
62 <    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
118 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
60 >    do_forces_initialized = .true.    
61 >    
62 >  end subroutine init_FF
63  
120 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
121 –    if (err_stat /= 0 ) then
122 –       status = -1
123 –       return
124 –    endif
64  
126 –    n_atypes = n_atypes + 1
65  
66 +  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
67 +  !------------------------------------------------------------->
68 +  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
69 +       error)
70 +    !! Position array provided by C, dimensioned by getNlocal
71 +    real ( kind = dp ), dimension(3,getNlocal()) :: q
72 +    !! Rotation Matrix for each long range particle in simulation.
73 +    real( kind = dp), dimension(9,getNlocal()) :: A    
74 +    !! Unit vectors for dipoles (lab frame)
75 +    real( kind = dp ), dimension(3,getNlocal()) :: u_l
76 +    !! Force array provided by C, dimensioned by getNlocal
77 +    real ( kind = dp ), dimension(3,getNlocal()) :: f
78 +    !! Torsion array provided by C, dimensioned by getNlocal
79 +    real( kind = dp ), dimension(3,getNlocal()) :: t    
80 +    !! Stress Tensor
81 +    real( kind = dp), dimension(9) :: tau  
82 +    real ( kind = dp ) :: pot
83 +    logical ( kind = 2) :: do_pot_c, do_stress_c
84 +    logical :: do_pot
85 +    logical :: do_stress
86 + #ifdef IS_MPI
87 +    real( kind = DP ) :: pot_local
88 +    integer :: nrow
89 +    integer :: ncol
90 + #endif
91 +    integer :: nlocal
92 +    integer :: natoms    
93 +    logical :: update_nlist  
94 +    integer :: i, j, jbeg, jend, jnab
95 +    integer :: nlist
96 +    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
97 +    real(kind=dp),dimension(3) :: d
98 +    real(kind=dp) :: rfpot, mu_i, virial
99 +    integer :: me_i
100 +    logical :: is_dp_i
101 +    integer :: neighborListSize
102 +    integer :: listerror, error
103 +    integer :: localError
104  
105 <  end subroutine new_atype
105 >    !! initialize local variables  
106  
107 + #ifdef IS_MPI
108 +    nlocal = getNlocal()
109 +    nrow   = getNrow(plan_row)
110 +    ncol   = getNcol(plan_col)
111 + #else
112 +    nlocal = getNlocal()
113 +    natoms = nlocal
114 + #endif
115  
116 <  subroutine init_FF(nComponents,ident, status)
117 < !! Number of components in ident array
118 <    integer, intent(inout) :: nComponents
119 < !! Array of identities nComponents long corresponding to
120 < !! ljatype ident.
121 <    integer, dimension(nComponents),intent(inout) :: ident
122 < !!  Result status, success = 0, error = -1
123 <    integer, intent(out) :: Status
116 >    call getRcut(rcut,rc2=rcutsq)
117 >    call getRlist(rlist,rlistsq)
118 >    
119 >    call check_initialization(localError)
120 >    if ( localError .ne. 0 ) then
121 >       error = -1
122 >       return
123 >    end if
124 >    call zero_work_arrays()
125  
126 <    integer :: alloc_stat
126 >    do_pot = do_pot_c
127 >    do_stress = do_stress_c
128  
129 <    integer :: thisStat
130 <    integer :: i
129 >    ! Gather all information needed by all force loops:
130 >    
131 > #ifdef IS_MPI    
132  
133 <    integer :: myNode
134 < #ifdef IS_MPI
135 <    integer, allocatable, dimension(:) :: identRow
136 <    integer, allocatable, dimension(:) :: identCol
137 <    integer :: nrow
138 <    integer :: ncol
133 >    call gather(q,q_Row,plan_row3d)
134 >    call gather(q,q_Col,plan_col3d)
135 >        
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 >    
144   #endif
153 –    status = 0
154 –  
155 –
145      
146 +    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 + #ifdef IS_MPI
159 +    
160 +    if (update_nlist) then
161 +      
162 +       !! save current configuration, construct neighbor list,
163 +       !! and calculate forces
164 +       call save_neighborList(q)
165 +      
166 +       neighborListSize = getNeighborListSize()
167 +       nlist = 0      
168 +      
169 +       do i = 1, nrow
170 +          point(i) = nlist + 1
171 +          
172 +          inner: do j = 1, ncol
173 +            
174 +             if (checkExcludes(i,j)) cycle inner
175 +            
176 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
177 +            
178 +             if (rijsq <  rlistsq) then            
179 +                
180 +                nlist = nlist + 1
181 +                
182 +                if (nlist > neighborListSize) then
183 +                   call expandNeighborList(nlocal, listerror)
184 +                   if (listerror /= 0) then
185 +                      error = -1
186 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
187 +                      return
188 +                   end if
189 +                endif
190 +                
191 +                list(nlist) = j
192 +                                
193 +                if (rijsq <  rcutsq) then
194 +                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
195 +                endif
196 +             endif
197 +          enddo inner
198 +       enddo
199  
200 < !! if were're not in MPI, we just update ljatypePtrList
201 < #ifndef IS_MPI
202 <    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
203 <    if ( thisStat /= 0 ) then
204 <       status = -1
205 <       return
200 >       point(nrow + 1) = nlist + 1
201 >      
202 >    else  !! (of update_check)
203 >
204 >       ! 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 >            
211 >             do jnab = jbeg, jend
212 >                j = list(jnab)
213 >
214 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
215 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
216 >
217 >             enddo
218 >          endif
219 >       enddo
220      endif
165 –
221      
167 – ! if were're in MPI, we also have to worry about row and col lists    
222   #else
223 <  
224 < ! We can only set up forces if mpiSimulation has been setup.
225 <    if (.not. isMPISimSet()) then
226 <       write(default_error,*) "MPI is not set"
227 <       status = -1
228 <       return
229 <    endif
230 <    nrow = getNrow(plan_row)
231 <    ncol = getNcol(plan_col)
232 <    mynode = getMyNode()
233 < !! Allocate temperary arrays to hold gather information
234 <    allocate(identRow(nrow),stat=alloc_stat)
235 <    if (alloc_stat /= 0 ) then
236 <       status = -1
237 <       return
238 <    endif
223 >    
224 >    if (update_nlist) then
225 >      
226 >       ! save current configuration, contruct neighbor list,
227 >       ! and calculate forces
228 >       call save_neighborList(q)
229 >      
230 >       neighborListSize = getNeighborListSize()
231 >       nlist = 0
232 >      
233 >       do i = 1, natoms-1
234 >          point(i) = nlist + 1
235 >          
236 >          inner: do j = i+1, natoms
237 >            
238 >             if (checkExcludes(i,j)) cycle inner
239 >            
240 >             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
241 >          
242 >             if (rijsq <  rlistsq) then
243 >                
244 >                nlist = nlist + 1
245 >                
246 >                if (nlist > neighborListSize) then
247 >                   call expandList(natoms, listerror)
248 >                   if (listerror /= 0) then
249 >                      error = -1
250 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
251 >                      return
252 >                   end if
253 >                endif
254 >                
255 >                list(nlist) = j
256 >                
257 >                if (rijsq <  rcutsq) then
258 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
259 >                endif
260 >             endif
261 >          enddo inner
262 >       enddo
263 >      
264 >       point(natoms) = nlist + 1
265 >      
266 >    else !! (update)
267 >      
268 >       ! use the list to find the neighbors
269 >       do i = 1, natoms-1
270 >          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  
278 <    allocate(identCol(ncol),stat=alloc_stat)
279 <    if (alloc_stat /= 0 ) then
188 <       status = -1
189 <       return
190 <    endif
278 >                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
279 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
280  
281 < !! Gather idents into row and column idents
282 <
283 <    call gather(ident,identRow,plan_row)
284 <    call gather(ident,identCol,plan_col)
281 >             enddo
282 >          endif
283 >       enddo
284 >    endif
285      
286 <  
287 < !! Create row and col pointer lists
288 <  
289 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
290 <    if (thisStat /= 0 ) then
291 <       status = -1
292 <       return
286 > #endif
287 >    
288 >    ! phew, done with main loop.
289 >    
290 > #ifdef IS_MPI
291 >    !!distribute forces
292 >    
293 >    call scatter(f_Row,f,plan_row3d)
294 >    call scatter(f_Col,f_temp,plan_col3d)
295 >    do i = 1,nlocal
296 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
297 >    end do
298 >    
299 >    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
300 >       call scatter(t_Row,t,plan_row3d)
301 >       call scatter(t_Col,t_temp,plan_col3d)
302 >      
303 >       do i = 1,nlocal
304 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
305 >       end do
306      endif
307 <  
308 <    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
309 <    if (thisStat /= 0 ) then
310 <       status = -1
311 <       return
312 <    endif
307 >    
308 >    if (do_pot) then
309 >       ! scatter/gather pot_row into the members of my column
310 >       call scatter(pot_Row, pot_Temp, plan_row)
311 >      
312 >       ! scatter/gather pot_local into all other procs
313 >       ! add resultant to get total pot
314 >       do i = 1, nlocal
315 >          pot_local = pot_local + pot_Temp(i)
316 >       enddo
317  
318 < !! free temporary ident arrays
213 <    if (allocated(identCol)) then
214 <       deallocate(identCol)
215 <    end if
216 <    if (allocated(identCol)) then
217 <       deallocate(identRow)
218 <    endif
318 >       pot_Temp = 0.0_DP
319  
320 +       call scatter(pot_Col, pot_Temp, plan_col)
321 +       do i = 1, nlocal
322 +          pot_local = pot_local + pot_Temp(i)
323 +       enddo
324 +      
325 +    endif    
326   #endif
221 –    
222 –    call initForce_Modules(thisStat)
223 –    if (thisStat /= 0) then
224 –       status = -1
225 –       return
226 –    endif
327  
328 < !! Create neighbor lists
329 <    call expandList(thisStat)
330 <    if (thisStat /= 0) then
331 <       status = -1
332 <       return
328 >    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 >             rfpot = 0.0_DP
343 > #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 >                call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
357 > #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  
235 –    isFFinit = .true.
367  
368 + #ifdef IS_MPI
369  
370 <  end subroutine init_FF
370 >    if (do_pot) then
371 >       pot = pot_local
372 >       !! 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 >    endif
375  
376 +    if (do_stress) then
377 +       call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
378 +            mpi_comm_world,mpi_err)
379 +       call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
380 +            mpi_comm_world,mpi_err)
381 +    endif
382  
383 + #else
384  
385 +    if (do_stress) then
386 +       tau = tau_Temp
387 +       virial = virial_Temp
388 +    endif
389  
390 <  subroutine initForce_Modules(thisStat)
244 <    integer, intent(out) :: thisStat
245 <    integer :: my_status
390 > #endif
391      
392 <    thisStat = 0
248 <    call init_lj_FF(ListHead,my_status)
249 <    if (my_status /= 0) then
250 <       thisStat = -1
251 <       return
252 <    end if
392 >  end subroutine do_force_loop
393  
254 –  end subroutine initForce_Modules
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 < !! FORCE routine Calculates Lennard Jones forces.
260 < !------------------------------------------------------------->
261 <  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
262 < !! Position array provided by C, dimensioned by getNlocal
263 <    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
401 >  end subroutine do_preForce
402  
403 <  !! Magnitude dipole moment
404 <    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
403 >  !! Calculate any post force loop components, i.e. reaction field, etc.
404 >  subroutine do_postForce()
405  
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
406  
283 –    
284 –    type(atype), pointer :: Atype_i
285 –    type(atype), pointer :: Atype_j
407  
408 +  end subroutine do_postForce
409  
410 +  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress)
411  
412 <
413 <  
414 <
415 < #ifdef IS_MPI
416 <  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
412 >    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 <  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
418 >    logical, intent(inout) :: do_pot, do_stress
419 >    integer, intent(in) :: i, j
420 >    real ( kind = dp ), intent(inout)    :: rijsq
421 >    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  
428 < #endif
429 <
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.
428 >    r = sqrt(rijsq)
429 >    
430   #ifdef IS_MPI
346 –  logical :: newtons_thrd = .true.
347 –  real( kind = dp ) :: pe_local
348 –  integer :: nlocal
349 – #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.
431  
432 +    me_i = atid_row(i)
433 +    me_j = atid_col(j)
434  
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
435   #else
386 –  nrow = getNrow(plan_row)
387 –  ncol = getNcol(plan_col)
388 –  nlocal = natoms
389 –  j_start = 1
390 – #endif
436  
437 <  
438 < !! 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
437 >    me_i = atid(i)
438 >    me_j = atid(j)
439  
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
440   #endif
441  
420 – ! communicate MPI positions
421 – #ifdef IS_MPI    
422 –    call gather(q,qRow,plan_row3d)
423 –    call gather(q,qCol,plan_col3d)
442  
443 <    call gather(mu,muRow,plan_row3d)
444 <    call gather(mu,muCol,plan_col3d)
443 >    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  
452 <    call gather(u_l,u_lRow,plan_row3d)
453 <    call gather(u_l,u_lCol,plan_col3d)
452 >    if (FF_uses_dipoles .and. SimUsesDipoles()) then
453 >       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 >       endif
468 >    endif
469  
470 <    call gather(A,ARow,plan_row_rotation)
432 <    call gather(A,ACol,plan_col_rotation)
470 >    if (FF_uses_Sticky .and. SimUsesSticky()) then
471  
472 < #endif
472 >       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 >    endif
480 >      
481 >  end subroutine do_pair
482  
483  
484 <  if (update_nlist) then
484 >  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
485 >    
486 >    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 <     ! save current configuration, contruct neighbor list,
440 <     ! and calculate forces
441 <     call save_neighborList(q)
442 <    
443 <     neighborListSize = getNeighborListSize()
444 <     nlist = 0
445 <    
491 >    d(1:3) = q_i(1:3) - q_j(1:3)
492      
493 +    ! Wrap back into periodic box if necessary
494 +    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 +    endif
498 +    
499 +    r_sq = dot_product(d,d)
500 +        
501 +  end subroutine get_interatomic_vector
502  
503 <     do i = 1, nrow
504 <        point(i) = nlist + 1
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 <        ljAtype_i => identPtrListRow(i)%this
515 <        tag_i = tagRow(i)
516 <        rxi = qRow(1,i)
517 <        ryi = qRow(2,i)
518 <        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)
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 <        inner: do j = j_start, ncol
522 < #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
521 >    return
522 >  end subroutine check_initialization
523  
524 <           rxij = wrap(rxi - qCol(1,j), 1)
525 <           ryij = wrap(ryi - qCol(2,j), 2)
526 <           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 <
524 >  
525 >  subroutine zero_work_arrays()
526 >    
527   #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
528  
529 +    q_Row = 0.0_dp
530 +    q_Col = 0.0_dp  
531      
532 <              if (rijsq <  rcutsq) then
533 <                
534 <                 r = dsqrt(rijsq)
535 <      
536 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
537 <      
538 < #ifdef IS_MPI
539 <                eRow(i) = eRow(i) + pot*0.5
540 <                eCol(i) = eCol(i) + pot*0.5
541 < #else
542 <                    pe = pe + pot
543 < #endif                
544 <            
545 <                drdx = -rxij / r
546 <                drdy = -ryij / r
547 <                drdz = -rzij / r
548 <                
526 <                fx = dudr * drdx
527 <                fy = dudr * drdy
528 <                fz = dudr * drdz
529 <                
530 < #ifdef IS_MPI
531 <                fCol(1,j) = fCol(1,j) - fx
532 <                fCol(2,j) = fCol(2,j) - fy
533 <                fCol(3,j) = fCol(3,j) - fz
534 <                
535 <                fRow(1,j) = fRow(1,j) + fx
536 <                fRow(2,j) = fRow(2,j) + fy
537 <                fRow(3,j) = fRow(3,j) + fz
538 < #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
557 <                endif
558 <             endif
559 <          enddo inner
560 <     enddo
532 >    u_l_Row = 0.0_dp
533 >    u_l_Col = 0.0_dp
534 >    
535 >    A_Row = 0.0_dp
536 >    A_Col = 0.0_dp
537 >    
538 >    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 >
546 >    pot_Row = 0.0_dp
547 >    pot_Col = 0.0_dp
548 >    pot_Temp = 0.0_dp
549  
562 – #ifdef IS_MPI
563 –     point(nrow + 1) = nlist + 1
564 – #else
565 –     point(natoms) = nlist + 1
550   #endif
551  
552 <  else
552 >    tau_Temp = 0.0_dp
553 >    virial_Temp = 0.0_dp
554 >    
555 >  end subroutine zero_work_arrays
556 >  
557  
558 <     ! use the list to find the neighbors
559 <     do i = 1, nrow
560 <        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
558 >  !! 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 <                 r = dsqrt(rijsq)
563 <                
564 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
562 >  function checkExcludes(atom1,atom2) result(do_cycle)
563 >    !--------------- Arguments--------------------------
564 >    ! Index i
565 >    integer,intent(in) :: atom1
566 >    ! Index j
567 >    integer,intent(in), optional :: atom2
568 >    ! Result do_cycle
569 >    logical :: do_cycle
570 >    !--------------- Local variables--------------------
571 >    integer :: tag_i
572 >    integer :: tag_j
573 >    integer :: i, j
574 >    !--------------- END DECLARATIONS------------------  
575 >    do_cycle = .false.
576 >    
577   #ifdef IS_MPI
578 <                eRow(i) = eRow(i) + pot*0.5
609 <                eCol(i) = eCol(i) + pot*0.5
578 >    tag_i = tagRow(atom1)
579   #else
580 <                pe = pe + pot
612 < #endif                
613 <  
614 <                drdx = -rxij / r
615 <                drdy = -ryij / r
616 <                drdz = -rzij / r
617 <                
618 <                fx = dudr * drdx
619 <                fy = dudr * drdy
620 <                fz = dudr * drdz
621 <                
622 < #ifdef IS_MPI
623 <                fCol(1,j) = fCol(1,j) - fx
624 <                fCol(2,j) = fCol(2,j) - fy
625 <                fCol(3,j) = fCol(3,j) - fz
626 <                
627 <                fRow(1,j) = fRow(1,j) + fx
628 <                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
580 >    tag_i = tag(atom1)
581   #endif
582 <                
583 <                if (do_stress) then
584 <                   tauTemp(1) = tauTemp(1) + fx * rxij
585 <                   tauTemp(2) = tauTemp(2) + fx * ryij
586 <                   tauTemp(3) = tauTemp(3) + fx * rzij
587 <                   tauTemp(4) = tauTemp(4) + fy * rxij
588 <                   tauTemp(5) = tauTemp(5) + fy * ryij
589 <                   tauTemp(6) = tauTemp(6) + fy * rzij
590 <                   tauTemp(7) = tauTemp(7) + fz * rxij
591 <                   tauTemp(8) = tauTemp(8) + fz * ryij
592 <                   tauTemp(9) = tauTemp(9) + fz * rzij
649 <                endif
650 <                
651 <                
652 <             endif
653 <          enddo
654 <       endif
655 <    enddo
656 < endif
657 <
582 >    
583 >    !! Check global excludes first
584 >    if (.not. present(atom2)) then
585 >       do i = 1, nExcludes_global
586 >          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 <
595 < #ifdef IS_MPI
661 <    !!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)
594 >    !! we return if atom2 not present here.
595 >    tag_j = tagColumn(atom2)
596      
597 <    do i = 1,nlocal
598 <       t(1:3,i) = t(1:3,i) + tMPITemp(1:3,i)
597 >    if (tag_i == tag_j) then
598 >       do_cycle = .true.
599 >       return
600 >    end if
601 >    
602 >    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 >            
609 >    do i = 1, nExcludes_local
610 >       if ((tag_i == excludesLocal(1,i)) .and. (excludesLocal(2,i) < 0)) then
611 >          do_cycle = .true.
612 >          return
613 >       end if
614      end do
615 +    
616 +    
617 +  end function checkExcludes
618  
619 <
620 <    if (do_pot) then
621 <       ! scatter/gather pot_row into the members of my column
622 <       call scatter(eRow,eTemp,plan_row)
623 <      
624 <       ! scatter/gather pot_local into all other procs
625 <       ! add resultant to get total pot
626 <       do i = 1, nlocal
627 <          pe_local = pe_local + eTemp(i)
628 <       enddo
629 <       if (newtons_thrd) then
630 <          eTemp = 0.0E0_DP
631 <          call scatter(eCol,eTemp,plan_col)
632 <          do i = 1, nlocal
633 <             pe_local = pe_local + eTemp(i)
634 <          enddo
694 <       endif
695 <       pe = pe_local
696 <    endif
697 <
698 < #endif
699 <
700 <    potE = pe
701 <
702 <
703 <    if (do_stress) then
704 < #ifdef IS_MPI
705 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
706 <            mpi_comm_world,mpi_err)
707 < #else
708 <       tau = tauTemp
709 < #endif      
710 <    endif
711 <
712 <
713 <  end subroutine do_force_loop
714 <
715 <
716 <
619 >  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   end module do_Forces

Diff Legend

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