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 898 by chuckv, Mon Jan 5 22:49:14 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.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.44 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
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 154 | Line 178 | contains
178      endif
179      if (.not. do_forces_initialized) then
180         !! Create neighbor lists
181 <       call expandNeighborList(getNlocal(), my_status)
181 >       call expandNeighborList(nLocal, my_status)
182         if (my_Status /= 0) then
183            write(default_error,*) "SimSetup: ExpandNeighborList returned error."
184            thisStat = -1
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 173 | Line 198 | contains
198    subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
199         error)
200      !! Position array provided by C, dimensioned by getNlocal
201 <    real ( kind = dp ), dimension(3,getNlocal()) :: q
201 >    real ( kind = dp ), dimension(3,nLocal) :: q
202      !! Rotation Matrix for each long range particle in simulation.
203 <    real( kind = dp), dimension(9,getNlocal()) :: A    
203 >    real( kind = dp), dimension(9,nLocal) :: A    
204      !! Unit vectors for dipoles (lab frame)
205 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
205 >    real( kind = dp ), dimension(3,nLocal) :: u_l
206      !! Force array provided by C, dimensioned by getNlocal
207 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
207 >    real ( kind = dp ), dimension(3,nLocal) :: f
208      !! Torsion array provided by C, dimensioned by getNlocal
209 <    real( kind = dp ), dimension(3,getNlocal()) :: t    
209 >    real( kind = dp ), dimension(3,nLocal) :: 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
196 –    integer :: nlocal
223      integer :: natoms    
224      logical :: update_nlist  
225      integer :: i, j, jbeg, jend, jnab
# Line 201 | Line 227 | contains
227      real( kind = DP ) ::  rijsq
228      real(kind=dp),dimension(3) :: d
229      real(kind=dp) :: rfpot, mu_i, virial
230 <    integer :: me_i
230 >    integer :: me_i, me_j
231      logical :: is_dp_i
232      integer :: neighborListSize
233      integer :: listerror, error
234      integer :: localError
235 +    integer :: propPack_i, propPack_j
236  
237 <    real(kind=dp) :: listSkin = 1.0
211 <    
237 >    real(kind=dp) :: listSkin = 1.0  
238  
239      !! initialize local variables  
240  
241   #ifdef IS_MPI
242      pot_local = 0.0_dp
217 –    nlocal = getNlocal()
243      nrow   = getNrow(plan_row)
244      ncol   = getNcol(plan_col)
245   #else
221 –    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 231 | Line 256 | contains
256  
257      do_pot = do_pot_c
258      do_stress = do_stress_c
259 +
260 +
261 + #ifdef IS_MPI
262 +    if (.not.allocated(propertyMapI)) then
263 +       allocate(propertyMapI(5,nrow))
264 +    endif
265 +
266 +    do i = 1, nrow
267 +       me_i = atid_row(i)
268 + #else
269 +    if (.not.allocated(propertyMapI)) then
270 +       allocate(propertyMapI(5,nlocal))
271 +    endif
272 +
273 +    do i = 1, natoms
274 +       me_i = atid(i)
275 + #endif
276 +      
277 +       propertyMapI(1:5,i) = .false.
278 +
279 +       call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
280 +    
281 +       ! unpack the properties
282 +      
283 +       if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
284 +            propertyMapI(1, i) = .true.
285 +       if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
286 +            propertyMapI(2, i) = .true.
287 +       if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
288 +            propertyMapI(3, i) = .true.
289 +       if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
290 +            propertyMapI(4, i) = .true.
291 +       if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
292 +            propertyMapI(5, i) = .true.
293 +
294 +    end do
295 +
296 + #ifdef IS_MPI
297 +    if (.not.allocated(propertyMapJ)) then
298 +       allocate(propertyMapJ(5,ncol))
299 +    endif
300 +
301 +    do j = 1, ncol
302 +       me_j = atid_col(j)
303 + #else
304 +    if (.not.allocated(propertyMapJ)) then
305 +       allocate(propertyMapJ(5,nlocal))
306 +    endif
307 +
308 +    do j = 1, natoms
309 +       me_j = atid(j)
310 + #endif
311 +      
312 +       propertyMapJ(1:5,j) = .false.
313 +
314 +       call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
315 +    
316 +       ! unpack the properties
317 +      
318 +       if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
319 +            propertyMapJ(1, j) = .true.
320 +       if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
321 +            propertyMapJ(2, j) = .true.
322 +       if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
323 +            propertyMapJ(3, j) = .true.
324 +       if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
325 +            propertyMapJ(4, j) = .true.
326 +       if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
327 +            propertyMapJ(5, j) = .true.
328  
329 +    end do
330 +
331      ! Gather all information needed by all force loops:
332      
333   #ifdef IS_MPI    
# Line 248 | Line 344 | contains
344      endif
345      
346   #endif
347 <    
347 >
348 > !! Begin force loop timing:
349 > #ifdef PROFILE
350 >    call cpu_time(forceTimeInitial)
351 >    nloops = nloops + 1
352 > #endif
353 >  
354      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
355         !! See if we need to update neighbor lists
356         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
# Line 256 | Line 358 | contains
358         !! do_prepair_loop_if_needed
359         !! if_mpi_scatter_stuff_from_prepair
360         !! if_mpi_gather_stuff_from_prepair_to_main_loop
361 +    
362 + !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
363 + #ifdef IS_MPI
364 +    
365 +    if (update_nlist) then
366 +      
367 +       !! save current configuration, construct neighbor list,
368 +       !! and calculate forces
369 +       call saveNeighborList(nlocal, q)
370 +      
371 +       neighborListSize = size(list)
372 +       nlist = 0      
373 +      
374 +       do i = 1, nrow
375 +          point(i) = nlist + 1
376 +          
377 +          prepair_inner: do j = 1, ncol
378 +            
379 +             if (skipThisPair(i,j)) cycle prepair_inner
380 +            
381 +             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
382 +            
383 +             if (rijsq < rlistsq) then            
384 +                
385 +                nlist = nlist + 1
386 +                
387 +                if (nlist > neighborListSize) then
388 +                   call expandNeighborList(nlocal, listerror)
389 +                   if (listerror /= 0) then
390 +                      error = -1
391 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
392 +                      return
393 +                   end if
394 +                   neighborListSize = size(list)
395 +                endif
396 +                
397 +                list(nlist) = j
398 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
399 +             endif
400 +          enddo prepair_inner
401 +       enddo
402 +
403 +       point(nrow + 1) = nlist + 1
404 +      
405 +    else  !! (of update_check)
406 +
407 +       ! use the list to find the neighbors
408 +       do i = 1, nrow
409 +          JBEG = POINT(i)
410 +          JEND = POINT(i+1) - 1
411 +          ! check thiat molecule i has neighbors
412 +          if (jbeg .le. jend) then
413 +            
414 +             do jnab = jbeg, jend
415 +                j = list(jnab)
416 +
417 +                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
418 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
419 +                     u_l, A, f, t, pot_local)
420 +
421 +             enddo
422 +          endif
423 +       enddo
424 +    endif
425 +    
426 + #else
427 +    
428 +    if (update_nlist) then
429 +      
430 +       ! save current configuration, contruct neighbor list,
431 +       ! and calculate forces
432 +       call saveNeighborList(natoms, q)
433 +      
434 +       neighborListSize = size(list)
435 +  
436 +       nlist = 0
437 +
438 +       do i = 1, natoms-1
439 +          point(i) = nlist + 1
440 +          
441 +          prepair_inner: do j = i+1, natoms
442 +            
443 +             if (skipThisPair(i,j))  cycle prepair_inner
444 +                          
445 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
446 +          
447 +
448 +             if (rijsq < rlistsq) then
449 +
450 +          
451 +                nlist = nlist + 1
452 +              
453 +                if (nlist > neighborListSize) then
454 +                   call expandNeighborList(natoms, listerror)
455 +                   if (listerror /= 0) then
456 +                      error = -1
457 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
458 +                      return
459 +                   end if
460 +                   neighborListSize = size(list)
461 +                endif
462 +                
463 +                list(nlist) = j
464 +                
465 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
466 +                        u_l, A, f, t, pot)
467 +                
468 +             endif
469 +          enddo prepair_inner
470 +       enddo
471 +      
472 +       point(natoms) = nlist + 1
473 +      
474 +    else !! (update)
475 +  
476 +       ! use the list to find the neighbors
477 +       do i = 1, natoms-1
478 +          JBEG = POINT(i)
479 +          JEND = POINT(i+1) - 1
480 +          ! check thiat molecule i has neighbors
481 +          if (jbeg .le. jend) then
482 +            
483 +             do jnab = jbeg, jend
484 +                j = list(jnab)
485 +
486 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
487 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
488 +                     u_l, A, f, t, pot)
489 +
490 +             enddo
491 +          endif
492 +       enddo
493 +    endif    
494 + #endif
495 +    !! Do rest of preforce calculations
496 +    !! do necessary preforce calculations  
497 +    call do_preforce(nlocal,pot)
498 +   ! we have already updated the neighbor list set it to false...
499 +   update_nlist = .false.
500      else
501 <       !! See if we need to update neighbor lists
501 >       !! See if we need to update neighbor lists for non pre-pair
502         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
503      endif
504 <    
504 >
505 >
506 >
507 >
508 >
509 > !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
510 >
511 >
512 >
513 >
514 >  
515   #ifdef IS_MPI
516      
517      if (update_nlist) then
267 –      
518         !! save current configuration, construct neighbor list,
519         !! and calculate forces
520         call saveNeighborList(nlocal, q)
# Line 273 | Line 523 | contains
523         nlist = 0      
524        
525         do i = 1, nrow
526 +
527            point(i) = nlist + 1
528            
529            inner: do j = 1, ncol
# Line 330 | Line 581 | contains
581   #else
582      
583      if (update_nlist) then
584 <      
584 >
585         ! save current configuration, contruct neighbor list,
586         ! and calculate forces
587         call saveNeighborList(natoms, q)
# Line 398 | Line 649 | contains
649   #endif
650      
651      ! phew, done with main loop.
652 <    
652 >
653 > !! Do timing
654 > #ifdef PROFILE
655 >    call cpu_time(forceTimeFinal)
656 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
657 > #endif
658 >
659 >
660   #ifdef IS_MPI
661      !!distribute forces
662    
# Line 460 | Line 718 | contains
718            end do
719   #endif
720            
721 <          do i = 1, getNlocal()
721 >          do i = 1, nLocal
722  
723               rfpot = 0.0_DP
724   #ifdef IS_MPI
# Line 510 | Line 768 | contains
768         tau = tau_Temp
769         virial = virial_Temp
770      endif
771 <
771 >    
772   #endif
773      
774 +    
775 +    
776    end subroutine do_force_loop
777  
778    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
779  
780      real( kind = dp ) :: pot
781 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
782 <    real (kind=dp), dimension(9,getNlocal()) :: A
783 <    real (kind=dp), dimension(3,getNlocal()) :: f
784 <    real (kind=dp), dimension(3,getNlocal()) :: t
781 >    real( kind = dp ), dimension(3,nLocal) :: u_l
782 >    real (kind=dp), dimension(9,nLocal) :: A
783 >    real (kind=dp), dimension(3,nLocal) :: f
784 >    real (kind=dp), dimension(3,nLocal) :: t
785  
786      logical, intent(inout) :: do_pot, do_stress
787      integer, intent(in) :: i, j
# Line 531 | Line 791 | contains
791      logical :: is_LJ_i, is_LJ_j
792      logical :: is_DP_i, is_DP_j
793      logical :: is_GB_i, is_GB_j
794 +    logical :: is_EAM_i,is_EAM_j
795      logical :: is_Sticky_i, is_Sticky_j
796      integer :: me_i, me_j
797 <
797 >    integer :: propPack_i
798 >    integer :: propPack_j
799      r = sqrt(rijsq)
800  
801   #ifdef IS_MPI
# Line 550 | Line 812 | contains
812      me_j = atid(j)
813  
814   #endif
815 <
815 >    
816      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)
817  
818 <       if ( is_LJ_i .and. is_LJ_j ) &
818 >       if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
819              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
820 +
821      endif
822  
823      if (FF_uses_dipoles .and. SimUsesDipoles()) then
824 <       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
825 <       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
565 <      
566 <       if ( is_DP_i .and. is_DP_j ) then
567 <          
824 >      
825 >       if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
826            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
827                 do_pot, do_stress)
828            if (FF_uses_RF .and. SimUsesRF()) then
# Line 577 | Line 835 | contains
835  
836      if (FF_uses_Sticky .and. SimUsesSticky()) then
837  
838 <       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
838 >       if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
839            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
840                 do_pot, do_stress)
841         endif
# Line 588 | Line 843 | contains
843  
844  
845      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)
846        
847 <       if ( is_GB_i .and. is_GB_j ) then
847 >       if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
848            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
849                 do_pot, do_stress)          
850         endif
851 +
852      endif
853      
854  
855 +  
856 +   if (FF_uses_EAM .and. SimUsesEAM()) then
857 +      
858 +      if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
859 +         call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
860 +      endif
861  
862 +   endif
863 +
864    end subroutine do_pair
865 +
866 +
867 +
868 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
869 +   real( kind = dp ) :: pot
870 +   real( kind = dp ), dimension(3,nLocal) :: u_l
871 +   real (kind=dp), dimension(9,nLocal) :: A
872 +   real (kind=dp), dimension(3,nLocal) :: f
873 +   real (kind=dp), dimension(3,nLocal) :: t
874 +  
875 +   logical, intent(inout) :: do_pot, do_stress
876 +   integer, intent(in) :: i, j
877 +   real ( kind = dp ), intent(inout)    :: rijsq
878 +   real ( kind = dp )                :: r
879 +   real ( kind = dp ), intent(inout) :: d(3)
880 +  
881 +   logical :: is_EAM_i, is_EAM_j
882 +  
883 +   integer :: me_i, me_j
884 +  
885 +   r = sqrt(rijsq)
886 +  
887 +
888 + #ifdef IS_MPI
889 +   if (tagRow(i) .eq. tagColumn(j)) then
890 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
891 +   endif
892 +  
893 +   me_i = atid_row(i)
894 +   me_j = atid_col(j)
895 +  
896 + #else
897 +  
898 +   me_i = atid(i)
899 +   me_j = atid(j)
900 +  
901 + #endif
902 +    
903 +   if (FF_uses_EAM .and. SimUsesEAM()) then
904 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
905 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
906 +      
907 +      if ( is_EAM_i .and. is_EAM_j ) &
908 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
909 +   endif
910  
911 + end subroutine do_prepair
912  
913 +
914 +
915 +
916 +  subroutine do_preforce(nlocal,pot)
917 +    integer :: nlocal
918 +    real( kind = dp ) :: pot
919 +
920 +    if (FF_uses_EAM .and. SimUsesEAM()) then
921 +       call calc_EAM_preforce_Frho(nlocal,pot)
922 +    endif
923 +
924 +
925 +  end subroutine do_preforce
926 +  
927 +  
928    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
929      
930      real (kind = dp), dimension(3) :: q_i
# Line 660 | Line 982 | contains
982      error = 0
983      ! Make sure we are properly initialized.
984      if (.not. do_forces_initialized) then
985 +       write(*,*) "Forces not initialized"
986         error = -1
987         return
988      endif
# Line 707 | Line 1030 | contains
1030  
1031   #endif
1032  
1033 +
1034 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1035 +       call clean_EAM()
1036 +    endif
1037 +
1038 +
1039 +
1040 +
1041 +
1042      rf = 0.0_dp
1043      tau_Temp = 0.0_dp
1044      virial_Temp = 0.0_dp
# Line 819 | Line 1151 | end module do_Forces
1151      doesit = FF_uses_RF
1152    end function FF_RequiresPostpairCalc
1153    
1154 + #ifdef PROFILE
1155 +  function getforcetime() result(totalforcetime)
1156 +    real(kind=dp) :: totalforcetime
1157 +    totalforcetime = forcetime
1158 +  end function getforcetime
1159 + #endif
1160 +
1161 + !! This cleans componets of force arrays belonging only to fortran
1162 +
1163   end module do_Forces

Diff Legend

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