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 597 by mmeineke, Mon Jul 14 21:28:54 2003 UTC vs.
Revision 801 by chuckv, Sat Oct 4 18:46:12 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.18 2003-07-14 21:28:54 mmeineke Exp $, $Date: 2003-07-14 21:28:54 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
7 > !! @version $Id: do_Forces.F90,v 1.31 2003-10-04 18:46:12 chuckv Exp $, $Date: 2003-10-04 18:46:12 $, $Name: not supported by cvs2svn $, $Revision: 1.31 $
8  
9   module do_Forces
10    use force_globals
# Line 17 | Line 17 | module do_Forces
17    use dipole_dipole
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 27 | Line 30 | module do_Forces
30   #define __FORTRAN90
31   #include "fForceField.h"
32  
33 <  logical, save :: do_forces_initialized = .false.
33 >  logical, save :: do_forces_initialized = .false., haveRlist = .false.
34 >  logical, save :: havePolicies = .false.
35    logical, save :: FF_uses_LJ
36    logical, save :: FF_uses_sticky
37    logical, save :: FF_uses_dipoles
# Line 35 | Line 39 | module do_Forces
39    logical, save :: FF_uses_GB
40    logical, save :: FF_uses_EAM
41  
42 +  real(kind=dp), save :: rlist, rlistsq
43 +
44    public :: init_FF
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 )
59 +    
60 +    real(kind=dp) :: this_rlist
61 +
62 +    rlist = this_rlist
63 +    rlistsq = rlist * rlist
64 +    
65 +    haveRlist = .true.
66 +    if( havePolicies ) do_forces_initialized = .true.
67 +
68 +  end subroutine setRlistDF    
69 +
70    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
71  
72      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 114 | contains
114      !! check to make sure the FF_uses_RF setting makes sense
115      
116      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
117         if (FF_uses_RF) then
118            dielect = getDielect()
119 <          call initialize_rf(rrf, rt, dielect)
119 >          call initialize_rf(dielect)
120         endif
121      else
122         if (FF_uses_RF) then          
# Line 100 | Line 124 | contains
124            thisStat = -1
125            return
126         endif
127 <    endif
127 >    endif
128  
129      if (FF_uses_LJ) then
130        
107 –       call getRcut(rcut)
108 –
131         select case (LJMIXPOLICY)
132         case (LB_MIXING_RULE)
133 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
133 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
134         case (EXPLICIT_MIXING_RULE)
135 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
135 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
136         case default
137            write(default_error,*) 'unknown LJ Mixing Policy!'
138            thisStat = -1
# Line 129 | 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 149 | Line 183 | contains
183            return
184         endif
185      endif
186 +    
187  
188 <    do_forces_initialized = .true.    
189 <
188 >    havePolicies = .true.
189 >    if( haveRlist ) do_forces_initialized = .true.
190 >
191    end subroutine init_FF
192    
193  
# Line 179 | 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    
222      logical :: update_nlist  
223      integer :: i, j, jbeg, jend, jnab
224      integer :: nlist
225 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
225 >    real( kind = DP ) ::  rijsq
226      real(kind=dp),dimension(3) :: d
227      real(kind=dp) :: rfpot, mu_i, virial
228      integer :: me_i
# Line 194 | Line 231 | contains
231      integer :: listerror, error
232      integer :: localError
233  
234 +    real(kind=dp) :: listSkin = 1.0
235 +    
236 +
237      !! initialize local variables  
238  
239   #ifdef IS_MPI
# Line 205 | Line 245 | contains
245      nlocal = getNlocal()
246      natoms = nlocal
247   #endif
248 <  
209 <    call getRcut(rcut,rc2=rcutsq)
210 <    call getRlist(rlist,rlistsq)
211 <    
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 219 | 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 235 | 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, rcut, rlist, update_nlist)  
286 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
287         !! if_mpi_gather_stuff_for_prepair
288         !! do_prepair_loop_if_needed
289         !! if_mpi_scatter_stuff_from_prepair
290         !! if_mpi_gather_stuff_from_prepair_to_main_loop
291 <    else
292 <       !! See if we need to update neighbor lists
293 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
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 for non pre-pair
432 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
433 +    endif
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 268 | Line 462 | contains
462              
463               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
464              
465 <             if (rijsq <  rlistsq) then            
465 >             if (rijsq < rlistsq) then            
466                  
467                  nlist = nlist + 1
468                  
# Line 284 | Line 478 | contains
478                  
479                  list(nlist) = j
480                                  
481 <                if (rijsq <  rcutsq) then
482 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
483 <                        u_l, A, f, t, pot_local)
290 <                endif
481 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
482 >                     u_l, A, f, t, pot_local)
483 >                
484               endif
485            enddo inner
486         enddo
# Line 337 | Line 530 | contains
530               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
531            
532  
533 <             if (rijsq <  rlistsq) then
533 >             if (rijsq < rlistsq) then
534                  
535                  nlist = nlist + 1
536                
# Line 353 | Line 546 | contains
546                  
547                  list(nlist) = j
548                  
549 <                if (rijsq <  rcutsq) then
357 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
549 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
550                          u_l, A, f, t, pot)
551 <                endif
551 >                
552               endif
553            enddo inner
554         enddo
# Line 387 | 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 501 | Line 700 | contains
700      endif
701  
702   #endif
504 –    
505 –    write(*,*) 'T(1) = '
506 –    write(*,'(3es12.3)') t(1,1), t(1,2), t(1,3)
507 –    write(*,*)
703  
704 <    write(*,*) 'T(2) = '
705 <    write(*,'(3es12.3)') t(2,1), t(2,2), t(2,3)
511 <    write(*,*)
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 +
736    end subroutine do_force_loop
737  
738    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 528 | Line 751 | contains
751      logical :: is_LJ_i, is_LJ_j
752      logical :: is_DP_i, is_DP_j
753      logical :: is_GB_i, is_GB_j
754 +    logical :: is_EAM_i,is_EAM_j
755      logical :: is_Sticky_i, is_Sticky_j
756      integer :: me_i, me_j
757  
758      r = sqrt(rijsq)
535 –
536 –    write(*,*) 'ul(1) = '
537 –    write(*,'(3es12.3)') u_l(1,1), u_l(1,2), u_l(1,3)
538 –    write(*,*)
539 –
540 –    write(*,*) 'ul(2) = '
541 –    write(*,'(3es12.3)') u_l(2,1), u_l(2,2), u_l(2,3)
542 –    write(*,*)
759  
544 –
545 –    write(*,*) 'A(1) = '
546 –    write(*,'(3es12.3)') A(1,1), A(2,1), A(3,1)
547 –    write(*,'(3es12.3)') A(4,1), A(5,1), A(6,1)
548 –    write(*,'(3es12.3)') A(7,1), A(8,1), A(9,1)
549 –    write(*,*)
550 –    write(*,*) 'A(2) = '
551 –    write(*,'(3es12.3)') A(1,2), A(2,2), A(3,2)
552 –    write(*,'(3es12.3)') A(4,2), A(5,2), A(6,2)
553 –    write(*,'(3es12.3)') A(7,2), A(8,2), A(9,2)
554 –    write(*,*)
555 –
556 –
760   #ifdef IS_MPI
761      if (tagRow(i) .eq. tagColumn(j)) then
762         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 582 | Line 785 | contains
785         call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
786        
787         if ( is_DP_i .and. is_DP_j ) then
585 –          
788            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
789                 do_pot, do_stress)
790            if (FF_uses_RF .and. SimUsesRF()) then
# Line 607 | Line 809 | contains
809  
810      if (FF_uses_GB .and. SimUsesGB()) then
811  
812 +
813         call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
814         call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
815        
# Line 617 | Line 820 | contains
820      endif
821      
822  
823 +  
824 +   if (FF_uses_EAM .and. SimUsesEAM()) then
825 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
826 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
827 +      
828 +      if ( is_EAM_i .and. is_EAM_j ) &
829 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
830 +   endif
831 +
832  
833 +
834 +
835    end subroutine do_pair
836  
837  
838 +
839 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
840 +   real( kind = dp ) :: pot
841 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
842 +   real (kind=dp), dimension(9,getNlocal()) :: A
843 +   real (kind=dp), dimension(3,getNlocal()) :: f
844 +   real (kind=dp), dimension(3,getNlocal()) :: t
845 +  
846 +   logical, intent(inout) :: do_pot, do_stress
847 +   integer, intent(in) :: i, j
848 +   real ( kind = dp ), intent(inout)    :: rijsq
849 +   real ( kind = dp )                :: r
850 +   real ( kind = dp ), intent(inout) :: d(3)
851 +  
852 +   logical :: is_EAM_i, is_EAM_j
853 +  
854 +   integer :: me_i, me_j
855 +  
856 +   r = sqrt(rijsq)
857 +  
858 +
859 + #ifdef IS_MPI
860 +   if (tagRow(i) .eq. tagColumn(j)) then
861 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
862 +   endif
863 +  
864 +   me_i = atid_row(i)
865 +   me_j = atid_col(j)
866 +  
867 + #else
868 +  
869 +   me_i = atid(i)
870 +   me_j = atid(j)
871 +  
872 + #endif
873 +    
874 +   if (FF_uses_EAM .and. SimUsesEAM()) then
875 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
876 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
877 +      
878 +      if ( is_EAM_i .and. is_EAM_j ) &
879 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
880 +   endif
881 +
882 + end subroutine do_prepair
883 +
884 +
885 +
886 +
887 +  subroutine do_preforce(nlocal,pot)
888 +    integer :: nlocal
889 +    real( kind = dp ) :: pot
890 +
891 +    if (FF_uses_EAM .and. SimUsesEAM()) then
892 +       call calc_EAM_preforce_Frho(nlocal,pot)
893 +    endif
894 +
895 +
896 +  end subroutine do_preforce
897 +  
898 +  
899    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
900      
901      real (kind = dp), dimension(3) :: q_i
# Line 678 | Line 953 | contains
953      error = 0
954      ! Make sure we are properly initialized.
955      if (.not. do_forces_initialized) then
956 +       write(*,*) "Forces not initialized"
957         error = -1
958         return
959      endif
# Line 725 | Line 1001 | contains
1001  
1002   #endif
1003  
1004 +
1005 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1006 +       call clean_EAM()
1007 +    endif
1008 +
1009 +
1010 +
1011 +
1012 +
1013      rf = 0.0_dp
1014      tau_Temp = 0.0_dp
1015      virial_Temp = 0.0_dp
# Line 837 | Line 1122 | end module do_Forces
1122      doesit = FF_uses_RF
1123    end function FF_RequiresPostpairCalc
1124    
1125 + !! This cleans componets of force arrays belonging only to fortran
1126 +
1127   end module do_Forces

Diff Legend

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