ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 631 by chuckv, Thu Jul 17 19:25:51 2003 UTC vs.
Revision 673 by chuckv, Fri Aug 8 21:22:37 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.22 2003-07-17 19:25:51 chuckv Exp $, $Date: 2003-07-17 19:25:51 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $
7 > !! @version $Id: do_Forces.F90,v 1.28 2003-08-08 21:22:37 chuckv Exp $, $Date: 2003-08-08 21:22:37 $, $Name: not supported by cvs2svn $, $Revision: 1.28 $
8  
9   module do_Forces
10    use force_globals
# Line 18 | Line 18 | module do_Forces
18    use reaction_field
19    use gb_pair
20    use vector_class
21 +  use eam
22 +  use status
23   #ifdef IS_MPI
24    use mpiSimulation
25   #endif
# Line 141 | Line 143 | contains
143            return
144         end if
145      endif
146 +
147 +
148 +    if (FF_uses_EAM) then
149 +       call init_EAM_FF(my_status)
150 +       if (my_status /= 0) then
151 +          thisStat = -1
152 +          return
153 +       end if
154 +    endif
155 +
156 +
157      
158      if (FF_uses_GB) then
159         call check_gb_pair_FF(my_status)
# Line 161 | Line 174 | contains
174            return
175         endif
176      endif
177 +    
178  
179      havePolicies = .true.
180      if( haveRlist ) do_forces_initialized = .true.
181 <    
181 >
182    end subroutine init_FF
183    
184  
# Line 221 | Line 235 | contains
235      nlocal = getNlocal()
236      natoms = nlocal
237   #endif
238 <  
238 >
239      call check_initialization(localError)
240      if ( localError .ne. 0 ) then
241 +       call handleError("do_force_loop","Not Initialized")
242         error = -1
243         return
244      end if
# Line 232 | Line 247 | contains
247      do_pot = do_pot_c
248      do_stress = do_stress_c
249  
250 +
251      ! Gather all information needed by all force loops:
252      
253   #ifdef IS_MPI    
# Line 248 | Line 264 | contains
264      endif
265      
266   #endif
267 <    
267 >  
268      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
269         !! See if we need to update neighbor lists
270         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
# Line 256 | Line 272 | contains
272         !! do_prepair_loop_if_needed
273         !! if_mpi_scatter_stuff_from_prepair
274         !! if_mpi_gather_stuff_from_prepair_to_main_loop
275 +    
276 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
277 + #ifdef IS_MPI
278 +    
279 +    if (update_nlist) then
280 +      
281 +       !! save current configuration, construct neighbor list,
282 +       !! and calculate forces
283 +       call saveNeighborList(nlocal, q)
284 +      
285 +       neighborListSize = size(list)
286 +       nlist = 0      
287 +      
288 +       do i = 1, nrow
289 +          point(i) = nlist + 1
290 +          
291 +          prepair_inner: do j = 1, ncol
292 +            
293 +             if (skipThisPair(i,j)) cycle prepair_inner
294 +            
295 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
296 +            
297 +             if (rijsq < rlistsq) then            
298 +                
299 +                nlist = nlist + 1
300 +                
301 +                if (nlist > neighborListSize) then
302 +                   call expandNeighborList(nlocal, listerror)
303 +                   if (listerror /= 0) then
304 +                      error = -1
305 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
306 +                      return
307 +                   end if
308 +                   neighborListSize = size(list)
309 +                endif
310 +                
311 +                list(nlist) = j
312 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
313 +             endif
314 +          enddo prepair_inner
315 +       enddo
316 +
317 +       point(nrow + 1) = nlist + 1
318 +      
319 +    else  !! (of update_check)
320 +
321 +       ! use the list to find the neighbors
322 +       do i = 1, nrow
323 +          JBEG = POINT(i)
324 +          JEND = POINT(i+1) - 1
325 +          ! check thiat molecule i has neighbors
326 +          if (jbeg .le. jend) then
327 +            
328 +             do jnab = jbeg, jend
329 +                j = list(jnab)
330 +
331 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
332 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
333 +                     u_l, A, f, t, pot_local)
334 +
335 +             enddo
336 +          endif
337 +       enddo
338 +    endif
339 +    
340 + #else
341 +    
342 +    if (update_nlist) then
343 +      
344 +       ! save current configuration, contruct neighbor list,
345 +       ! and calculate forces
346 +       call saveNeighborList(natoms, q)
347 +      
348 +       neighborListSize = size(list)
349 +  
350 +       nlist = 0
351 +
352 +       do i = 1, natoms-1
353 +          point(i) = nlist + 1
354 +          
355 +          prepair_inner: do j = i+1, natoms
356 +            
357 +             if (skipThisPair(i,j))  cycle prepair_inner
358 +                          
359 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
360 +          
361 +
362 +             if (rijsq < rlistsq) then
363 +
364 +          
365 +                nlist = nlist + 1
366 +              
367 +                if (nlist > neighborListSize) then
368 +                   call expandNeighborList(natoms, listerror)
369 +                   if (listerror /= 0) then
370 +                      error = -1
371 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
372 +                      return
373 +                   end if
374 +                   neighborListSize = size(list)
375 +                endif
376 +                
377 +                list(nlist) = j
378 +                
379 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
380 +                        u_l, A, f, t, pot)
381 +                
382 +             endif
383 +          enddo prepair_inner
384 +       enddo
385 +      
386 +       point(natoms) = nlist + 1
387 +      
388 +    else !! (update)
389 +  
390 +       ! use the list to find the neighbors
391 +       do i = 1, natoms-1
392 +          JBEG = POINT(i)
393 +          JEND = POINT(i+1) - 1
394 +          ! check thiat molecule i has neighbors
395 +          if (jbeg .le. jend) then
396 +            
397 +             do jnab = jbeg, jend
398 +                j = list(jnab)
399 +
400 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
401 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
402 +                     u_l, A, f, t, pot)
403 +
404 +             enddo
405 +          endif
406 +       enddo
407 +    endif    
408 + #endif
409 +    !! Do rest of preforce calculations
410 +    !! do necessary preforce calculations  
411 +    call do_preforce(nlocal,pot)
412 +   ! we have already updated the neighbor list set it to false...
413 +   update_nlist = .false.
414      else
415 <       !! See if we need to update neighbor lists
415 >       !! See if we need to update neighbor lists for non pre-pair
416         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
417      endif
418 <    
418 >
419 >
420 >
421 >
422 >
423 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
424 >
425 >
426 >
427 >
428 >  
429   #ifdef IS_MPI
430      
431      if (update_nlist) then
# Line 531 | Line 696 | contains
696      logical :: is_LJ_i, is_LJ_j
697      logical :: is_DP_i, is_DP_j
698      logical :: is_GB_i, is_GB_j
699 +    logical :: is_EAM_i,is_EAM_j
700      logical :: is_Sticky_i, is_Sticky_j
701      integer :: me_i, me_j
702  
# Line 598 | Line 764 | contains
764         endif
765      endif
766      
767 +
768 +  
769 +   if (FF_uses_EAM .and. SimUsesEAM()) then
770 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
771 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
772 +      
773 +      if ( is_EAM_i .and. is_EAM_j ) &
774 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
775 +   endif
776 +
777  
778  
779 +
780    end subroutine do_pair
781  
782  
783  
784 <  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
784 >  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
785     real( kind = dp ) :: pot
786     real( kind = dp ), dimension(3,getNlocal()) :: u_l
787     real (kind=dp), dimension(9,getNlocal()) :: A
# Line 623 | Line 800 | contains
800    
801     r = sqrt(rijsq)
802    
803 +
804   #ifdef IS_MPI
805     if (tagRow(i) .eq. tagColumn(j)) then
806        write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 637 | Line 815 | contains
815     me_j = atid(j)
816    
817   #endif
818 <  
818 >    
819     if (FF_uses_EAM .and. SimUsesEAM()) then
820        call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
821        call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
822        
823        if ( is_EAM_i .and. is_EAM_j ) &
824 <           call calc_EAM_prepair(i, j, d, r, rijsq )
824 >           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
825     endif
826 +
827 + end subroutine do_prepair
828 +
829 +
830 +
831 +
832 +  subroutine do_preforce(nlocal,pot)
833 +    integer :: nlocal
834 +    real( kind = dp ) :: pot
835 +
836 +    if (FF_uses_EAM .and. SimUsesEAM()) then
837 +       call calc_EAM_preforce_Frho(nlocal,pot)
838 +    endif
839 +
840 +
841    end subroutine do_preforce
842    
843    
# Line 705 | Line 898 | contains
898      error = 0
899      ! Make sure we are properly initialized.
900      if (.not. do_forces_initialized) then
901 +       write(*,*) "Forces not initialized"
902         error = -1
903         return
904      endif
# Line 752 | Line 946 | contains
946  
947   #endif
948  
949 +
950 +    if (FF_uses_EAM .and. SimUsesEAM()) then
951 +       call clean_EAM()
952 +    endif
953 +
954 +
955 +
956 +
957 +
958      rf = 0.0_dp
959      tau_Temp = 0.0_dp
960      virial_Temp = 0.0_dp
# Line 864 | Line 1067 | end module do_Forces
1067      doesit = FF_uses_RF
1068    end function FF_RequiresPostpairCalc
1069    
1070 + !! This cleans componets of force arrays belonging only to fortran
1071 +
1072   end module do_Forces

Diff Legend

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