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 304 by gezelter, Mon Mar 10 18:56:37 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.7 2003-03-10 18:56:37 gezelter Exp $, $Date: 2003-03-10 18:56:37 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
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
317 –
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
68  
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  
437 –  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 –    
446 –    
447 –
448 –     do i = 1, nrow
449 –        point(i) = nlist + 1
158   #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
463 –
464 –        inner: do j = j_start, ncol
465 – #ifdef IS_MPI
466 – ! Assign identity pointers and tags
467 –           ljAtype_j => identPtrListColumn(j)%this
468 –           tag_j = tagColumn(j)
469 –           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
476 –
477 –           rxij = wrap(rxi - qCol(1,j), 1)
478 –           ryij = wrap(ryi - qCol(2,j), 2)
479 –           rzij = wrap(rzi - qCol(3,j), 3)
480 – #else          
481 –           ljAtype_j   => identPtrList(j)%this
482 –           rxij = wrap(rxi - q(1,j), 1)
483 –           ryij = wrap(ryi - q(2,j), 2)
484 –           rzij = wrap(rzi - q(3,j), 3)
485 –      
486 – #endif          
487 –           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
488 –
489 – #ifdef IS_MPI
490 –             if (rijsq <=  rlistsq .AND. &
491 –                  tag_j /= tag_i) then
492 – #else
493 –          
494 –             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
507 –
159      
160 <              if (rijsq <  rcutsq) then
510 <                
511 <                 r = dsqrt(rijsq)
160 >    if (update_nlist) then
161        
162 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
162 >       ! save current configuration, contruct neighbor list,
163 >       ! and calculate forces
164 >       call save_neighborList(q)
165        
166 < #ifdef IS_MPI
167 <                eRow(i) = eRow(i) + pot*0.5
168 <                eCol(i) = eCol(i) + pot*0.5
169 < #else
170 <                    pe = pe + pot
171 < #endif                
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 <                drdx = -rxij / r
181 <                drdy = -ryij / r
182 <                drdz = -rzij / r
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 <                fx = dudr * drdx
527 <                fy = dudr * drdy
528 <                fz = dudr * drdz
188 >                nlist = nlist + 1
189                  
190 < #ifdef IS_MPI
191 <                fCol(1,j) = fCol(1,j) - fx
192 <                fCol(2,j) = fCol(2,j) - fy
193 <                fCol(3,j) = fCol(3,j) - fz
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 <                fRow(1,j) = fRow(1,j) + fx
536 <                fRow(2,j) = fRow(2,j) + fy
537 <                fRow(3,j) = fRow(3,j) + fz
538 < #else
539 <                f(1,j) = f(1,j) - fx
540 <                f(2,j) = f(2,j) - fy
541 <                f(3,j) = f(3,j) - fz
542 <                f(1,i) = f(1,i) + fx
543 <                f(2,i) = f(2,i) + fy
544 <                f(3,i) = f(3,i) + fz
545 < #endif
199 >                list(nlist) = j
200                  
201 <                if (do_stress) then
202 <                   tauTemp(1) = tauTemp(1) + fx * rxij
203 <                   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
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
207 >       enddo
208  
209 < #ifdef IS_MPI
210 <     point(nrow + 1) = nlist + 1
211 < #else
565 <     point(natoms) = nlist + 1
566 < #endif
209 >       point(nrow + 1) = nlist + 1
210 >      
211 >    else !! (update)
212  
213 <  else
214 <
215 <     ! use the list to find the neighbors
216 <     do i = 1, nrow
217 <        JBEG = POINT(i)
218 <        JEND = POINT(i+1) - 1
574 <        ! check thiat molecule i has neighbors
575 <        if (jbeg .le. jend) then
576 < #ifdef IS_MPI
577 <           ljAtype_i => identPtrListRow(i)%this
578 <           rxi = qRow(1,i)
579 <           ryi = qRow(2,i)
580 <           rzi = qRow(3,i)
581 < #else
582 <           ljAtype_i   => identPtrList(i)%this
583 <           rxi = q(1,i)
584 <           ryi = q(2,i)
585 <           rzi = q(3,i)
586 < #endif
587 <           do jnab = jbeg, jend
588 <              j = list(jnab)
589 < #ifdef IS_MPI
590 <              ljAtype_j = identPtrListColumn(j)%this
591 <              rxij = wrap(rxi - qCol(1,j), 1)
592 <              ryij = wrap(ryi - qCol(2,j), 2)
593 <              rzij = wrap(rzi - qCol(3,j), 3)
594 < #else
595 <              ljAtype_j = identPtrList(j)%this
596 <              rxij = wrap(rxi - q(1,j), 1)
597 <              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
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 <                 r = dsqrt(rijsq)
221 <                
222 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
223 < #ifdef IS_MPI
224 <                eRow(i) = eRow(i) + pot*0.5
225 <                eCol(i) = eCol(i) + pot*0.5
610 < #else
611 <                pe = pe + pot
612 < #endif                
613 <  
614 <                drdx = -rxij / r
615 <                drdy = -ryij / r
616 <                drdz = -rzij / r
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 <                fx = dudr * drdx
228 <                fy = dudr * drdy
229 <                fz = dudr * drdz
230 <                
231 < #ifdef IS_MPI
232 <                fCol(1,j) = fCol(1,j) - fx
624 <                fCol(2,j) = fCol(2,j) - fy
625 <                fCol(3,j) = fCol(3,j) - fz
626 <                
627 <                fRow(1,j) = fRow(1,j) + fx
628 <                fRow(2,j) = fRow(2,j) + fy
629 <                fRow(3,j) = fRow(3,j) + fz
227 >                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
228 >             enddo
229 >          endif
230 >       enddo
231 >    endif
232 >
233   #else
234 <                f(1,j) = f(1,j) - fx
235 <                f(2,j) = f(2,j) - fy
236 <                f(3,j) = f(3,j) - fz
237 <                f(1,i) = f(1,i) + fx
238 <                f(2,i) = f(2,i) + fy
239 <                f(3,i) = f(3,i) + fz
240 < #endif
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
255 >
256 >                nlist = nlist + 1
257                  
258 <                if (do_stress) then
259 <                   tauTemp(1) = tauTemp(1) + fx * rxij
260 <                   tauTemp(2) = tauTemp(2) + fx * ryij
261 <                   tauTemp(3) = tauTemp(3) + fx * rzij
262 <                   tauTemp(4) = tauTemp(4) + fy * rxij
263 <                   tauTemp(5) = tauTemp(5) + fy * ryij
264 <                   tauTemp(6) = tauTemp(6) + fy * rzij
646 <                   tauTemp(7) = tauTemp(7) + fz * rxij
647 <                   tauTemp(8) = tauTemp(8) + fz * ryij
648 <                   tauTemp(9) = tauTemp(9) + fz * rzij
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 <                
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
275 <       endif
276 <    enddo
277 < endif
278 <
274 >          enddo inner
275 >       enddo
276 >      
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 +    !! 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 +
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 +
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
332 +
333 +
334 +
335 +
336 + #ifdef IS_MPI
337      !!distribute forces
338  
339 <    call scatter(fRow,f,plan_row3d)
340 <    call scatter(fCol,fMPITemp,plan_col3d)
339 >    call scatter(fRow,fTemp1,plan_row3d)
340 >    call scatter(fCol,fTemp2,plan_col3d)
341  
342 +
343      do i = 1,nlocal
344 <       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
344 >       fTemp(1:3,i) = fTemp1(1:3,i) + fTemp2(1:3,i)
345      end do
346  
347 <
348 <    call scatter(tRow,t,plan_row3d)
349 <    call scatter(tCol,tMPITemp,plan_col3d)
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 <       t(1:3,i) = t(1:3,i) + tMPITemp(1:3,i)
353 <    end do
351 >       do i = 1,nlocal
352 >          tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(1:3,i)
353 >       end do
354 >    endif
355  
678 –
356      if (do_pot) then
357         ! scatter/gather pot_row into the members of my column
358         call scatter(eRow,eTemp,plan_row)
# Line 685 | Line 362 | contains
362         do i = 1, nlocal
363            pe_local = pe_local + eTemp(i)
364         enddo
365 <       if (newtons_thrd) then
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 <       endif
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 <
374 > #else
375 > ! Copy local array into return array for c
376 >    f = f+fTemp
377 >    t = t+tTemp
378   #endif
379  
380      potE = pe
# Line 709 | Line 389 | contains
389   #endif      
390      endif
391  
712 –
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 + #ifdef IS_MPI
466 +
467 +    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
468 +       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
469 +    endif
470 +
471 +    if (Atype_i%is_dp .and. Atype_j%is_dp) then
472 +
473 +       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
474 +            ulRow(:,i), ulCol(:,j), rt, rrf, pot)
475 +
476 +       if (do_reaction_field) then
477 +          call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
478 +               ulRow(:i), ulCol(:,j), rt, rrf)
479 +       endif
480 +
481 +    endif
482 +
483 +    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
484 +       call getstickyforce(r, pot, dudr, Atype_i, Atype_j)
485 +    endif
486 +
487 + #else
488 +
489 +    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
490 +       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
491 +    endif
492 +
493 +    if (Atype_i%is_dp .and. Atype_j%is_dp) then
494 +       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
495 +            ul(:,i), ul(:,j), rt, rrf, pot)
496 +
497 +       if (do_reaction_field) then
498 +          call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
499 +               ul(:,i), ul(:,j), rt, rrf)
500 +       endif
501 +
502 +    endif
503 +
504 +    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
505 +       call getstickyforce(r,pot,dudr, Atype_i, Atype_j)
506 +    endif
507 +
508 + #endif
509 +
510 +      
511 + #ifdef IS_MPI
512 +                eRow(i) = eRow(i) + pot*0.5
513 +                eCol(i) = eCol(i) + pot*0.5
514 + #else
515 +                    pe = pe + pot
516 + #endif                
517 +            
518 +                drdx = -rxij / r
519 +                drdy = -ryij / r
520 +                drdz = -rzij / r
521 +                
522 +                fx = dudr * drdx
523 +                fy = dudr * drdy
524 +                fz = dudr * drdz
525 +                
526 + #ifdef IS_MPI
527 +                fCol(1,j) = fCol(1,j) - fx
528 +                fCol(2,j) = fCol(2,j) - fy
529 +                fCol(3,j) = fCol(3,j) - fz
530 +                
531 +                fRow(1,j) = fRow(1,j) + fx
532 +                fRow(2,j) = fRow(2,j) + fy
533 +                fRow(3,j) = fRow(3,j) + fz
534 + #else
535 +                fTemp(1,j) = fTemp(1,j) - fx
536 +                fTemp(2,j) = fTemp(2,j) - fy
537 +                fTemp(3,j) = fTemp(3,j) - fz
538 +                fTemp(1,i) = fTemp(1,i) + fx
539 +                fTemp(2,i) = fTemp(2,i) + fy
540 +                fTemp(3,i) = fTemp(3,i) + fz
541 + #endif
542 +                
543 +                if (do_stress) then
544 +                   tauTemp(1) = tauTemp(1) + fx * rxij
545 +                   tauTemp(2) = tauTemp(2) + fx * ryij
546 +                   tauTemp(3) = tauTemp(3) + fx * rzij
547 +                   tauTemp(4) = tauTemp(4) + fy * rxij
548 +                   tauTemp(5) = tauTemp(5) + fy * ryij
549 +                   tauTemp(6) = tauTemp(6) + fy * rzij
550 +                   tauTemp(7) = tauTemp(7) + fz * rxij
551 +                   tauTemp(8) = tauTemp(8) + fz * ryij
552 +                   tauTemp(9) = tauTemp(9) + fz * rzij
553 +                endif
554 +
555 +
556 +
557 +  end subroutine do_pair
558 +
559 +
560 +
561 +
562 +
563 +
564 +
565 +
566 +
567 +
568 +
569 +
570 +
571 +
572 +
573 +
574 +  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
575 + !---------------- Arguments-------------------------------
576 +   !! index i
577 +
578 +    !! Position array
579 +    real (kind = dp), dimension(3) :: q_i
580 +    real (kind = dp), dimension(3) :: q_j
581 +    !! x component of vector between i and j
582 +    real ( kind = dp ), intent(out)  :: rx_ij
583 +    !! y component of vector between i and j
584 +    real ( kind = dp ), intent(out)  :: ry_ij
585 +    !! z component of vector between i and j
586 +    real ( kind = dp ), intent(out)  :: rz_ij
587 +    !! magnitude of r squared
588 +    real ( kind = dp ), intent(out) :: r_sq
589 +    !! magnitude of vector r between atoms i and j.
590 +    real ( kind = dp ), intent(out) :: r_ij
591 +    !! wrap into periodic box.
592 +    logical, intent(in) :: wrap
593 +
594 + !--------------- Local Variables---------------------------
595 +    !! Distance between i and j
596 +    real( kind = dp ) :: d(3)
597 + !---------------- END DECLARATIONS-------------------------
598 +
599 +
600 + ! Find distance between i and j
601 +    d(1:3) = q_i(1:3) - q_j(1:3)
602 +
603 + ! Wrap back into periodic box if necessary
604 +    if ( wrap ) then
605 +       d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
606 +            int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
607 +    end if
608 +    
609 + !   Find Magnitude of the vector
610 +    r_sq = dot_product(d,d)
611 +    r_ij = sqrt(r_sq)
612 +
613 + !   Set each component for force calculation
614 +    rx_ij = d(1)
615 +    ry_ij = d(2)
616 +    rz_ij = d(3)
617 +
618 +
619 +  end subroutine get_interatomic_vector
620 +
621 +  subroutine zero_module_variables()
622 +
623 + #ifndef IS_MPI
624 +
625 +    pe = 0.0E0_DP
626 +    tauTemp = 0.0_dp
627 +    fTemp = 0.0_dp
628 +    tTemp = 0.0_dp
629 + #else
630 +    qRow = 0.0_dp
631 +    qCol = 0.0_dp
632 +    
633 +    muRow = 0.0_dp
634 +    muCol = 0.0_dp
635 +    
636 +    u_lRow = 0.0_dp
637 +    u_lCol = 0.0_dp
638 +    
639 +    ARow = 0.0_dp
640 +    ACol = 0.0_dp
641 +    
642 +    fRow = 0.0_dp
643 +    fCol = 0.0_dp
644 +    
645 +  
646 +    tRow = 0.0_dp
647 +    tCol = 0.0_dp
648 +
649 +  
650 +
651 +    eRow = 0.0_dp
652 +    eCol = 0.0_dp
653 +    eTemp = 0.0_dp
654 + #endif
655 +
656 +  end subroutine zero_module_variables
657 +
658 + #ifdef IS_MPI
659 + !! Function to properly build neighbor lists in MPI using newtons 3rd law.
660 + !! We don't want 2 processors doing the same i j pair twice.
661 + !! Also checks to see if i and j are the same particle.
662 +  function mpi_cycle_jLoop(i,j) result(do_cycle)
663 + !--------------- Arguments--------------------------
664 + ! Index i
665 +    integer,intent(in) :: i
666 + ! Index j
667 +    integer,intent(in) :: j
668 + ! Result do_cycle
669 +    logical :: do_cycle
670 + !--------------- Local variables--------------------
671 +    integer :: tag_i
672 +    integer :: tag_j
673 + !--------------- END DECLARATIONS------------------    
674 +    tag_i = tagRow(i)
675 +    tag_j = tagColumn(j)
676 +
677 +    do_cycle = .false.
678 +
679 +    if (tag_i == tag_j) then
680 +       do_cycle = .true.
681 +       return
682 +    end if
683 +
684 +    if (tag_i < tag_j) then
685 +       if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
686 +       return
687 +    else                
688 +       if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
689 +    endif
690 +  end function mpi_cycle_jLoop
691 + #endif
692 +
693   end module do_Forces

Diff Legend

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