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 841 by mmeineke, Wed Oct 29 17:55:28 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.34 2003-10-29 17:55:28 mmeineke Exp $, $Date: 2003-10-29 17:55:28 $, $Name: not supported by cvs2svn $, $Revision: 1.34 $
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 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 >    write(*,*) "In F, after: Tau = "
736 >    write(*,*) tau(1), tau(2), tau(3)
737 >    write(*,*) tau(4), tau(5), tau(6)
738 >    write(*,*) tau(7), tau(8), tau(9)
739 >    write(*,*) " "
740 >
741 >
742    end subroutine do_force_loop
743  
744    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 531 | Line 757 | contains
757      logical :: is_LJ_i, is_LJ_j
758      logical :: is_DP_i, is_DP_j
759      logical :: is_GB_i, is_GB_j
760 +    logical :: is_EAM_i,is_EAM_j
761      logical :: is_Sticky_i, is_Sticky_j
762      integer :: me_i, me_j
763  
# Line 564 | Line 791 | contains
791         call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
792        
793         if ( is_DP_i .and. is_DP_j ) then
567 –          
794            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
795                 do_pot, do_stress)
796            if (FF_uses_RF .and. SimUsesRF()) then
# Line 589 | Line 815 | contains
815  
816      if (FF_uses_GB .and. SimUsesGB()) then
817  
818 +
819         call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
820         call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
821        
# Line 598 | Line 825 | contains
825         endif
826      endif
827      
828 +
829 +  
830 +   if (FF_uses_EAM .and. SimUsesEAM()) then
831 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
832 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
833 +      
834 +      if ( is_EAM_i .and. is_EAM_j ) &
835 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
836 +   endif
837 +
838  
839  
840 +
841    end subroutine do_pair
842  
843  
844 +
845 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
846 +   real( kind = dp ) :: pot
847 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
848 +   real (kind=dp), dimension(9,getNlocal()) :: A
849 +   real (kind=dp), dimension(3,getNlocal()) :: f
850 +   real (kind=dp), dimension(3,getNlocal()) :: t
851 +  
852 +   logical, intent(inout) :: do_pot, do_stress
853 +   integer, intent(in) :: i, j
854 +   real ( kind = dp ), intent(inout)    :: rijsq
855 +   real ( kind = dp )                :: r
856 +   real ( kind = dp ), intent(inout) :: d(3)
857 +  
858 +   logical :: is_EAM_i, is_EAM_j
859 +  
860 +   integer :: me_i, me_j
861 +  
862 +   r = sqrt(rijsq)
863 +  
864 +
865 + #ifdef IS_MPI
866 +   if (tagRow(i) .eq. tagColumn(j)) then
867 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
868 +   endif
869 +  
870 +   me_i = atid_row(i)
871 +   me_j = atid_col(j)
872 +  
873 + #else
874 +  
875 +   me_i = atid(i)
876 +   me_j = atid(j)
877 +  
878 + #endif
879 +    
880 +   if (FF_uses_EAM .and. SimUsesEAM()) then
881 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
882 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
883 +      
884 +      if ( is_EAM_i .and. is_EAM_j ) &
885 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
886 +   endif
887 +
888 + end subroutine do_prepair
889 +
890 +
891 +
892 +
893 +  subroutine do_preforce(nlocal,pot)
894 +    integer :: nlocal
895 +    real( kind = dp ) :: pot
896 +
897 +    if (FF_uses_EAM .and. SimUsesEAM()) then
898 +       call calc_EAM_preforce_Frho(nlocal,pot)
899 +    endif
900 +
901 +
902 +  end subroutine do_preforce
903 +  
904 +  
905    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
906      
907      real (kind = dp), dimension(3) :: q_i
# Line 660 | Line 959 | contains
959      error = 0
960      ! Make sure we are properly initialized.
961      if (.not. do_forces_initialized) then
962 +       write(*,*) "Forces not initialized"
963         error = -1
964         return
965      endif
# Line 707 | Line 1007 | contains
1007  
1008   #endif
1009  
1010 +
1011 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1012 +       call clean_EAM()
1013 +    endif
1014 +
1015 +
1016 +
1017 +
1018 +
1019      rf = 0.0_dp
1020      tau_Temp = 0.0_dp
1021      virial_Temp = 0.0_dp
# Line 819 | Line 1128 | end module do_Forces
1128      doesit = FF_uses_RF
1129    end function FF_RequiresPostpairCalc
1130    
1131 + !! This cleans componets of force arrays belonging only to fortran
1132 +
1133   end module do_Forces

Diff Legend

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