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

Diff Legend

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