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 297 by gezelter, Thu Mar 6 22:08:29 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.4 2003-03-06 22:08:29 gezelter Exp $, $Date: 2003-03-06 22:08:29 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
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(9,getNrow(plan_row)) :: ARow = 0.0_dp
66 –  real(kind = dp), dimension(9,getNcol(plan_col)) :: ACol = 0.0_dp  
67 –
68 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
69 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
70 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp1 = 0.0_dp
71 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp1 = 0.0_dp
72 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp2 = 0.0_dp
73 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp2 = 0.0_dp
74 –  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
75 –  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp
76 –
77 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
78 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
79 –
80 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: rflRow = 0.0_dp
81 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: rflCol = 0.0_dp
82 –  real(kind = dp), dimension(3,getNlocal()) :: rflTemp = 0.0_dp
83 –
84 –  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
85 –  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
86 –
87 –  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
88 – #endif
89 –  real(kind = dp) :: pe = 0.0_dp
90 –  real(kind = dp), dimension(3,natoms) :: fTemp = 0.0_dp
91 –  real(kind = dp), dimension(3,natoms) :: tTemp = 0.0_dp
92 –  real(kind = dp), dimension(3,natoms) :: rflTemp = 0.0_dp
93 –  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp
94 –
95 –  logical :: do_preForce  = .false.
96 –  logical :: do_postForce = .false.
97 –
98 –
99 –
100 – !! Public methods and data
101 –  public :: new_atype
102 –  public :: do_forceLoop
36    public :: init_FF
37 +  public :: do_forces
38  
105 –  
106 –
107 –
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
114 <    real( kind = dp ), intent(in) :: epsilon
115 <    real( kind = dp ), intent(in) :: sigma
116 <    real( kind = dp ), intent(in) :: w0
117 <    real( kind = dp ), intent(in) :: v0
118 <    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
121 <    integer, intent(out) :: status
122 <    integer, intent(in) :: is_Sticky
123 <    integer, intent(in) :: is_DP
124 <    integer, intent(in) :: is_GB
125 <    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
131 <    integer :: alloc_size
132 <    integer :: err_stat
133 <    status = 0
134 <
135 <
136 <
137 < ! allocate a new atype    
138 <    allocate(the_new_atype,stat=alloc_error)
139 <    if (alloc_error /= 0 ) then
140 <       status = -1
48 >    
49 >    call init_lj_FF(my_status)
50 >    if (my_status /= 0) then
51 >       thisStat = -1
52         return
53      end if
143 –
144 – ! assign our new atype information
145 –    the_new_atype%mass        = mass
146 –    the_new_atype%epsilon     = epsilon
147 –    the_new_atype%sigma       = sigma
148 –    the_new_atype%sigma2      = sigma * sigma
149 –    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
150 –         * the_new_atype%sigma2
151 –    the_new_atype%w0       = w0
152 –    the_new_atype%v0       = v0
153 –    the_new_atype%dipoleMoment       = dipoleMoment
154 –
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.
163 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
61 >    do_forces_initialized = .true.    
62 >    
63 >  end subroutine init_FF
64  
165 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
166 –    if (err_stat /= 0 ) then
167 –       status = -1
168 –       return
169 –    endif
65  
171 –    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
214 <  
215 < ! We can only set up forces if mpiSimulation has been setup.
216 <    if (.not. isMPISimSet()) then
217 <       write(default_error,*) "MPI is not set"
218 <       status = -1
219 <       return
220 <    endif
221 <    nrow = getNrow(plan_row)
222 <    ncol = getNcol(plan_col)
223 <    mynode = getMyNode()
224 < !! Allocate temperary arrays to hold gather information
225 <    allocate(identRow(nrow),stat=alloc_stat)
226 <    if (alloc_stat /= 0 ) then
227 <       status = -1
228 <       return
229 <    endif
115 >    call check_initialization()
116 >    call zero_work_arrays()
117  
118 <    allocate(identCol(ncol),stat=alloc_stat)
119 <    if (alloc_stat /= 0 ) then
233 <       status = -1
234 <       return
235 <    endif
118 >    do_pot = do_pot_c
119 >    do_stress = do_stress_c
120  
121 < !! Gather idents into row and column idents
238 <
239 <    call gather(ident,identRow,plan_row)
240 <    call gather(ident,identCol,plan_col)
121 >    ! Gather all information needed by all force loops:
122      
123 <  
124 < !! Create row and col pointer lists
125 <  
126 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
127 <    if (thisStat /= 0 ) then
128 <       status = -1
129 <       return
123 > #ifdef IS_MPI    
124 >
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 <  
251 <    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
252 <    if (thisStat /= 0 ) then
253 <       status = -1
254 <       return
255 <    endif
256 <
257 < !! free temporary ident arrays
258 <    if (allocated(identCol)) then
259 <       deallocate(identCol)
260 <    end if
261 <    if (allocated(identCol)) then
262 <       deallocate(identRow)
263 <    endif
264 <
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
272 –
273 – !! Create neighbor lists
274 –    call expandList(thisStat)
275 –    if (thisStat /= 0) then
276 –       status = -1
277 –       return
278 –    endif
279 –
280 –    isFFinit = .true.
281 –
282 –
283 –  end subroutine init_FF
284 –
285 –
286 –
287 –
288 –  subroutine initForce_Modules(thisStat)
289 –    integer, intent(out) :: thisStat
290 –    integer :: my_status
149      
292 –    thisStat = 0
293 –    call init_lj_FF(ListHead,my_status)
294 –    if (my_status /= 0) then
295 –       thisStat = -1
296 –       return
297 –    end if
298 –
299 –  end subroutine initForce_Modules
300 –
301 –
302 –
303 –
304 – !! FORCE routine Calculates Lennard Jones forces.
305 – !------------------------------------------------------------->
306 –  subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
307 – !! Position array provided by C, dimensioned by getNlocal
308 –    real ( kind = dp ), dimension(3,getNlocal()) :: q
309 –  !! Rotation Matrix for each long range particle in simulation.
310 –    real( kind = dp), dimension(9,getNlocal()) :: A
311 –
312 –  !! Magnitude dipole moment
313 –    real( kind = dp ), dimension(3,getNlocal()) :: mu
314 –  !! Unit vectors for dipoles (lab frame)
315 –    real( kind = dp ), dimension(3,getNlocal()) :: u_l
316 – !! Force array provided by C, dimensioned by getNlocal
317 –    real ( kind = dp ), dimension(3,getNlocal()) :: f
318 – !! Torsion array provided by C, dimensioned by getNlocal
319 –    real( kind = dp ), dimension(3,getNlocal()) :: t
320 –
321 – !! Stress Tensor
322 –    real( kind = dp), dimension(9) :: tau
323 –
324 –    real ( kind = dp ) :: potE
325 –    logical ( kind = 2) :: do_pot
326 –    integer :: FFerror
327 –
328 –    
329 –    type(atype), pointer :: Atype_i
330 –    type(atype), pointer :: Atype_j
331 –
332 –
333 –
334 –
335 –  
336 –
150   #ifdef IS_MPI
338 –  real( kind = DP ) :: pot_local
339 –
340 – !! Local arrays needed for MPI
341 –
342 – #endif
343 –
344 –
345 –
346 –  real( kind = DP )   :: pe
347 –  logical             :: update_nlist
348 –
349 –
350 –  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
351 –  integer :: nlist
352 –  integer :: j_start
353 –
354 –  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
355 –
356 –  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
357 –  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
358 –
359 –  real( kind = DP ) :: dielectric = 0.0_dp
360 –
361 – ! a rig that need to be fixed.
362 – #ifdef IS_MPI
363 –  real( kind = dp ) :: pe_local
364 –  integer :: nlocal
365 – #endif
366 –  integer :: nrow
367 –  integer :: ncol
368 –  integer :: natoms
369 –  integer :: neighborListSize
370 –  integer :: listerror
371 – !! should we calculate the stress tensor
372 –  logical  :: do_stress = .false.
373 –
374 –
375 –  FFerror = 0
376 –
377 – ! Make sure we are properly initialized.
378 –  if (.not. isFFInit) then
379 –     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
380 –     FFerror = -1
381 –     return
382 –  endif
383 – #ifdef IS_MPI
384 –    if (.not. isMPISimSet()) then
385 –     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
386 –     FFerror = -1
387 –     return
388 –  endif
389 – #endif
390 –
391 – !! initialize local variables  
392 –  natoms = getNlocal()
393 –  call getRcut(rcut,rcut2=rcutsq)
394 –  call getRlist(rlist,rlistsq)
395 –
396 – !! Find ensemble
397 –  if (isEnsemble("NPT")) do_stress = .true.
398 – !! set to wrap
399 –  if (isPBC()) wrap = .true.
400 –
401 –
402 –
403 –  
404 – !! See if we need to update neighbor lists
405 –  call check(q,update_nlist)
406 –
407 – !--------------WARNING...........................
408 – ! Zero variables, NOTE:::: Forces are zeroed in C
409 – ! Zeroing them here could delete previously computed
410 – ! Forces.
411 – !------------------------------------------------
412 –  call zero_module_variables()
413 –
414 –
415 – ! communicate MPI positions
416 – #ifdef IS_MPI    
417 –    call gather(q,qRow,plan_row3d)
418 –    call gather(q,qCol,plan_col3d)
419 –
420 –    call gather(mu,muRow,plan_row3d)
421 –    call gather(mu,muCol,plan_col3d)
422 –
423 –    call gather(u_l,u_lRow,plan_row3d)
424 –    call gather(u_l,u_lCol,plan_col3d)
425 –
426 –    call gather(A,ARow,plan_row_rotation)
427 –    call gather(A,ACol,plan_col_rotation)
428 – #endif
429 –
430 –
431 – #ifdef IS_MPI
151      
152      if (update_nlist) then
153        
154 <       ! save current configuration, contruct neighbor list,
155 <       ! and calculate forces
154 >       !! save current configuration, construct neighbor list,
155 >       !! and calculate forces
156         call save_neighborList(q)
157        
158         neighborListSize = getNeighborListSize()
159 <       nlist = 0
159 >       nlist = 0      
160        
442 –       nrow = getNrow(plan_row)
443 –       ncol = getNcol(plan_col)
444 –       nlocal = getNlocal()
445 –      
161         do i = 1, nrow
162            point(i) = nlist + 1
448 –          Atype_i => identPtrListRow(i)%this
163            
164            inner: do j = 1, ncol
451 –             Atype_j => identPtrListColumn(j)%this
165              
166 <             call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
454 <                  rxij,ryij,rzij,rijsq,r)
166 >             if (check_exclude(i,j)) cycle inner:
167              
168 <             ! skip the loop if the atoms are identical
169 <             if (mpi_cycle_jLoop(i,j)) cycle inner:
458 <            
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 expandList(listerror)
175 >                   call expandNeighborList(nlocal, listerror)
176                     if (listerror /= 0) then
177 <                      FFerror = -1
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 <                
474 <                
184 >                                
185                  if (rijsq <  rcutsq) then
186 <                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
186 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
187                  endif
188               endif
189            enddo inner
# Line 481 | Line 191 | contains
191  
192         point(nrow + 1) = nlist + 1
193        
194 <    else !! (update)
194 >    else  !! (of update_check)
195  
196         ! use the list to find the neighbors
197         do i = 1, nrow
# Line 489 | Line 199 | contains
199            JEND = POINT(i+1) - 1
200            ! check thiat molecule i has neighbors
201            if (jbeg .le. jend) then
202 <
493 <             Atype_i => identPtrListRow(i)%this
202 >            
203               do jnab = jbeg, jend
204                  j = list(jnab)
205 <                Atype_j = identPtrListColumn(j)%this
206 <                call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
207 <                     rxij,ryij,rzij,rijsq,r)
208 <                
500 <                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
205 >
206 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
207 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
208 >
209               enddo
210            endif
211         enddo
212      endif
213 <
213 >    
214   #else
215      
216      if (update_nlist) then
# Line 514 | Line 222 | contains
222         neighborListSize = getNeighborListSize()
223         nlist = 0
224        
517 –    
225         do i = 1, natoms-1
226            point(i) = nlist + 1
227 <          Atype_i   => identPtrList(i)%this
521 <
227 >          
228            inner: do j = i+1, natoms
229 <             Atype_j   => identPtrList(j)%this
230 <             call get_interatomic_vector(i,j,q(:,i),q(:,j),&
231 <                  rxij,ryij,rzij,rijsq,r)
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 <
235 >                
236                  nlist = nlist + 1
237                  
238                  if (nlist > neighborListSize) then
239 <                   call expandList(listerror)
239 >                   call expandList(natoms, listerror)
240                     if (listerror /= 0) then
241 <                      FFerror = -1
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 <
542 <    
248 >                
249                  if (rijsq <  rcutsq) then
250 <                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
250 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
251                  endif
252               endif
253            enddo inner
# Line 558 | Line 264 | contains
264            ! check thiat molecule i has neighbors
265            if (jbeg .le. jend) then
266              
561 –             Atype_i => identPtrList(i)%this
267               do jnab = jbeg, jend
268                  j = list(jnab)
269 <                Atype_j = identPtrList(j)%this
270 <                call get_interatomic_vector(i,j,q(:,i),q(:,j),&
271 <                     rxij,ryij,rzij,rijsq,r)
272 <                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
269 >
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               enddo
274            endif
275         enddo
276      endif
277 <
277 >    
278   #endif
279 <
280 <
279 >    
280 >    ! phew, done with main loop.
281 >    
282   #ifdef IS_MPI
577 –    !! distribute all reaction field stuff (or anything for post-pair):
578 –    call scatter(rflRow,rflTemp1,plan_row3d)
579 –    call scatter(rflCol,rflTemp2,plan_col3d)
580 –    do i = 1,nlocal
581 –       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
582 –    end do
583 – #endif
584 –
585 – ! This is the post-pair loop:
586 – #ifdef IS_MPI
587 –
588 –    if (system_has_postpair_atoms) then
589 –       do i = 1, nlocal
590 –          Atype_i => identPtrListRow(i)%this
591 –          call do_postpair(i, Atype_i)
592 –       enddo
593 –    endif
594 –
595 – #else
596 –
597 –    if (system_has_postpair_atoms) then
598 –       do i = 1, natoms
599 –          Atype_i => identPtr(i)%this
600 –          call do_postpair(i, Atype_i)
601 –       enddo
602 –    endif
603 –
604 – #endif
605 –
606 –
607 –
608 –
609 – #ifdef IS_MPI
283      !!distribute forces
284 <
285 <    call scatter(fRow,fTemp1,plan_row3d)
286 <    call scatter(fCol,fTemp2,plan_col3d)
614 <
615 <
284 >    
285 >    call scatter(f_Row,f,plan_row3d)
286 >    call scatter(f_Col,f_temp,plan_col3d)
287      do i = 1,nlocal
288 <       fTemp(1:3,i) = fTemp1(1:3,i) + fTemp2(1:3,i)
288 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
289      end do
619 –
620 –    if (do_torque) then
621 –       call scatter(tRow,tTemp1,plan_row3d)
622 –       call scatter(tCol,tTemp2,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 <          tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(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        
645 –       pe = pe_local
317      endif
318 < #else
319 < ! Copy local array into return array for c
320 <    f = f+fTemp
321 <    t = 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  
653 –    potE = pe
654 –
655 –
656 –    if (do_stress) then
333   #ifdef IS_MPI
334 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
659 <            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  
665 –  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  
674 –
675 –
385   !! Calculate any pre-force loop components and update nlist if necessary.
386    subroutine do_preForce(updateNlist)
387      logical, intent(inout) :: updateNlist
# Line 681 | Line 390 | contains
390  
391    end subroutine do_preForce
392  
684 –
685 –
686 –
687 –
688 –
689 –
690 –
691 –
692 –
693 –
694 –
695 –
393   !! Calculate any post force loop components, i.e. reaction field, etc.
394    subroutine do_postForce()
395  
# Line 700 | 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 <
706 <
707 <
708 <
709 <
710 <
711 <
712 <
713 <
714 <
715 <
716 <
717 <
718 <
719 <  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
720 <    type (atype ), pointer, intent(inout) :: atype_i
721 <    type (atype ), pointer, intent(inout) :: atype_j
722 <    integer :: i
723 <    integer :: j
724 <    real ( kind = dp ), intent(inout) :: rx_ij
725 <    real ( kind = dp ), intent(inout) :: ry_ij
726 <    real ( kind = dp ), intent(inout) :: rz_ij
727 <
728 <
729 <    real( kind = dp ) :: fx = 0.0_dp
730 <    real( kind = dp ) :: fy = 0.0_dp
731 <    real( kind = dp ) :: fz = 0.0_dp  
732 <  
733 <    real( kind = dp ) ::  drdx = 0.0_dp
734 <    real( kind = dp ) ::  drdy = 0.0_dp
735 <    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 <    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
739 <       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
740 <    endif
415 > #ifdef IS_MPI
416  
417 <    if (Atype_i%is_dp .and. Atype_j%is_dp) then
417 >    me_i = atid_row(i)
418 >    me_j = atid_col(j)
419  
744 – #ifdef IS_MPI
745 –       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
746 –            ulRow(:,i), ulCol(:,j), rt, rrf, pot)
420   #else
748 –       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
749 –            ul(:,i), ul(:,j), rt, rrf, pot)
750 – #endif
421  
422 <       if (do_reaction_field) then
423 < #ifdef IS_MPI
424 <          call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
755 <               ulRow(:i), ulCol(:,j), rt, rrf)
756 < #else
757 <          call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
758 <               ul(:,i), ul(:,j), rt, rrf)
422 >    me_i = atid(i)
423 >    me_j = atid(j)
424 >
425   #endif
760 –       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 (Atype_i%is_sticky .and. Atype_j%is_sticky) then
438 <       call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
439 <    endif
768 <
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 < #ifdef IS_MPI
442 <                eRow(i) = eRow(i) + pot*0.5
443 <                eCol(i) = eCol(i) + pot*0.5
444 < #else
445 <                    pe = pe + pot
775 < #endif                
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 <                drdx = -rxij / r
448 <                drdy = -ryij / r
449 <                drdz = -rzij / r
450 <                
451 <                fx = dudr * drdx
452 <                fy = dudr * drdy
453 <                fz = dudr * drdz
784 <
785 <
786 <
787 <
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 +    if (FF_uses_Sticky .and. SimUsesSticky()) then
456  
457 <                
458 < #ifdef IS_MPI
459 <                fCol(1,j) = fCol(1,j) - fx
460 <                fCol(2,j) = fCol(2,j) - fy
461 <                fCol(3,j) = fCol(3,j) - fz
462 <                
463 <                fRow(1,j) = fRow(1,j) + fx
464 <                fRow(2,j) = fRow(2,j) + fy
465 <                fRow(3,j) = fRow(3,j) + fz
799 < #else
800 <                fTemp(1,j) = fTemp(1,j) - fx
801 <                fTemp(2,j) = fTemp(2,j) - fy
802 <                fTemp(3,j) = fTemp(3,j) - fz
803 <                fTemp(1,i) = fTemp(1,i) + fx
804 <                fTemp(2,i) = fTemp(2,i) + fy
805 <                fTemp(3,i) = fTemp(3,i) + fz
806 < #endif
807 <                
808 <                if (do_stress) then
809 <                   tauTemp(1) = tauTemp(1) + fx * rxij
810 <                   tauTemp(2) = tauTemp(2) + fx * ryij
811 <                   tauTemp(3) = tauTemp(3) + fx * rzij
812 <                   tauTemp(4) = tauTemp(4) + fy * rxij
813 <                   tauTemp(5) = tauTemp(5) + fy * ryij
814 <                   tauTemp(6) = tauTemp(6) + fy * rzij
815 <                   tauTemp(7) = tauTemp(7) + fz * rxij
816 <                   tauTemp(8) = tauTemp(8) + fz * ryij
817 <                   tauTemp(9) = tauTemp(9) + fz * rzij
818 <                endif
819 <
820 <
821 <
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 <
827 <
828 <
829 <
830 <
831 <
832 <
833 <
834 <
835 <
836 <
837 <
838 <
839 <  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
840 < !---------------- Arguments-------------------------------
841 <   !! index i
842 <
843 <    !! 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
846 –    !! x component of vector between i and j
847 –    real ( kind = dp ), intent(out)  :: rx_ij
848 –    !! y component of vector between i and j
849 –    real ( kind = dp ), intent(out)  :: ry_ij
850 –    !! z component of vector between i and j
851 –    real ( kind = dp ), intent(out)  :: rz_ij
852 –    !! magnitude of r squared
473      real ( kind = dp ), intent(out) :: r_sq
854 –    !! magnitude of vector r between atoms i and j.
855 –    real ( kind = dp ), intent(out) :: r_ij
856 –    !! wrap into periodic box.
857 –    logical, intent(in) :: wrap
858 –
859 – !--------------- Local Variables---------------------------
860 –    !! Distance between i and j
474      real( kind = dp ) :: d(3)
862 – !---------------- END DECLARATIONS-------------------------
475  
864 –
865 – ! 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      
874 – !   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)
880 <    ry_ij = d(2)
881 <    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 <
890 <    pe = 0.0E0_DP
891 <    tauTemp = 0.0_dp
892 <    fTemp = 0.0_dp
893 <    tTemp = 0.0_dp
894 < #else
895 <    qRow = 0.0_dp
896 <    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 <  
911 <    tRow = 0.0_dp
912 <    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  
916 –    eRow = 0.0_dp
917 –    eCol = 0.0_dp
918 –    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  
923 – #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 952 | Line 591 | contains
591      else                
592         if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
593      endif
955 –  end function mpi_cycle_jLoop
956 – #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