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 829 by gezelter, Tue Oct 28 16:03: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.32 2003-10-28 16:03:35 gezelter Exp $, $Date: 2003-10-28 16:03:35 $, $Name: not supported by cvs2svn $, $Revision: 1.32 $
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 43 | Line 45 | contains
45    public :: do_force_loop
46    public :: setRlistDF
47  
48 + #ifdef PROFILE
49 +  real(kind = dp) :: forceTime
50 +  real(kind = dp) :: forceTimeInitial, forceTimeFinal
51 +  real(kind = dp) :: globalForceTime
52 +  real(kind = dp) :: maxForceTime
53 +  integer, save :: nloops = 0
54 + #endif
55 +
56   contains
57  
58    subroutine setRlistDF( this_rlist )
# Line 119 | Line 129 | contains
129      if (FF_uses_LJ) then
130        
131         select case (LJMIXPOLICY)
132 <       case (LB_MIXING_RULE)
133 <          call init_lj_FF(LB_MIXING_RULE, my_status)            
134 <       case (EXPLICIT_MIXING_RULE)
135 <          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
132 >       case (1)
133 >          call init_lj_FF(1, my_status)            
134 >       case (2)
135 >          call init_lj_FF(2, my_status)
136         case default
137            write(default_error,*) 'unknown LJ Mixing Policy!'
138            thisStat = -1
# Line 141 | Line 151 | contains
151            return
152         end if
153      endif
154 +
155 +
156 +    if (FF_uses_EAM) then
157 +         call init_EAM_FF(my_status)
158 +       if (my_status /= 0) then
159 +          write(*,*) "init_EAM_FF returned a bad status"
160 +          thisStat = -1
161 +          return
162 +       end if
163 +    endif
164 +
165 +
166      
167      if (FF_uses_GB) then
168         call check_gb_pair_FF(my_status)
# Line 161 | Line 183 | contains
183            return
184         endif
185      endif
186 +    
187  
188      havePolicies = .true.
189      if( haveRlist ) do_forces_initialized = .true.
190 <    
190 >
191    end subroutine init_FF
192    
193  
# Line 192 | Line 215 | contains
215      real( kind = DP ) :: pot_local
216      integer :: nrow
217      integer :: ncol
218 +    integer :: nprocs
219   #endif
220      integer :: nlocal
221      integer :: natoms    
# Line 221 | Line 245 | contains
245      nlocal = getNlocal()
246      natoms = nlocal
247   #endif
248 <  
248 >
249      call check_initialization(localError)
250      if ( localError .ne. 0 ) then
251 +       call handleError("do_force_loop","Not Initialized")
252         error = -1
253         return
254      end if
# Line 232 | Line 257 | contains
257      do_pot = do_pot_c
258      do_stress = do_stress_c
259  
260 +
261      ! Gather all information needed by all force loops:
262      
263   #ifdef IS_MPI    
# Line 248 | Line 274 | contains
274      endif
275      
276   #endif
277 <    
277 >
278 > !! Begin force loop timing:
279 > #ifdef PROFILE
280 >    call cpu_time(forceTimeInitial)
281 >    nloops = nloops + 1
282 > #endif
283 >  
284      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
285         !! See if we need to update neighbor lists
286         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
# Line 256 | Line 288 | contains
288         !! do_prepair_loop_if_needed
289         !! if_mpi_scatter_stuff_from_prepair
290         !! if_mpi_gather_stuff_from_prepair_to_main_loop
291 +    
292 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
293 + #ifdef IS_MPI
294 +    
295 +    if (update_nlist) then
296 +      
297 +       !! save current configuration, construct neighbor list,
298 +       !! and calculate forces
299 +       call saveNeighborList(nlocal, q)
300 +      
301 +       neighborListSize = size(list)
302 +       nlist = 0      
303 +      
304 +       do i = 1, nrow
305 +          point(i) = nlist + 1
306 +          
307 +          prepair_inner: do j = 1, ncol
308 +            
309 +             if (skipThisPair(i,j)) cycle prepair_inner
310 +            
311 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
312 +            
313 +             if (rijsq < rlistsq) then            
314 +                
315 +                nlist = nlist + 1
316 +                
317 +                if (nlist > neighborListSize) then
318 +                   call expandNeighborList(nlocal, listerror)
319 +                   if (listerror /= 0) then
320 +                      error = -1
321 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
322 +                      return
323 +                   end if
324 +                   neighborListSize = size(list)
325 +                endif
326 +                
327 +                list(nlist) = j
328 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
329 +             endif
330 +          enddo prepair_inner
331 +       enddo
332 +
333 +       point(nrow + 1) = nlist + 1
334 +      
335 +    else  !! (of update_check)
336 +
337 +       ! use the list to find the neighbors
338 +       do i = 1, nrow
339 +          JBEG = POINT(i)
340 +          JEND = POINT(i+1) - 1
341 +          ! check thiat molecule i has neighbors
342 +          if (jbeg .le. jend) then
343 +            
344 +             do jnab = jbeg, jend
345 +                j = list(jnab)
346 +
347 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
348 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
349 +                     u_l, A, f, t, pot_local)
350 +
351 +             enddo
352 +          endif
353 +       enddo
354 +    endif
355 +    
356 + #else
357 +    
358 +    if (update_nlist) then
359 +      
360 +       ! save current configuration, contruct neighbor list,
361 +       ! and calculate forces
362 +       call saveNeighborList(natoms, q)
363 +      
364 +       neighborListSize = size(list)
365 +  
366 +       nlist = 0
367 +
368 +       do i = 1, natoms-1
369 +          point(i) = nlist + 1
370 +          
371 +          prepair_inner: do j = i+1, natoms
372 +            
373 +             if (skipThisPair(i,j))  cycle prepair_inner
374 +                          
375 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
376 +          
377 +
378 +             if (rijsq < rlistsq) then
379 +
380 +          
381 +                nlist = nlist + 1
382 +              
383 +                if (nlist > neighborListSize) then
384 +                   call expandNeighborList(natoms, listerror)
385 +                   if (listerror /= 0) then
386 +                      error = -1
387 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
388 +                      return
389 +                   end if
390 +                   neighborListSize = size(list)
391 +                endif
392 +                
393 +                list(nlist) = j
394 +                
395 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
396 +                        u_l, A, f, t, pot)
397 +                
398 +             endif
399 +          enddo prepair_inner
400 +       enddo
401 +      
402 +       point(natoms) = nlist + 1
403 +      
404 +    else !! (update)
405 +  
406 +       ! use the list to find the neighbors
407 +       do i = 1, natoms-1
408 +          JBEG = POINT(i)
409 +          JEND = POINT(i+1) - 1
410 +          ! check thiat molecule i has neighbors
411 +          if (jbeg .le. jend) then
412 +            
413 +             do jnab = jbeg, jend
414 +                j = list(jnab)
415 +
416 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
417 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
418 +                     u_l, A, f, t, pot)
419 +
420 +             enddo
421 +          endif
422 +       enddo
423 +    endif    
424 + #endif
425 +    !! Do rest of preforce calculations
426 +    !! do necessary preforce calculations  
427 +    call do_preforce(nlocal,pot)
428 +   ! we have already updated the neighbor list set it to false...
429 +   update_nlist = .false.
430      else
431 <       !! See if we need to update neighbor lists
431 >       !! See if we need to update neighbor lists for non pre-pair
432         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
433      endif
434 <    
434 >
435 >
436 >
437 >
438 >
439 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
440 >
441 >
442 >
443 >
444 >  
445   #ifdef IS_MPI
446      
447      if (update_nlist) then
# Line 398 | Line 579 | contains
579   #endif
580      
581      ! phew, done with main loop.
582 <    
582 >
583 > !! Do timing
584 > #ifdef PROFILE
585 >    call cpu_time(forceTimeFinal)
586 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
587 > #endif
588 >
589 >
590   #ifdef IS_MPI
591      !!distribute forces
592    
# Line 512 | Line 700 | contains
700      endif
701  
702   #endif
703 <    
703 >
704 > #ifdef PROFILE
705 >    if (do_pot) then
706 >
707 > #ifdef IS_MPI
708 >
709 >      
710 >       call printCommTime()
711 >
712 >       call mpi_allreduce(forceTime,globalForceTime,1,MPI_DOUBLE_PRECISION, &
713 >            mpi_sum,mpi_comm_world,mpi_err)
714 >
715 >       call mpi_allreduce(forceTime,maxForceTime,1,MPI_DOUBLE_PRECISION, &
716 >            MPI_MAX,mpi_comm_world,mpi_err)
717 >      
718 >       call mpi_comm_size( MPI_COMM_WORLD, nprocs,mpi_err)
719 >      
720 >       if (getMyNode() == 0) then
721 >          write(*,*) "Total processor time spent in force calculations is: ", globalForceTime
722 >          write(*,*) "Total Time spent in force loop per processor is: ", globalforceTime/nprocs
723 >          write(*,*) "Maximum force time on any processor is: ", maxForceTime
724 >       end if
725 > #else
726 >       write(*,*) "Time spent in force loop is: ", forceTime
727 > #endif
728 >
729 >    
730 >    endif
731 >
732 > #endif
733 >
734 >
735 >
736    end subroutine do_force_loop
737  
738    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 531 | Line 751 | contains
751      logical :: is_LJ_i, is_LJ_j
752      logical :: is_DP_i, is_DP_j
753      logical :: is_GB_i, is_GB_j
754 +    logical :: is_EAM_i,is_EAM_j
755      logical :: is_Sticky_i, is_Sticky_j
756      integer :: me_i, me_j
757  
# Line 564 | Line 785 | contains
785         call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
786        
787         if ( is_DP_i .and. is_DP_j ) then
567 –          
788            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
789                 do_pot, do_stress)
790            if (FF_uses_RF .and. SimUsesRF()) then
# Line 589 | Line 809 | contains
809  
810      if (FF_uses_GB .and. SimUsesGB()) then
811  
812 +
813         call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
814         call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
815        
# Line 598 | Line 819 | contains
819         endif
820      endif
821      
822 +
823 +  
824 +   if (FF_uses_EAM .and. SimUsesEAM()) then
825 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
826 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
827 +      
828 +      if ( is_EAM_i .and. is_EAM_j ) &
829 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
830 +   endif
831 +
832 +
833  
834  
835    end subroutine do_pair
836  
837  
838  
839 <  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
839 >  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
840     real( kind = dp ) :: pot
841     real( kind = dp ), dimension(3,getNlocal()) :: u_l
842     real (kind=dp), dimension(9,getNlocal()) :: A
# Line 623 | Line 855 | contains
855    
856     r = sqrt(rijsq)
857    
858 +
859   #ifdef IS_MPI
860     if (tagRow(i) .eq. tagColumn(j)) then
861        write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 637 | Line 870 | contains
870     me_j = atid(j)
871    
872   #endif
873 <  
873 >    
874     if (FF_uses_EAM .and. SimUsesEAM()) then
875        call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
876        call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
877        
878        if ( is_EAM_i .and. is_EAM_j ) &
879 <           call calc_EAM_prepair(i, j, d, r, rijsq )
879 >           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
880     endif
881 +
882 + end subroutine do_prepair
883 +
884 +
885 +
886 +
887 +  subroutine do_preforce(nlocal,pot)
888 +    integer :: nlocal
889 +    real( kind = dp ) :: pot
890 +
891 +    if (FF_uses_EAM .and. SimUsesEAM()) then
892 +       call calc_EAM_preforce_Frho(nlocal,pot)
893 +    endif
894 +
895 +
896    end subroutine do_preforce
897    
898    
# Line 705 | Line 953 | contains
953      error = 0
954      ! Make sure we are properly initialized.
955      if (.not. do_forces_initialized) then
956 +       write(*,*) "Forces not initialized"
957         error = -1
958         return
959      endif
# Line 752 | Line 1001 | contains
1001  
1002   #endif
1003  
1004 +
1005 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1006 +       call clean_EAM()
1007 +    endif
1008 +
1009 +
1010 +
1011 +
1012 +
1013      rf = 0.0_dp
1014      tau_Temp = 0.0_dp
1015      virial_Temp = 0.0_dp
# Line 864 | Line 1122 | end module do_Forces
1122      doesit = FF_uses_RF
1123    end function FF_RequiresPostpairCalc
1124    
1125 + !! This cleans componets of force arrays belonging only to fortran
1126 +
1127   end module do_Forces

Diff Legend

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