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 657 by chuckv, Wed Jul 30 21:17:01 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.26 2003-07-30 21:17:01 chuckv Exp $, $Date: 2003-07-30 21:17:01 $, $Name: not supported by cvs2svn $, $Revision: 1.26 $
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 >    write(*,*) "Inside do_Force Loop"
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 256 | Line 271 | contains
271         !! do_prepair_loop_if_needed
272         !! if_mpi_scatter_stuff_from_prepair
273         !! if_mpi_gather_stuff_from_prepair_to_main_loop
274 +
275 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
276 + #ifdef IS_MPI
277 +    
278 +    if (update_nlist) then
279 +      
280 +       !! save current configuration, construct neighbor list,
281 +       !! and calculate forces
282 +       call saveNeighborList(nlocal, q)
283 +      
284 +       neighborListSize = size(list)
285 +       nlist = 0      
286 +      
287 +       do i = 1, nrow
288 +          point(i) = nlist + 1
289 +          
290 +          prepair_inner: do j = 1, ncol
291 +            
292 +             if (skipThisPair(i,j)) cycle prepair_inner
293 +            
294 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
295 +            
296 +             if (rijsq < rlistsq) then            
297 +                
298 +                nlist = nlist + 1
299 +                
300 +                if (nlist > neighborListSize) then
301 +                   call expandNeighborList(nlocal, listerror)
302 +                   if (listerror /= 0) then
303 +                      error = -1
304 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
305 +                      return
306 +                   end if
307 +                   neighborListSize = size(list)
308 +                endif
309 +                
310 +                list(nlist) = j
311 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
312 +             endif
313 +          enddo prepair_inner
314 +       enddo
315 +
316 +       point(nrow + 1) = nlist + 1
317 +      
318 +    else  !! (of update_check)
319 +
320 +       ! use the list to find the neighbors
321 +       do i = 1, nrow
322 +          JBEG = POINT(i)
323 +          JEND = POINT(i+1) - 1
324 +          ! check thiat molecule i has neighbors
325 +          if (jbeg .le. jend) then
326 +            
327 +             do jnab = jbeg, jend
328 +                j = list(jnab)
329 +
330 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
331 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
332 +                     u_l, A, f, t, pot_local)
333 +
334 +             enddo
335 +          endif
336 +       enddo
337 +    endif
338 +    
339 + #else
340 +    
341 +    if (update_nlist) then
342 +      
343 +       ! save current configuration, contruct neighbor list,
344 +       ! and calculate forces
345 +       call saveNeighborList(natoms, q)
346 +      
347 +       neighborListSize = size(list)
348 +  
349 +       nlist = 0
350 +      
351 +       do i = 1, natoms-1
352 +          point(i) = nlist + 1
353 +          
354 +          prepair_inner: do j = i+1, natoms
355 +            
356 +             if (skipThisPair(i,j))  cycle prepair_inner
357 +                          
358 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
359 +          
360 +
361 +             if (rijsq < rlistsq) then
362 +                
363 +                nlist = nlist + 1
364 +              
365 +                if (nlist > neighborListSize) then
366 +                   call expandNeighborList(natoms, listerror)
367 +                   if (listerror /= 0) then
368 +                      error = -1
369 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
370 +                      return
371 +                   end if
372 +                   neighborListSize = size(list)
373 +                endif
374 +                
375 +                list(nlist) = j
376 +                
377 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
378 +                        u_l, A, f, t, pot)
379 +                
380 +             endif
381 +          enddo prepair_inner
382 +       enddo
383 +      
384 +       point(natoms) = nlist + 1
385 +      
386 +    else !! (update)
387 +      
388 +       ! use the list to find the neighbors
389 +       do i = 1, natoms-1
390 +          JBEG = POINT(i)
391 +          JEND = POINT(i+1) - 1
392 +          ! check thiat molecule i has neighbors
393 +          if (jbeg .le. jend) then
394 +            
395 +             do jnab = jbeg, jend
396 +                j = list(jnab)
397 +
398 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
399 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
400 +                     u_l, A, f, t, pot)
401 +
402 +             enddo
403 +          endif
404 +       enddo
405 +    endif    
406 + #endif
407 +    !! Do rest of preforce calculations
408 +   call do_preforce(nlocal,pot)
409      else
410 <       !! See if we need to update neighbor lists
410 >       !! See if we need to update neighbor lists for non pre-pair
411         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
412      endif
413 <    
413 >
414 >
415 >
416 >
417 >
418 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
419 >
420 >
421 >
422 >
423 >  
424   #ifdef IS_MPI
425      
426      if (update_nlist) then
# Line 531 | Line 691 | contains
691      logical :: is_LJ_i, is_LJ_j
692      logical :: is_DP_i, is_DP_j
693      logical :: is_GB_i, is_GB_j
694 +    logical :: is_EAM_i,is_EAM_j
695      logical :: is_Sticky_i, is_Sticky_j
696      integer :: me_i, me_j
697  
# Line 598 | Line 759 | contains
759         endif
760      endif
761      
762 +
763 +  
764 +   if (FF_uses_EAM .and. SimUsesEAM()) then
765 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
766 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
767 +      
768 +      if ( is_EAM_i .and. is_EAM_j ) &
769 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
770 +   endif
771 +
772  
773  
774 +
775    end subroutine do_pair
776  
777  
778  
779 <  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
779 >  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
780     real( kind = dp ) :: pot
781     real( kind = dp ), dimension(3,getNlocal()) :: u_l
782     real (kind=dp), dimension(9,getNlocal()) :: A
# Line 643 | Line 815 | contains
815        call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
816        
817        if ( is_EAM_i .and. is_EAM_j ) &
818 <           call calc_EAM_prepair(i, j, d, r, rijsq )
818 >           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
819 >   endif
820 >  end subroutine do_prepair
821 >
822 >
823 >
824 >
825 >  subroutine do_preforce(nlocal,pot)
826 >    integer :: nlocal
827 >    real( kind = dp ) :: pot
828 >
829 >   if (FF_uses_EAM .and. SimUsesEAM()) then
830 >      call calc_EAM_preforce_Frho(nlocal,pot)
831     endif
832 +
833 +
834    end subroutine do_preforce
835    
836    
# Line 705 | Line 891 | contains
891      error = 0
892      ! Make sure we are properly initialized.
893      if (.not. do_forces_initialized) then
894 +       write(*,*) "Forces not initialized"
895         error = -1
896         return
897      endif

Diff Legend

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