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

Diff Legend

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