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 571 by gezelter, Tue Jul 1 22:39:53 2003 UTC vs.
Revision 843 by mmeineke, Wed Oct 29 20:41:39 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.16 2003-07-01 22:39:53 gezelter Exp $, $Date: 2003-07-01 22:39:53 $, $Name: not supported by cvs2svn $, $Revision: 1.16 $
7 > !! @version $Id: do_Forces.F90,v 1.35 2003-10-29 20:41:39 mmeineke Exp $, $Date: 2003-10-29 20:41:39 $, $Name: not supported by cvs2svn $, $Revision: 1.35 $
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
703 <    
703 >
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    end subroutine do_force_loop
735  
736    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 520 | Line 749 | contains
749      logical :: is_LJ_i, is_LJ_j
750      logical :: is_DP_i, is_DP_j
751      logical :: is_GB_i, is_GB_j
752 +    logical :: is_EAM_i,is_EAM_j
753      logical :: is_Sticky_i, is_Sticky_j
754      integer :: me_i, me_j
755  
756      r = sqrt(rijsq)
757  
528 –
529 –
758   #ifdef IS_MPI
759      if (tagRow(i) .eq. tagColumn(j)) then
760         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 555 | Line 783 | contains
783         call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
784        
785         if ( is_DP_i .and. is_DP_j ) then
558 –          
786            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
787                 do_pot, do_stress)
788            if (FF_uses_RF .and. SimUsesRF()) then
# Line 580 | Line 807 | contains
807  
808      if (FF_uses_GB .and. SimUsesGB()) then
809  
810 +
811         call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
812         call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
813        
# Line 589 | Line 817 | contains
817         endif
818      endif
819      
820 +
821 +  
822 +   if (FF_uses_EAM .and. SimUsesEAM()) then
823 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
824 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
825 +      
826 +      if ( is_EAM_i .and. is_EAM_j ) &
827 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
828 +   endif
829 +
830 +
831 +
832 +
833    end subroutine do_pair
834 +
835 +
836 +
837 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
838 +   real( kind = dp ) :: pot
839 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
840 +   real (kind=dp), dimension(9,getNlocal()) :: A
841 +   real (kind=dp), dimension(3,getNlocal()) :: f
842 +   real (kind=dp), dimension(3,getNlocal()) :: t
843 +  
844 +   logical, intent(inout) :: do_pot, do_stress
845 +   integer, intent(in) :: i, j
846 +   real ( kind = dp ), intent(inout)    :: rijsq
847 +   real ( kind = dp )                :: r
848 +   real ( kind = dp ), intent(inout) :: d(3)
849 +  
850 +   logical :: is_EAM_i, is_EAM_j
851 +  
852 +   integer :: me_i, me_j
853 +  
854 +   r = sqrt(rijsq)
855 +  
856 +
857 + #ifdef IS_MPI
858 +   if (tagRow(i) .eq. tagColumn(j)) then
859 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
860 +   endif
861 +  
862 +   me_i = atid_row(i)
863 +   me_j = atid_col(j)
864 +  
865 + #else
866 +  
867 +   me_i = atid(i)
868 +   me_j = atid(j)
869 +  
870 + #endif
871 +    
872 +   if (FF_uses_EAM .and. SimUsesEAM()) then
873 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
874 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
875 +      
876 +      if ( is_EAM_i .and. is_EAM_j ) &
877 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
878 +   endif
879  
880 + end subroutine do_prepair
881  
882 +
883 +
884 +
885 +  subroutine do_preforce(nlocal,pot)
886 +    integer :: nlocal
887 +    real( kind = dp ) :: pot
888 +
889 +    if (FF_uses_EAM .and. SimUsesEAM()) then
890 +       call calc_EAM_preforce_Frho(nlocal,pot)
891 +    endif
892 +
893 +
894 +  end subroutine do_preforce
895 +  
896 +  
897    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
898      
899      real (kind = dp), dimension(3) :: q_i
# Line 608 | Line 910 | contains
910         if( .not.boxIsOrthorhombic ) then
911            ! calc the scaled coordinates.
912            
913 <          scaled = matmul(d, HmatInv)
913 >          scaled = matmul(HmatInv, d)
914            
915            ! wrap the scaled coordinates
916  
917 <          do i = 1, 3
918 <             scaled(i) = scaled(i) - anint(scaled(i))
617 <          enddo
917 >          scaled = scaled  - anint(scaled)
918 >          
919  
920            ! calc the wrapped real coordinates from the wrapped scaled
921            ! coordinates
922  
923 <          d = matmul(scaled,Hmat)
923 >          d = matmul(Hmat,scaled)
924  
925         else
926            ! calc the scaled coordinates.
# Line 650 | Line 951 | contains
951      error = 0
952      ! Make sure we are properly initialized.
953      if (.not. do_forces_initialized) then
954 +       write(*,*) "Forces not initialized"
955         error = -1
956         return
957      endif
# Line 697 | Line 999 | contains
999  
1000   #endif
1001  
1002 +
1003 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1004 +       call clean_EAM()
1005 +    endif
1006 +
1007 +
1008 +
1009 +
1010 +
1011      rf = 0.0_dp
1012      tau_Temp = 0.0_dp
1013      virial_Temp = 0.0_dp
# Line 809 | Line 1120 | end module do_Forces
1120      doesit = FF_uses_RF
1121    end function FF_RequiresPostpairCalc
1122    
1123 + !! This cleans componets of force arrays belonging only to fortran
1124 +
1125   end module do_Forces

Diff Legend

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