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 292 by chuckv, Thu Mar 6 14:52:44 2003 UTC vs.
Revision 298 by chuckv, Fri Mar 7 18:26:30 2003 UTC

# Line 1 | Line 1
1   !! Calculates Long Range forces.
2   !! @author Charles F. Vardeman II
3   !! @author Matthew Meineke
4 < !! @version $Id: do_Forces.F90,v 1.1 2003-03-06 14:52:44 chuckv Exp $, $Date: 2003-03-06 14:52:44 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
4 > !! @version $Id: do_Forces.F90,v 1.5 2003-03-07 18:26:30 chuckv Exp $, $Date: 2003-03-07 18:26:30 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $
5  
6  
7  
8   module do_Forces
9    use simulation
10    use definitions
11 <  use generic_atypes
11 >  use forceGlobals
12 >  use atype_typedefs
13    use neighborLists
14 +
15    
16    use lj_FF
17    use sticky_FF
# Line 22 | Line 24 | module do_Forces
24    implicit none
25    PRIVATE
26  
27 < !! Number of lj_atypes in lj_atype_list
26 <  integer, save :: n_atypes = 0
27 > public :: do_force_loop
28  
28 – !! Global list of lj atypes in simulation
29 –  type (atype), pointer :: ListHead => null()
30 –  type (atype), pointer :: ListTail => null()
31 –
32 –
33 – !! LJ mixing array  
34 – !  type (lj_atype), dimension(:,:), pointer :: ljMixed => null()
35 –
36 –
37 –  logical, save :: firstTime = .True.
38 –
39 – !! Atype identity pointer lists
40 – #ifdef IS_MPI
41 – !! Row lj_atype pointer list
42 –  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
43 – !! Column lj_atype pointer list
44 –  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
45 – #else
46 –  type(identPtrList ), dimension(:), pointer :: identPtrList => null()
47 – #endif
48 –
49 –
50 – !! Logical has lj force field module been initialized?
51 –  logical, save :: isFFinit = .false.
52 –
53 –
54 – !! Public methods and data
55 –  public :: new_atype
56 –  public :: do_lj_ff
57 –  public :: getLjPot
58 –  public :: init_FF
59 –
60 –  
61 –
62 –
29   contains
30  
65 – !! Adds a new lj_atype to the list.
66 –  subroutine new_atype(ident,mass,epsilon,sigma, &
67 –       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
68 –    real( kind = dp ), intent(in) :: mass
69 –    real( kind = dp ), intent(in) :: epsilon
70 –    real( kind = dp ), intent(in) :: sigma
71 –    real( kind = dp ), intent(in) :: w0
72 –    real( kind = dp ), intent(in) :: v0
73 –    real( kind = dp ), intent(in) :: dipoleMoment
74 –
75 –    integer, intent(in) :: ident
76 –    integer, intent(out) :: status
77 –    integer, intent(in) :: is_Sticky
78 –    integer, intent(in) :: is_DP
79 –    integer, intent(in) :: is_GB
80 –    integer, intent(in) :: is_LJ
81 –
82 –
83 –    type (atype), pointer :: the_new_atype
84 –    integer :: alloc_error
85 –    integer :: atype_counter = 0
86 –    integer :: alloc_size
87 –    integer :: err_stat
88 –    status = 0
89 –
90 –
91 –
92 – ! allocate a new atype    
93 –    allocate(the_new_atype,stat=alloc_error)
94 –    if (alloc_error /= 0 ) then
95 –       status = -1
96 –       return
97 –    end if
98 –
99 – ! assign our new atype information
100 –    the_new_atype%mass        = mass
101 –    the_new_atype%epsilon     = epsilon
102 –    the_new_atype%sigma       = sigma
103 –    the_new_atype%sigma2      = sigma * sigma
104 –    the_new_atype%sigma6      = the_new_atype%sigma2 * the_new_atype%sigma2 &
105 –         * the_new_atype%sigma2
106 –    the_new_atype%w0       = w0
107 –    the_new_atype%v0       = v0
108 –    the_new_atype%dipoleMoment       = dipoleMoment
109 –
110 –    
111 – ! assume that this atype will be successfully added
112 –    the_new_atype%atype_ident = ident
113 –    the_new_atype%atype_number = n_lj_atypes + 1
114 –    
115 –    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
116 –    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
117 –    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
118 –    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
119 –
120 –    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
121 –    if (err_stat /= 0 ) then
122 –       status = -1
123 –       return
124 –    endif
125 –
126 –    n_atypes = n_atypes + 1
127 –
128 –
129 –  end subroutine new_atype
130 –
131 –
132 –  subroutine init_FF(nComponents,ident, status)
133 – !! Number of components in ident array
134 –    integer, intent(inout) :: nComponents
135 – !! Array of identities nComponents long corresponding to
136 – !! ljatype ident.
137 –    integer, dimension(nComponents),intent(inout) :: ident
138 – !!  Result status, success = 0, error = -1
139 –    integer, intent(out) :: Status
140 –
141 –    integer :: alloc_stat
142 –
143 –    integer :: thisStat
144 –    integer :: i
145 –
146 –    integer :: myNode
147 – #ifdef IS_MPI
148 –    integer, allocatable, dimension(:) :: identRow
149 –    integer, allocatable, dimension(:) :: identCol
150 –    integer :: nrow
151 –    integer :: ncol
152 – #endif
153 –    status = 0
154 –  
155 –
156 –    
157 –
158 – !! if were're not in MPI, we just update ljatypePtrList
159 – #ifndef IS_MPI
160 –    call create_IdentPtrlst(ident,ListHead,identPtrList,thisStat)
161 –    if ( thisStat /= 0 ) then
162 –       status = -1
163 –       return
164 –    endif
165 –
166 –    
167 – ! if were're in MPI, we also have to worry about row and col lists    
168 – #else
169 –  
170 – ! We can only set up forces if mpiSimulation has been setup.
171 –    if (.not. isMPISimSet()) then
172 –       write(default_error,*) "MPI is not set"
173 –       status = -1
174 –       return
175 –    endif
176 –    nrow = getNrow(plan_row)
177 –    ncol = getNcol(plan_col)
178 –    mynode = getMyNode()
179 – !! Allocate temperary arrays to hold gather information
180 –    allocate(identRow(nrow),stat=alloc_stat)
181 –    if (alloc_stat /= 0 ) then
182 –       status = -1
183 –       return
184 –    endif
185 –
186 –    allocate(identCol(ncol),stat=alloc_stat)
187 –    if (alloc_stat /= 0 ) then
188 –       status = -1
189 –       return
190 –    endif
191 –
192 – !! Gather idents into row and column idents
193 –
194 –    call gather(ident,identRow,plan_row)
195 –    call gather(ident,identCol,plan_col)
196 –    
197 –  
198 – !! Create row and col pointer lists
199 –  
200 –    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
201 –    if (thisStat /= 0 ) then
202 –       status = -1
203 –       return
204 –    endif
205 –  
206 –    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
207 –    if (thisStat /= 0 ) then
208 –       status = -1
209 –       return
210 –    endif
211 –
212 – !! free temporary ident arrays
213 –    if (allocated(identCol)) then
214 –       deallocate(identCol)
215 –    end if
216 –    if (allocated(identCol)) then
217 –       deallocate(identRow)
218 –    endif
219 –
220 – #endif
221 –    
222 –    call initForce_Modules(thisStat)
223 –    if (thisStat /= 0) then
224 –       status = -1
225 –       return
226 –    endif
227 –
228 – !! Create neighbor lists
229 –    call expandList(thisStat)
230 –    if (thisStat /= 0) then
231 –       status = -1
232 –       return
233 –    endif
234 –
235 –    isFFinit = .true.
236 –
237 –
238 –  end subroutine init_FF
239 –
240 –
241 –
242 –
243 –  subroutine initForce_Modules(thisStat)
244 –    integer, intent(out) :: thisStat
245 –    integer :: my_status
246 –    
247 –    thisStat = 0
248 –    call init_lj_FF(ListHead,my_status)
249 –    if (my_status /= 0) then
250 –       thisStat = -1
251 –       return
252 –    end if
253 –
254 –  end subroutine initForce_Modules
255 –
256 –
257 –
258 –
31   !! FORCE routine Calculates Lennard Jones forces.
32   !------------------------------------------------------------->
33 <  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
33 >  subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
34   !! Position array provided by C, dimensioned by getNlocal
35      real ( kind = dp ), dimension(3,getNlocal()) :: q
36    !! Rotation Matrix for each long range particle in simulation.
# Line 275 | Line 47 | contains
47  
48   !! Stress Tensor
49      real( kind = dp), dimension(9) :: tau
50 <    real( kind = dp), dimension(9) :: tauTemp
50 >
51      real ( kind = dp ) :: potE
52      logical ( kind = 2) :: do_pot
53      integer :: FFerror
# Line 293 | Line 65 | contains
65    real( kind = DP ) :: pot_local
66  
67   !! Local arrays needed for MPI
296 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
297 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
298 –
299 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
300 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
301 –
302 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
303 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
304 –
305 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: ARow = 0.0_dp
306 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp
307 –
308 –  
309 –
310 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
311 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
312 –  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp = 0.0_dp
313 –
314 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
315 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
316 –  real(kind = dp), dimension(3,getNlocal()) :: tMPITemp = 0.0_dp
68  
318 –
319 –  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
320 –  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
321 –
322 –  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
323 –
69   #endif
70  
71  
# Line 332 | Line 77 | contains
77    integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
78    integer :: nlist
79    integer :: j_start
80 <  integer :: tag_i,tag_j
81 <  real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
82 <  real( kind = dp ) :: fx,fy,fz
83 <  real( kind = DP ) ::  drdx, drdy, drdz
339 <  real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
80 >
81 >  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
82 >
83 >  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
84    real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
85  
86    real( kind = DP ) :: dielectric = 0.0_dp
87  
88   ! a rig that need to be fixed.
89   #ifdef IS_MPI
346 –  logical :: newtons_thrd = .true.
90    real( kind = dp ) :: pe_local
91    integer :: nlocal
92   #endif
# Line 376 | Line 119 | contains
119    natoms = getNlocal()
120    call getRcut(rcut,rcut2=rcutsq)
121    call getRlist(rlist,rlistsq)
122 +
123   !! Find ensemble
124    if (isEnsemble("NPT")) do_stress = .true.
125 + !! set to wrap
126 +  if (isPBC()) wrap = .true.
127  
382 – #ifndef IS_MPI
383 –  nrow = natoms - 1
384 –  ncol = natoms
385 – #else
386 –  nrow = getNrow(plan_row)
387 –  ncol = getNcol(plan_col)
388 –  nlocal = natoms
389 –  j_start = 1
390 – #endif
128  
129 +
130    
131   !! See if we need to update neighbor lists
132    call check(q,update_nlist)
395 – !  if (firstTime) then
396 – !     update_nlist = .true.
397 – !     firstTime = .false.
398 – !  endif
133  
134   !--------------WARNING...........................
135   ! Zero variables, NOTE:::: Forces are zeroed in C
136   ! Zeroing them here could delete previously computed
137   ! Forces.
138   !------------------------------------------------
139 < #ifndef IS_MPI
406 < !  nloops = nloops + 1
407 <  pe = 0.0E0_DP
408 <
409 < #else
410 <    fRow = 0.0E0_DP
411 <    fCol = 0.0E0_DP
139 >  call zero_module_variables()
140  
413 –    pe_local = 0.0E0_DP
141  
415 –    eRow = 0.0E0_DP
416 –    eCol = 0.0E0_DP
417 –    eTemp = 0.0E0_DP
418 – #endif
419 –
142   ! communicate MPI positions
143   #ifdef IS_MPI    
144      call gather(q,qRow,plan_row3d)
# Line 430 | Line 152 | contains
152  
153      call gather(A,ARow,plan_row_rotation)
154      call gather(A,ACol,plan_col_rotation)
433 –
155   #endif
156  
157  
158 <  if (update_nlist) then
438 <
439 <     ! save current configuration, contruct neighbor list,
440 <     ! and calculate forces
441 <     call save_neighborList(q)
442 <    
443 <     neighborListSize = getNeighborListSize()
444 <     nlist = 0
445 <    
158 > #ifdef IS_MPI
159      
160 +    if (update_nlist) then
161 +      
162 +       ! save current configuration, contruct neighbor list,
163 +       ! and calculate forces
164 +       call save_neighborList(q)
165 +      
166 +       neighborListSize = getNeighborListSize()
167 +       nlist = 0
168 +      
169 +       nrow = getNrow(plan_row)
170 +       ncol = getNcol(plan_col)
171 +       nlocal = getNlocal()
172 +      
173 +       do i = 1, nrow
174 +          point(i) = nlist + 1
175 +          Atype_i => identPtrListRow(i)%this
176 +          
177 +          inner: do j = 1, ncol
178 +             Atype_j => identPtrListColumn(j)%this
179 +            
180 +             call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
181 +                  rxij,ryij,rzij,rijsq,r)
182 +            
183 +             ! skip the loop if the atoms are identical
184 +             if (mpi_cycle_jLoop(i,j)) cycle inner:
185 +            
186 +             if (rijsq <  rlistsq) then            
187 +                
188 +                nlist = nlist + 1
189 +                
190 +                if (nlist > neighborListSize) then
191 +                   call expandList(listerror)
192 +                   if (listerror /= 0) then
193 +                      FFerror = -1
194 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
195 +                      return
196 +                   end if
197 +                endif
198 +                
199 +                list(nlist) = j
200 +                
201 +                
202 +                if (rijsq <  rcutsq) then
203 +                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
204 +                endif
205 +             endif
206 +          enddo inner
207 +       enddo
208  
209 <     do i = 1, nrow
210 <        point(i) = nlist + 1
211 < #ifdef IS_MPI
451 <        ljAtype_i => identPtrListRow(i)%this
452 <        tag_i = tagRow(i)
453 <        rxi = qRow(1,i)
454 <        ryi = qRow(2,i)
455 <        rzi = qRow(3,i)
456 < #else
457 <        ljAtype_i   => identPtrList(i)%this
458 <        j_start = i + 1
459 <        rxi = q(1,i)
460 <        ryi = q(2,i)
461 <        rzi = q(3,i)
462 < #endif
209 >       point(nrow + 1) = nlist + 1
210 >      
211 >    else !! (update)
212  
213 <        inner: do j = j_start, ncol
214 < #ifdef IS_MPI
215 < ! Assign identity pointers and tags
216 <           ljAtype_j => identPtrListColumn(j)%this
217 <           tag_j = tagColumn(j)
218 <           if (newtons_thrd) then
470 <              if (tag_i <= tag_j) then
471 <                 if (mod(tag_i + tag_j,2) == 0) cycle inner
472 <              else                
473 <                 if (mod(tag_i + tag_j,2) == 1) cycle inner
474 <              endif
475 <           endif
213 >       ! use the list to find the neighbors
214 >       do i = 1, nrow
215 >          JBEG = POINT(i)
216 >          JEND = POINT(i+1) - 1
217 >          ! check thiat molecule i has neighbors
218 >          if (jbeg .le. jend) then
219  
220 <           rxij = wrap(rxi - qCol(1,j), 1)
221 <           ryij = wrap(ryi - qCol(2,j), 2)
222 <           rzij = wrap(rzi - qCol(3,j), 3)
223 < #else          
224 <           ljAtype_j   => identPtrList(j)%this
225 <           rxij = wrap(rxi - q(1,j), 1)
226 <           ryij = wrap(ryi - q(2,j), 2)
227 <           rzij = wrap(rzi - q(3,j), 3)
228 <      
229 < #endif          
230 <           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
220 >             Atype_i => identPtrListRow(i)%this
221 >             do jnab = jbeg, jend
222 >                j = list(jnab)
223 >                Atype_j = identPtrListColumn(j)%this
224 >                call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
225 >                     rxij,ryij,rzij,rijsq,r)
226 >                
227 >                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
228 >             enddo
229 >          endif
230 >       enddo
231 >    endif
232  
489 – #ifdef IS_MPI
490 –             if (rijsq <=  rlistsq .AND. &
491 –                  tag_j /= tag_i) then
233   #else
234 +    
235 +    if (update_nlist) then
236 +      
237 +       ! save current configuration, contruct neighbor list,
238 +       ! and calculate forces
239 +       call save_neighborList(q)
240 +      
241 +       neighborListSize = getNeighborListSize()
242 +       nlist = 0
243 +      
244 +    
245 +       do i = 1, natoms-1
246 +          point(i) = nlist + 1
247 +          Atype_i   => identPtrList(i)%this
248 +
249 +          inner: do j = i+1, natoms
250 +             Atype_j   => identPtrList(j)%this
251 +             call get_interatomic_vector(i,j,q(:,i),q(:,j),&
252 +                  rxij,ryij,rzij,rijsq,r)
253            
254               if (rijsq <  rlistsq) then
495 – #endif
496 –            
497 –              nlist = nlist + 1
498 –              if (nlist > neighborListSize) then
499 –                 call expandList(listerror)
500 –                 if (listerror /= 0) then
501 –                    FFerror = -1
502 –                    write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
503 –                    return
504 –                 end if
505 –              endif
506 –              list(nlist) = j
255  
256 <    
509 <              if (rijsq <  rcutsq) then
256 >                nlist = nlist + 1
257                  
258 <                 r = dsqrt(rijsq)
258 >                if (nlist > neighborListSize) then
259 >                   call expandList(listerror)
260 >                   if (listerror /= 0) then
261 >                      FFerror = -1
262 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
263 >                      return
264 >                   end if
265 >                endif
266 >                
267 >                list(nlist) = j
268 >
269 >    
270 >                if (rijsq <  rcutsq) then
271 >                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
272 >                endif
273 >             endif
274 >          enddo inner
275 >       enddo
276        
277 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
278 <      
277 >       point(natoms) = nlist + 1
278 >      
279 >    else !! (update)
280 >      
281 >       ! use the list to find the neighbors
282 >       do i = 1, nrow
283 >          JBEG = POINT(i)
284 >          JEND = POINT(i+1) - 1
285 >          ! check thiat molecule i has neighbors
286 >          if (jbeg .le. jend) then
287 >            
288 >             Atype_i => identPtrList(i)%this
289 >             do jnab = jbeg, jend
290 >                j = list(jnab)
291 >                Atype_j = identPtrList(j)%this
292 >                call get_interatomic_vector(i,j,q(:,i),q(:,j),&
293 >                     rxij,ryij,rzij,rijsq,r)
294 >                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
295 >             enddo
296 >          endif
297 >       enddo
298 >    endif
299 >
300 > #endif
301 >
302 >
303   #ifdef IS_MPI
304 <                eRow(i) = eRow(i) + pot*0.5
305 <                eCol(i) = eCol(i) + pot*0.5
306 < #else
307 <                    pe = pe + pot
308 < #endif                
309 <            
310 <                drdx = -rxij / r
311 <                drdy = -ryij / r
312 <                drdz = -rzij / r
525 <                
526 <                fx = dudr * drdx
527 <                fy = dudr * drdy
528 <                fz = dudr * drdz
529 <                
304 >    !! distribute all reaction field stuff (or anything for post-pair):
305 >    call scatter(rflRow,rflTemp1,plan_row3d)
306 >    call scatter(rflCol,rflTemp2,plan_col3d)
307 >    do i = 1,nlocal
308 >       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
309 >    end do
310 > #endif
311 >
312 > ! This is the post-pair loop:
313   #ifdef IS_MPI
314 <                fCol(1,j) = fCol(1,j) - fx
315 <                fCol(2,j) = fCol(2,j) - fy
316 <                fCol(3,j) = fCol(3,j) - fz
317 <                
318 <                fRow(1,j) = fRow(1,j) + fx
319 <                fRow(2,j) = fRow(2,j) + fy
320 <                fRow(3,j) = fRow(3,j) + fz
314 >
315 >    if (system_has_postpair_atoms) then
316 >       do i = 1, nlocal
317 >          Atype_i => identPtrListRow(i)%this
318 >          call do_postpair(i, Atype_i)
319 >       enddo
320 >    endif
321 >
322   #else
323 <                f(1,j) = f(1,j) - fx
324 <                f(2,j) = f(2,j) - fy
325 <                f(3,j) = f(3,j) - fz
326 <                f(1,i) = f(1,i) + fx
327 <                f(2,i) = f(2,i) + fy
328 <                f(3,i) = f(3,i) + fz
323 >
324 >    if (system_has_postpair_atoms) then
325 >       do i = 1, natoms
326 >          Atype_i => identPtr(i)%this
327 >          call do_postpair(i, Atype_i)
328 >       enddo
329 >    endif
330 >
331   #endif
546 –                
547 –                if (do_stress) then
548 –                   tauTemp(1) = tauTemp(1) + fx * rxij
549 –                   tauTemp(2) = tauTemp(2) + fx * ryij
550 –                   tauTemp(3) = tauTemp(3) + fx * rzij
551 –                   tauTemp(4) = tauTemp(4) + fy * rxij
552 –                   tauTemp(5) = tauTemp(5) + fy * ryij
553 –                   tauTemp(6) = tauTemp(6) + fy * rzij
554 –                   tauTemp(7) = tauTemp(7) + fz * rxij
555 –                   tauTemp(8) = tauTemp(8) + fz * ryij
556 –                   tauTemp(9) = tauTemp(9) + fz * rzij
557 –                endif
558 –             endif
559 –          enddo inner
560 –     enddo
332  
333 +
334 +
335 +
336   #ifdef IS_MPI
337 <     point(nrow + 1) = nlist + 1
337 >    !!distribute forces
338 >
339 >    call scatter(fRow,fTemp1,plan_row3d)
340 >    call scatter(fCol,fTemp2,plan_col3d)
341 >
342 >
343 >    do i = 1,nlocal
344 >       fTemp(1:3,i) = fTemp1(1:3,i) + fTemp2(1:3,i)
345 >    end do
346 >
347 >    if (do_torque) then
348 >       call scatter(tRow,tTemp1,plan_row3d)
349 >       call scatter(tCol,tTemp2,plan_col3d)
350 >    
351 >       do i = 1,nlocal
352 >          tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(1:3,i)
353 >       end do
354 >    endif
355 >
356 >    if (do_pot) then
357 >       ! scatter/gather pot_row into the members of my column
358 >       call scatter(eRow,eTemp,plan_row)
359 >      
360 >       ! scatter/gather pot_local into all other procs
361 >       ! add resultant to get total pot
362 >       do i = 1, nlocal
363 >          pe_local = pe_local + eTemp(i)
364 >       enddo
365 >
366 >       eTemp = 0.0E0_DP
367 >       call scatter(eCol,eTemp,plan_col)
368 >       do i = 1, nlocal
369 >          pe_local = pe_local + eTemp(i)
370 >       enddo
371 >      
372 >       pe = pe_local
373 >    endif
374   #else
375 <     point(natoms) = nlist + 1
375 > ! Copy local array into return array for c
376 >    f = f+fTemp
377 >    t = t+tTemp
378   #endif
379  
380 <  else
380 >    potE = pe
381  
382 <     ! use the list to find the neighbors
383 <     do i = 1, nrow
572 <        JBEG = POINT(i)
573 <        JEND = POINT(i+1) - 1
574 <        ! check thiat molecule i has neighbors
575 <        if (jbeg .le. jend) then
382 >
383 >    if (do_stress) then
384   #ifdef IS_MPI
385 <           ljAtype_i => identPtrListRow(i)%this
386 <           rxi = qRow(1,i)
579 <           ryi = qRow(2,i)
580 <           rzi = qRow(3,i)
385 >       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
386 >            mpi_comm_world,mpi_err)
387   #else
388 <           ljAtype_i   => identPtrList(i)%this
389 <           rxi = q(1,i)
390 <           ryi = q(2,i)
391 <           rzi = q(3,i)
392 < #endif
393 <           do jnab = jbeg, jend
394 <              j = list(jnab)
388 >       tau = tauTemp
389 > #endif      
390 >    endif
391 >
392 >  end subroutine do_force_loop
393 >
394 >
395 >
396 >
397 >
398 >
399 >
400 >
401 >
402 >
403 > !! Calculate any pre-force loop components and update nlist if necessary.
404 >  subroutine do_preForce(updateNlist)
405 >    logical, intent(inout) :: updateNlist
406 >
407 >
408 >
409 >  end subroutine do_preForce
410 >
411 >
412 >
413 >
414 >
415 >
416 >
417 >
418 >
419 >
420 >
421 >
422 >
423 > !! Calculate any post force loop components, i.e. reaction field, etc.
424 >  subroutine do_postForce()
425 >
426 >
427 >
428 >  end subroutine do_postForce
429 >
430 >
431 >
432 >
433 >
434 >
435 >
436 >
437 >
438 >
439 >
440 >
441 >
442 >
443 >
444 >
445 >
446 >  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
447 >    type (atype ), pointer, intent(inout) :: atype_i
448 >    type (atype ), pointer, intent(inout) :: atype_j
449 >    integer :: i
450 >    integer :: j
451 >    real ( kind = dp ), intent(inout) :: rx_ij
452 >    real ( kind = dp ), intent(inout) :: ry_ij
453 >    real ( kind = dp ), intent(inout) :: rz_ij
454 >
455 >
456 >    real( kind = dp ) :: fx = 0.0_dp
457 >    real( kind = dp ) :: fy = 0.0_dp
458 >    real( kind = dp ) :: fz = 0.0_dp  
459 >  
460 >    real( kind = dp ) ::  drdx = 0.0_dp
461 >    real( kind = dp ) ::  drdy = 0.0_dp
462 >    real( kind = dp ) ::  drdz = 0.0_dp
463 >    
464 >
465 >    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
466 >       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
467 >    endif
468 >
469 >    if (Atype_i%is_dp .and. Atype_j%is_dp) then
470 >
471   #ifdef IS_MPI
472 <              ljAtype_j = identPtrListColumn(j)%this
473 <              rxij = wrap(rxi - qCol(1,j), 1)
592 <              ryij = wrap(ryi - qCol(2,j), 2)
593 <              rzij = wrap(rzi - qCol(3,j), 3)
472 >       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
473 >            ulRow(:,i), ulCol(:,j), rt, rrf, pot)
474   #else
475 <              ljAtype_j = identPtrList(j)%this
476 <              rxij = wrap(rxi - q(1,j), 1)
477 <              ryij = wrap(ryi - q(2,j), 2)
598 <              rzij = wrap(rzi - q(3,j), 3)
599 < #endif
600 <              rijsq = rxij*rxij + ryij*ryij + rzij*rzij
601 <              
602 <              if (rijsq <  rcutsq) then
475 >       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
476 >            ul(:,i), ul(:,j), rt, rrf, pot)
477 > #endif
478  
479 <                 r = dsqrt(rijsq)
605 <                
606 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
479 >       if (do_reaction_field) then
480   #ifdef IS_MPI
481 +          call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
482 +               ulRow(:i), ulCol(:,j), rt, rrf)
483 + #else
484 +          call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
485 +               ul(:,i), ul(:,j), rt, rrf)
486 + #endif
487 +       endif
488 +
489 +
490 +    endif
491 +
492 +    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
493 +       call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
494 +    endif
495 +
496 +      
497 + #ifdef IS_MPI
498                  eRow(i) = eRow(i) + pot*0.5
499                  eCol(i) = eCol(i) + pot*0.5
500   #else
501 <                pe = pe + pot
501 >                    pe = pe + pot
502   #endif                
503 <  
503 >            
504                  drdx = -rxij / r
505                  drdy = -ryij / r
506                  drdz = -rzij / r
# Line 618 | Line 508 | contains
508                  fx = dudr * drdx
509                  fy = dudr * drdy
510                  fz = dudr * drdz
511 +
512 +
513 +
514 +
515 +
516 +
517                  
518   #ifdef IS_MPI
519                  fCol(1,j) = fCol(1,j) - fx
# Line 628 | Line 524 | contains
524                  fRow(2,j) = fRow(2,j) + fy
525                  fRow(3,j) = fRow(3,j) + fz
526   #else
527 <                f(1,j) = f(1,j) - fx
528 <                f(2,j) = f(2,j) - fy
529 <                f(3,j) = f(3,j) - fz
530 <                f(1,i) = f(1,i) + fx
531 <                f(2,i) = f(2,i) + fy
532 <                f(3,i) = f(3,i) + fz
527 >                fTemp(1,j) = fTemp(1,j) - fx
528 >                fTemp(2,j) = fTemp(2,j) - fy
529 >                fTemp(3,j) = fTemp(3,j) - fz
530 >                fTemp(1,i) = fTemp(1,i) + fx
531 >                fTemp(2,i) = fTemp(2,i) + fy
532 >                fTemp(3,i) = fTemp(3,i) + fz
533   #endif
534                  
535                  if (do_stress) then
# Line 647 | Line 543 | contains
543                     tauTemp(8) = tauTemp(8) + fz * ryij
544                     tauTemp(9) = tauTemp(9) + fz * rzij
545                  endif
650 –                
651 –                
652 –             endif
653 –          enddo
654 –       endif
655 –    enddo
656 – endif
657 –
546  
547  
660 – #ifdef IS_MPI
661 –    !!distribute forces
548  
549 <    call scatter(fRow,f,plan_row3d)
664 <    call scatter(fCol,fMPITemp,plan_col3d)
549 >  end subroutine do_pair
550  
666 –    do i = 1,nlocal
667 –       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
668 –    end do
551  
552  
553 <    call scatter(tRow,t,plan_row3d)
554 <    call scatter(tCol,tMPITemp,plan_col3d)
553 >
554 >
555 >
556 >
557 >
558 >
559 >
560 >
561 >
562 >
563 >
564 >
565 >
566 >  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
567 > !---------------- Arguments-------------------------------
568 >   !! index i
569 >
570 >    !! Position array
571 >    real (kind = dp), dimension(3) :: q_i
572 >    real (kind = dp), dimension(3) :: q_j
573 >    !! x component of vector between i and j
574 >    real ( kind = dp ), intent(out)  :: rx_ij
575 >    !! y component of vector between i and j
576 >    real ( kind = dp ), intent(out)  :: ry_ij
577 >    !! z component of vector between i and j
578 >    real ( kind = dp ), intent(out)  :: rz_ij
579 >    !! magnitude of r squared
580 >    real ( kind = dp ), intent(out) :: r_sq
581 >    !! magnitude of vector r between atoms i and j.
582 >    real ( kind = dp ), intent(out) :: r_ij
583 >    !! wrap into periodic box.
584 >    logical, intent(in) :: wrap
585 >
586 > !--------------- Local Variables---------------------------
587 >    !! Distance between i and j
588 >    real( kind = dp ) :: d(3)
589 > !---------------- END DECLARATIONS-------------------------
590 >
591 >
592 > ! Find distance between i and j
593 >    d(1:3) = q_i(1:3) - q_j(1:3)
594 >
595 > ! Wrap back into periodic box if necessary
596 >    if ( wrap ) then
597 >       d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
598 >            int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
599 >    end if
600      
601 <    do i = 1,nlocal
602 <       t(1:3,i) = t(1:3,i) + tMPITemp(1:3,i)
603 <    end do
601 > !   Find Magnitude of the vector
602 >    r_sq = dot_product(d,d)
603 >    r_ij = sqrt(r_sq)
604  
605 + !   Set each component for force calculation
606 +    rx_ij = d(1)
607 +    ry_ij = d(2)
608 +    rz_ij = d(3)
609  
679 –    if (do_pot) then
680 –       ! scatter/gather pot_row into the members of my column
681 –       call scatter(eRow,eTemp,plan_row)
682 –      
683 –       ! scatter/gather pot_local into all other procs
684 –       ! add resultant to get total pot
685 –       do i = 1, nlocal
686 –          pe_local = pe_local + eTemp(i)
687 –       enddo
688 –       if (newtons_thrd) then
689 –          eTemp = 0.0E0_DP
690 –          call scatter(eCol,eTemp,plan_col)
691 –          do i = 1, nlocal
692 –             pe_local = pe_local + eTemp(i)
693 –          enddo
694 –       endif
695 –       pe = pe_local
696 –    endif
610  
611 < #endif
611 >  end subroutine get_interatomic_vector
612  
613 <    potE = pe
613 >  subroutine zero_module_variables()
614  
615 + #ifndef IS_MPI
616  
617 <    if (do_stress) then
618 < #ifdef IS_MPI
619 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
620 <            mpi_comm_world,mpi_err)
617 >    pe = 0.0E0_DP
618 >    tauTemp = 0.0_dp
619 >    fTemp = 0.0_dp
620 >    tTemp = 0.0_dp
621   #else
622 <       tau = tauTemp
623 < #endif      
624 <    endif
622 >    qRow = 0.0_dp
623 >    qCol = 0.0_dp
624 >    
625 >    muRow = 0.0_dp
626 >    muCol = 0.0_dp
627 >    
628 >    u_lRow = 0.0_dp
629 >    u_lCol = 0.0_dp
630 >    
631 >    ARow = 0.0_dp
632 >    ACol = 0.0_dp
633 >    
634 >    fRow = 0.0_dp
635 >    fCol = 0.0_dp
636 >    
637 >  
638 >    tRow = 0.0_dp
639 >    tCol = 0.0_dp
640 >
641 >  
642  
643 +    eRow = 0.0_dp
644 +    eCol = 0.0_dp
645 +    eTemp = 0.0_dp
646 + #endif
647  
648 <  end subroutine do_force_loop
648 >  end subroutine zero_module_variables
649  
650 + #ifdef IS_MPI
651 + !! Function to properly build neighbor lists in MPI using newtons 3rd law.
652 + !! We don't want 2 processors doing the same i j pair twice.
653 + !! Also checks to see if i and j are the same particle.
654 +  function mpi_cycle_jLoop(i,j) result(do_cycle)
655 + !--------------- Arguments--------------------------
656 + ! Index i
657 +    integer,intent(in) :: i
658 + ! Index j
659 +    integer,intent(in) :: j
660 + ! Result do_cycle
661 +    logical :: do_cycle
662 + !--------------- Local variables--------------------
663 +    integer :: tag_i
664 +    integer :: tag_j
665 + !--------------- END DECLARATIONS------------------    
666 +    tag_i = tagRow(i)
667 +    tag_j = tagColumn(j)
668  
669 +    do_cycle = .false.
670  
671 +    if (tag_i == tag_j) then
672 +       do_cycle = .true.
673 +       return
674 +    end if
675 +
676 +    if (tag_i < tag_j) then
677 +       if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
678 +       return
679 +    else                
680 +       if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
681 +    endif
682 +  end function mpi_cycle_jLoop
683 + #endif
684 +
685   end module do_Forces

Diff Legend

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