ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90 (file contents):
Revision 295 by chuckv, Thu Mar 6 19:57:03 2003 UTC vs.
Revision 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.3 2003-03-06 19:57:03 chuckv Exp $, $Date: 2003-03-06 19:57:03 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
7 > !! @version $Id: do_Forces.F90,v 1.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 –
34 –
35 –  logical, save :: firstTime = .True.
36 –
37 – !! Atype identity pointer lists
38 – #ifdef IS_MPI
39 – !! Row lj_atype pointer list
40 –  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
41 – !! Column lj_atype pointer list
42 –  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
43 – #else
44 –  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
45 – #endif
46 –
47 –
48 – !! Logical has lj force field module been initialized?
49 –  logical, save :: isFFinit = .false.
50 –
51 – !! Use periodic boundry conditions
52 –  logical :: wrap = .false.
53 –
54 – !! Potential energy global module variables
55 – #ifdef IS_MPI
56 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
57 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
58 –
59 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
60 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
61 –
62 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
63 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
64 –
65 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: ARow = 0.0_dp
66 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp
67 –
68 –  
69 –
70 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
71 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
72 –
73 –
74 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
75 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
76 –
77 –
78 –
79 –  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
80 –  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
81 –
82 –  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
83 – #endif
84 –  real(kind = dp) :: pe = 0.0_dp
85 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
86 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp
87 –  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp
88 –
89 –  logical :: do_preForce  = .false.
90 –  logical :: do_postForce = .false.
91 –
92 –
93 –
94 – !! Public methods and data
95 –  public :: new_atype
96 –  public :: do_forceLoop
36    public :: init_FF
37 +  public :: do_forces
38  
99 –  
100 –
101 –
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
108 <    real( kind = dp ), intent(in) :: epsilon
109 <    real( kind = dp ), intent(in) :: sigma
110 <    real( kind = dp ), intent(in) :: w0
111 <    real( kind = dp ), intent(in) :: v0
112 <    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
115 <    integer, intent(out) :: status
116 <    integer, intent(in) :: is_Sticky
117 <    integer, intent(in) :: is_DP
118 <    integer, intent(in) :: is_GB
119 <    integer, intent(in) :: is_LJ
46 >    ! be a smarter subroutine.
47  
48 <
49 <    type (atype), pointer :: the_new_atype
50 <    integer :: alloc_error
51 <    integer :: atype_counter = 0
125 <    integer :: alloc_size
126 <    integer :: err_stat
127 <    status = 0
128 <
129 <
130 <
131 < ! allocate a new atype    
132 <    allocate(the_new_atype,stat=alloc_error)
133 <    if (alloc_error /= 0 ) then
134 <       status = -1
48 >    
49 >    call init_lj_FF(my_status)
50 >    if (my_status /= 0) then
51 >       thisStat = -1
52         return
53      end if
137 –
138 – ! assign our new atype information
139 –    the_new_atype%mass        = mass
140 –    the_new_atype%epsilon     = epsilon
141 –    the_new_atype%sigma       = sigma
142 –    the_new_atype%sigma2      = sigma * sigma
143 –    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
144 –         * the_new_atype%sigma2
145 –    the_new_atype%w0       = w0
146 –    the_new_atype%v0       = v0
147 –    the_new_atype%dipoleMoment       = dipoleMoment
148 –
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.
157 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
61 >    do_forces_initialized = .true.    
62 >    
63 >  end subroutine init_FF
64  
159 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
160 –    if (err_stat /= 0 ) then
161 –       status = -1
162 –       return
163 –    endif
65  
165 –    n_atypes = n_atypes + 1
66  
67 <
68 <  end subroutine new_atype
69 <
70 <
71 <  subroutine init_FF(nComponents,ident, status)
72 < !! Number of components in ident array
73 <    integer, intent(inout) :: nComponents
74 < !! Array of identities nComponents long corresponding to
75 < !! ljatype ident.
76 <    integer, dimension(nComponents),intent(inout) :: ident
77 < !!  Result status, success = 0, error = -1
78 <    integer, intent(out) :: Status
79 <
80 <    integer :: alloc_stat
81 <
82 <    integer :: thisStat
83 <    integer :: i
84 <
85 <    integer :: myNode
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 <    integer, allocatable, dimension(:) :: identRow
89 <    integer, allocatable, dimension(:) :: identCol
88 >    real( kind = DP ) :: pot_local
89 >    integer :: nlocal
90      integer :: nrow
91      integer :: ncol
92   #endif
93 <    status = 0
94 <  
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 <    
101 >    !! initialize local variables  
102  
103 < !! if were're not in MPI, we just update ljatypePtrList
104 < #ifndef IS_MPI
105 <    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
106 <    if ( thisStat /= 0 ) then
107 <       status = -1
108 <       return
109 <    endif
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 +    call getRcut(rcut,rcut2=rcutsq)
113 +    call getRlist(rlist,rlistsq)
114      
115 < ! if were're in MPI, we also have to worry about row and col lists    
116 < #else
208 <  
209 < ! We can only set up forces if mpiSimulation has been setup.
210 <    if (.not. isMPISimSet()) then
211 <       write(default_error,*) "MPI is not set"
212 <       status = -1
213 <       return
214 <    endif
215 <    nrow = getNrow(plan_row)
216 <    ncol = getNcol(plan_col)
217 <    mynode = getMyNode()
218 < !! Allocate temperary arrays to hold gather information
219 <    allocate(identRow(nrow),stat=alloc_stat)
220 <    if (alloc_stat /= 0 ) then
221 <       status = -1
222 <       return
223 <    endif
115 >    call check_initialization()
116 >    call zero_work_arrays()
117  
118 <    allocate(identCol(ncol),stat=alloc_stat)
119 <    if (alloc_stat /= 0 ) then
227 <       status = -1
228 <       return
229 <    endif
118 >    do_pot = do_pot_c
119 >    do_stress = do_stress_c
120  
121 < !! Gather idents into row and column idents
232 <
233 <    call gather(ident,identRow,plan_row)
234 <    call gather(ident,identCol,plan_col)
121 >    ! Gather all information needed by all force loops:
122      
123 <  
237 < !! Create row and col pointer lists
238 <  
239 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
240 <    if (thisStat /= 0 ) then
241 <       status = -1
242 <       return
243 <    endif
244 <  
245 <    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
246 <    if (thisStat /= 0 ) then
247 <       status = -1
248 <       return
249 <    endif
123 > #ifdef IS_MPI    
124  
125 < !! free temporary ident arrays
126 <    if (allocated(identCol)) then
127 <       deallocate(identCol)
128 <    end if
129 <    if (allocated(identCol)) then
130 <       deallocate(identRow)
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 <
135 >    
136   #endif
137      
138 <    call initForce_Modules(thisStat)
139 <    if (thisStat /= 0) then
140 <       status = -1
141 <       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
266 –
267 – !! Create neighbor lists
268 –    call expandList(thisStat)
269 –    if (thisStat /= 0) then
270 –       status = -1
271 –       return
272 –    endif
273 –
274 –    isFFinit = .true.
275 –
276 –
277 –  end subroutine init_FF
278 –
279 –
280 –
281 –
282 –  subroutine initForce_Modules(thisStat)
283 –    integer, intent(out) :: thisStat
284 –    integer :: my_status
149      
286 –    thisStat = 0
287 –    call init_lj_FF(ListHead,my_status)
288 –    if (my_status /= 0) then
289 –       thisStat = -1
290 –       return
291 –    end if
292 –
293 –  end subroutine initForce_Modules
294 –
295 –
296 –
297 –
298 – !! FORCE routine Calculates Lennard Jones forces.
299 – !------------------------------------------------------------->
300 –  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
301 – !! Position array provided by C, dimensioned by getNlocal
302 –    real ( kind = dp ), dimension(3,getNlocal()) :: q
303 –  !! Rotation Matrix for each long range particle in simulation.
304 –    real( kind = dp), dimension(9,getNlocal()) :: A
305 –
306 –  !! Magnitude dipole moment
307 –    real( kind = dp ), dimension(3,getNlocal()) :: mu
308 –  !! Unit vectors for dipoles (lab frame)
309 –    real( kind = dp ), dimension(3,getNlocal()) :: u_l
310 – !! Force array provided by C, dimensioned by getNlocal
311 –    real ( kind = dp ), dimension(3,getNlocal()) :: f
312 – !! Torsion array provided by C, dimensioned by getNlocal
313 –    real( kind = dp ), dimension(3,getNlocal()) :: t
314 –
315 – !! Stress Tensor
316 –    real( kind = dp), dimension(9) :: tau
317 –
318 –    real ( kind = dp ) :: potE
319 –    logical ( kind = 2) :: do_pot
320 –    integer :: FFerror
321 –
322 –    
323 –    type(atype), pointer :: Atype_i
324 –    type(atype), pointer :: Atype_j
325 –
326 –
327 –
328 –
329 –  
330 –
150   #ifdef IS_MPI
332 –  real( kind = DP ) :: pot_local
333 –
334 – !! Local arrays needed for MPI
335 –
336 – #endif
337 –
338 –
339 –
340 –  real( kind = DP )   :: pe
341 –  logical             :: update_nlist
342 –
343 –
344 –  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
345 –  integer :: nlist
346 –  integer :: j_start
347 –
348 –  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
349 –
350 –  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
351 –  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
352 –
353 –  real( kind = DP ) :: dielectric = 0.0_dp
354 –
355 – ! a rig that need to be fixed.
356 – #ifdef IS_MPI
357 –  real( kind = dp ) :: pe_local
358 –  integer :: nlocal
359 – #endif
360 –  integer :: nrow
361 –  integer :: ncol
362 –  integer :: natoms
363 –  integer :: neighborListSize
364 –  integer :: listerror
365 – !! should we calculate the stress tensor
366 –  logical  :: do_stress = .false.
367 –
368 –
369 –  FFerror = 0
370 –
371 – ! Make sure we are properly initialized.
372 –  if (.not. isFFInit) then
373 –     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
374 –     FFerror = -1
375 –     return
376 –  endif
377 – #ifdef IS_MPI
378 –    if (.not. isMPISimSet()) then
379 –     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
380 –     FFerror = -1
381 –     return
382 –  endif
383 – #endif
384 –
385 – !! initialize local variables  
386 –  natoms = getNlocal()
387 –  call getRcut(rcut,rcut2=rcutsq)
388 –  call getRlist(rlist,rlistsq)
389 –
390 – !! Find ensemble
391 –  if (isEnsemble("NPT")) do_stress = .true.
392 – !! set to wrap
393 –  if (isPBC()) wrap = .true.
394 –
395 –
396 – #ifndef IS_MPI
397 –  nrow = natoms - 1
398 –  ncol = natoms
399 – #else
400 –  nrow = getNrow(plan_row)
401 –  ncol = getNcol(plan_col)
402 –  nlocal = natoms
403 –  j_start = 1
404 – #endif
405 –
406 –  
407 – !! See if we need to update neighbor lists
408 –  call check(q,update_nlist)
409 –
410 – !--------------WARNING...........................
411 – ! Zero variables, NOTE:::: Forces are zeroed in C
412 – ! Zeroing them here could delete previously computed
413 – ! Forces.
414 – !------------------------------------------------
415 –  call zero_module_variables()
416 –
417 –
418 – ! communicate MPI positions
419 – #ifdef IS_MPI    
420 –    call gather(q,qRow,plan_row3d)
421 –    call gather(q,qCol,plan_col3d)
422 –
423 –    call gather(mu,muRow,plan_row3d)
424 –    call gather(mu,muCol,plan_col3d)
425 –
426 –    call gather(u_l,u_lRow,plan_row3d)
427 –    call gather(u_l,u_lCol,plan_col3d)
428 –
429 –    call gather(A,ARow,plan_row_rotation)
430 –    call gather(A,ACol,plan_col_rotation)
431 – #endif
432 –
433 –
434 –  if (update_nlist) then
435 –
436 –     ! save current configuration, contruct neighbor list,
437 –     ! and calculate forces
438 –     call save_neighborList(q)
439 –    
440 –     neighborListSize = getNeighborListSize()
441 –     nlist = 0
442 –    
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 <     do i = 1, nrow
446 <        point(i) = nlist + 1
447 < #ifdef IS_MPI
448 <        Atype_i => identPtrListRow(i)%this
449 <        tag_i = tagRow(i)
450 < #else
451 <        Atype_i   => identPtrList(i)%this
452 <        j_start = i + 1
453 < #endif
454 <
455 <        inner: do j = j_start, ncol
456 < ! Assign identity pointers and tags
457 < #ifdef IS_MPI
458 <           Atype_j => identPtrListColumn(j)%this
459 <          
460 <           call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
461 <                rxij,ryij,rzij,rijsq,r)
462 < !! For mpi, use newtons 3rd law when building neigbor list
463 < !! Also check to see the particle i != j.
464 <           if (mpi_cycle_jLoop(i,j)) cycle inner:
465 <
466 < #else          
467 <           Atype_j   => identPtrList(j)%this
468 <           call get_interatomic_vector(i,j,q(:,i),q(:,j),&
469 <                rxij,ryij,rzij,rijsq,r)
192 >       point(nrow + 1) = nlist + 1
193        
194 < #endif          
472 <          
473 <           if (rijsq <  rlistsq) then
194 >    else  !! (of update_check)
195  
196 <              nlist = nlist + 1
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 <              if (nlist > neighborListSize) then
207 <                 call expandList(listerror)
479 <                 if (listerror /= 0) then
480 <                    FFerror = -1
481 <                    write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
482 <                    return
483 <                 end if
484 <              endif
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 <              list(nlist) = j
210 <
209 >             enddo
210 >          endif
211 >       enddo
212 >    endif
213      
214 <              if (rijsq <  rcutsq) then
215 <                 call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
216 <              endif
217 <          enddo inner
218 <     enddo
214 > #else
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 < #ifdef IS_MPI
271 <     point(nrow + 1) = nlist + 1
497 < #else
498 <     point(natoms) = nlist + 1
499 < #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 <  else !! (update)
274 <
275 <     ! use the list to find the neighbors
276 <     do i = 1, nrow
277 <        JBEG = POINT(i)
278 <        JEND = POINT(i+1) - 1
279 <        ! check thiat molecule i has neighbors
280 <        if (jbeg .le. jend) then
281 <
273 >             enddo
274 >          endif
275 >       enddo
276 >    endif
277 >    
278 > #endif
279 >    
280 >    ! phew, done with main loop.
281 >    
282   #ifdef IS_MPI
511 –           ljAtype_i => identPtrListRow(i)%this
512 – #else
513 –           ljAtype_i => identPtrList(i)%this
514 – #endif
515 –           do jnab = jbeg, jend
516 –              j = list(jnab)
517 – #ifdef IS_MPI
518 –              ljAtype_j = identPtrListColumn(j)%this
519 –              call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
520 –                   rxij,ryij,rzij,rijsq,r)
521 –              
522 – #else
523 –              ljAtype_j = identPtrList(j)%this
524 –              call get_interatomic_vector(i,j,q(:,i),q(:,j),&
525 –                   rxij,ryij,rzij,rijsq,r)
526 – #endif
527 –              call do_pair(i,j,r,rxij,ryij,rzij)
528 –          enddo
529 –       endif
530 –    enddo
531 – endif
532 –
533 –
534 –
535 – #ifdef IS_MPI
283      !!distribute forces
284 <
285 <    call scatter(fRow,f,plan_row3d)
286 <    call scatter(fCol,fTemp,plan_col3d)
540 <
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) + fTemp(1:3,i)
288 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
289      end do
544 –
545 –    if (do_torque) then
546 –       call scatter(tRow,t,plan_row3d)
547 –       call scatter(tCol,tTemp,plan_col3d)
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) + tTemp(1:3,i)
296 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
297         end do
298      endif
299 <
299 >    
300      if (do_pot) then
301         ! scatter/gather pot_row into the members of my column
302 <       call scatter(eRow,eTemp,plan_row)
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 <          pe_local = pe_local + eTemp(i)
307 >          pot_local = pot_local + pot_Temp(i)
308         enddo
309  
310 <       eTemp = 0.0E0_DP
311 <       call scatter(eCol,eTemp,plan_col)
310 >       pot_Temp = 0.0_DP
311 >
312 >       call scatter(pot_Col, pot_Temp, plan_col)
313         do i = 1, nlocal
314 <          pe_local = pe_local + eTemp(i)
314 >          pot_local = pot_local + pot_Temp(i)
315         enddo
316        
570 –       pe = pe_local
317      endif
318 < #else
319 < ! Copy local array into return array for c
320 <    f = fTemp
321 <    t = tTemp
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  
578 –    potE = pe
579 –
580 –
581 –    if (do_stress) then
333   #ifdef IS_MPI
334 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
584 <            mpi_comm_world,mpi_err)
334 >             me_i = atid_row(i)
335   #else
336 <       tau = tauTemp
337 < #endif      
336 >             me_i = atid(i)
337 > #endif
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  
590 –  end subroutine do_force_loop
357  
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 +    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 + #endif
381 +    
382 +  end subroutine do_force_loop
383  
384  
599 –
600 –
385   !! Calculate any pre-force loop components and update nlist if necessary.
386    subroutine do_preForce(updateNlist)
387      logical, intent(inout) :: updateNlist
# Line 606 | Line 390 | contains
390  
391    end subroutine do_preForce
392  
609 –
610 –
611 –
612 –
613 –
614 –
615 –
616 –
617 –
618 –
619 –
620 –
393   !! Calculate any post force loop components, i.e. reaction field, etc.
394    subroutine do_postForce()
395  
# Line 625 | Line 397 | contains
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 <
631 <
632 <
633 <
634 <
635 <
636 <
637 <
638 <
639 <
640 <
641 <
642 <
643 <
644 <  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
645 <    type (atype ), pointer, intent(inout) :: atype_i
646 <    type (atype ), pointer, intent(inout) :: atype_j
647 <    integer :: i
648 <    integer :: j
649 <    real ( kind = dp ), intent(inout) :: rx_ij
650 <    real ( kind = dp ), intent(inout) :: ry_ij
651 <    real ( kind = dp ), intent(inout) :: rz_ij
652 <
653 <
654 <    real( kind = dp ) :: fx = 0.0_dp
655 <    real( kind = dp ) :: fy = 0.0_dp
656 <    real( kind = dp ) :: fz = 0.0_dp  
657 <
658 <    real( kind = dp ) ::  drdx = 0.0_dp
659 <    real( kind = dp ) ::  drdy = 0.0_dp
660 <    real( kind = dp ) ::  drdz = 0.0_dp
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 +    me_i = atid_row(i)
418 +    me_j = atid_col(j)
419  
665 –
666 –
667 –    call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j)
668 –      
669 – #ifdef IS_MPI
670 –                eRow(i) = eRow(i) + pot*0.5
671 –                eCol(i) = eCol(i) + pot*0.5
420   #else
673 –                    pe = pe + pot
674 – #endif                
675 –            
676 –                drdx = -rxij / r
677 –                drdy = -ryij / r
678 –                drdz = -rzij / r
679 –                
680 –                fx = dudr * drdx
681 –                fy = dudr * drdy
682 –                fz = dudr * drdz
421  
422 +    me_i = atid(i)
423 +    me_j = atid(j)
424  
425 + #endif
426  
427  
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 <                
690 < #ifdef IS_MPI
691 <                fCol(1,j) = fCol(1,j) - fx
692 <                fCol(2,j) = fCol(2,j) - fy
693 <                fCol(3,j) = fCol(3,j) - fz
694 <                
695 <                fRow(1,j) = fRow(1,j) + fx
696 <                fRow(2,j) = fRow(2,j) + fy
697 <                fRow(3,j) = fRow(3,j) + fz
698 < #else
699 <                fTemp(1,j) = fTemp(1,j) - fx
700 <                fTemp(2,j) = fTemp(2,j) - fy
701 <                fTemp(3,j) = fTemp(3,j) - fz
702 <                fTemp(1,i) = fTemp(1,i) + fx
703 <                fTemp(2,i) = fTemp(2,i) + fy
704 <                fTemp(3,i) = fTemp(3,i) + fz
705 < #endif
706 <                
707 <                if (do_stress) then
708 <                   tauTemp(1) = tauTemp(1) + fx * rxij
709 <                   tauTemp(2) = tauTemp(2) + fx * ryij
710 <                   tauTemp(3) = tauTemp(3) + fx * rzij
711 <                   tauTemp(4) = tauTemp(4) + fy * rxij
712 <                   tauTemp(5) = tauTemp(5) + fy * ryij
713 <                   tauTemp(6) = tauTemp(6) + fy * rzij
714 <                   tauTemp(7) = tauTemp(7) + fz * rxij
715 <                   tauTemp(8) = tauTemp(8) + fz * ryij
716 <                   tauTemp(9) = tauTemp(9) + fz * rzij
717 <                endif
455 >    if (FF_uses_Sticky .and. SimUsesSticky()) then
456  
457 <
458 <
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  
468  
469 <
470 <
726 <
727 <
728 <
729 <
730 <
731 <
732 <
733 <
734 <
735 <
736 <
737 <
738 <  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
739 < !---------------- Arguments-------------------------------
740 <   !! index i
741 <
742 <    !! Position array
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
745 –    !! x component of vector between i and j
746 –    real ( kind = dp ), intent(out)  :: rx_ij
747 –    !! y component of vector between i and j
748 –    real ( kind = dp ), intent(out)  :: ry_ij
749 –    !! z component of vector between i and j
750 –    real ( kind = dp ), intent(out)  :: rz_ij
751 –    !! magnitude of r squared
473      real ( kind = dp ), intent(out) :: r_sq
753 –    !! magnitude of vector r between atoms i and j.
754 –    real ( kind = dp ), intent(out) :: r_ij
755 –    !! wrap into periodic box.
756 –    logical, intent(in) :: wrap
757 –
758 – !--------------- Local Variables---------------------------
759 –    !! Distance between i and j
474      real( kind = dp ) :: d(3)
761 – !---------------- END DECLARATIONS-------------------------
475  
763 –
764 – ! Find distance between i and j
476      d(1:3) = q_i(1:3) - q_j(1:3)
477 <
478 < ! Wrap back into periodic box if necessary
479 <    if ( wrap ) then
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 <    end if
482 >    endif
483      
773 – !   Find Magnitude of the vector
484      r_sq = dot_product(d,d)
485 <    r_ij = sqrt(r_sq)
485 >        
486 >  end subroutine get_interatomic_vector
487  
488 < !   Set each component for force calculation
489 <    rx_ij = d(1)
779 <    ry_ij = d(2)
780 <    rz_ij = d(3)
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 initialized!"
501 +       error = -1
502 +       return
503 +    endif
504 + #endif
505  
506 <  end subroutine get_interatomic_vector
506 >    return
507 >  end subroutine check_initialization
508  
509 <  subroutine zero_module_variables()
509 >  
510 >  subroutine zero_work_arrays()
511 >    
512 > #ifdef IS_MPI
513  
514 < #ifndef IS_MPI
515 <
789 <    pe = 0.0E0_DP
790 <    tauTemp = 0.0_dp
791 <    fTemp = 0.0_dp
792 <    tTemp = 0.0_dp
793 < #else
794 <    qRow = 0.0_dp
795 <    qCol = 0.0_dp
514 >    q_Row = 0.0_dp
515 >    q_Col = 0.0_dp  
516      
517 <    muRow = 0.0_dp
518 <    muCol = 0.0_dp
517 >    u_l_Row = 0.0_dp
518 >    u_l_Col = 0.0_dp
519      
520 <    u_lRow = 0.0_dp
521 <    u_lCol = 0.0_dp
520 >    A_Row = 0.0_dp
521 >    A_Col = 0.0_dp
522      
523 <    ARow = 0.0_dp
524 <    ACol = 0.0_dp
525 <    
526 <    fRow = 0.0_dp
527 <    fCol = 0.0_dp
528 <    
529 <  
810 <    tRow = 0.0_dp
811 <    tCol = 0.0_dp
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 <  
531 >    pot_Row = 0.0_dp
532 >    pot_Col = 0.0_dp
533 >    pot_Temp = 0.0_dp
534  
815 –    eRow = 0.0_dp
816 –    eCol = 0.0_dp
817 –    eTemp = 0.0_dp
535   #endif
536  
537 <  end subroutine zero_module_variables
537 >    tau_Temp = 0.0_dp
538 >    virial_Temp = 0.0_dp
539 >    
540 >  end subroutine zero_work_arrays
541 >  
542  
822 – #ifdef IS_MPI
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 mpi_cycle_jLoop(i,j) result(do_cycle)
546 >  function checkExcludes(atom1,atom2) result(do_cycle)
547   !--------------- Arguments--------------------------
548   ! Index i
549 <    integer,intent(in) :: i
549 >    integer,intent(in) :: atom1
550   ! Index j
551 <    integer,intent(in) :: 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 < !--------------- END DECLARATIONS------------------    
558 <    tag_i = tagRow(i)
557 >    integer :: i
558 > !--------------- END DECLARATIONS------------------  
559 >    do_cycle = .false.
560 >
561 > #ifdef IS_MPI
562 >    tag_i = tagRow(atom1)
563 > #else
564 >    tag_i = tag(atom1)
565 > #endif
566 >
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 > !! we return if j not present here.
579      tag_j = tagColumn(j)
580  
581 <    do_cycle = .false.
581 >
582  
583      if (tag_i == tag_j) then
584         do_cycle = .true.
# Line 851 | Line 591 | contains
591      else                
592         if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
593      endif
854 –  end function mpi_cycle_jLoop
855 – #endif
594  
595 +
596 +
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
603 +      
604 +
605 +  end function checkExcludes
606 +
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