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 634 by mmeineke, Thu Jul 17 19:38:23 2003 UTC vs.
Revision 895 by chuckv, Mon Jan 5 22:12:11 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.23 2003-07-17 19:38:23 mmeineke Exp $, $Date: 2003-07-17 19:38:23 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
7 > !! @version $Id: do_Forces.F90,v 1.42 2004-01-05 22:12:11 chuckv Exp $, $Date: 2004-01-05 22:12:11 $, $Name: not supported by cvs2svn $, $Revision: 1.42 $
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 +  logical, allocatable :: propertyMapI(:,:)
56 +  logical, allocatable :: propertyMapJ(:,:)
57 +
58   contains
59  
60    subroutine setRlistDF( this_rlist )
# Line 141 | Line 153 | contains
153            return
154         end if
155      endif
156 +
157 +
158 +    if (FF_uses_EAM) then
159 +         call init_EAM_FF(my_status)
160 +       if (my_status /= 0) then
161 +          write(*,*) "init_EAM_FF returned a bad status"
162 +          thisStat = -1
163 +          return
164 +       end if
165 +    endif
166 +
167 +
168      
169      if (FF_uses_GB) then
170         call check_gb_pair_FF(my_status)
# Line 161 | Line 185 | contains
185            return
186         endif
187      endif
188 +    
189  
190      havePolicies = .true.
191      if( haveRlist ) do_forces_initialized = .true.
192 <    
192 >
193    end subroutine init_FF
194    
195  
# Line 182 | Line 207 | contains
207      real ( kind = dp ), dimension(3,getNlocal()) :: f
208      !! Torsion array provided by C, dimensioned by getNlocal
209      real( kind = dp ), dimension(3,getNlocal()) :: t    
210 +
211      !! Stress Tensor
212      real( kind = dp), dimension(9) :: tau  
213      real ( kind = dp ) :: pot
# Line 192 | Line 218 | contains
218      real( kind = DP ) :: pot_local
219      integer :: nrow
220      integer :: ncol
221 +    integer :: nprocs
222   #endif
223      integer :: nlocal
224      integer :: natoms    
# Line 207 | Line 234 | contains
234      integer :: listerror, error
235      integer :: localError
236  
237 <    real(kind=dp) :: listSkin = 1.0
211 <    
237 >    real(kind=dp) :: listSkin = 1.0  
238  
239      !! initialize local variables  
240  
# Line 221 | Line 247 | contains
247      nlocal = getNlocal()
248      natoms = nlocal
249   #endif
250 <  
250 >
251      call check_initialization(localError)
252      if ( localError .ne. 0 ) then
253 +       call handleError("do_force_loop","Not Initialized")
254         error = -1
255         return
256      end if
# Line 232 | Line 259 | contains
259      do_pot = do_pot_c
260      do_stress = do_stress_c
261  
262 +
263 + #ifdef IS_MPI
264 +    if (.not.allocated(propertyMapI)) then
265 +       allocate(propertyMapI(5,getNrow())
266 +    endif
267 +
268 +    do i = 1, nrow
269 +       me_i = atid_row(i)
270 + #else
271 +    if (.not.allocated(propertyMapI)) then
272 +       allocate(propertyMapI(5,getNlocal())
273 +    endif
274 +
275 +    do i = 1, natoms
276 +       me_i = atid(i)
277 + #endif
278 +      
279 +       propertyMapI(1:5,i) = .false.
280 +
281 +       call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
282 +    
283 +       ! unpack the properties
284 +      
285 +       if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
286 +            propertyMapI(1, i) = .true.
287 +       if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
288 +            propertyMapI(2, i) = .true.
289 +       if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
290 +            propertyMapI(3, i) = .true.
291 +       if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
292 +            propertyMapI(4, i) = .true.
293 +       if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
294 +            propertyMapI(5, i) = .true.
295 +
296 +    end do
297 +
298 + #ifdef IS_MPI
299 +    if (.not.allocated(propertyMapJ)) then
300 +       allocate(propertyMapJ(5,getNcol())
301 +    endif
302 +
303 +    do j = 1, ncol
304 +       me_j = atid_col(j)
305 + #else
306 +    if (.not.allocated(propertyMapJ)) then
307 +       allocate(propertyMapJ(5,getNlocal())
308 +    endif
309 +
310 +    do j = 1, natoms
311 +       me_j = atid(j)
312 + #endif
313 +      
314 +       propertyMapJ(1:5,j) = .false.
315 +
316 +       call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
317 +    
318 +       ! unpack the properties
319 +      
320 +       if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
321 +            propertyMapJ(1, j) = .true.
322 +       if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
323 +            propertyMapJ(2, j) = .true.
324 +       if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
325 +            propertyMapJ(3, j) = .true.
326 +       if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
327 +            propertyMapJ(4, j) = .true.
328 +       if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
329 +            propertyMapJ(5, j) = .true.
330 +
331 +    end do
332 +
333      ! Gather all information needed by all force loops:
334      
335   #ifdef IS_MPI    
# Line 248 | Line 346 | contains
346      endif
347      
348   #endif
349 <    
349 >
350 > !! Begin force loop timing:
351 > #ifdef PROFILE
352 >    call cpu_time(forceTimeInitial)
353 >    nloops = nloops + 1
354 > #endif
355 >  
356      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
357         !! See if we need to update neighbor lists
358         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
# Line 256 | Line 360 | contains
360         !! do_prepair_loop_if_needed
361         !! if_mpi_scatter_stuff_from_prepair
362         !! if_mpi_gather_stuff_from_prepair_to_main_loop
363 +    
364 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
365 + #ifdef IS_MPI
366 +    
367 +    if (update_nlist) then
368 +      
369 +       !! save current configuration, construct neighbor list,
370 +       !! and calculate forces
371 +       call saveNeighborList(nlocal, q)
372 +      
373 +       neighborListSize = size(list)
374 +       nlist = 0      
375 +      
376 +       do i = 1, nrow
377 +          point(i) = nlist + 1
378 +          
379 +          prepair_inner: do j = 1, ncol
380 +            
381 +             if (skipThisPair(i,j)) cycle prepair_inner
382 +            
383 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
384 +            
385 +             if (rijsq < rlistsq) then            
386 +                
387 +                nlist = nlist + 1
388 +                
389 +                if (nlist > neighborListSize) then
390 +                   call expandNeighborList(nlocal, listerror)
391 +                   if (listerror /= 0) then
392 +                      error = -1
393 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
394 +                      return
395 +                   end if
396 +                   neighborListSize = size(list)
397 +                endif
398 +                
399 +                list(nlist) = j
400 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
401 +             endif
402 +          enddo prepair_inner
403 +       enddo
404 +
405 +       point(nrow + 1) = nlist + 1
406 +      
407 +    else  !! (of update_check)
408 +
409 +       ! use the list to find the neighbors
410 +       do i = 1, nrow
411 +          JBEG = POINT(i)
412 +          JEND = POINT(i+1) - 1
413 +          ! check thiat molecule i has neighbors
414 +          if (jbeg .le. jend) then
415 +            
416 +             do jnab = jbeg, jend
417 +                j = list(jnab)
418 +
419 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
420 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
421 +                     u_l, A, f, t, pot_local)
422 +
423 +             enddo
424 +          endif
425 +       enddo
426 +    endif
427 +    
428 + #else
429 +    
430 +    if (update_nlist) then
431 +      
432 +       ! save current configuration, contruct neighbor list,
433 +       ! and calculate forces
434 +       call saveNeighborList(natoms, q)
435 +      
436 +       neighborListSize = size(list)
437 +  
438 +       nlist = 0
439 +
440 +       do i = 1, natoms-1
441 +          point(i) = nlist + 1
442 +          
443 +          prepair_inner: do j = i+1, natoms
444 +            
445 +             if (skipThisPair(i,j))  cycle prepair_inner
446 +                          
447 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
448 +          
449 +
450 +             if (rijsq < rlistsq) then
451 +
452 +          
453 +                nlist = nlist + 1
454 +              
455 +                if (nlist > neighborListSize) then
456 +                   call expandNeighborList(natoms, listerror)
457 +                   if (listerror /= 0) then
458 +                      error = -1
459 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
460 +                      return
461 +                   end if
462 +                   neighborListSize = size(list)
463 +                endif
464 +                
465 +                list(nlist) = j
466 +                
467 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
468 +                        u_l, A, f, t, pot)
469 +                
470 +             endif
471 +          enddo prepair_inner
472 +       enddo
473 +      
474 +       point(natoms) = nlist + 1
475 +      
476 +    else !! (update)
477 +  
478 +       ! use the list to find the neighbors
479 +       do i = 1, natoms-1
480 +          JBEG = POINT(i)
481 +          JEND = POINT(i+1) - 1
482 +          ! check thiat molecule i has neighbors
483 +          if (jbeg .le. jend) then
484 +            
485 +             do jnab = jbeg, jend
486 +                j = list(jnab)
487 +
488 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
489 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
490 +                     u_l, A, f, t, pot)
491 +
492 +             enddo
493 +          endif
494 +       enddo
495 +    endif    
496 + #endif
497 +    !! Do rest of preforce calculations
498 +    !! do necessary preforce calculations  
499 +    call do_preforce(nlocal,pot)
500 +   ! we have already updated the neighbor list set it to false...
501 +   update_nlist = .false.
502      else
503 <       !! See if we need to update neighbor lists
503 >       !! See if we need to update neighbor lists for non pre-pair
504         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
505      endif
506 <    
506 >
507 >
508 >
509 >
510 >
511 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
512 >
513 >
514 >
515 >
516 >  
517   #ifdef IS_MPI
518      
519      if (update_nlist) then
267 –      
520         !! save current configuration, construct neighbor list,
521         !! and calculate forces
522         call saveNeighborList(nlocal, q)
# Line 273 | Line 525 | contains
525         nlist = 0      
526        
527         do i = 1, nrow
528 +
529            point(i) = nlist + 1
530            
531            inner: do j = 1, ncol
# Line 330 | Line 583 | contains
583   #else
584      
585      if (update_nlist) then
586 <      
586 >
587         ! save current configuration, contruct neighbor list,
588         ! and calculate forces
589         call saveNeighborList(natoms, q)
# Line 398 | Line 651 | contains
651   #endif
652      
653      ! phew, done with main loop.
654 <    
654 >
655 > !! Do timing
656 > #ifdef PROFILE
657 >    call cpu_time(forceTimeFinal)
658 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
659 > #endif
660 >
661 >
662   #ifdef IS_MPI
663      !!distribute forces
664    
# Line 510 | Line 770 | contains
770         tau = tau_Temp
771         virial = virial_Temp
772      endif
773 <
773 >    
774   #endif
775      
776 +    
777 +    
778    end subroutine do_force_loop
779  
780    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 531 | Line 793 | contains
793      logical :: is_LJ_i, is_LJ_j
794      logical :: is_DP_i, is_DP_j
795      logical :: is_GB_i, is_GB_j
796 +    logical :: is_EAM_i,is_EAM_j
797      logical :: is_Sticky_i, is_Sticky_j
798      integer :: me_i, me_j
799 <
799 >    integer :: propPack_i
800 >    integer :: propPack_j
801      r = sqrt(rijsq)
802  
803   #ifdef IS_MPI
# Line 550 | Line 814 | contains
814      me_j = atid(j)
815  
816   #endif
817 <
817 >    
818      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)
819  
820 <       if ( is_LJ_i .and. is_LJ_j ) &
820 >       if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
821              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
822 +
823      endif
824  
825      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)
826        
827 <       if ( is_DP_i .and. is_DP_j ) then
567 <          
827 >       if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
828            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
829                 do_pot, do_stress)
830            if (FF_uses_RF .and. SimUsesRF()) then
# Line 577 | Line 837 | contains
837  
838      if (FF_uses_Sticky .and. SimUsesSticky()) then
839  
840 <       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
581 <       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
582 <
583 <       if ( is_Sticky_i .and. is_Sticky_j ) then
840 >       if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
841            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
842                 do_pot, do_stress)
843         endif
# Line 588 | Line 845 | contains
845  
846  
847      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)
848        
849 <       if ( is_GB_i .and. is_GB_j ) then
849 >       if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
850            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
851                 do_pot, do_stress)          
852         endif
853 +
854      endif
855      
856  
857 +  
858 +   if (FF_uses_EAM .and. SimUsesEAM()) then
859 +      
860 +      if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
861 +         call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
862 +      endif
863  
864 +   endif
865 +
866    end subroutine do_pair
867  
868  
869  
870 <  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
870 >  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
871     real( kind = dp ) :: pot
872     real( kind = dp ), dimension(3,getNlocal()) :: u_l
873     real (kind=dp), dimension(9,getNlocal()) :: A
# Line 623 | Line 886 | contains
886    
887     r = sqrt(rijsq)
888    
889 +
890   #ifdef IS_MPI
891     if (tagRow(i) .eq. tagColumn(j)) then
892        write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 637 | Line 901 | contains
901     me_j = atid(j)
902    
903   #endif
904 <  
904 >    
905     if (FF_uses_EAM .and. SimUsesEAM()) then
906        call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
907        call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
908        
909 < !!$      if ( is_EAM_i .and. is_EAM_j ) &
910 < !!$           call calc_EAM_prepair(i, j, d, r, rijsq )
909 >      if ( is_EAM_i .and. is_EAM_j ) &
910 >           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
911     endif
912 +
913 + end subroutine do_prepair
914 +
915 +
916 +
917 +
918 +  subroutine do_preforce(nlocal,pot)
919 +    integer :: nlocal
920 +    real( kind = dp ) :: pot
921 +
922 +    if (FF_uses_EAM .and. SimUsesEAM()) then
923 +       call calc_EAM_preforce_Frho(nlocal,pot)
924 +    endif
925 +
926 +
927    end subroutine do_preforce
928    
929    
# Line 705 | Line 984 | contains
984      error = 0
985      ! Make sure we are properly initialized.
986      if (.not. do_forces_initialized) then
987 +       write(*,*) "Forces not initialized"
988         error = -1
989         return
990      endif
# Line 752 | Line 1032 | contains
1032  
1033   #endif
1034  
1035 +
1036 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1037 +       call clean_EAM()
1038 +    endif
1039 +
1040 +
1041 +
1042 +
1043 +
1044      rf = 0.0_dp
1045      tau_Temp = 0.0_dp
1046      virial_Temp = 0.0_dp
# Line 864 | Line 1153 | end module do_Forces
1153      doesit = FF_uses_RF
1154    end function FF_RequiresPostpairCalc
1155    
1156 + #ifdef PROFILE
1157 +  function getforcetime() result(totalforcetime)
1158 +    real(kind=dp) :: totalforcetime
1159 +    totalforcetime = forcetime
1160 +  end function getforcetime
1161 + #endif
1162 +
1163 + !! This cleans componets of force arrays belonging only to fortran
1164 +
1165   end module do_Forces

Diff Legend

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