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 295 by chuckv, Thu Mar 6 19:57:03 2003 UTC vs.
Revision 328 by gezelter, Wed Mar 12 20:00:58 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.3 2003-03-06 19:57:03 chuckv Exp $, $Date: 2003-03-06 19:57:03 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
7 > !! @version $Id: do_Forces.F90,v 1.14 2003-03-12 20:00:58 gezelter Exp $, $Date: 2003-03-12 20:00:58 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $
8  
9  
10  
11   module do_Forces
12    use simulation
13    use definitions
14 <  use generic_atypes
15 <  use neighborLists
16 <  
14 >  use forceGlobals
15 >  use atype_module
16 >  use neighborLists  
17    use lj_FF
18    use sticky_FF
19 <  use dp_FF
19 >  use dipole_dipole
20    use gb_FF
21  
22   #ifdef IS_MPI
# Line 22 | Line 25 | module do_Forces
25    implicit none
26    PRIVATE
27  
28 < !! Number of lj_atypes in lj_atype_list
26 <  integer, save :: n_atypes = 0
28 >  public :: do_force_loop
29  
30 < !! Global list of lj atypes in simulation
31 <  type (atype), pointer :: ListHead => null()
30 <  type (atype), pointer :: ListTail => null()
30 >  logical :: do_pot
31 >  logical :: do_stress
32  
33 + contains
34  
35 +  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
36 +  !------------------------------------------------------------->
37 +  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
38 +       FFerror)
39 +    !! Position array provided by C, dimensioned by getNlocal
40 +    real ( kind = dp ), dimension(3,getNlocal()) :: q
41 +    !! Rotation Matrix for each long range particle in simulation.
42 +    real( kind = dp), dimension(9,getNlocal()) :: A    
43 +    !! Unit vectors for dipoles (lab frame)
44 +    real( kind = dp ), dimension(3,getNlocal()) :: u_l
45 +    !! Force array provided by C, dimensioned by getNlocal
46 +    real ( kind = dp ), dimension(3,getNlocal()) :: f
47 +    !! Torsion array provided by C, dimensioned by getNlocal
48 +    real( kind = dp ), dimension(3,getNlocal()) :: t    
49 +    !! Stress Tensor
50 +    real( kind = dp), dimension(9) :: tau  
51 +    real ( kind = dp ) :: pot
52 +    logical ( kind = 2) :: do_pot_c, do_stress_c
53 +    integer :: FFerror
54  
55 + #ifdef IS_MPI
56 +    real( kind = DP ) :: pot_local
57 + #endif
58 +    
59 +    logical :: update_nlist  
60 +    integer :: i, j, jbeg, jend, jnab
61 +    integer :: nlist
62 +    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
63  
35 –  logical, save :: firstTime = .True.
36 –
37 – !! Atype identity pointer lists
64   #ifdef IS_MPI
65 < !! Row lj_atype pointer list
40 <  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
41 < !! Column lj_atype pointer list
42 <  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
43 < #else
44 <  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
65 >    integer :: nlocal
66   #endif
67 +    integer :: nrow
68 +    integer :: ncol
69 +    integer :: natoms
70 +    integer :: neighborListSize
71 +    integer :: listerror
72 +    FFerror = 0
73  
74 +    do_pot = do_pot_c
75 +    do_stress = do_stress_c
76  
77 < !! Logical has lj force field module been initialized?
78 <  logical, save :: isFFinit = .false.
79 <
80 < !! Use periodic boundry conditions
81 <  logical :: wrap = .false.
82 <
54 < !! Potential energy global module variables
77 >    ! Make sure we are properly initialized.
78 >    if (.not. isFFInit) then
79 >       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
80 >       FFerror = -1
81 >       return
82 >    endif
83   #ifdef IS_MPI
84 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
85 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
84 >    if (.not. isMPISimSet()) then
85 >       write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
86 >       FFerror = -1
87 >       return
88 >    endif
89 > #endif
90 >    
91 >    !! initialize local variables  
92 >    natoms = getNlocal()
93 >    call getRcut(rcut,rcut2=rcutsq)
94 >    call getRlist(rlist,rlistsq)
95 >    
96 >    !! See if we need to update neighbor lists
97 >    call checkNeighborList(natoms, q, rcut, rlist, update_nlist)
98 >    
99 >    !--------------WARNING...........................
100 >    ! Zero variables, NOTE:::: Forces are zeroed in C
101 >    ! Zeroing them here could delete previously computed
102 >    ! Forces.
103 >    !------------------------------------------------
104 >    call zero_module_variables()
105  
106 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
107 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
108 <
109 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
110 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
111 <
112 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: ARow = 0.0_dp
113 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp
114 <
115 <  
69 <
70 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
71 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
72 <
73 <
74 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
75 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
76 <
77 <
78 <
79 <  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
80 <  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
81 <
82 <  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
106 >    ! communicate MPI positions
107 > #ifdef IS_MPI    
108 >    call gather(q,q_Row,plan_row3d)
109 >    call gather(q,q_Col,plan_col3d)
110 >    
111 >    call gather(u_l,u_l_Row,plan_row3d)
112 >    call gather(u_l,u_l_Col,plan_col3d)
113 >    
114 >    call gather(A,A_Row,plan_row_rotation)
115 >    call gather(A,A_Col,plan_col_rotation)
116   #endif
84 –  real(kind = dp) :: pe = 0.0_dp
85 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
86 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp
87 –  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp
117  
89 –  logical :: do_preForce  = .false.
90 –  logical :: do_postForce = .false.
118  
119 + #ifdef IS_MPI
120 +    
121 +    if (update_nlist) then
122 +      
123 +       ! save current configuration, contruct neighbor list,
124 +       ! and calculate forces
125 +       call save_neighborList(q)
126 +      
127 +       neighborListSize = getNeighborListSize()
128 +       nlist = 0
129 +      
130 +       nrow = getNrow(plan_row)
131 +       ncol = getNcol(plan_col)
132 +       nlocal = getNlocal()
133 +      
134 +       do i = 1, nrow
135 +          point(i) = nlist + 1
136 +          
137 +          inner: do j = 1, ncol
138 +            
139 +             if (check_exclude(i,j)) cycle inner:
140  
141 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
142 +            
143 +             if (rijsq <  rlistsq) then            
144 +                
145 +                nlist = nlist + 1
146 +                
147 +                if (nlist > neighborListSize) then
148 +                   call expandNeighborList(nlocal, listerror)
149 +                   if (listerror /= 0) then
150 +                      FFerror = -1
151 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
152 +                      return
153 +                   end if
154 +                endif
155 +                
156 +                list(nlist) = j
157 +                                
158 +                if (rijsq <  rcutsq) then
159 +                   call do_pair(i, j, rijsq, d)
160 +                endif
161 +             endif
162 +          enddo inner
163 +       enddo
164  
165 < !! Public methods and data
166 <  public :: new_atype
167 <  public :: do_forceLoop
97 <  public :: init_FF
165 >       point(nrow + 1) = nlist + 1
166 >      
167 >    else !! (update)
168  
169 <  
169 >       ! use the list to find the neighbors
170 >       do i = 1, nrow
171 >          JBEG = POINT(i)
172 >          JEND = POINT(i+1) - 1
173 >          ! check thiat molecule i has neighbors
174 >          if (jbeg .le. jend) then
175 >            
176 >             do jnab = jbeg, jend
177 >                j = list(jnab)
178  
179 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
180 +                call do_pair(i, j, rijsq, d)
181  
182 < contains
182 >             enddo
183 >          endif
184 >       enddo
185 >    endif
186  
187 < !! Adds a new lj_atype to the list.
188 <  subroutine new_atype(ident,mass,epsilon,sigma, &
189 <       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
190 <    real( kind = dp ), intent(in) :: mass
191 <    real( kind = dp ), intent(in) :: epsilon
192 <    real( kind = dp ), intent(in) :: sigma
193 <    real( kind = dp ), intent(in) :: w0
194 <    real( kind = dp ), intent(in) :: v0
195 <    real( kind = dp ), intent(in) :: dipoleMoment
187 > #else
188 >    
189 >    if (update_nlist) then
190 >      
191 >       ! save current configuration, contruct neighbor list,
192 >       ! and calculate forces
193 >       call save_neighborList(q)
194 >      
195 >       neighborListSize = getNeighborListSize()
196 >       nlist = 0
197 >          
198 >       do i = 1, natoms-1
199 >          point(i) = nlist + 1
200  
201 <    integer, intent(in) :: ident
115 <    integer, intent(out) :: status
116 <    integer, intent(in) :: is_Sticky
117 <    integer, intent(in) :: is_DP
118 <    integer, intent(in) :: is_GB
119 <    integer, intent(in) :: is_LJ
201 >          inner: do j = i+1, natoms
202  
203 +             if (check_exclude(i,j)) cycle inner:
204  
205 <    type (atype), pointer :: the_new_atype
206 <    integer :: alloc_error
207 <    integer :: atype_counter = 0
208 <    integer :: alloc_size
209 <    integer :: err_stat
210 <    status = 0
211 <
205 >             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
206 >          
207 >             if (rijsq <  rlistsq) then
208 >                
209 >                nlist = nlist + 1
210 >                
211 >                if (nlist > neighborListSize) then
212 >                   call expandList(natoms, listerror)
213 >                   if (listerror /= 0) then
214 >                      FFerror = -1
215 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
216 >                      return
217 >                   end if
218 >                endif
219 >                
220 >                list(nlist) = j
221 >                    
222 >                if (rijsq <  rcutsq) then
223 >                   call do_pair(i, j, rijsq, d)
224 >                endif
225 >             endif
226 >          enddo inner
227 >       enddo
228 >      
229 >       point(natoms) = nlist + 1
230 >      
231 >    else !! (update)
232 >      
233 >       ! use the list to find the neighbors
234 >       do i = 1, nrow
235 >          JBEG = POINT(i)
236 >          JEND = POINT(i+1) - 1
237 >          ! check thiat molecule i has neighbors
238 >          if (jbeg .le. jend) then
239 >            
240 >             do jnab = jbeg, jend
241 >                j = list(jnab)
242  
243 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
244 +                call do_pair(i, j, rijsq, d)
245  
246 < ! allocate a new atype    
247 <    allocate(the_new_atype,stat=alloc_error)
248 <    if (alloc_error /= 0 ) then
134 <       status = -1
135 <       return
136 <    end if
137 <
138 < ! assign our new atype information
139 <    the_new_atype%mass        = mass
140 <    the_new_atype%epsilon     = epsilon
141 <    the_new_atype%sigma       = sigma
142 <    the_new_atype%sigma2      = sigma * sigma
143 <    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
144 <         * the_new_atype%sigma2
145 <    the_new_atype%w0       = w0
146 <    the_new_atype%v0       = v0
147 <    the_new_atype%dipoleMoment       = dipoleMoment
148 <
149 <    
150 < ! assume that this atype will be successfully added
151 <    the_new_atype%atype_ident = ident
152 <    the_new_atype%atype_number = n_lj_atypes + 1
153 <    
154 <    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
155 <    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
156 <    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
157 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
158 <
159 <    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
160 <    if (err_stat /= 0 ) then
161 <       status = -1
162 <       return
246 >             enddo
247 >          endif
248 >       enddo
249      endif
164 –
165 –    n_atypes = n_atypes + 1
166 –
167 –
168 –  end subroutine new_atype
169 –
170 –
171 –  subroutine init_FF(nComponents,ident, status)
172 – !! Number of components in ident array
173 –    integer, intent(inout) :: nComponents
174 – !! Array of identities nComponents long corresponding to
175 – !! ljatype ident.
176 –    integer, dimension(nComponents),intent(inout) :: ident
177 – !!  Result status, success = 0, error = -1
178 –    integer, intent(out) :: Status
179 –
180 –    integer :: alloc_stat
181 –
182 –    integer :: thisStat
183 –    integer :: i
184 –
185 –    integer :: myNode
186 – #ifdef IS_MPI
187 –    integer, allocatable, dimension(:) :: identRow
188 –    integer, allocatable, dimension(:) :: identCol
189 –    integer :: nrow
190 –    integer :: ncol
191 – #endif
192 –    status = 0
193 –  
194 –
250      
251 <
197 < !! if were're not in MPI, we just update ljatypePtrList
198 < #ifndef IS_MPI
199 <    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
200 <    if ( thisStat /= 0 ) then
201 <       status = -1
202 <       return
203 <    endif
204 <
251 > #endif
252      
206 – ! if were're in MPI, we also have to worry about row and col lists    
207 – #else
208 –  
209 – ! We can only set up forces if mpiSimulation has been setup.
210 –    if (.not. isMPISimSet()) then
211 –       write(default_error,*) "MPI is not set"
212 –       status = -1
213 –       return
214 –    endif
215 –    nrow = getNrow(plan_row)
216 –    ncol = getNcol(plan_col)
217 –    mynode = getMyNode()
218 – !! Allocate temperary arrays to hold gather information
219 –    allocate(identRow(nrow),stat=alloc_stat)
220 –    if (alloc_stat /= 0 ) then
221 –       status = -1
222 –       return
223 –    endif
224 –
225 –    allocate(identCol(ncol),stat=alloc_stat)
226 –    if (alloc_stat /= 0 ) then
227 –       status = -1
228 –       return
229 –    endif
230 –
231 – !! Gather idents into row and column idents
232 –
233 –    call gather(ident,identRow,plan_row)
234 –    call gather(ident,identCol,plan_col)
253      
236 –  
237 – !! Create row and col pointer lists
238 –  
239 –    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
240 –    if (thisStat /= 0 ) then
241 –       status = -1
242 –       return
243 –    endif
244 –  
245 –    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
246 –    if (thisStat /= 0 ) then
247 –       status = -1
248 –       return
249 –    endif
250 –
251 – !! free temporary ident arrays
252 –    if (allocated(identCol)) then
253 –       deallocate(identCol)
254 –    end if
255 –    if (allocated(identCol)) then
256 –       deallocate(identRow)
257 –    endif
258 –
259 – #endif
260 –    
261 –    call initForce_Modules(thisStat)
262 –    if (thisStat /= 0) then
263 –       status = -1
264 –       return
265 –    endif
266 –
267 – !! Create neighbor lists
268 –    call expandList(thisStat)
269 –    if (thisStat /= 0) then
270 –       status = -1
271 –       return
272 –    endif
273 –
274 –    isFFinit = .true.
275 –
276 –
277 –  end subroutine init_FF
278 –
279 –
280 –
281 –
282 –  subroutine initForce_Modules(thisStat)
283 –    integer, intent(out) :: thisStat
284 –    integer :: my_status
285 –    
286 –    thisStat = 0
287 –    call init_lj_FF(ListHead,my_status)
288 –    if (my_status /= 0) then
289 –       thisStat = -1
290 –       return
291 –    end if
292 –
293 –  end subroutine initForce_Modules
294 –
295 –
296 –
297 –
298 – !! FORCE routine Calculates Lennard Jones forces.
299 – !------------------------------------------------------------->
300 –  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
301 – !! Position array provided by C, dimensioned by getNlocal
302 –    real ( kind = dp ), dimension(3,getNlocal()) :: q
303 –  !! Rotation Matrix for each long range particle in simulation.
304 –    real( kind = dp), dimension(9,getNlocal()) :: A
305 –
306 –  !! Magnitude dipole moment
307 –    real( kind = dp ), dimension(3,getNlocal()) :: mu
308 –  !! Unit vectors for dipoles (lab frame)
309 –    real( kind = dp ), dimension(3,getNlocal()) :: u_l
310 – !! Force array provided by C, dimensioned by getNlocal
311 –    real ( kind = dp ), dimension(3,getNlocal()) :: f
312 – !! Torsion array provided by C, dimensioned by getNlocal
313 –    real( kind = dp ), dimension(3,getNlocal()) :: t
314 –
315 – !! Stress Tensor
316 –    real( kind = dp), dimension(9) :: tau
317 –
318 –    real ( kind = dp ) :: potE
319 –    logical ( kind = 2) :: do_pot
320 –    integer :: FFerror
321 –
322 –    
323 –    type(atype), pointer :: Atype_i
324 –    type(atype), pointer :: Atype_j
325 –
326 –
327 –
328 –
329 –  
330 –
254   #ifdef IS_MPI
255 <  real( kind = DP ) :: pot_local
256 <
257 < !! Local arrays needed for MPI
258 <
259 < #endif
260 <
338 <
339 <
340 <  real( kind = DP )   :: pe
341 <  logical             :: update_nlist
342 <
343 <
344 <  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
345 <  integer :: nlist
346 <  integer :: j_start
347 <
348 <  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
349 <
350 <  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
351 <  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
352 <
353 <  real( kind = DP ) :: dielectric = 0.0_dp
354 <
355 < ! a rig that need to be fixed.
356 < #ifdef IS_MPI
357 <  real( kind = dp ) :: pe_local
358 <  integer :: nlocal
255 >    !! distribute all reaction field stuff (or anything for post-pair):
256 >    call scatter(rflRow,rflTemp1,plan_row3d)
257 >    call scatter(rflCol,rflTemp2,plan_col3d)
258 >    do i = 1,nlocal
259 >       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
260 >    end do
261   #endif
360 –  integer :: nrow
361 –  integer :: ncol
362 –  integer :: natoms
363 –  integer :: neighborListSize
364 –  integer :: listerror
365 – !! should we calculate the stress tensor
366 –  logical  :: do_stress = .false.
367 –
368 –
369 –  FFerror = 0
370 –
371 – ! Make sure we are properly initialized.
372 –  if (.not. isFFInit) then
373 –     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
374 –     FFerror = -1
375 –     return
376 –  endif
377 – #ifdef IS_MPI
378 –    if (.not. isMPISimSet()) then
379 –     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
380 –     FFerror = -1
381 –     return
382 –  endif
383 – #endif
384 –
385 – !! initialize local variables  
386 –  natoms = getNlocal()
387 –  call getRcut(rcut,rcut2=rcutsq)
388 –  call getRlist(rlist,rlistsq)
389 –
390 – !! Find ensemble
391 –  if (isEnsemble("NPT")) do_stress = .true.
392 – !! set to wrap
393 –  if (isPBC()) wrap = .true.
394 –
395 –
396 – #ifndef IS_MPI
397 –  nrow = natoms - 1
398 –  ncol = natoms
399 – #else
400 –  nrow = getNrow(plan_row)
401 –  ncol = getNcol(plan_col)
402 –  nlocal = natoms
403 –  j_start = 1
404 – #endif
405 –
406 –  
407 – !! See if we need to update neighbor lists
408 –  call check(q,update_nlist)
409 –
410 – !--------------WARNING...........................
411 – ! Zero variables, NOTE:::: Forces are zeroed in C
412 – ! Zeroing them here could delete previously computed
413 – ! Forces.
414 – !------------------------------------------------
415 –  call zero_module_variables()
416 –
417 –
418 – ! communicate MPI positions
419 – #ifdef IS_MPI    
420 –    call gather(q,qRow,plan_row3d)
421 –    call gather(q,qCol,plan_col3d)
422 –
423 –    call gather(mu,muRow,plan_row3d)
424 –    call gather(mu,muCol,plan_col3d)
425 –
426 –    call gather(u_l,u_lRow,plan_row3d)
427 –    call gather(u_l,u_lCol,plan_col3d)
428 –
429 –    call gather(A,ARow,plan_row_rotation)
430 –    call gather(A,ACol,plan_col_rotation)
431 – #endif
432 –
433 –
434 –  if (update_nlist) then
435 –
436 –     ! save current configuration, contruct neighbor list,
437 –     ! and calculate forces
438 –     call save_neighborList(q)
439 –    
440 –     neighborListSize = getNeighborListSize()
441 –     nlist = 0
442 –    
262      
263 <
445 <     do i = 1, nrow
446 <        point(i) = nlist + 1
263 > ! This is the post-pair loop:
264   #ifdef IS_MPI
265 <        Atype_i => identPtrListRow(i)%this
266 <        tag_i = tagRow(i)
265 >    
266 >    if (system_has_postpair_atoms) then
267 >       do i = 1, nlocal
268 >          Atype_i => identPtrListRow(i)%this
269 >          call do_postpair(i, Atype_i)
270 >       enddo
271 >    endif
272 >    
273   #else
274 <        Atype_i   => identPtrList(i)%this
275 <        j_start = i + 1
274 >    
275 >    if (system_has_postpair_atoms) then
276 >       do i = 1, natoms
277 >          Atype_i => identPtr(i)%this
278 >          call do_postpair(i, Atype_i)
279 >       enddo
280 >    endif
281 >    
282   #endif
454 –
455 –        inner: do j = j_start, ncol
456 – ! Assign identity pointers and tags
457 – #ifdef IS_MPI
458 –           Atype_j => identPtrListColumn(j)%this
459 –          
460 –           call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
461 –                rxij,ryij,rzij,rijsq,r)
462 – !! For mpi, use newtons 3rd law when building neigbor list
463 – !! Also check to see the particle i != j.
464 –           if (mpi_cycle_jLoop(i,j)) cycle inner:
465 –
466 – #else          
467 –           Atype_j   => identPtrList(j)%this
468 –           call get_interatomic_vector(i,j,q(:,i),q(:,j),&
469 –                rxij,ryij,rzij,rijsq,r)
470 –      
471 – #endif          
472 –          
473 –           if (rijsq <  rlistsq) then
474 –
475 –              nlist = nlist + 1
476 –
477 –              if (nlist > neighborListSize) then
478 –                 call expandList(listerror)
479 –                 if (listerror /= 0) then
480 –                    FFerror = -1
481 –                    write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
482 –                    return
483 –                 end if
484 –              endif
485 –
486 –              list(nlist) = j
487 –
283      
489 –              if (rijsq <  rcutsq) then
490 –                 call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
491 –              endif
492 –          enddo inner
493 –     enddo
284  
285   #ifdef IS_MPI
496 –     point(nrow + 1) = nlist + 1
497 – #else
498 –     point(natoms) = nlist + 1
499 – #endif
500 –
501 –  else !! (update)
502 –
503 –     ! use the list to find the neighbors
504 –     do i = 1, nrow
505 –        JBEG = POINT(i)
506 –        JEND = POINT(i+1) - 1
507 –        ! check thiat molecule i has neighbors
508 –        if (jbeg .le. jend) then
509 –
510 – #ifdef IS_MPI
511 –           ljAtype_i => identPtrListRow(i)%this
512 – #else
513 –           ljAtype_i => identPtrList(i)%this
514 – #endif
515 –           do jnab = jbeg, jend
516 –              j = list(jnab)
517 – #ifdef IS_MPI
518 –              ljAtype_j = identPtrListColumn(j)%this
519 –              call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
520 –                   rxij,ryij,rzij,rijsq,r)
521 –              
522 – #else
523 –              ljAtype_j = identPtrList(j)%this
524 –              call get_interatomic_vector(i,j,q(:,i),q(:,j),&
525 –                   rxij,ryij,rzij,rijsq,r)
526 – #endif
527 –              call do_pair(i,j,r,rxij,ryij,rzij)
528 –          enddo
529 –       endif
530 –    enddo
531 – endif
532 –
533 –
534 –
535 – #ifdef IS_MPI
286      !!distribute forces
287  
288 <    call scatter(fRow,f,plan_row3d)
289 <    call scatter(fCol,fTemp,plan_col3d)
540 <
288 >    call scatter(f_Row,f,plan_row3d)
289 >    call scatter(f_Col,f_temp,plan_col3d)
290      do i = 1,nlocal
291 <       f(1:3,i) = f(1:3,i) + fTemp(1:3,i)
291 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
292      end do
293  
294 <    if (do_torque) then
295 <       call scatter(tRow,t,plan_row3d)
296 <       call scatter(tCol,tTemp,plan_col3d)
294 >    if (doTorque()) then
295 >       call scatter(t_Row,t,plan_row3d)
296 >       call scatter(t_Col,t_temp,plan_col3d)
297      
298         do i = 1,nlocal
299 <          t(1:3,i) = t(1:3,i) + tTemp(1:3,i)
299 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
300         end do
301      endif
302 <
302 >    
303      if (do_pot) then
304         ! scatter/gather pot_row into the members of my column
305 <       call scatter(eRow,eTemp,plan_row)
305 >       call scatter(pot_Row, pot_Temp, plan_row)
306        
307         ! scatter/gather pot_local into all other procs
308         ! add resultant to get total pot
309         do i = 1, nlocal
310 <          pe_local = pe_local + eTemp(i)
310 >          pot_local = pot_local + pot_Temp(i)
311         enddo
312  
313 <       eTemp = 0.0E0_DP
314 <       call scatter(eCol,eTemp,plan_col)
313 >       pot_Temp = 0.0_DP
314 >
315 >       call scatter(pot_Col, pot_Temp, plan_col)
316         do i = 1, nlocal
317 <          pe_local = pe_local + eTemp(i)
317 >          pot_local = pot_local + pot_Temp(i)
318         enddo
319        
320 <       pe = pe_local
320 >       pot = pot_local
321      endif
572 – #else
573 – ! Copy local array into return array for c
574 –    f = fTemp
575 –    t = tTemp
576 – #endif
322  
323 <    potE = pe
323 >    if (doStress()) then
324 >       mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
325 >            mpi_comm_world,mpi_err)
326 >       mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
327 >            mpi_comm_world,mpi_err)
328 >    endif
329  
330 + #endif
331  
332 <    if (do_stress) then
333 < #ifdef IS_MPI
334 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
584 <            mpi_comm_world,mpi_err)
585 < #else
586 <       tau = tauTemp
587 < #endif      
332 >    if (doStress()) then
333 >       tau = tau_Temp
334 >       virial = virial_Temp
335      endif
336  
337    end subroutine do_force_loop
338  
339  
593 –
594 –
595 –
596 –
597 –
598 –
599 –
600 –
340   !! Calculate any pre-force loop components and update nlist if necessary.
341    subroutine do_preForce(updateNlist)
342      logical, intent(inout) :: updateNlist
# Line 606 | Line 345 | contains
345  
346    end subroutine do_preForce
347  
609 –
610 –
611 –
612 –
613 –
614 –
615 –
616 –
617 –
618 –
619 –
620 –
348   !! Calculate any post force loop components, i.e. reaction field, etc.
349    subroutine do_postForce()
350  
# Line 625 | Line 352 | contains
352  
353    end subroutine do_postForce
354  
355 +  subroutine do_pair(i, j, rijsq, d)
356  
357 +    integer, intent(in) :: i, j
358 +    real ( kind = dp ), intent(in)    :: rijsq
359 +    real ( kind = dp )                :: r
360 +    real ( kind = dp ), intent(inout) :: d(3)
361  
362 +    r = sqrt(rijsq)
363 +    
364 +    logical :: is_LJ_i, is_LJ_j
365 +    logical :: is_DP_i, is_DP_j
366 +    logical :: is_Sticky_i, is_Sticky_j
367 +    integer :: me_i, me_j
368  
369 + #ifdef IS_MPI
370  
371 +    me_i = atid_row(i)
372 +    me_j = atid_col(j)
373  
374 + #else
375  
376 +    me_i = atid(i)
377 +    me_j = atid(j)
378  
379 + #endif
380  
381 +    call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
382 +    call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
383  
384 +    if ( is_LJ_i .and. is_LJ_j ) &
385 +         call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
386  
387 +    call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
388 +    call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
389  
390 +    if ( is_DP_i .and. is_DP_j ) then
391  
392 +       call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
393  
394 +       if (do_reaction_field) then
395 +          call accumulate_rf(i, j, r)
396 +       endif
397  
398 +    endif
399  
400 +    call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
401 +    call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
402  
403 <  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
404 <    type (atype ), pointer, intent(inout) :: atype_i
405 <    type (atype ), pointer, intent(inout) :: atype_j
647 <    integer :: i
648 <    integer :: j
649 <    real ( kind = dp ), intent(inout) :: rx_ij
650 <    real ( kind = dp ), intent(inout) :: ry_ij
651 <    real ( kind = dp ), intent(inout) :: rz_ij
403 >    if ( is_Sticky_i .and. is_Sticky_j ) then
404 >       call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, do_pot, do_stress)
405 >    endif
406  
407 <
654 <    real( kind = dp ) :: fx = 0.0_dp
655 <    real( kind = dp ) :: fy = 0.0_dp
656 <    real( kind = dp ) :: fz = 0.0_dp  
657 <
658 <    real( kind = dp ) ::  drdx = 0.0_dp
659 <    real( kind = dp ) ::  drdy = 0.0_dp
660 <    real( kind = dp ) ::  drdz = 0.0_dp
661 <    
662 <
663 <
664 <
665 <
666 <
667 <    call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j)
668 <      
669 < #ifdef IS_MPI
670 <                eRow(i) = eRow(i) + pot*0.5
671 <                eCol(i) = eCol(i) + pot*0.5
672 < #else
673 <                    pe = pe + pot
674 < #endif                
675 <            
676 <                drdx = -rxij / r
677 <                drdy = -ryij / r
678 <                drdz = -rzij / r
679 <                
680 <                fx = dudr * drdx
681 <                fy = dudr * drdy
682 <                fz = dudr * drdz
683 <
684 <
685 <
686 <
687 <
688 <
689 <                
690 < #ifdef IS_MPI
691 <                fCol(1,j) = fCol(1,j) - fx
692 <                fCol(2,j) = fCol(2,j) - fy
693 <                fCol(3,j) = fCol(3,j) - fz
694 <                
695 <                fRow(1,j) = fRow(1,j) + fx
696 <                fRow(2,j) = fRow(2,j) + fy
697 <                fRow(3,j) = fRow(3,j) + fz
698 < #else
699 <                fTemp(1,j) = fTemp(1,j) - fx
700 <                fTemp(2,j) = fTemp(2,j) - fy
701 <                fTemp(3,j) = fTemp(3,j) - fz
702 <                fTemp(1,i) = fTemp(1,i) + fx
703 <                fTemp(2,i) = fTemp(2,i) + fy
704 <                fTemp(3,i) = fTemp(3,i) + fz
705 < #endif
706 <                
707 <                if (do_stress) then
708 <                   tauTemp(1) = tauTemp(1) + fx * rxij
709 <                   tauTemp(2) = tauTemp(2) + fx * ryij
710 <                   tauTemp(3) = tauTemp(3) + fx * rzij
711 <                   tauTemp(4) = tauTemp(4) + fy * rxij
712 <                   tauTemp(5) = tauTemp(5) + fy * ryij
713 <                   tauTemp(6) = tauTemp(6) + fy * rzij
714 <                   tauTemp(7) = tauTemp(7) + fz * rxij
715 <                   tauTemp(8) = tauTemp(8) + fz * ryij
716 <                   tauTemp(9) = tauTemp(9) + fz * rzij
717 <                endif
718 <
719 <
720 <
407 >      
408    end subroutine do_pair
409  
410  
411 <
412 <
726 <
727 <
728 <
729 <
730 <
731 <
732 <
733 <
734 <
735 <
736 <
737 <
738 <  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
739 < !---------------- Arguments-------------------------------
740 <   !! index i
741 <
742 <    !! Position array
411 >  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
412 >    
413      real (kind = dp), dimension(3) :: q_i
414      real (kind = dp), dimension(3) :: q_j
745 –    !! x component of vector between i and j
746 –    real ( kind = dp ), intent(out)  :: rx_ij
747 –    !! y component of vector between i and j
748 –    real ( kind = dp ), intent(out)  :: ry_ij
749 –    !! z component of vector between i and j
750 –    real ( kind = dp ), intent(out)  :: rz_ij
751 –    !! magnitude of r squared
415      real ( kind = dp ), intent(out) :: r_sq
753 –    !! magnitude of vector r between atoms i and j.
754 –    real ( kind = dp ), intent(out) :: r_ij
755 –    !! wrap into periodic box.
756 –    logical, intent(in) :: wrap
757 –
758 – !--------------- Local Variables---------------------------
759 –    !! Distance between i and j
416      real( kind = dp ) :: d(3)
761 – !---------------- END DECLARATIONS-------------------------
417  
763 –
764 – ! Find distance between i and j
418      d(1:3) = q_i(1:3) - q_j(1:3)
419 <
420 < ! Wrap back into periodic box if necessary
421 <    if ( wrap ) then
419 >    
420 >    ! Wrap back into periodic box if necessary
421 >    if ( isPBC() ) then
422         d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
423              int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
424 <    end if
424 >    endif
425      
773 – !   Find Magnitude of the vector
426      r_sq = dot_product(d,d)
427 <    r_ij = sqrt(r_sq)
776 <
777 < !   Set each component for force calculation
778 <    rx_ij = d(1)
779 <    ry_ij = d(2)
780 <    rz_ij = d(3)
781 <
782 <
427 >        
428    end subroutine get_interatomic_vector
429 +  
430 +  subroutine zero_work_arrays()
431 +    
432 + #ifdef IS_MPI
433  
434 <  subroutine zero_module_variables()
435 <
787 < #ifndef IS_MPI
788 <
789 <    pe = 0.0E0_DP
790 <    tauTemp = 0.0_dp
791 <    fTemp = 0.0_dp
792 <    tTemp = 0.0_dp
793 < #else
794 <    qRow = 0.0_dp
795 <    qCol = 0.0_dp
434 >    q_Row = 0.0_dp
435 >    q_Col = 0.0_dp  
436      
437 <    muRow = 0.0_dp
438 <    muCol = 0.0_dp
437 >    u_l_Row = 0.0_dp
438 >    u_l_Col = 0.0_dp
439      
440 <    u_lRow = 0.0_dp
441 <    u_lCol = 0.0_dp
440 >    A_Row = 0.0_dp
441 >    A_Col = 0.0_dp
442      
443 <    ARow = 0.0_dp
444 <    ACol = 0.0_dp
445 <    
446 <    fRow = 0.0_dp
447 <    fCol = 0.0_dp
448 <    
449 <  
810 <    tRow = 0.0_dp
811 <    tCol = 0.0_dp
443 >    f_Row = 0.0_dp
444 >    f_Col = 0.0_dp
445 >    f_Temp = 0.0_dp
446 >      
447 >    t_Row = 0.0_dp
448 >    t_Col = 0.0_dp
449 >    t_Temp = 0.0_dp
450  
451 <  
451 >    pot_Row = 0.0_dp
452 >    pot_Col = 0.0_dp
453 >    pot_Temp = 0.0_dp
454  
815 –    eRow = 0.0_dp
816 –    eCol = 0.0_dp
817 –    eTemp = 0.0_dp
455   #endif
456  
457 <  end subroutine zero_module_variables
457 >    tau_Temp = 0.0_dp
458 >    virial_Temp = 0.0_dp
459 >    
460 >  end subroutine zero_work_arrays
461 >  
462  
822 – #ifdef IS_MPI
463   !! Function to properly build neighbor lists in MPI using newtons 3rd law.
464   !! We don't want 2 processors doing the same i j pair twice.
465   !! Also checks to see if i and j are the same particle.
466 <  function mpi_cycle_jLoop(i,j) result(do_cycle)
466 >  function checkExcludes(atom1,atom2) result(do_cycle)
467   !--------------- Arguments--------------------------
468   ! Index i
469 <    integer,intent(in) :: i
469 >    integer,intent(in) :: atom1
470   ! Index j
471 <    integer,intent(in) :: j
471 >    integer,intent(in), optional :: atom2
472   ! Result do_cycle
473      logical :: do_cycle
474   !--------------- Local variables--------------------
475      integer :: tag_i
476      integer :: tag_j
477 < !--------------- END DECLARATIONS------------------    
478 <    tag_i = tagRow(i)
477 >    integer :: i
478 > !--------------- END DECLARATIONS------------------  
479 >    do_cycle = .false.
480 >
481 > #ifdef IS_MPI
482 >    tag_i = tagRow(atom1)
483 > #else
484 >    tag_i = tag(atom1)
485 > #endif
486 >
487 > !! Check global excludes first
488 >    if (.not. present(atom2)) then
489 >       do i = 1,nGlobalExcludes
490 >          if (excludeGlobal(i) == tag_i) then
491 >             do_cycle = .true.
492 >             return
493 >          end if
494 >       end do
495 >       return !! return after checking globals
496 >    end if
497 >
498 > !! we return if j not present here.
499      tag_j = tagColumn(j)
500  
501 <    do_cycle = .false.
501 >
502  
503      if (tag_i == tag_j) then
504         do_cycle = .true.
# Line 851 | Line 511 | contains
511      else                
512         if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
513      endif
854 –  end function mpi_cycle_jLoop
855 – #endif
514  
515 +
516 +
517 +    do i = 1, nLocalExcludes
518 +       if (tag_i = excludes(1,i) .and. excludes(2,i) < 0) then
519 +          do_cycle = .true.
520 +          return
521 +       end if
522 +    end do
523 +      
524 +
525 +  end function checkExcludes
526 +
527 +
528   end module do_Forces

Diff Legend

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