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 325 by gezelter, Wed Mar 12 19:10:54 2003 UTC

# Line 1 | Line 1
1 + !! do_Forces.F90
2 + !! module do_Forces
3   !! Calculates Long Range forces.
4 +
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.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.13 2003-03-12 19:10:54 gezelter Exp $, $Date: 2003-03-12 19:10:54 $, $Name: not supported by cvs2svn $, $Revision: 1.13 $
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_typedefs
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 < !! LJ mixing array  
36 < !  type (lj_atype), dimension(:,:), pointer :: ljMixed => null()
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  
37 –  logical, save :: firstTime = .True.
38 –
39 – !! Atype identity pointer lists
64   #ifdef IS_MPI
65 < !! 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()
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 <
54 < !! Public methods and data
55 <  public :: new_atype
56 <  public :: do_forceLoop
57 <  public :: getLjPot
58 <  public :: init_FF
59 <
60 <  
61 <
62 <
63 < contains
64 <
65 < !! Adds a new lj_atype to the list.
66 <  subroutine new_atype(ident,mass,epsilon,sigma, &
67 <       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
68 <    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
74 <
75 <    integer, intent(in) :: ident
76 <    integer, intent(out) :: status
77 <    integer, intent(in) :: is_Sticky
78 <    integer, intent(in) :: is_DP
79 <    integer, intent(in) :: is_GB
80 <    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
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 <    end if
82 >    endif
83 > #ifdef IS_MPI
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 < ! assign our new atype information
107 <    the_new_atype%mass        = mass
108 <    the_new_atype%epsilon     = epsilon
109 <    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 <
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 < ! assume that this atype will be successfully added
112 <    the_new_atype%atype_ident = ident
113 <    the_new_atype%atype_number = n_lj_atypes + 1
111 >    call gather(u_l,u_l_Row,plan_row3d)
112 >    call gather(u_l,u_l_Col,plan_col3d)
113      
114 <    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
115 <    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
116 <    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
118 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
114 >    call gather(A,A_Row,plan_row_rotation)
115 >    call gather(A,A_Col,plan_col_rotation)
116 > #endif
117  
120 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
121 –    if (err_stat /= 0 ) then
122 –       status = -1
123 –       return
124 –    endif
118  
119 <    n_atypes = n_atypes + 1
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 <  end subroutine new_atype
165 >       point(nrow + 1) = nlist + 1
166 >      
167 >    else !! (update)
168  
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 <  subroutine init_FF(nComponents,ident, status)
180 < !! Number of components in ident array
134 <    integer, intent(inout) :: nComponents
135 < !! Array of identities nComponents long corresponding to
136 < !! ljatype ident.
137 <    integer, dimension(nComponents),intent(inout) :: ident
138 < !!  Result status, success = 0, error = -1
139 <    integer, intent(out) :: Status
179 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
180 >                call do_pair(i, j, rijsq, d)
181  
182 <    integer :: alloc_stat
182 >             enddo
183 >          endif
184 >       enddo
185 >    endif
186  
187 <    integer :: thisStat
188 <    integer :: i
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 :: myNode
147 < #ifdef IS_MPI
148 <    integer, allocatable, dimension(:) :: identRow
149 <    integer, allocatable, dimension(:) :: identCol
150 <    integer :: nrow
151 <    integer :: ncol
152 < #endif
153 <    status = 0
154 <  
201 >          inner: do j = i+1, natoms
202  
203 <    
203 >             if (check_exclude(i,j)) cycle inner:
204  
205 < !! if were're not in MPI, we just update ljatypePtrList
206 < #ifndef IS_MPI
207 <    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
208 <    if ( thisStat /= 0 ) then
209 <       status = -1
210 <       return
211 <    endif
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 +             enddo
247 +          endif
248 +       enddo
249 +    endif
250      
251 < ! if were're in MPI, we also have to worry about row and col lists    
251 > #endif
252 >    
253 >    
254 > #ifdef IS_MPI
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
262 >    
263 > ! This is the post-pair loop:
264 > #ifdef IS_MPI
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 <  
275 < ! We can only set up forces if mpiSimulation has been setup.
276 <    if (.not. isMPISimSet()) then
277 <       write(default_error,*) "MPI is not set"
278 <       status = -1
279 <       return
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 <    nrow = getNrow(plan_row)
282 <    ncol = getNcol(plan_col)
283 <    mynode = getMyNode()
179 < !! Allocate temperary arrays to hold gather information
180 <    allocate(identRow(nrow),stat=alloc_stat)
181 <    if (alloc_stat /= 0 ) then
182 <       status = -1
183 <       return
184 <    endif
281 >    
282 > #endif
283 >    
284  
285 <    allocate(identCol(ncol),stat=alloc_stat)
286 <    if (alloc_stat /= 0 ) then
188 <       status = -1
189 <       return
190 <    endif
285 > #ifdef IS_MPI
286 >    !!distribute forces
287  
288 < !! Gather idents into row and column idents
289 <
290 <    call gather(ident,identRow,plan_row)
291 <    call gather(ident,identCol,plan_col)
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) + f_temp(1:3,i)
292 >    end do
293 >
294 >    if (doTorque()) then
295 >       call scatter(t_Row,t,plan_row3d)
296 >       call scatter(t_Col,t_temp,plan_col3d)
297      
298 <  
299 < !! Create row and col pointer lists
300 <  
200 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
201 <    if (thisStat /= 0 ) then
202 <       status = -1
203 <       return
298 >       do i = 1,nlocal
299 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
300 >       end do
301      endif
302 <  
303 <    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
304 <    if (thisStat /= 0 ) then
305 <       status = -1
306 <       return
302 >    
303 >    if (do_pot) then
304 >       ! scatter/gather pot_row into the members of my column
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 >          pot_local = pot_local + pot_Temp(i)
311 >       enddo
312 >
313 >       pot_Temp = 0.0_DP
314 >
315 >       call scatter(pot_Col, pot_Temp, plan_col)
316 >       do i = 1, nlocal
317 >          pot_local = pot_local + pot_Temp(i)
318 >       enddo
319 >      
320 >       pot = pot_local
321      endif
322  
323 < !! free temporary ident arrays
324 <    if (allocated(identCol)) then
325 <       deallocate(identCol)
326 <    end if
327 <    if (allocated(identCol)) then
217 <       deallocate(identRow)
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
221 –    
222 –    call initForce_Modules(thisStat)
223 –    if (thisStat /= 0) then
224 –       status = -1
225 –       return
226 –    endif
331  
332 < !! Create neighbor lists
333 <    call expandList(thisStat)
334 <    if (thisStat /= 0) then
231 <       status = -1
232 <       return
332 >    if (doStress()) then
333 >       tau = tau_Temp
334 >       virial = virial_Temp
335      endif
336  
337 <    isFFinit = .true.
337 >  end subroutine do_force_loop
338  
339  
340 <  end subroutine init_FF
340 > !! Calculate any pre-force loop components and update nlist if necessary.
341 >  subroutine do_preForce(updateNlist)
342 >    logical, intent(inout) :: updateNlist
343  
344  
345  
346 +  end subroutine do_preForce
347  
348 <  subroutine initForce_Modules(thisStat)
349 <    integer, intent(out) :: thisStat
245 <    integer :: my_status
246 <    
247 <    thisStat = 0
248 <    call init_lj_FF(ListHead,my_status)
249 <    if (my_status /= 0) then
250 <       thisStat = -1
251 <       return
252 <    end if
348 > !! Calculate any post force loop components, i.e. reaction field, etc.
349 >  subroutine do_postForce()
350  
254 –  end subroutine initForce_Modules
351  
352  
353 +  end subroutine do_postForce
354  
355 +  subroutine do_pair(i, j, rijsq, d)
356  
357 < !! FORCE routine Calculates Lennard Jones forces.
358 < !------------------------------------------------------------->
359 <  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
360 < !! 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
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 <  !! Magnitude dipole moment
268 <    real( kind = dp ), dimension(3,getNlocal()) :: mu
269 <  !! Unit vectors for dipoles (lab frame)
270 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
271 < !! Force array provided by C, dimensioned by getNlocal
272 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
273 < !! Torsion array provided by C, dimensioned by getNlocal
274 <    real( kind = dp ), dimension(3,getNlocal()) :: t
275 <
276 < !! Stress Tensor
277 <    real( kind = dp), dimension(9) :: tau
278 <    real( kind = dp), dimension(9) :: tauTemp
279 <    real ( kind = dp ) :: potE
280 <    logical ( kind = 2) :: do_pot
281 <    integer :: FFerror
282 <
362 >    r = sqrt(rijsq)
363      
364 <    type(atype), pointer :: Atype_i
365 <    type(atype), pointer :: Atype_j
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  
287 –
288 –
289 –
290 –  
291 –
369   #ifdef IS_MPI
293 –  real( kind = DP ) :: pot_local
370  
371 < !! Local arrays needed for MPI
372 <  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
371 >    me_i = atid_row(i)
372 >    me_j = atid_col(j)
373  
374 <  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
374 > #else
375  
376 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
377 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
376 >    me_i = atid(i)
377 >    me_j = atid(j)
378  
379 <  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
379 > #endif
380  
381 <  
381 >    call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
382 >    call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
383  
384 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
385 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
312 <  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp = 0.0_dp
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 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
388 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
316 <  real(kind = dp), dimension(3,getNlocal()) :: tMPITemp = 0.0_dp
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 <  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
320 <  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
392 >       call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
393  
394 <  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
394 >       if (do_reaction_field) then
395 >          call accumulate_rf(i, j, r)
396 >       endif
397  
398 < #endif
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 +    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 <  real( kind = DP )   :: pe
408 <  logical             :: update_nlist
407 >      
408 >  end subroutine do_pair
409  
410  
411 <  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
412 <  integer :: nlist
413 <  integer :: j_start
414 <  integer :: tag_i,tag_j
415 <  real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
416 <  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
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
415 >    real ( kind = dp ), intent(out) :: r_sq
416 >    real( kind = dp ) :: d(3)
417  
418 <  real( kind = DP ) :: dielectric = 0.0_dp
419 <
420 < ! a rig that need to be fixed.
418 >    d(1:3) = q_i(1:3) - q_j(1:3)
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 >    endif
425 >    
426 >    r_sq = dot_product(d,d)
427 >        
428 >  end subroutine get_interatomic_vector
429 >  
430 >  subroutine zero_work_arrays()
431 >    
432   #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.
433  
434 <
435 <  FFerror = 0
436 <
437 < ! Make sure we are properly initialized.
438 <  if (.not. isFFInit) then
439 <     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
440 <     FFerror = -1
441 <     return
442 <  endif
443 < #ifdef IS_MPI
444 <    if (.not. isMPISimSet()) then
445 <     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
446 <     FFerror = -1
447 <     return
448 <  endif
449 < #endif
374 <
375 < !! initialize local variables  
376 <  natoms = getNlocal()
377 <  call getRcut(rcut,rcut2=rcutsq)
378 <  call getRlist(rlist,rlistsq)
379 < !! Find ensemble
380 <  if (isEnsemble("NPT")) do_stress = .true.
381 <
382 < #ifndef IS_MPI
383 <  nrow = natoms - 1
384 <  ncol = natoms
385 < #else
386 <  nrow = getNrow(plan_row)
387 <  ncol = getNcol(plan_col)
388 <  nlocal = natoms
389 <  j_start = 1
390 < #endif
391 <
392 <  
393 < !! See if we need to update neighbor lists
394 <  call check(q,update_nlist)
395 < !  if (firstTime) then
396 < !     update_nlist = .true.
397 < !     firstTime = .false.
398 < !  endif
399 <
400 < !--------------WARNING...........................
401 < ! Zero variables, NOTE:::: Forces are zeroed in C
402 < ! Zeroing them here could delete previously computed
403 < ! Forces.
404 < !------------------------------------------------
405 < #ifndef IS_MPI
406 < !  nloops = nloops + 1
407 <  pe = 0.0E0_DP
434 >    q_Row = 0.0_dp
435 >    q_Col = 0.0_dp  
436 >    
437 >    u_l_Row = 0.0_dp
438 >    u_l_Col = 0.0_dp
439 >    
440 >    A_Row = 0.0_dp
441 >    A_Col = 0.0_dp
442 >    
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 < #else
452 <    fRow = 0.0E0_DP
453 <    fCol = 0.0E0_DP
451 >    pot_Row = 0.0_dp
452 >    pot_Col = 0.0_dp
453 >    pot_Temp = 0.0_dp
454  
413 –    pe_local = 0.0E0_DP
414 –
415 –    eRow = 0.0E0_DP
416 –    eCol = 0.0E0_DP
417 –    eTemp = 0.0E0_DP
455   #endif
456  
457 < ! communicate MPI positions
458 < #ifdef IS_MPI    
422 <    call gather(q,qRow,plan_row3d)
423 <    call gather(q,qCol,plan_col3d)
424 <
425 <    call gather(mu,muRow,plan_row3d)
426 <    call gather(mu,muCol,plan_col3d)
427 <
428 <    call gather(u_l,u_lRow,plan_row3d)
429 <    call gather(u_l,u_lCol,plan_col3d)
430 <
431 <    call gather(A,ARow,plan_row_rotation)
432 <    call gather(A,ACol,plan_col_rotation)
433 <
434 < #endif
435 <
436 <
437 <  if (update_nlist) then
438 <
439 <     ! save current configuration, contruct neighbor list,
440 <     ! and calculate forces
441 <     call save_neighborList(q)
442 <    
443 <     neighborListSize = getNeighborListSize()
444 <     nlist = 0
445 <    
457 >    tau_Temp = 0.0_dp
458 >    virial_Temp = 0.0_dp
459      
460 +  end subroutine zero_work_arrays
461 +  
462  
463 <     do i = 1, nrow
464 <        point(i) = nlist + 1
465 < #ifdef IS_MPI
466 <        ljAtype_i => identPtrListRow(i)%this
467 <        tag_i = tagRow(i)
468 <        rxi = qRow(1,i)
469 <        ryi = qRow(2,i)
470 <        rzi = qRow(3,i)
471 < #else
472 <        ljAtype_i   => identPtrList(i)%this
473 <        j_start = i + 1
474 <        rxi = q(1,i)
475 <        ryi = q(2,i)
476 <        rzi = q(3,i)
477 < #endif
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 checkExcludes(atom1,atom2) result(do_cycle)
467 > !--------------- Arguments--------------------------
468 > ! Index i
469 >    integer,intent(in) :: atom1
470 > ! Index 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 >    integer :: i
478 > !--------------- END DECLARATIONS------------------  
479 >    do_cycle = .false.
480  
464 –        inner: do j = j_start, ncol
481   #ifdef IS_MPI
482 < ! Assign identity pointers and tags
467 <           ljAtype_j => identPtrListColumn(j)%this
468 <           tag_j = tagColumn(j)
469 <           if (newtons_thrd) then
470 <              if (tag_i <= tag_j) then
471 <                 if (mod(tag_i + tag_j,2) == 0) cycle inner
472 <              else                
473 <                 if (mod(tag_i + tag_j,2) == 1) cycle inner
474 <              endif
475 <           endif
476 <
477 <           rxij = wrap(rxi - qCol(1,j), 1)
478 <           ryij = wrap(ryi - qCol(2,j), 2)
479 <           rzij = wrap(rzi - qCol(3,j), 3)
480 < #else          
481 <           ljAtype_j   => identPtrList(j)%this
482 <           rxij = wrap(rxi - q(1,j), 1)
483 <           ryij = wrap(ryi - q(2,j), 2)
484 <           rzij = wrap(rzi - q(3,j), 3)
485 <      
486 < #endif          
487 <           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
488 <
489 < #ifdef IS_MPI
490 <             if (rijsq <=  rlistsq .AND. &
491 <                  tag_j /= tag_i) then
482 >    tag_i = tagRow(atom1)
483   #else
484 <          
494 <             if (rijsq <  rlistsq) then
484 >    tag_i = tag(atom1)
485   #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
486  
487 <    
488 <              if (rijsq <  rcutsq) then
489 <                
490 <                 r = dsqrt(rijsq)
491 <      
492 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
493 <      
494 < #ifdef IS_MPI
495 <                eRow(i) = eRow(i) + pot*0.5
496 <                eCol(i) = eCol(i) + pot*0.5
518 < #else
519 <                    pe = pe + pot
520 < #endif                
521 <            
522 <                drdx = -rxij / r
523 <                drdy = -ryij / r
524 <                drdz = -rzij / r
525 <                
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
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 < #ifdef IS_MPI
499 <     point(nrow + 1) = nlist + 1
564 < #else
565 <     point(natoms) = nlist + 1
566 < #endif
498 > !! we return if j not present here.
499 >    tag_j = tagColumn(j)
500  
568 –  else
569 –
570 –     ! use the list to find the neighbors
571 –     do i = 1, nrow
572 –        JBEG = POINT(i)
573 –        JEND = POINT(i+1) - 1
574 –        ! check thiat molecule i has neighbors
575 –        if (jbeg .le. jend) then
576 – #ifdef IS_MPI
577 –           ljAtype_i => identPtrListRow(i)%this
578 –           rxi = qRow(1,i)
579 –           ryi = qRow(2,i)
580 –           rzi = qRow(3,i)
581 – #else
582 –           ljAtype_i   => identPtrList(i)%this
583 –           rxi = q(1,i)
584 –           ryi = q(2,i)
585 –           rzi = q(3,i)
586 – #endif
587 –           do jnab = jbeg, jend
588 –              j = list(jnab)
589 – #ifdef IS_MPI
590 –              ljAtype_j = identPtrListColumn(j)%this
591 –              rxij = wrap(rxi - qCol(1,j), 1)
592 –              ryij = wrap(ryi - qCol(2,j), 2)
593 –              rzij = wrap(rzi - qCol(3,j), 3)
594 – #else
595 –              ljAtype_j = identPtrList(j)%this
596 –              rxij = wrap(rxi - q(1,j), 1)
597 –              ryij = wrap(ryi - q(2,j), 2)
598 –              rzij = wrap(rzi - q(3,j), 3)
599 – #endif
600 –              rijsq = rxij*rxij + ryij*ryij + rzij*rzij
601 –              
602 –              if (rijsq <  rcutsq) then
603 –
604 –                 r = dsqrt(rijsq)
605 –                
606 –                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
607 – #ifdef IS_MPI
608 –                eRow(i) = eRow(i) + pot*0.5
609 –                eCol(i) = eCol(i) + pot*0.5
610 – #else
611 –                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
637 – #endif
638 –                
639 –                if (do_stress) then
640 –                   tauTemp(1) = tauTemp(1) + fx * rxij
641 –                   tauTemp(2) = tauTemp(2) + fx * ryij
642 –                   tauTemp(3) = tauTemp(3) + fx * rzij
643 –                   tauTemp(4) = tauTemp(4) + fy * rxij
644 –                   tauTemp(5) = tauTemp(5) + fy * ryij
645 –                   tauTemp(6) = tauTemp(6) + fy * rzij
646 –                   tauTemp(7) = tauTemp(7) + fz * rxij
647 –                   tauTemp(8) = tauTemp(8) + fz * ryij
648 –                   tauTemp(9) = tauTemp(9) + fz * rzij
649 –                endif
650 –                
651 –                
652 –             endif
653 –          enddo
654 –       endif
655 –    enddo
656 – endif
501  
502  
503 +    if (tag_i == tag_j) then
504 +       do_cycle = .true.
505 +       return
506 +    end if
507  
508 < #ifdef IS_MPI
509 <    !!distribute forces
508 >    if (tag_i < tag_j) then
509 >       if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
510 >       return
511 >    else                
512 >       if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
513 >    endif
514  
663 –    call scatter(fRow,f,plan_row3d)
664 –    call scatter(fCol,fMPITemp,plan_col3d)
515  
666 –    do i = 1,nlocal
667 –       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
668 –    end do
516  
517 <
518 <    call scatter(tRow,t,plan_row3d)
519 <    call scatter(tCol,tMPITemp,plan_col3d)
520 <    
521 <    do i = 1,nlocal
675 <       t(1:3,i) = t(1:3,i) + tMPITemp(1:3,i)
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
677 –
678 –
679 –    if (do_pot) then
680 –       ! scatter/gather pot_row into the members of my column
681 –       call scatter(eRow,eTemp,plan_row)
523        
683 –       ! scatter/gather pot_local into all other procs
684 –       ! add resultant to get total pot
685 –       do i = 1, nlocal
686 –          pe_local = pe_local + eTemp(i)
687 –       enddo
688 –       if (newtons_thrd) then
689 –          eTemp = 0.0E0_DP
690 –          call scatter(eCol,eTemp,plan_col)
691 –          do i = 1, nlocal
692 –             pe_local = pe_local + eTemp(i)
693 –          enddo
694 –       endif
695 –       pe = pe_local
696 –    endif
524  
525 < #endif
525 >  end function checkExcludes
526  
700 –    potE = pe
527  
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 –
528   end module do_Forces

Diff Legend

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