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 328 by gezelter, Wed Mar 12 20:00:58 2003 UTC vs.
Revision 329 by gezelter, Wed Mar 12 22:27:59 2003 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.14 2003-03-12 20:00:58 gezelter Exp $, $Date: 2003-03-12 20:00:58 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $
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 forceGlobals
14    use atype_module
15    use neighborLists  
16    use lj_FF
# Line 25 | Line 24 | module do_Forces
24    implicit none
25    PRIVATE
26  
27 <  public :: do_force_loop
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  
30  logical :: do_pot
31  logical :: do_stress
35  
36 +  public :: init_FF
37 +  public :: do_forces
38 +
39   contains
40  
41 +  subroutine init_FF(thisStat)
42 +  
43 +    integer, intent(out) :: my_status
44 +    integer :: thisStat = 0
45 +
46 +    ! be a smarter subroutine.
47 +
48 +    
49 +    call init_lj_FF(my_status)
50 +    if (my_status /= 0) then
51 +       thisStat = -1
52 +       return
53 +    end if
54 +    
55 +    call check_sticky_FF(my_status)
56 +    if (my_status /= 0) then
57 +       thisStat = -1
58 +       return
59 +    end if
60 +    
61 +    do_forces_initialized = .true.    
62 +    
63 +  end subroutine init_FF
64 +
65 +
66 +
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 <       FFerror)
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.
# Line 50 | Line 82 | contains
82      real( kind = dp), dimension(9) :: tau  
83      real ( kind = dp ) :: pot
84      logical ( kind = 2) :: do_pot_c, do_stress_c
85 <    integer :: FFerror
86 <
85 >    logical :: do_pot
86 >    logical :: do_stress
87   #ifdef IS_MPI
88      real( kind = DP ) :: pot_local
89 +    integer :: nlocal
90 +    integer :: nrow
91 +    integer :: ncol
92   #endif
93 <    
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
63
64 #ifdef IS_MPI
65    integer :: nlocal
66 #endif
67    integer :: nrow
68    integer :: ncol
69    integer :: natoms
98      integer :: neighborListSize
99      integer :: listerror
72    FFerror = 0
100  
101 <    do_pot = do_pot_c
75 <    do_stress = do_stress_c
101 >    !! initialize local variables  
102  
77    ! Make sure we are properly initialized.
78    if (.not. isFFInit) then
79       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
80       FFerror = -1
81       return
82    endif
103   #ifdef IS_MPI
104 <    if (.not. isMPISimSet()) then
105 <       write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
106 <       FFerror = -1
107 <       return
108 <    endif
104 >    nlocal = getNlocal()
105 >    nrow   = getNrow(plan_row)
106 >    ncol   = getNcol(plan_col)
107 > #else
108 >    nlocal = getNlocal()
109 >    natoms = nlocal
110   #endif
111 <    
91 <    !! initialize local variables  
92 <    natoms = getNlocal()
111 >
112      call getRcut(rcut,rcut2=rcutsq)
113      call getRlist(rlist,rlistsq)
114      
115 <    !! See if we need to update neighbor lists
116 <    call checkNeighborList(natoms, q, rcut, rlist, update_nlist)
115 >    call check_initialization()
116 >    call zero_work_arrays()
117 >
118 >    do_pot = do_pot_c
119 >    do_stress = do_stress_c
120 >
121 >    ! Gather all information needed by all force loops:
122      
123 <    !--------------WARNING...........................
100 <    ! Zero variables, NOTE:::: Forces are zeroed in C
101 <    ! Zeroing them here could delete previously computed
102 <    ! Forces.
103 <    !------------------------------------------------
104 <    call zero_module_variables()
123 > #ifdef IS_MPI    
124  
106    ! communicate MPI positions
107 #ifdef IS_MPI    
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      
111    call gather(u_l,u_l_Row,plan_row3d)
112    call gather(u_l,u_l_Col,plan_col3d)
113    
114    call gather(A,A_Row,plan_row_rotation)
115    call gather(A,A_Col,plan_col_rotation)
136   #endif
137 <
138 <
137 >    
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
149 >    
150   #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        
130       nrow = getNrow(plan_row)
131       ncol = getNcol(plan_col)
132       nlocal = getNlocal()
133      
161         do i = 1, nrow
162            point(i) = nlist + 1
163            
164            inner: do j = 1, ncol
165              
166               if (check_exclude(i,j)) cycle inner:
167 <
167 >            
168               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
169              
170               if (rijsq <  rlistsq) then            
# Line 147 | Line 174 | contains
174                  if (nlist > neighborListSize) then
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
# Line 156 | Line 183 | contains
183                  list(nlist) = j
184                                  
185                  if (rijsq <  rcutsq) then
186 <                   call do_pair(i, j, rijsq, d)
186 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
187                  endif
188               endif
189            enddo inner
# Line 164 | 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 177 | Line 204 | contains
204                  j = list(jnab)
205  
206                  call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
207 <                call do_pair(i, j, rijsq, d)
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 194 | Line 221 | contains
221        
222         neighborListSize = getNeighborListSize()
223         nlist = 0
224 <          
224 >      
225         do i = 1, natoms-1
226            point(i) = nlist + 1
227 <
227 >          
228            inner: do j = i+1, natoms
229 <
229 >            
230               if (check_exclude(i,j)) cycle inner:
231 <
231 >            
232               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
233            
234               if (rijsq <  rlistsq) then
# Line 211 | Line 238 | contains
238                  if (nlist > neighborListSize) then
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 <                    
248 >                
249                  if (rijsq <  rcutsq) then
250 <                   call do_pair(i, j, rijsq, d)
250 >                   call do_pair(i, j, rijsq, d, do_pot, do_stress)
251                  endif
252               endif
253            enddo inner
# Line 241 | Line 268 | contains
268                  j = list(jnab)
269  
270                  call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
271 <                call do_pair(i, j, rijsq, d)
271 >                call do_pair(i, j, rijsq, d, do_pot, do_stress)
272  
273               enddo
274            endif
# Line 250 | Line 277 | contains
277      
278   #endif
279      
280 <    
254 < #ifdef IS_MPI
255 <    !! distribute all reaction field stuff (or anything for post-pair):
256 <    call scatter(rflRow,rflTemp1,plan_row3d)
257 <    call scatter(rflCol,rflTemp2,plan_col3d)
258 <    do i = 1,nlocal
259 <       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
260 <    end do
261 < #endif
262 <    
263 < ! This is the post-pair loop:
264 < #ifdef IS_MPI
265 <    
266 <    if (system_has_postpair_atoms) then
267 <       do i = 1, nlocal
268 <          Atype_i => identPtrListRow(i)%this
269 <          call do_postpair(i, Atype_i)
270 <       enddo
271 <    endif
280 >    ! phew, done with main loop.
281      
273 #else
274    
275    if (system_has_postpair_atoms) then
276       do i = 1, natoms
277          Atype_i => identPtr(i)%this
278          call do_postpair(i, Atype_i)
279       enddo
280    endif
281    
282 #endif
283    
284
282   #ifdef IS_MPI
283      !!distribute forces
284 <
284 >    
285      call scatter(f_Row,f,plan_row3d)
286      call scatter(f_Col,f_temp,plan_col3d)
287      do i = 1,nlocal
288         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
289      end do
290 <
291 <    if (doTorque()) then
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 <    
294 >      
295         do i = 1,nlocal
296            t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
297         end do
# Line 317 | Line 314 | contains
314            pot_local = pot_local + pot_Temp(i)
315         enddo
316        
317 +    endif
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 +
333 + #ifdef IS_MPI
334 +             me_i = atid_row(i)
335 + #else
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 +
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 (doStress()) then
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 < #endif
373 > #else
374  
375 <    if (doStress()) then
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  
# Line 352 | Line 397 | contains
397  
398    end subroutine do_postForce
399  
400 <  subroutine do_pair(i, j, rijsq, d)
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
# Line 378 | Line 424 | contains
424  
425   #endif
426  
381    call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
382    call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
427  
428 <    if ( is_LJ_i .and. is_LJ_j ) &
429 <         call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
430 <
431 <    call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
432 <    call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
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 ( is_DP_i .and. is_DP_j ) then
438 <
439 <       call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
440 <
441 <       if (do_reaction_field) then
442 <          call accumulate_rf(i, j, r)
437 >    if (FF_uses_DP .and. SimUsesDP()) then
438 >       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
439 >       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
440 >      
441 >       if ( is_DP_i .and. is_DP_j ) then
442 >          
443 >          call do_dipole_pair(i, j, d, r, pot, u_l, f, t, do_pot, do_stress)
444 >          
445 >          if (FF_uses_RF .and. SimUsesRF()) then
446 >            
447 >             call accumulate_rf(i, j, r, u_l)
448 >             call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
449 >            
450 >          endif
451 >          
452         endif
397
453      endif
454  
455 <    call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
401 <    call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
455 >    if (FF_uses_Sticky .and. SimUsesSticky()) then
456  
457 <    if ( is_Sticky_i .and. is_Sticky_j ) then
458 <       call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, do_pot, do_stress)
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
406
465        
466    end subroutine do_pair
467  
# Line 426 | Line 484 | contains
484      r_sq = dot_product(d,d)
485          
486    end subroutine get_interatomic_vector
487 +
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 +    return
507 +  end subroutine check_initialization
508 +
509    
510    subroutine zero_work_arrays()
511      
# Line 524 | Line 604 | contains
604  
605    end function checkExcludes
606  
607 <
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