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 317 by gezelter, Tue Mar 11 23:13:06 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.11 2003-03-11 23:13:06 gezelter Exp $, $Date: 2003-03-11 23:13:06 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
8  
9  
10  
11   module do_Forces
12    use simulation
13    use definitions
14 <  use generic_atypes
14 >  use forceGlobals
15 >  use atype_typedefs
16    use neighborLists
17 +
18    
19 <  use lj_FF
19 >  use lj
20    use sticky_FF
21 <  use dp_FF
21 >  use dipole_dipole
22    use gb_FF
23  
24   #ifdef IS_MPI
# Line 22 | Line 27 | module do_Forces
27    implicit none
28    PRIVATE
29  
30 < !! Number of lj_atypes in lj_atype_list
26 <  integer, save :: n_atypes = 0
30 > public :: do_force_loop
31  
32 < !! Global list of lj atypes in simulation
29 <  type (atype), pointer :: ListHead => null()
30 <  type (atype), pointer :: ListTail => null()
32 > contains
33  
34 +  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
35 +  !------------------------------------------------------------->
36 +  subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
37 +    !! Position array provided by C, dimensioned by getNlocal
38 +    real ( kind = dp ), dimension(3,getNlocal()) :: q
39 +    !! Rotation Matrix for each long range particle in simulation.
40 +    real( kind = dp), dimension(9,getNlocal()) :: A
41 +    
42 +    !! Magnitude dipole moment
43 +    real( kind = dp ), dimension(3,getNlocal()) :: mu
44 +    !! Unit vectors for dipoles (lab frame)
45 +    real( kind = dp ), dimension(3,getNlocal()) :: u_l
46 +    !! Force array provided by C, dimensioned by getNlocal
47 +    real ( kind = dp ), dimension(3,getNlocal()) :: f
48 +    !! Torsion array provided by C, dimensioned by getNlocal
49 +    real( kind = dp ), dimension(3,getNlocal()) :: t
50 +    
51 +    !! Stress Tensor
52 +    real( kind = dp), dimension(9) :: tau
53 +    
54 +    real ( kind = dp ) :: potE
55 +    logical ( kind = 2) :: do_pot
56 +    integer :: FFerror
57  
58 + #ifdef IS_MPI
59 +    real( kind = DP ) :: pot_local
60 + #endif
61 +    
62 +    real( kind = DP )   :: pe
63 +    logical             :: update_nlist
64 +    
65  
66 +    integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
67 +    integer :: nlist
68 +    integer :: j_start
69 +    
70 +    real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
71 +    
72 +    real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
73 +    real( kind = DP ) ::  rlistsq, rcutsq, rlist, rcut
74  
75 <  logical, save :: firstTime = .True.
36 <
37 < !! Atype identity pointer lists
75 >    ! a rig that need to be fixed.
76   #ifdef IS_MPI
77 < !! Row lj_atype pointer list
78 <  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()
77 >    real( kind = dp ) :: pe_local
78 >    integer :: nlocal
79   #endif
80 +    integer :: nrow
81 +    integer :: ncol
82 +    integer :: natoms
83 +    integer :: neighborListSize
84 +    integer :: listerror
85 +    !! should we calculate the stress tensor
86 +    logical  :: do_stress = .false.
87 +    FFerror = 0
88  
89 <
90 < !! Logical has lj force field module been initialized?
91 <  logical, save :: isFFinit = .false.
92 <
93 < !! Use periodic boundry conditions
94 <  logical :: wrap = .false.
53 <
54 < !! Potential energy global module variables
89 >    ! Make sure we are properly initialized.
90 >    if (.not. isFFInit) then
91 >       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
92 >       FFerror = -1
93 >       return
94 >    endif
95   #ifdef IS_MPI
96 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
97 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
96 >    if (.not. isMPISimSet()) then
97 >       write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
98 >       FFerror = -1
99 >       return
100 >    endif
101 > #endif
102 >    
103 >    !! initialize local variables  
104 >    natoms = getNlocal()
105 >    call getRcut(rcut,rcut2=rcutsq)
106 >    call getRlist(rlist,rlistsq)
107 >    
108 >    !! See if we need to update neighbor lists
109 >    call check(q, update_nlist)
110 >    
111 >    !--------------WARNING...........................
112 >    ! Zero variables, NOTE:::: Forces are zeroed in C
113 >    ! Zeroing them here could delete previously computed
114 >    ! Forces.
115 >    !------------------------------------------------
116 >    call zero_module_variables()
117  
118 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
119 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
120 <
121 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
122 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
123 <
124 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: ARow = 0.0_dp
125 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp
126 <
127 <  
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
118 >    ! communicate MPI positions
119 > #ifdef IS_MPI    
120 >    call gather(q,q_Row,plan_row3d)
121 >    call gather(q,q_Col,plan_col3d)
122 >    
123 >    call gather(u_l,u_l_Row,plan_row3d)
124 >    call gather(u_l,u_l_Col,plan_col3d)
125 >    
126 >    call gather(A,A_Row,plan_row_rotation)
127 >    call gather(A,A_Col,plan_col_rotation)
128   #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
129  
89 –  logical :: do_preForce  = .false.
90 –  logical :: do_postForce = .false.
130  
131 + #ifdef IS_MPI
132 +    
133 +    if (update_nlist) then
134 +      
135 +       ! save current configuration, contruct neighbor list,
136 +       ! and calculate forces
137 +       call save_neighborList(q)
138 +      
139 +       neighborListSize = getNeighborListSize()
140 +       nlist = 0
141 +      
142 +       nrow = getNrow(plan_row)
143 +       ncol = getNcol(plan_col)
144 +       nlocal = getNlocal()
145 +      
146 +       do i = 1, nrow
147 +          point(i) = nlist + 1
148 +          
149 +          inner: do j = 1, ncol
150 +            
151 +             if (check_exclude(i,j)) cycle inner:
152  
153 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
154 +            
155 +             if (rijsq <  rlistsq) then            
156 +                
157 +                nlist = nlist + 1
158 +                
159 +                if (nlist > neighborListSize) then
160 +                   call expandList(listerror)
161 +                   if (listerror /= 0) then
162 +                      FFerror = -1
163 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
164 +                      return
165 +                   end if
166 +                endif
167 +                
168 +                list(nlist) = j
169 +                                
170 +                if (rijsq <  rcutsq) then
171 +                   call do_pair(i, j, rijsq, d)
172 +                endif
173 +             endif
174 +          enddo inner
175 +       enddo
176  
177 < !! Public methods and data
178 <  public :: new_atype
179 <  public :: do_forceLoop
97 <  public :: init_FF
177 >       point(nrow + 1) = nlist + 1
178 >      
179 >    else !! (update)
180  
181 <  
181 >       ! use the list to find the neighbors
182 >       do i = 1, nrow
183 >          JBEG = POINT(i)
184 >          JEND = POINT(i+1) - 1
185 >          ! check thiat molecule i has neighbors
186 >          if (jbeg .le. jend) then
187 >            
188 >             do jnab = jbeg, jend
189 >                j = list(jnab)
190  
191 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
192 +                call do_pair(i, j, rijsq, d)
193  
194 < contains
194 >             enddo
195 >          endif
196 >       enddo
197 >    endif
198  
199 < !! Adds a new lj_atype to the list.
200 <  subroutine new_atype(ident,mass,epsilon,sigma, &
201 <       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
202 <    real( kind = dp ), intent(in) :: mass
203 <    real( kind = dp ), intent(in) :: epsilon
204 <    real( kind = dp ), intent(in) :: sigma
205 <    real( kind = dp ), intent(in) :: w0
206 <    real( kind = dp ), intent(in) :: v0
207 <    real( kind = dp ), intent(in) :: dipoleMoment
199 > #else
200 >    
201 >    if (update_nlist) then
202 >      
203 >       ! save current configuration, contruct neighbor list,
204 >       ! and calculate forces
205 >       call save_neighborList(q)
206 >      
207 >       neighborListSize = getNeighborListSize()
208 >       nlist = 0
209 >          
210 >       do i = 1, natoms-1
211 >          point(i) = nlist + 1
212  
213 <    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
213 >          inner: do j = i+1, natoms
214  
215 +             if (check_exclude(i,j)) cycle inner:
216  
217 <    type (atype), pointer :: the_new_atype
218 <    integer :: alloc_error
219 <    integer :: atype_counter = 0
220 <    integer :: alloc_size
221 <    integer :: err_stat
222 <    status = 0
223 <
217 >             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
218 >          
219 >             if (rijsq <  rlistsq) then
220 >                
221 >                nlist = nlist + 1
222 >                
223 >                if (nlist > neighborListSize) then
224 >                   call expandList(listerror)
225 >                   if (listerror /= 0) then
226 >                      FFerror = -1
227 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
228 >                      return
229 >                   end if
230 >                endif
231 >                
232 >                list(nlist) = j
233 >                    
234 >                if (rijsq <  rcutsq) then
235 >                   call do_pair(i, j, rijsq, d)
236 >                endif
237 >             endif
238 >          enddo inner
239 >       enddo
240 >      
241 >       point(natoms) = nlist + 1
242 >      
243 >    else !! (update)
244 >      
245 >       ! use the list to find the neighbors
246 >       do i = 1, nrow
247 >          JBEG = POINT(i)
248 >          JEND = POINT(i+1) - 1
249 >          ! check thiat molecule i has neighbors
250 >          if (jbeg .le. jend) then
251 >            
252 >             do jnab = jbeg, jend
253 >                j = list(jnab)
254  
255 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
256 +                call do_pair(i, j, rijsq, d)
257  
258 < ! allocate a new atype    
259 <    allocate(the_new_atype,stat=alloc_error)
260 <    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
258 >             enddo
259 >          endif
260 >       enddo
261      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 –
262      
263 <
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 <
263 > #endif
264      
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)
265      
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 –
266   #ifdef IS_MPI
267 <  real( kind = DP ) :: pot_local
268 <
269 < !! Local arrays needed for MPI
270 <
271 < #endif
272 <
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
267 >    !! distribute all reaction field stuff (or anything for post-pair):
268 >    call scatter(rflRow,rflTemp1,plan_row3d)
269 >    call scatter(rflCol,rflTemp2,plan_col3d)
270 >    do i = 1,nlocal
271 >       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
272 >    end do
273   #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 –    
274      
275 <
445 <     do i = 1, nrow
446 <        point(i) = nlist + 1
275 > ! This is the post-pair loop:
276   #ifdef IS_MPI
277 <        Atype_i => identPtrListRow(i)%this
278 <        tag_i = tagRow(i)
277 >    
278 >    if (system_has_postpair_atoms) then
279 >       do i = 1, nlocal
280 >          Atype_i => identPtrListRow(i)%this
281 >          call do_postpair(i, Atype_i)
282 >       enddo
283 >    endif
284 >    
285   #else
286 <        Atype_i   => identPtrList(i)%this
287 <        j_start = i + 1
286 >    
287 >    if (system_has_postpair_atoms) then
288 >       do i = 1, natoms
289 >          Atype_i => identPtr(i)%this
290 >          call do_postpair(i, Atype_i)
291 >       enddo
292 >    endif
293 >    
294   #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 –
295      
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
296  
297   #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
298      !!distribute forces
299  
300 <    call scatter(fRow,f,plan_row3d)
301 <    call scatter(fCol,fTemp,plan_col3d)
540 <
300 >    call scatter(f_Row,f,plan_row3d)
301 >    call scatter(f_Col,f_temp,plan_col3d)
302      do i = 1,nlocal
303 <       f(1:3,i) = f(1:3,i) + fTemp(1:3,i)
303 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
304      end do
305  
306 <    if (do_torque) then
307 <       call scatter(tRow,t,plan_row3d)
308 <       call scatter(tCol,tTemp,plan_col3d)
306 >    if (doTorque()) then
307 >       call scatter(t_Row,t,plan_row3d)
308 >       call scatter(t_Col,t_temp,plan_col3d)
309      
310         do i = 1,nlocal
311 <          t(1:3,i) = t(1:3,i) + tTemp(1:3,i)
311 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
312         end do
313      endif
314 <
314 >    
315      if (do_pot) then
316         ! scatter/gather pot_row into the members of my column
317 <       call scatter(eRow,eTemp,plan_row)
317 >       call scatter(pot_Row, pot_Temp, plan_row)
318        
319         ! scatter/gather pot_local into all other procs
320         ! add resultant to get total pot
321         do i = 1, nlocal
322 <          pe_local = pe_local + eTemp(i)
322 >          pot_local = pot_local + pot_Temp(i)
323         enddo
324  
325 <       eTemp = 0.0E0_DP
326 <       call scatter(eCol,eTemp,plan_col)
325 >       pot_Temp = 0.0_DP
326 >
327 >       call scatter(pot_Col, pot_Temp, plan_col)
328         do i = 1, nlocal
329 <          pe_local = pe_local + eTemp(i)
329 >          pot_local = pot_local + pot_Temp(i)
330         enddo
331        
332 <       pe = pe_local
332 >       pot = pot_local
333      endif
572 – #else
573 – ! Copy local array into return array for c
574 –    f = fTemp
575 –    t = tTemp
576 – #endif
334  
335 <    potE = pe
335 >    if (doStress()) then
336 >       mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
337 >            mpi_comm_world,mpi_err)
338 >       mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
339 >            mpi_comm_world,mpi_err)
340 >    endif
341  
342 + #endif
343  
344 <    if (do_stress) then
345 < #ifdef IS_MPI
346 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
584 <            mpi_comm_world,mpi_err)
585 < #else
586 <       tau = tauTemp
587 < #endif      
344 >    if (doStress()) then
345 >       tau = tau_Temp
346 >       virial = virial_Temp
347      endif
348  
349    end subroutine do_force_loop
350  
351  
593 –
594 –
595 –
596 –
597 –
598 –
599 –
600 –
352   !! Calculate any pre-force loop components and update nlist if necessary.
353    subroutine do_preForce(updateNlist)
354      logical, intent(inout) :: updateNlist
# Line 606 | Line 357 | contains
357  
358    end subroutine do_preForce
359  
609 –
610 –
611 –
612 –
613 –
614 –
615 –
616 –
617 –
618 –
619 –
620 –
360   !! Calculate any post force loop components, i.e. reaction field, etc.
361    subroutine do_postForce()
362  
# Line 625 | Line 364 | contains
364  
365    end subroutine do_postForce
366  
367 +  subroutine do_pair(i, j, rijsq, d)
368  
369 +    integer, intent(in) :: i, j
370 +    real ( kind = dp ), intent(in)    :: rijsq
371 +    real ( kind = dp )                :: r
372 +    real ( kind = dp ), intent(inout) :: d(3)
373  
374 <
631 <
632 <
633 <
634 <
635 <
636 <
637 <
638 <
639 <
640 <
641 <
642 <
643 <
644 <  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
645 <    type (atype ), pointer, intent(inout) :: atype_i
646 <    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
652 <
653 <
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
374 >    r = sqrt(rijsq)
375      
376 +    logical :: is_LJ_i, is_LJ_j
377 +    logical :: is_DP_i, is_DP_j
378 +    logical :: is_Sticky_i, is_Sticky_j
379 +    integer :: me_i, me_j
380  
381 + #ifdef IS_MPI
382  
383 +    me_i = atid_row(i)
384 +    me_j = atid_col(j)
385  
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
386   #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
387  
388 +    me_i = atid(i)
389 +    me_j = atid(j)
390  
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
391   #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
392  
393 +    call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
394 +    call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
395  
396 +    if ( is_LJ_i .and. is_LJ_j ) call do_lj_pair(i, j, d, r, pot, f)
397  
398 <  end subroutine do_pair
398 >    call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
399 >    call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
400  
401 +    if ( is_DP_i .and. is_DP_j ) then
402  
403 +       call do_dipole_pair(i, j, d, r, pot, u_l, f, t)
404  
405 +       if (do_reaction_field) then
406 +          call accumulate_rf(i, j, r_ij)
407 +       endif
408  
409 +    endif
410  
411 +    call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
412 +    call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
413  
414 +    if ( is_Sticky_i .and. is_Sticky_j ) then
415 +       call do_sticky_pair(i, j, d, r, pot, u_l, f, t)
416 +    endif
417  
418 +      
419 +  end subroutine do_pair
420  
421  
422 <
423 <
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
422 >  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
423 >    
424      real (kind = dp), dimension(3) :: q_i
425      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
426      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
427      real( kind = dp ) :: d(3)
761 – !---------------- END DECLARATIONS-------------------------
428  
763 –
764 – ! Find distance between i and j
429      d(1:3) = q_i(1:3) - q_j(1:3)
430 <
431 < ! Wrap back into periodic box if necessary
432 <    if ( wrap ) then
430 >    
431 >    ! Wrap back into periodic box if necessary
432 >    if ( isPBC() ) then
433         d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
434              int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
435 <    end if
435 >    endif
436      
773 – !   Find Magnitude of the vector
437      r_sq = dot_product(d,d)
438 <    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 <
438 >        
439    end subroutine get_interatomic_vector
440 <
440 >  
441    subroutine zero_module_variables()
442  
443   #ifndef IS_MPI
# Line 819 | Line 475 | contains
475  
476    end subroutine zero_module_variables
477  
478 < #ifdef IS_MPI
478 >
479   !! Function to properly build neighbor lists in MPI using newtons 3rd law.
480   !! We don't want 2 processors doing the same i j pair twice.
481   !! Also checks to see if i and j are the same particle.
482 <  function mpi_cycle_jLoop(i,j) result(do_cycle)
482 >  function checkExcludes(atom1,atom2) result(do_cycle)
483   !--------------- Arguments--------------------------
484   ! Index i
485 <    integer,intent(in) :: i
485 >    integer,intent(in) :: atom1
486   ! Index j
487 <    integer,intent(in) :: j
487 >    integer,intent(in), optional :: atom2
488   ! Result do_cycle
489      logical :: do_cycle
490   !--------------- Local variables--------------------
491      integer :: tag_i
492      integer :: tag_j
493 < !--------------- END DECLARATIONS------------------    
494 <    tag_i = tagRow(i)
493 >    integer :: i
494 > !--------------- END DECLARATIONS------------------  
495 >    do_cycle = .false.
496 >
497 > #ifdef IS_MPI
498 >    tag_i = tagRow(atom1)
499 > #else
500 >    tag_i = tag(atom1)
501 > #endif
502 >
503 > !! Check global excludes first
504 >    if (.not. present(atom2)) then
505 >       do i = 1,nGlobalExcludes
506 >          if (excludeGlobal(i) == tag_i) then
507 >             do_cycle = .true.
508 >             return
509 >          end if
510 >       end do
511 >       return !! return after checking globals
512 >    end if
513 >
514 > !! we return if j not present here.
515      tag_j = tagColumn(j)
516  
517 <    do_cycle = .false.
517 >
518  
519      if (tag_i == tag_j) then
520         do_cycle = .true.
# Line 851 | Line 527 | contains
527      else                
528         if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
529      endif
854 –  end function mpi_cycle_jLoop
855 – #endif
530  
531 +
532 +
533 +    do i = 1, nLocalExcludes
534 +       if (tag_i = excludes(1,i) .and. excludes(2,i) < 0) then
535 +          do_cycle = .true.
536 +          return
537 +       end if
538 +    end do
539 +      
540 +
541 +  end function checkExcludes
542 +
543 +
544   end module do_Forces

Diff Legend

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