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

Diff Legend

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