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

Diff Legend

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