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 306 by chuckv, Mon Mar 10 19:26:45 2003 UTC vs.
Revision 332 by gezelter, Thu Mar 13 15:28:43 2003 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.8 2003-03-10 19:26:45 chuckv Exp $, $Date: 2003-03-10 19:26:45 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
7 > !! @version $Id: do_Forces.F90,v 1.18 2003-03-13 15:28:43 gezelter Exp $, $Date: 2003-03-13 15:28:43 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
8  
9  
10  
11   module do_Forces
12    use simulation
13    use definitions
14 <  use forceGlobals
15 <  use atype_typedefs
16 <  use neighborLists
14 >  use atype_module
15 >  use neighborLists  
16 >  use lj
17 >  use sticky_pair
18 >  use dipole_dipole
19 >  use reaction_field
20  
18 –  
19 –  use lj_FF
20 –  use sticky_FF
21 –  use dp_FF
22 –  use gb_FF
23 –
21   #ifdef IS_MPI
22    use mpiSimulation
23   #endif
24    implicit none
25    PRIVATE
26  
27 < public :: do_force_loop
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  
35 +  public :: init_FF
36 +  public :: do_force_loop
37 +
38   contains
39  
40 < !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
41 < !------------------------------------------------------------->
42 <  subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
43 < !! Position array provided by C, dimensioned by getNlocal
44 <    real ( kind = dp ), dimension(3,getNlocal()) :: q
45 <  !! Rotation Matrix for each long range particle in simulation.
46 <    real( kind = dp), dimension(9,getNlocal()) :: A
40 >  subroutine init_FF(LJ_mix_policy, use_RF_c, thisStat)
41 >    logical(kind = 2), intent(in) :: use_RF_c
42 >    logical :: use_RF_f
43 >    integer, intent(out) :: thisStat  
44 >    integer :: my_status, nMatches
45 >    character(len = 100) :: LJ_mix_Policy
46 >    integer, pointer :: MatchList(:)
47 >    
48 >    !! Fortran's version of a cast:
49 >    use_RF_f = use_RF_c
50  
51 <  !! Magnitude dipole moment
52 <    real( kind = dp ), dimension(3,getNlocal()) :: mu
53 <  !! Unit vectors for dipoles (lab frame)
54 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
55 < !! Force array provided by C, dimensioned by getNlocal
56 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
57 < !! Torsion array provided by C, dimensioned by getNlocal
58 <    real( kind = dp ), dimension(3,getNlocal()) :: t
51 >    !! assume things are copacetic, unless they aren't
52 >    thisStat = 0
53 >    
54 >    !! init_FF is called *after* all of the atom types have been
55 >    !! defined in atype_module using the new_atype subroutine.
56 >    !!
57 >    !! this will scan through the known atypes and figure out what
58 >    !! interactions are used by the force field.    
59  
60 < !! Stress Tensor
61 <    real( kind = dp), dimension(9) :: tau
60 >    FF_uses_LJ = .false.
61 >    FF_uses_sticky = .false.
62 >    FF_uses_dipoles = .false.
63 >    FF_uses_GB = .false.
64 >    FF_uses_EAM = .false.
65 >    
66 >    call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
67 >    deallocate(MatchList)
68 >    if (nMatches .gt. 0) FF_uses_LJ = .true.
69  
70 <    real ( kind = dp ) :: potE
71 <    logical ( kind = 2) :: do_pot
72 <    integer :: FFerror
70 >    call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
71 >    deallocate(MatchList)
72 >    if (nMatches .gt. 0) FF_uses_dipoles = .true.
73  
74 +    call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
75 +         MatchList)
76 +    deallocate(MatchList)
77 +    if (nMatches .gt. 0) FF_uses_Sticky = .true.
78      
79 <    type(atype), pointer :: Atype_i
80 <    type(atype), pointer :: Atype_j
79 >    call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
80 >    deallocate(MatchList)
81 >    if (nMatches .gt. 0) FF_uses_GB = .true.
82  
83 +    call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
84 +    deallocate(MatchList)
85 +    if (nMatches .gt. 0) FF_uses_EAM = .true.
86  
87 +    !! check to make sure the use_RF setting makes sense
88 +    if (use_RF_f) then
89 +       if (FF_uses_dipoles) then
90 +          FF_uses_RF = .true.
91 +          call initialize_rf()
92 +       else
93 +          write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
94 +          thisStat = -1
95 +          return
96 +       endif
97 +    endif
98 +    
99 +    call init_lj_FF(LJ_mix_Policy, my_status)
100 +    if (my_status /= 0) then
101 +       thisStat = -1
102 +       return
103 +    end if
104 +    
105 +    call check_sticky_FF(my_status)
106 +    if (my_status /= 0) then
107 +       thisStat = -1
108 +       return
109 +    end if
110 +    
111 +    do_forces_initialized = .true.    
112 +    
113 +  end subroutine init_FF
114  
115  
65 –  
116  
117 +  !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
118 +  !------------------------------------------------------------->
119 +  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
120 +       error)
121 +    !! Position array provided by C, dimensioned by getNlocal
122 +    real ( kind = dp ), dimension(3,getNlocal()) :: q
123 +    !! Rotation Matrix for each long range particle in simulation.
124 +    real( kind = dp), dimension(9,getNlocal()) :: A    
125 +    !! Unit vectors for dipoles (lab frame)
126 +    real( kind = dp ), dimension(3,getNlocal()) :: u_l
127 +    !! Force array provided by C, dimensioned by getNlocal
128 +    real ( kind = dp ), dimension(3,getNlocal()) :: f
129 +    !! Torsion array provided by C, dimensioned by getNlocal
130 +    real( kind = dp ), dimension(3,getNlocal()) :: t    
131 +    !! Stress Tensor
132 +    real( kind = dp), dimension(9) :: tau  
133 +    real ( kind = dp ) :: pot
134 +    logical ( kind = 2) :: do_pot_c, do_stress_c
135 +    logical :: do_pot
136 +    logical :: do_stress
137   #ifdef IS_MPI
138 <  real( kind = DP ) :: pot_local
139 <
140 < !! Local arrays needed for MPI
71 <
72 < #endif
73 <
74 <
75 <
76 <  real( kind = DP )   :: pe
77 <  logical             :: update_nlist
78 <
79 <
80 <  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
81 <  integer :: nlist
82 <  integer :: j_start
83 <
84 <  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
85 <
86 <  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
87 <  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
88 <
89 <  
90 <
91 < ! a rig that need to be fixed.
92 < #ifdef IS_MPI
93 <  real( kind = dp ) :: pe_local
94 <  integer :: nlocal
138 >    real( kind = DP ) :: pot_local
139 >    integer :: nrow
140 >    integer :: ncol
141   #endif
142 <  integer :: nrow
143 <  integer :: ncol
144 <  integer :: natoms
145 <  integer :: neighborListSize
146 <  integer :: listerror
147 < !! should we calculate the stress tensor
148 <  logical  :: do_stress = .false.
142 >    integer :: nlocal
143 >    integer :: natoms    
144 >    logical :: update_nlist  
145 >    integer :: i, j, jbeg, jend, jnab
146 >    integer :: nlist
147 >    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
148 >    real(kind=dp),dimension(3) :: d
149 >    real(kind=dp) :: rfpot, mu_i, virial
150 >    integer :: me_i
151 >    logical :: is_dp_i
152 >    integer :: neighborListSize
153 >    integer :: listerror, error
154 >    integer :: localError
155  
156 +    !! initialize local variables  
157  
105 –  FFerror = 0
106 –
107 – ! Make sure we are properly initialized.
108 –  if (.not. isFFInit) then
109 –     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
110 –     FFerror = -1
111 –     return
112 –  endif
158   #ifdef IS_MPI
159 <    if (.not. isMPISimSet()) then
160 <     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
161 <     FFerror = -1
162 <     return
163 <  endif
159 >    nlocal = getNlocal()
160 >    nrow   = getNrow(plan_row)
161 >    ncol   = getNcol(plan_col)
162 > #else
163 >    nlocal = getNlocal()
164 >    natoms = nlocal
165   #endif
166  
167 < !! initialize local variables  
168 <  natoms = getNlocal()
169 <  call getRcut(rcut,rcut2=rcutsq)
170 <  call getRlist(rlist,rlistsq)
167 >    call getRcut(rcut,rc2=rcutsq)
168 >    call getRlist(rlist,rlistsq)
169 >    
170 >    call check_initialization(localError)
171 >    if ( localError .ne. 0 ) then
172 >       error = -1
173 >       return
174 >    end if
175 >    call zero_work_arrays()
176  
177 < !! Find ensemble
178 <  if (isEnsemble("NPT")) do_stress = .true.
128 < !! set to wrap
129 <  if (isPBC()) wrap = .true.
177 >    do_pot = do_pot_c
178 >    do_stress = do_stress_c
179  
180 <
181 <
133 <  
134 < !! See if we need to update neighbor lists
135 <  call check(q,update_nlist)
136 <
137 < !--------------WARNING...........................
138 < ! Zero variables, NOTE:::: Forces are zeroed in C
139 < ! Zeroing them here could delete previously computed
140 < ! Forces.
141 < !------------------------------------------------
142 <  call zero_module_variables()
143 <
144 <
145 < ! communicate MPI positions
180 >    ! Gather all information needed by all force loops:
181 >    
182   #ifdef IS_MPI    
147 –    call gather(q,qRow,plan_row3d)
148 –    call gather(q,qCol,plan_col3d)
183  
184 <    call gather(mu,muRow,plan_row3d)
185 <    call gather(mu,muCol,plan_col3d)
186 <
187 <    call gather(u_l,u_lRow,plan_row3d)
188 <    call gather(u_l,u_lCol,plan_col3d)
189 <
190 <    call gather(A,ARow,plan_row_rotation)
191 <    call gather(A,ACol,plan_col_rotation)
184 >    call gather(q,q_Row,plan_row3d)
185 >    call gather(q,q_Col,plan_col3d)
186 >        
187 >    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
188 >       call gather(u_l,u_l_Row,plan_row3d)
189 >       call gather(u_l,u_l_Col,plan_col3d)
190 >      
191 >       call gather(A,A_Row,plan_row_rotation)
192 >       call gather(A,A_Col,plan_col_rotation)
193 >    endif
194 >    
195   #endif
196 <
197 <
196 >    
197 >    if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
198 >       !! See if we need to update neighbor lists
199 >       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
200 >       !! if_mpi_gather_stuff_for_prepair
201 >       !! do_prepair_loop_if_needed
202 >       !! if_mpi_scatter_stuff_from_prepair
203 >       !! if_mpi_gather_stuff_from_prepair_to_main_loop
204 >    else
205 >       !! See if we need to update neighbor lists
206 >       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
207 >    endif
208 >    
209   #ifdef IS_MPI
210      
211      if (update_nlist) then
212        
213 <       ! save current configuration, contruct neighbor list,
214 <       ! and calculate forces
213 >       !! save current configuration, construct neighbor list,
214 >       !! and calculate forces
215         call save_neighborList(q)
216        
217         neighborListSize = getNeighborListSize()
218 <       nlist = 0
218 >       nlist = 0      
219        
172 –       nrow = getNrow(plan_row)
173 –       ncol = getNcol(plan_col)
174 –       nlocal = getNlocal()
175 –      
220         do i = 1, nrow
221            point(i) = nlist + 1
178 –          Atype_i => identPtrListRow(i)%this
222            
223            inner: do j = 1, ncol
181 –             Atype_j => identPtrListColumn(j)%this
224              
225 <             call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
184 <                  rxij,ryij,rzij,rijsq,r)
225 >             if (checkExcludes(i,j)) cycle inner
226              
227 <             ! skip the loop if the atoms are identical
187 <             if (mpi_cycle_jLoop(i,j)) cycle inner:
227 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
228              
229               if (rijsq <  rlistsq) then            
230                  
231                  nlist = nlist + 1
232                  
233                  if (nlist > neighborListSize) then
234 <                   call expandList(listerror)
234 >                   call expandNeighborList(nlocal, listerror)
235                     if (listerror /= 0) then
236 <                      FFerror = -1
236 >                      error = -1
237                        write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
238                        return
239                     end if
240                  endif
241                  
242                  list(nlist) = j
243 <                
204 <                
243 >                                
244                  if (rijsq <  rcutsq) then
245 <                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
245 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
246                  endif
247               endif
248            enddo inner
# Line 211 | Line 250 | contains
250  
251         point(nrow + 1) = nlist + 1
252        
253 <    else !! (update)
253 >    else  !! (of update_check)
254  
255         ! use the list to find the neighbors
256         do i = 1, nrow
# Line 219 | Line 258 | contains
258            JEND = POINT(i+1) - 1
259            ! check thiat molecule i has neighbors
260            if (jbeg .le. jend) then
261 <
223 <             Atype_i => identPtrListRow(i)%this
261 >            
262               do jnab = jbeg, jend
263                  j = list(jnab)
264 <                Atype_j = identPtrListColumn(j)%this
265 <                call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
266 <                     rxij,ryij,rzij,rijsq,r)
267 <                
230 <                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
264 >
265 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
266 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
267 >
268               enddo
269            endif
270         enddo
271      endif
272 <
272 >    
273   #else
274      
275      if (update_nlist) then
# Line 244 | Line 281 | contains
281         neighborListSize = getNeighborListSize()
282         nlist = 0
283        
247 –    
284         do i = 1, natoms-1
285            point(i) = nlist + 1
286 <          Atype_i   => identPtrList(i)%this
251 <
286 >          
287            inner: do j = i+1, natoms
288 <             Atype_j   => identPtrList(j)%this
289 <             call get_interatomic_vector(i,j,q(:,i),q(:,j),&
290 <                  rxij,ryij,rzij,rijsq,r)
288 >            
289 >             if (checkExcludes(i,j)) cycle inner
290 >            
291 >             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
292            
293               if (rijsq <  rlistsq) then
294 <
294 >                
295                  nlist = nlist + 1
296                  
297                  if (nlist > neighborListSize) then
298 <                   call expandList(listerror)
298 >                   call expandList(natoms, listerror)
299                     if (listerror /= 0) then
300 <                      FFerror = -1
300 >                      error = -1
301                        write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
302                        return
303                     end if
304                  endif
305                  
306                  list(nlist) = j
307 <
272 <    
307 >                
308                  if (rijsq <  rcutsq) then
309 <                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
309 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
310                  endif
311               endif
312            enddo inner
# Line 282 | Line 317 | contains
317      else !! (update)
318        
319         ! use the list to find the neighbors
320 <       do i = 1, nrow
320 >       do i = 1, natoms-1
321            JBEG = POINT(i)
322            JEND = POINT(i+1) - 1
323            ! check thiat molecule i has neighbors
324            if (jbeg .le. jend) then
325              
291 –             Atype_i => identPtrList(i)%this
326               do jnab = jbeg, jend
327                  j = list(jnab)
328 <                Atype_j = identPtrList(j)%this
329 <                call get_interatomic_vector(i,j,q(:,i),q(:,j),&
330 <                     rxij,ryij,rzij,rijsq,r)
331 <                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
328 >
329 >                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
330 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
331 >
332               enddo
333            endif
334         enddo
335      endif
336 <
336 >    
337   #endif
338 <
339 <
338 >    
339 >    ! phew, done with main loop.
340 >    
341   #ifdef IS_MPI
307 –    !! distribute all reaction field stuff (or anything for post-pair):
308 –    call scatter(rflRow,rflTemp1,plan_row3d)
309 –    call scatter(rflCol,rflTemp2,plan_col3d)
310 –    do i = 1,nlocal
311 –       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
312 –    end do
313 – #endif
314 –
315 – ! This is the post-pair loop:
316 – #ifdef IS_MPI
317 –
318 –    if (system_has_postpair_atoms) then
319 –       do i = 1, nlocal
320 –          Atype_i => identPtrListRow(i)%this
321 –          call do_postpair(i, Atype_i)
322 –       enddo
323 –    endif
324 –
325 – #else
326 –
327 –    if (system_has_postpair_atoms) then
328 –       do i = 1, natoms
329 –          Atype_i => identPtr(i)%this
330 –          call do_postpair(i, Atype_i)
331 –       enddo
332 –    endif
333 –
334 – #endif
335 –
336 –
337 –
338 –
339 – #ifdef IS_MPI
342      !!distribute forces
343 <
344 <    call scatter(fRow,fTemp1,plan_row3d)
345 <    call scatter(fCol,fTemp2,plan_col3d)
344 <
345 <
343 >    
344 >    call scatter(f_Row,f,plan_row3d)
345 >    call scatter(f_Col,f_temp,plan_col3d)
346      do i = 1,nlocal
347 <       fTemp(1:3,i) = fTemp1(1:3,i) + fTemp2(1:3,i)
347 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
348      end do
349 –
350 –    if (do_torque) then
351 –       call scatter(tRow,tTemp1,plan_row3d)
352 –       call scatter(tCol,tTemp2,plan_col3d)
349      
350 +    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
351 +       call scatter(t_Row,t,plan_row3d)
352 +       call scatter(t_Col,t_temp,plan_col3d)
353 +      
354         do i = 1,nlocal
355 <          tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(1:3,i)
355 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
356         end do
357      endif
358 <
358 >    
359      if (do_pot) then
360         ! scatter/gather pot_row into the members of my column
361 <       call scatter(eRow,eTemp,plan_row)
361 >       call scatter(pot_Row, pot_Temp, plan_row)
362        
363         ! scatter/gather pot_local into all other procs
364         ! add resultant to get total pot
365         do i = 1, nlocal
366 <          pe_local = pe_local + eTemp(i)
366 >          pot_local = pot_local + pot_Temp(i)
367         enddo
368  
369 <       eTemp = 0.0E0_DP
370 <       call scatter(eCol,eTemp,plan_col)
369 >       pot_Temp = 0.0_DP
370 >
371 >       call scatter(pot_Col, pot_Temp, plan_col)
372         do i = 1, nlocal
373 <          pe_local = pe_local + eTemp(i)
373 >          pot_local = pot_local + pot_Temp(i)
374         enddo
375        
376 <       pe = pe_local
376 <    endif
377 < #else
378 < ! Copy local array into return array for c
379 <    f = f+fTemp
380 <    t = t+tTemp
376 >    endif    
377   #endif
378  
379 <    potE = pe
379 >    if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
380 >      
381 >       if (FF_uses_RF .and. SimUsesRF()) then
382 >          
383 > #ifdef IS_MPI
384 >          call scatter(rf_Row,rf,plan_row3d)
385 >          call scatter(rf_Col,rf_Temp,plan_col3d)
386 >          do i = 1,nlocal
387 >             rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
388 >          end do
389 > #endif
390 >          
391 >          do i = 1, getNlocal()
392  
393 <
386 <    if (do_stress) then
393 >             rfpot = 0.0_DP
394   #ifdef IS_MPI
395 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
389 <            mpi_comm_world,mpi_err)
395 >             me_i = atid_row(i)
396   #else
397 <       tau = tauTemp
398 < #endif      
397 >             me_i = atid(i)
398 > #endif
399 >             call getElementProperty(atypes, me_i, "is_DP", is_DP_i)      
400 >             if ( is_DP_i ) then
401 >                call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
402 >                !! The reaction field needs to include a self contribution
403 >                !! to the field:
404 >                call accumulate_self_rf(i, mu_i, u_l)            
405 >                !! Get the reaction field contribution to the
406 >                !! potential and torques:
407 >                call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
408 > #ifdef IS_MPI
409 >                pot_local = pot_local + rfpot
410 > #else
411 >                pot = pot + rfpot
412 > #endif
413 >             endif            
414 >          enddo
415 >       endif
416      endif
417  
395 –  end subroutine do_force_loop
418  
419 <
419 > #ifdef IS_MPI
420  
421 +    if (do_pot) then
422 +       pot = pot_local
423 +       !! we assume the c code will do the allreduce to get the total potential
424 +       !! we could do it right here if we needed to...
425 +    endif
426  
427 +    if (do_stress) then
428 +       call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
429 +            mpi_comm_world,mpi_err)
430 +       call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
431 +            mpi_comm_world,mpi_err)
432 +    endif
433  
434 + #else
435  
436 +    if (do_stress) then
437 +       tau = tau_Temp
438 +       virial = virial_Temp
439 +    endif
440  
441 + #endif
442 +    
443 +  end subroutine do_force_loop
444  
445 +  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress)
446  
447 +    real( kind = dp ) :: pot
448 +    real( kind = dp ), dimension(3,getNlocal()) :: u_l
449 +    real (kind=dp), dimension(9,getNlocal()) :: A
450 +    real (kind=dp), dimension(3,getNlocal()) :: f
451 +    real (kind=dp), dimension(3,getNlocal()) :: t
452  
453 < !! Calculate any pre-force loop components and update nlist if necessary.
454 <  subroutine do_preForce(updateNlist)
455 <    logical, intent(inout) :: updateNlist
453 >    logical, intent(inout) :: do_pot, do_stress
454 >    integer, intent(in) :: i, j
455 >    real ( kind = dp ), intent(inout)    :: rijsq
456 >    real ( kind = dp )                :: r
457 >    real ( kind = dp ), intent(inout) :: d(3)
458 >    logical :: is_LJ_i, is_LJ_j
459 >    logical :: is_DP_i, is_DP_j
460 >    logical :: is_Sticky_i, is_Sticky_j
461 >    integer :: me_i, me_j
462  
463 <
411 <
412 <  end subroutine do_preForce
413 <
414 <
415 <
416 <
417 <
418 <
419 <
420 <
421 <
422 <
423 <
424 <
425 <
426 < !! Calculate any post force loop components, i.e. reaction field, etc.
427 <  subroutine do_postForce()
428 <
429 <
430 <
431 <  end subroutine do_postForce
432 <
433 <
434 <
435 <
436 <
437 <
438 <
439 <
440 <
441 <
442 <
443 <
444 <
445 <
446 <
447 <
448 <
449 <  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
450 <    type (atype ), pointer, intent(inout) :: atype_i
451 <    type (atype ), pointer, intent(inout) :: atype_j
452 <    integer :: i
453 <    integer :: j
454 <    real ( kind = dp ), intent(inout) :: rx_ij
455 <    real ( kind = dp ), intent(inout) :: ry_ij
456 <    real ( kind = dp ), intent(inout) :: rz_ij
457 <
458 <
459 <    real( kind = dp ) :: fx = 0.0_dp
460 <    real( kind = dp ) :: fy = 0.0_dp
461 <    real( kind = dp ) :: fz = 0.0_dp  
462 <  
463 <    real( kind = dp ) ::  drdx = 0.0_dp
464 <    real( kind = dp ) ::  drdy = 0.0_dp
465 <    real( kind = dp ) ::  drdz = 0.0_dp
463 >    r = sqrt(rijsq)
464      
467 –
465   #ifdef IS_MPI
466  
467 <    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
468 <       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
472 <    endif
467 >    me_i = atid_row(i)
468 >    me_j = atid_col(j)
469  
474 –    if (Atype_i%is_dp .and. Atype_j%is_dp) then
475 –
476 –       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
477 –            ulRow(:,i), ulCol(:,j), rt, rrf, pot)
478 –
479 –       if (do_reaction_field) then
480 –          call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
481 –               ulRow(:i), ulCol(:,j), rt, rrf)
482 –       endif
483 –
484 –    endif
485 –
486 –    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
487 –       call getstickyforce(r, pot, dudr, Atype_i, Atype_j)
488 –    endif
489 –
470   #else
471  
472 <    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
473 <       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
494 <    endif
472 >    me_i = atid(i)
473 >    me_j = atid(j)
474  
475 <    if (Atype_i%is_dp .and. Atype_j%is_dp) then
497 <       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
498 <            ul(:,i), ul(:,j), rt, rrf, pot)
475 > #endif
476  
500 –       if (do_reaction_field) then
501 –          call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
502 –               ul(:,i), ul(:,j), rt, rrf)
503 –       endif
477  
478 +    if (FF_uses_LJ .and. SimUsesLJ()) then
479 +       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
480 +       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
481 +      
482 +       if ( is_LJ_i .and. is_LJ_j ) &
483 +            call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
484      endif
485 +      
486  
487 <    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
488 <       call getstickyforce(r,pot,dudr, Atype_i, Atype_j)
489 <    endif
510 <
511 < #endif
512 <
487 >    if (FF_uses_dipoles .and. SimUsesDipoles()) then
488 >       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
489 >       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
490        
491 < #ifdef IS_MPI
492 <                eRow(i) = eRow(i) + pot*0.5
493 <                eCol(i) = eCol(i) + pot*0.5
494 < #else
495 <                    pe = pe + pot
519 < #endif                
491 >       if ( is_DP_i .and. is_DP_j ) then
492 >          
493 >          call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
494 >          
495 >          if (FF_uses_RF .and. SimUsesRF()) then
496              
497 <                drdx = -rxij / r
498 <                drdy = -ryij / r
499 <                drdz = -rzij / r
500 <                
501 <                fx = dudr * drdx
502 <                fy = dudr * drdy
503 <                fz = dudr * drdz
528 <                
529 < #ifdef IS_MPI
530 <                fCol(1,j) = fCol(1,j) - fx
531 <                fCol(2,j) = fCol(2,j) - fy
532 <                fCol(3,j) = fCol(3,j) - fz
533 <                
534 <                fRow(1,j) = fRow(1,j) + fx
535 <                fRow(2,j) = fRow(2,j) + fy
536 <                fRow(3,j) = fRow(3,j) + fz
537 < #else
538 <                fTemp(1,j) = fTemp(1,j) - fx
539 <                fTemp(2,j) = fTemp(2,j) - fy
540 <                fTemp(3,j) = fTemp(3,j) - fz
541 <                fTemp(1,i) = fTemp(1,i) + fx
542 <                fTemp(2,i) = fTemp(2,i) + fy
543 <                fTemp(3,i) = fTemp(3,i) + fz
544 < #endif
545 <                
546 <                if (do_stress) then
547 <                   tauTemp(1) = tauTemp(1) + fx * rxij
548 <                   tauTemp(2) = tauTemp(2) + fx * ryij
549 <                   tauTemp(3) = tauTemp(3) + fx * rzij
550 <                   tauTemp(4) = tauTemp(4) + fy * rxij
551 <                   tauTemp(5) = tauTemp(5) + fy * ryij
552 <                   tauTemp(6) = tauTemp(6) + fy * rzij
553 <                   tauTemp(7) = tauTemp(7) + fz * rxij
554 <                   tauTemp(8) = tauTemp(8) + fz * ryij
555 <                   tauTemp(9) = tauTemp(9) + fz * rzij
556 <                endif
497 >             call accumulate_rf(i, j, r, u_l)
498 >             call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
499 >            
500 >          endif
501 >          
502 >       endif
503 >    endif
504  
505 +    if (FF_uses_Sticky .and. SimUsesSticky()) then
506  
507 <
507 >       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
508 >       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
509 >      
510 >       if ( is_Sticky_i .and. is_Sticky_j ) then
511 >          call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
512 >               do_pot, do_stress)
513 >       endif
514 >    endif
515 >      
516    end subroutine do_pair
517  
518  
519 +  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
520 +    
521 +    real (kind = dp), dimension(3) :: q_i
522 +    real (kind = dp), dimension(3) :: q_j
523 +    real ( kind = dp ), intent(out) :: r_sq
524 +    real( kind = dp ) :: d(3)
525  
526 +    d(1:3) = q_i(1:3) - q_j(1:3)
527 +    
528 +    ! Wrap back into periodic box if necessary
529 +    if ( SimUsesPBC() ) then
530 +       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,box(1:3)) * &
531 +            int(abs(d(1:3)/box(1:3) + 0.5_dp))
532 +    endif
533 +    
534 +    r_sq = dot_product(d,d)
535 +        
536 +  end subroutine get_interatomic_vector
537  
538 <
539 <
540 <
541 <
542 <
543 <
544 <
545 <
538 >  subroutine check_initialization(error)
539 >    integer, intent(out) :: error
540 >    
541 >    error = 0
542 >    ! Make sure we are properly initialized.
543 >    if (.not. do_forces_initialized) then
544 >       write(default_error,*) "ERROR: do_Forces has not been initialized!"
545 >       error = -1
546 >       return
547 >    endif
548  
549 + #ifdef IS_MPI
550 +    if (.not. isMPISimSet()) then
551 +       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
552 +       error = -1
553 +       return
554 +    endif
555 + #endif
556 +    
557 +    return
558 +  end subroutine check_initialization
559  
560 <
561 <
577 <  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
578 < !---------------- Arguments-------------------------------
579 <   !! index i
580 <
581 <    !! Position array
582 <    real (kind = dp), dimension(3) :: q_i
583 <    real (kind = dp), dimension(3) :: q_j
584 <    !! x component of vector between i and j
585 <    real ( kind = dp ), intent(out)  :: rx_ij
586 <    !! y component of vector between i and j
587 <    real ( kind = dp ), intent(out)  :: ry_ij
588 <    !! z component of vector between i and j
589 <    real ( kind = dp ), intent(out)  :: rz_ij
590 <    !! magnitude of r squared
591 <    real ( kind = dp ), intent(out) :: r_sq
592 <    !! magnitude of vector r between atoms i and j.
593 <    real ( kind = dp ), intent(out) :: r_ij
594 <    !! wrap into periodic box.
595 <    logical, intent(in) :: wrap
596 <
597 < !--------------- Local Variables---------------------------
598 <    !! Distance between i and j
599 <    real( kind = dp ) :: d(3)
600 < !---------------- END DECLARATIONS-------------------------
601 <
602 <
603 < ! Find distance between i and j
604 <    d(1:3) = q_i(1:3) - q_j(1:3)
605 <
606 < ! Wrap back into periodic box if necessary
607 <    if ( wrap ) then
608 <       d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
609 <            int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
610 <    end if
560 >  
561 >  subroutine zero_work_arrays()
562      
563 < !   Find Magnitude of the vector
613 <    r_sq = dot_product(d,d)
614 <    r_ij = sqrt(r_sq)
563 > #ifdef IS_MPI
564  
565 < !   Set each component for force calculation
566 <    rx_ij = d(1)
618 <    ry_ij = d(2)
619 <    rz_ij = d(3)
620 <
621 <
622 <  end subroutine get_interatomic_vector
623 <
624 <  subroutine zero_module_variables()
625 <
626 < #ifndef IS_MPI
627 <
628 <    pe = 0.0E0_DP
629 <    tauTemp = 0.0_dp
630 <    fTemp = 0.0_dp
631 <    tTemp = 0.0_dp
632 < #else
633 <    qRow = 0.0_dp
634 <    qCol = 0.0_dp
565 >    q_Row = 0.0_dp
566 >    q_Col = 0.0_dp  
567      
568 <    muRow = 0.0_dp
569 <    muCol = 0.0_dp
568 >    u_l_Row = 0.0_dp
569 >    u_l_Col = 0.0_dp
570      
571 <    u_lRow = 0.0_dp
572 <    u_lCol = 0.0_dp
571 >    A_Row = 0.0_dp
572 >    A_Col = 0.0_dp
573      
574 <    ARow = 0.0_dp
575 <    ACol = 0.0_dp
576 <    
577 <    fRow = 0.0_dp
578 <    fCol = 0.0_dp
579 <    
580 <  
649 <    tRow = 0.0_dp
650 <    tCol = 0.0_dp
574 >    f_Row = 0.0_dp
575 >    f_Col = 0.0_dp
576 >    f_Temp = 0.0_dp
577 >      
578 >    t_Row = 0.0_dp
579 >    t_Col = 0.0_dp
580 >    t_Temp = 0.0_dp
581  
582 <  
582 >    pot_Row = 0.0_dp
583 >    pot_Col = 0.0_dp
584 >    pot_Temp = 0.0_dp
585  
586 <    eRow = 0.0_dp
587 <    eCol = 0.0_dp
588 <    eTemp = 0.0_dp
586 >    rf_Row = 0.0_dp
587 >    rf_Col = 0.0_dp
588 >    rf_Temp = 0.0_dp
589 >
590   #endif
591  
592 <  end subroutine zero_module_variables
592 >    rf = 0.0_dp
593 >    tau_Temp = 0.0_dp
594 >    virial_Temp = 0.0_dp
595 >    
596 >  end subroutine zero_work_arrays
597 >  
598  
599 +  !! Function to properly build neighbor lists in MPI using newtons 3rd law.
600 +  !! We don't want 2 processors doing the same i j pair twice.
601 +  !! Also checks to see if i and j are the same particle.
602  
662 – !! Function to properly build neighbor lists in MPI using newtons 3rd law.
663 – !! We don't want 2 processors doing the same i j pair twice.
664 – !! Also checks to see if i and j are the same particle.
603    function checkExcludes(atom1,atom2) result(do_cycle)
604 < !--------------- Arguments--------------------------
667 < ! Index i
604 >    
605      integer,intent(in) :: atom1
669 – ! Index j
606      integer,intent(in), optional :: atom2
671 – ! Result do_cycle
607      logical :: do_cycle
608 < !--------------- Local variables--------------------
609 <    integer :: tag_i
675 <    integer :: tag_j
676 <    integer :: i
677 < !--------------- END DECLARATIONS------------------  
678 <    do_cycle = .false.
608 >    integer :: unique_id_1, unique_id_2
609 >    integer :: i, j
610  
611 +    do_cycle = .false.
612 +    
613   #ifdef IS_MPI
614 <    tag_i = tagRow(atom1)
614 >    unique_id_1 = tagRow(atom1)
615   #else
616 <    tag_i = tag(atom1)
616 >    unique_id_1 = tag(atom1)
617   #endif
618 <
619 < !! Check global excludes first
618 >    
619 >    !! Check global excludes first
620      if (.not. present(atom2)) then
621 <       do i = 1,nGlobalExcludes
622 <          if (excludeGlobal(i) == tag_i) then
621 >       do i = 1, nExcludes_global
622 >          if (excludesGlobal(i) == unique_id_1) then
623               do_cycle = .true.
624               return
625            end if
# Line 694 | Line 627 | contains
627         return !! return after checking globals
628      end if
629  
630 < !! we return if j not present here.
631 <    tag_j = tagColumn(j)
632 <
633 <
634 <
635 <    if (tag_i == tag_j) then
630 >    !! we return if atom2 not present here.
631 >    
632 > #ifdef IS_MPI
633 >    unique_id_2 = tagColumn(atom2)
634 > #else
635 >    unique_id_2 = tag(atom2)
636 > #endif
637 >    
638 >    if (unique_id_1 == unique_id_2) then
639         do_cycle = .true.
640         return
641      end if
642 <
643 <    if (tag_i < tag_j) then
644 <       if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
642 >    
643 >    if (unique_id_1 < unique_id_2) then
644 >       if (mod(unique_id_1 + unique_id_2,2) == 0) do_cycle = .true.
645         return
646      else                
647 <       if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
647 >       if (mod(unique_id_1 + unique_id_2,2) == 1) do_cycle = .true.
648      endif
649 <
650 <
651 <
652 <    do i = 1, nLocalExcludes
717 <       if (tag_i = excludes(1,i) .and. excludes(2,i) < 0) then
649 >    
650 >    do i = 1, nExcludes_local
651 >       if ((unique_id_1 == excludesLocal(1,i)) .and.  &
652 >            (excludesLocal(2,i) < 0)) then
653            do_cycle = .true.
654            return
655         end if
656      end do
657 <      
723 <
657 >    
658    end function checkExcludes
659  
660 <
660 >  function FF_UsesDirectionalAtoms() result(doesit)
661 >    logical :: doesit
662 >    doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
663 >         FF_uses_GB .or. FF_uses_RF
664 >  end function FF_UsesDirectionalAtoms
665 >  
666 >  function FF_RequiresPrepairCalc() result(doesit)
667 >    logical :: doesit
668 >    doesit = FF_uses_EAM
669 >  end function FF_RequiresPrepairCalc
670 >  
671 >  function FF_RequiresPostpairCalc() result(doesit)
672 >    logical :: doesit
673 >    doesit = FF_uses_RF
674 >  end function FF_RequiresPostpairCalc
675 >  
676   end module do_Forces

Diff Legend

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