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 894 by chuckv, Mon Jan 5 21:00:05 2004 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.41 2004-01-05 21:00:05 chuckv Exp $, $Date: 2004-01-05 21:00:05 $, $Name: not supported by cvs2svn $, $Revision: 1.41 $
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 +  public :: getforcetime
50 +  real, save :: forceTime = 0
51 +  real :: forceTimeInitial, forceTimeFinal
52 +  integer :: nLoops
53 + #endif
54 +
55   contains
56  
57    subroutine setRlistDF( this_rlist )
# Line 141 | Line 150 | contains
150            return
151         end if
152      endif
153 +
154 +
155 +    if (FF_uses_EAM) then
156 +         call init_EAM_FF(my_status)
157 +       if (my_status /= 0) then
158 +          write(*,*) "init_EAM_FF returned a bad status"
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 207 | Line 230 | contains
230      integer :: listerror, error
231      integer :: localError
232  
233 <    real(kind=dp) :: listSkin = 1.0
211 <    
233 >    real(kind=dp) :: listSkin = 1.0  
234  
235      !! initialize local variables  
236  
# Line 221 | Line 243 | contains
243      nlocal = getNlocal()
244      natoms = nlocal
245   #endif
246 <  
246 >
247      call check_initialization(localError)
248      if ( localError .ne. 0 ) then
249 +       call handleError("do_force_loop","Not Initialized")
250         error = -1
251         return
252      end if
# Line 232 | Line 255 | contains
255      do_pot = do_pot_c
256      do_stress = do_stress_c
257  
258 +
259      ! Gather all information needed by all force loops:
260      
261   #ifdef IS_MPI    
# Line 248 | Line 272 | contains
272      endif
273      
274   #endif
275 <    
275 >
276 > !! Begin force loop timing:
277 > #ifdef PROFILE
278 >    call cpu_time(forceTimeInitial)
279 >    nloops = nloops + 1
280 > #endif
281 >  
282      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
283         !! See if we need to update neighbor lists
284         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
# Line 256 | Line 286 | contains
286         !! do_prepair_loop_if_needed
287         !! if_mpi_scatter_stuff_from_prepair
288         !! if_mpi_gather_stuff_from_prepair_to_main_loop
289 +    
290 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
291 + #ifdef IS_MPI
292 +    
293 +    if (update_nlist) then
294 +      
295 +       !! save current configuration, construct neighbor list,
296 +       !! and calculate forces
297 +       call saveNeighborList(nlocal, q)
298 +      
299 +       neighborListSize = size(list)
300 +       nlist = 0      
301 +      
302 +       do i = 1, nrow
303 +          point(i) = nlist + 1
304 +          
305 +          prepair_inner: do j = 1, ncol
306 +            
307 +             if (skipThisPair(i,j)) cycle prepair_inner
308 +            
309 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
310 +            
311 +             if (rijsq < rlistsq) then            
312 +                
313 +                nlist = nlist + 1
314 +                
315 +                if (nlist > neighborListSize) then
316 +                   call expandNeighborList(nlocal, listerror)
317 +                   if (listerror /= 0) then
318 +                      error = -1
319 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
320 +                      return
321 +                   end if
322 +                   neighborListSize = size(list)
323 +                endif
324 +                
325 +                list(nlist) = j
326 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
327 +             endif
328 +          enddo prepair_inner
329 +       enddo
330 +
331 +       point(nrow + 1) = nlist + 1
332 +      
333 +    else  !! (of update_check)
334 +
335 +       ! use the list to find the neighbors
336 +       do i = 1, nrow
337 +          JBEG = POINT(i)
338 +          JEND = POINT(i+1) - 1
339 +          ! check thiat molecule i has neighbors
340 +          if (jbeg .le. jend) then
341 +            
342 +             do jnab = jbeg, jend
343 +                j = list(jnab)
344 +
345 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
346 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
347 +                     u_l, A, f, t, pot_local)
348 +
349 +             enddo
350 +          endif
351 +       enddo
352 +    endif
353 +    
354 + #else
355 +    
356 +    if (update_nlist) then
357 +      
358 +       ! save current configuration, contruct neighbor list,
359 +       ! and calculate forces
360 +       call saveNeighborList(natoms, q)
361 +      
362 +       neighborListSize = size(list)
363 +  
364 +       nlist = 0
365 +
366 +       do i = 1, natoms-1
367 +          point(i) = nlist + 1
368 +          
369 +          prepair_inner: do j = i+1, natoms
370 +            
371 +             if (skipThisPair(i,j))  cycle prepair_inner
372 +                          
373 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
374 +          
375 +
376 +             if (rijsq < rlistsq) then
377 +
378 +          
379 +                nlist = nlist + 1
380 +              
381 +                if (nlist > neighborListSize) then
382 +                   call expandNeighborList(natoms, listerror)
383 +                   if (listerror /= 0) then
384 +                      error = -1
385 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
386 +                      return
387 +                   end if
388 +                   neighborListSize = size(list)
389 +                endif
390 +                
391 +                list(nlist) = j
392 +                
393 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
394 +                        u_l, A, f, t, pot)
395 +                
396 +             endif
397 +          enddo prepair_inner
398 +       enddo
399 +      
400 +       point(natoms) = nlist + 1
401 +      
402 +    else !! (update)
403 +  
404 +       ! use the list to find the neighbors
405 +       do i = 1, natoms-1
406 +          JBEG = POINT(i)
407 +          JEND = POINT(i+1) - 1
408 +          ! check thiat molecule i has neighbors
409 +          if (jbeg .le. jend) then
410 +            
411 +             do jnab = jbeg, jend
412 +                j = list(jnab)
413 +
414 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
415 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
416 +                     u_l, A, f, t, pot)
417 +
418 +             enddo
419 +          endif
420 +       enddo
421 +    endif    
422 + #endif
423 +    !! Do rest of preforce calculations
424 +    !! do necessary preforce calculations  
425 +    call do_preforce(nlocal,pot)
426 +   ! we have already updated the neighbor list set it to false...
427 +   update_nlist = .false.
428      else
429 <       !! See if we need to update neighbor lists
429 >       !! See if we need to update neighbor lists for non pre-pair
430         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
431      endif
432 <    
432 >
433 >
434 >
435 >
436 >
437 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
438 >
439 >
440 >
441 >
442 >  
443   #ifdef IS_MPI
444      
445      if (update_nlist) then
267 –      
446         !! save current configuration, construct neighbor list,
447         !! and calculate forces
448         call saveNeighborList(nlocal, q)
# Line 330 | Line 508 | contains
508   #else
509      
510      if (update_nlist) then
511 <      
511 >
512         ! save current configuration, contruct neighbor list,
513         ! and calculate forces
514         call saveNeighborList(natoms, q)
# Line 398 | Line 576 | contains
576   #endif
577      
578      ! phew, done with main loop.
579 <    
579 >
580 > !! Do timing
581 > #ifdef PROFILE
582 >    call cpu_time(forceTimeFinal)
583 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
584 > #endif
585 >
586 >
587   #ifdef IS_MPI
588      !!distribute forces
589    
# Line 510 | Line 695 | contains
695         tau = tau_Temp
696         virial = virial_Temp
697      endif
698 <
698 >    
699   #endif
700      
701 +    
702 +    
703    end subroutine do_force_loop
704  
705    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 531 | Line 718 | contains
718      logical :: is_LJ_i, is_LJ_j
719      logical :: is_DP_i, is_DP_j
720      logical :: is_GB_i, is_GB_j
721 +    logical :: is_EAM_i,is_EAM_j
722      logical :: is_Sticky_i, is_Sticky_j
723      integer :: me_i, me_j
724 <
724 >    integer :: propPack_i
725 >    integer :: propPack_j
726      r = sqrt(rijsq)
727  
728   #ifdef IS_MPI
# Line 550 | Line 739 | contains
739      me_j = atid(j)
740  
741   #endif
742 +
743 +    call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
744 +    call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
745 +    
746 +    ! unpack the properties
747 +
748 +    if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) is_LJ_i = .true.
749 +    if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) is_DP_i = .true.
750 +    if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) is_Sticky_i = .true.
751 +    if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) is_GB_i = .true.
752 +    if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) is_EAM_i = .true.
753  
754 +    if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) is_LJ_j = .true.
755 +    if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) is_DP_j = .true.
756 +    if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) is_Sticky_j = .true.
757 +    if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) is_GB_j = .true.
758 +    if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) is_EAM_j = .true.
759 +
760 +
761      if (FF_uses_LJ .and. SimUsesLJ()) then
555 –       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
556 –       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
762  
763         if ( is_LJ_i .and. is_LJ_j ) &
764              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
765 +
766      endif
767  
768      if (FF_uses_dipoles .and. SimUsesDipoles()) then
563 –       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
564 –       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
769        
770         if ( is_DP_i .and. is_DP_j ) then
567 –          
771            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
772                 do_pot, do_stress)
773            if (FF_uses_RF .and. SimUsesRF()) then
# Line 577 | Line 780 | contains
780  
781      if (FF_uses_Sticky .and. SimUsesSticky()) then
782  
580 –       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
581 –       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
582 –
783         if ( is_Sticky_i .and. is_Sticky_j ) then
784            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
785                 do_pot, do_stress)
# Line 588 | Line 788 | contains
788  
789  
790      if (FF_uses_GB .and. SimUsesGB()) then
591 –
592 –       call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
593 –       call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
791        
792         if ( is_GB_i .and. is_GB_j ) then
793            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
794                 do_pot, do_stress)          
795         endif
796 +
797      endif
798      
799  
800 +  
801 +   if (FF_uses_EAM .and. SimUsesEAM()) then
802 +      
803 +      if ( is_EAM_i .and. is_EAM_j ) &
804 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
805  
806 +   endif
807 +
808    end subroutine do_pair
809  
810  
811  
812 <  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
812 >  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
813     real( kind = dp ) :: pot
814     real( kind = dp ), dimension(3,getNlocal()) :: u_l
815     real (kind=dp), dimension(9,getNlocal()) :: A
# Line 623 | Line 828 | contains
828    
829     r = sqrt(rijsq)
830    
831 +
832   #ifdef IS_MPI
833     if (tagRow(i) .eq. tagColumn(j)) then
834        write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 637 | Line 843 | contains
843     me_j = atid(j)
844    
845   #endif
846 <  
846 >    
847     if (FF_uses_EAM .and. SimUsesEAM()) then
848        call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
849        call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
850        
851        if ( is_EAM_i .and. is_EAM_j ) &
852 <           call calc_EAM_prepair(i, j, d, r, rijsq )
852 >           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
853     endif
854 +
855 + end subroutine do_prepair
856 +
857 +
858 +
859 +
860 +  subroutine do_preforce(nlocal,pot)
861 +    integer :: nlocal
862 +    real( kind = dp ) :: pot
863 +
864 +    if (FF_uses_EAM .and. SimUsesEAM()) then
865 +       call calc_EAM_preforce_Frho(nlocal,pot)
866 +    endif
867 +
868 +
869    end subroutine do_preforce
870    
871    
# Line 705 | Line 926 | contains
926      error = 0
927      ! Make sure we are properly initialized.
928      if (.not. do_forces_initialized) then
929 +       write(*,*) "Forces not initialized"
930         error = -1
931         return
932      endif
# Line 752 | Line 974 | contains
974  
975   #endif
976  
977 +
978 +    if (FF_uses_EAM .and. SimUsesEAM()) then
979 +       call clean_EAM()
980 +    endif
981 +
982 +
983 +
984 +
985 +
986      rf = 0.0_dp
987      tau_Temp = 0.0_dp
988      virial_Temp = 0.0_dp
# Line 864 | Line 1095 | end module do_Forces
1095      doesit = FF_uses_RF
1096    end function FF_RequiresPostpairCalc
1097    
1098 + #ifdef PROFILE
1099 +  function getforcetime() result(totalforcetime)
1100 +    real(kind=dp) :: totalforcetime
1101 +    totalforcetime = forcetime
1102 +  end function getforcetime
1103 + #endif
1104 +
1105 + !! This cleans componets of force arrays belonging only to fortran
1106 +
1107   end module do_Forces

Diff Legend

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