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 572 by mmeineke, Wed Jul 2 21:26:55 2003 UTC vs.
Revision 891 by mmeineke, Fri Dec 19 20:36:35 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.17 2003-07-02 21:26:55 mmeineke Exp $, $Date: 2003-07-02 21:26:55 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $
7 > !! @version $Id: do_Forces.F90,v 1.40 2003-12-19 20:36:35 mmeineke Exp $, $Date: 2003-12-19 20:36:35 $, $Name: not supported by cvs2svn $, $Revision: 1.40 $
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 +  public :: getforcetime
50 +  real, save :: forceTime = 0
51 +  real :: forceTimeInitial, forceTimeFinal
52 +  integer :: nLoops
53 + #endif
54 +
55   contains
56  
57 +  subroutine setRlistDF( this_rlist )
58 +    
59 +    real(kind=dp) :: this_rlist
60 +
61 +    rlist = this_rlist
62 +    rlistsq = rlist * rlist
63 +    
64 +    haveRlist = .true.
65 +    if( havePolicies ) do_forces_initialized = .true.
66 +
67 +  end subroutine setRlistDF    
68 +
69    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
70  
71      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 113 | contains
113      !! check to make sure the FF_uses_RF setting makes sense
114      
115      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
116         if (FF_uses_RF) then
117            dielect = getDielect()
118 <          call initialize_rf(rrf, rt, dielect)
118 >          call initialize_rf(dielect)
119         endif
120      else
121         if (FF_uses_RF) then          
# Line 100 | Line 123 | contains
123            thisStat = -1
124            return
125         endif
126 <    endif
126 >    endif
127  
128      if (FF_uses_LJ) then
129        
107 –       call getRcut(rcut)
108 –
130         select case (LJMIXPOLICY)
131         case (LB_MIXING_RULE)
132 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
132 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
133         case (EXPLICIT_MIXING_RULE)
134 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
134 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
135         case default
136            write(default_error,*) 'unknown LJ Mixing Policy!'
137            thisStat = -1
# Line 129 | Line 150 | contains
150            return
151         end if
152      endif
153 +
154 +
155 +    if (FF_uses_EAM) then
156 +         call init_EAM_FF(my_status)
157 +       if (my_status /= 0) then
158 +          write(*,*) "init_EAM_FF returned a bad status"
159 +          thisStat = -1
160 +          return
161 +       end if
162 +    endif
163 +
164 +
165      
166      if (FF_uses_GB) then
167         call check_gb_pair_FF(my_status)
# Line 149 | Line 182 | contains
182            return
183         endif
184      endif
185 +    
186  
187 <    do_forces_initialized = .true.    
188 <
187 >    havePolicies = .true.
188 >    if( haveRlist ) do_forces_initialized = .true.
189 >
190    end subroutine init_FF
191    
192  
# Line 179 | Line 214 | contains
214      real( kind = DP ) :: pot_local
215      integer :: nrow
216      integer :: ncol
217 +    integer :: nprocs
218   #endif
219      integer :: nlocal
220      integer :: natoms    
221      logical :: update_nlist  
222      integer :: i, j, jbeg, jend, jnab
223      integer :: nlist
224 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
224 >    real( kind = DP ) ::  rijsq
225      real(kind=dp),dimension(3) :: d
226      real(kind=dp) :: rfpot, mu_i, virial
227      integer :: me_i
# Line 194 | Line 230 | contains
230      integer :: listerror, error
231      integer :: localError
232  
233 +    real(kind=dp) :: listSkin = 1.0  
234 +
235      !! initialize local variables  
236  
237   #ifdef IS_MPI
# Line 205 | Line 243 | contains
243      nlocal = getNlocal()
244      natoms = nlocal
245   #endif
246 <  
209 <    call getRcut(rcut,rc2=rcutsq)
210 <    call getRlist(rlist,rlistsq)
211 <    
246 >
247      call check_initialization(localError)
248      if ( localError .ne. 0 ) then
249 +       call handleError("do_force_loop","Not Initialized")
250         error = -1
251         return
252      end if
# Line 219 | Line 255 | contains
255      do_pot = do_pot_c
256      do_stress = do_stress_c
257  
258 +
259      ! Gather all information needed by all force loops:
260      
261   #ifdef IS_MPI    
# Line 235 | Line 272 | contains
272      endif
273      
274   #endif
275 <    
275 >
276 > !! Begin force loop timing:
277 > #ifdef PROFILE
278 >    call cpu_time(forceTimeInitial)
279 >    nloops = nloops + 1
280 > #endif
281 >  
282      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
283         !! See if we need to update neighbor lists
284 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
284 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
285         !! if_mpi_gather_stuff_for_prepair
286         !! do_prepair_loop_if_needed
287         !! if_mpi_scatter_stuff_from_prepair
288         !! if_mpi_gather_stuff_from_prepair_to_main_loop
289 <    else
290 <       !! See if we need to update neighbor lists
291 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
289 >    
290 > !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
291 > #ifdef IS_MPI
292 >    
293 >    if (update_nlist) then
294 >      
295 >       !! save current configuration, construct neighbor list,
296 >       !! and calculate forces
297 >       call saveNeighborList(nlocal, q)
298 >      
299 >       neighborListSize = size(list)
300 >       nlist = 0      
301 >      
302 >       do i = 1, nrow
303 >          point(i) = nlist + 1
304 >          
305 >          prepair_inner: do j = 1, ncol
306 >            
307 >             if (skipThisPair(i,j)) cycle prepair_inner
308 >            
309 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
310 >            
311 >             if (rijsq < rlistsq) then            
312 >                
313 >                nlist = nlist + 1
314 >                
315 >                if (nlist > neighborListSize) then
316 >                   call expandNeighborList(nlocal, listerror)
317 >                   if (listerror /= 0) then
318 >                      error = -1
319 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
320 >                      return
321 >                   end if
322 >                   neighborListSize = size(list)
323 >                endif
324 >                
325 >                list(nlist) = j
326 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
327 >             endif
328 >          enddo prepair_inner
329 >       enddo
330 >
331 >       point(nrow + 1) = nlist + 1
332 >      
333 >    else  !! (of update_check)
334 >
335 >       ! use the list to find the neighbors
336 >       do i = 1, nrow
337 >          JBEG = POINT(i)
338 >          JEND = POINT(i+1) - 1
339 >          ! check thiat molecule i has neighbors
340 >          if (jbeg .le. jend) then
341 >            
342 >             do jnab = jbeg, jend
343 >                j = list(jnab)
344 >
345 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
346 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
347 >                     u_l, A, f, t, pot_local)
348 >
349 >             enddo
350 >          endif
351 >       enddo
352      endif
353      
354 < #ifdef IS_MPI
354 > #else
355      
356      if (update_nlist) then
357 +      
358 +       ! save current configuration, contruct neighbor list,
359 +       ! and calculate forces
360 +       call saveNeighborList(natoms, q)
361 +      
362 +       neighborListSize = size(list)
363 +  
364 +       nlist = 0
365 +
366 +       do i = 1, natoms-1
367 +          point(i) = nlist + 1
368 +          
369 +          prepair_inner: do j = i+1, natoms
370 +            
371 +             if (skipThisPair(i,j))  cycle prepair_inner
372 +                          
373 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
374 +          
375 +
376 +             if (rijsq < rlistsq) then
377 +
378 +          
379 +                nlist = nlist + 1
380 +              
381 +                if (nlist > neighborListSize) then
382 +                   call expandNeighborList(natoms, listerror)
383 +                   if (listerror /= 0) then
384 +                      error = -1
385 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
386 +                      return
387 +                   end if
388 +                   neighborListSize = size(list)
389 +                endif
390 +                
391 +                list(nlist) = j
392 +                
393 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
394 +                        u_l, A, f, t, pot)
395 +                
396 +             endif
397 +          enddo prepair_inner
398 +       enddo
399        
400 +       point(natoms) = nlist + 1
401 +      
402 +    else !! (update)
403 +  
404 +       ! use the list to find the neighbors
405 +       do i = 1, natoms-1
406 +          JBEG = POINT(i)
407 +          JEND = POINT(i+1) - 1
408 +          ! check thiat molecule i has neighbors
409 +          if (jbeg .le. jend) then
410 +            
411 +             do jnab = jbeg, jend
412 +                j = list(jnab)
413 +
414 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
415 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
416 +                     u_l, A, f, t, pot)
417 +
418 +             enddo
419 +          endif
420 +       enddo
421 +    endif    
422 + #endif
423 +    !! Do rest of preforce calculations
424 +    !! do necessary preforce calculations  
425 +    call do_preforce(nlocal,pot)
426 +   ! we have already updated the neighbor list set it to false...
427 +   update_nlist = .false.
428 +    else
429 +       !! See if we need to update neighbor lists for non pre-pair
430 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
431 +    endif
432 +
433 +
434 +
435 +
436 +
437 + !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
438 +
439 +
440 +
441 +
442 +  
443 + #ifdef IS_MPI
444 +    
445 +    if (update_nlist) then
446         !! save current configuration, construct neighbor list,
447         !! and calculate forces
448         call saveNeighborList(nlocal, q)
# Line 268 | Line 459 | contains
459              
460               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
461              
462 <             if (rijsq <  rlistsq) then            
462 >             if (rijsq < rlistsq) then            
463                  
464                  nlist = nlist + 1
465                  
# Line 284 | Line 475 | contains
475                  
476                  list(nlist) = j
477                                  
478 <                if (rijsq <  rcutsq) then
479 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
480 <                        u_l, A, f, t, pot_local)
290 <                endif
478 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
479 >                     u_l, A, f, t, pot_local)
480 >                
481               endif
482            enddo inner
483         enddo
# Line 318 | Line 508 | contains
508   #else
509      
510      if (update_nlist) then
511 <      
511 >
512         ! save current configuration, contruct neighbor list,
513         ! and calculate forces
514         call saveNeighborList(natoms, q)
# Line 337 | Line 527 | contains
527               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
528            
529  
530 <             if (rijsq <  rlistsq) then
530 >             if (rijsq < rlistsq) then
531                  
532                  nlist = nlist + 1
533                
# Line 353 | Line 543 | contains
543                  
544                  list(nlist) = j
545                  
546 <                if (rijsq <  rcutsq) then
357 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
546 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
547                          u_l, A, f, t, pot)
548 <                endif
548 >                
549               endif
550            enddo inner
551         enddo
# Line 387 | Line 576 | contains
576   #endif
577      
578      ! phew, done with main loop.
579 <    
579 >
580 > !! Do timing
581 > #ifdef PROFILE
582 >    call cpu_time(forceTimeFinal)
583 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
584 > #endif
585 >
586 >
587   #ifdef IS_MPI
588      !!distribute forces
589    
# Line 499 | Line 695 | contains
695         tau = tau_Temp
696         virial = virial_Temp
697      endif
698 <
698 >    
699   #endif
700      
701 +    
702 +    
703    end subroutine do_force_loop
704  
705    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 520 | Line 718 | contains
718      logical :: is_LJ_i, is_LJ_j
719      logical :: is_DP_i, is_DP_j
720      logical :: is_GB_i, is_GB_j
721 +    logical :: is_EAM_i,is_EAM_j
722      logical :: is_Sticky_i, is_Sticky_j
723      integer :: me_i, me_j
724  
725      r = sqrt(rijsq)
726  
528 –
529 –
727   #ifdef IS_MPI
728      if (tagRow(i) .eq. tagColumn(j)) then
729         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 555 | Line 752 | contains
752         call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
753        
754         if ( is_DP_i .and. is_DP_j ) then
558 –          
755            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
756                 do_pot, do_stress)
757            if (FF_uses_RF .and. SimUsesRF()) then
# Line 580 | Line 776 | contains
776  
777      if (FF_uses_GB .and. SimUsesGB()) then
778  
779 +
780         call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
781         call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
782        
# Line 589 | Line 786 | contains
786         endif
787      endif
788      
789 +
790 +  
791 +   if (FF_uses_EAM .and. SimUsesEAM()) then
792 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
793 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
794 +      
795 +      if ( is_EAM_i .and. is_EAM_j ) &
796 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
797 +   endif
798 +
799 +
800 +
801 +
802    end subroutine do_pair
803 +
804 +
805 +
806 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
807 +   real( kind = dp ) :: pot
808 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
809 +   real (kind=dp), dimension(9,getNlocal()) :: A
810 +   real (kind=dp), dimension(3,getNlocal()) :: f
811 +   real (kind=dp), dimension(3,getNlocal()) :: t
812 +  
813 +   logical, intent(inout) :: do_pot, do_stress
814 +   integer, intent(in) :: i, j
815 +   real ( kind = dp ), intent(inout)    :: rijsq
816 +   real ( kind = dp )                :: r
817 +   real ( kind = dp ), intent(inout) :: d(3)
818 +  
819 +   logical :: is_EAM_i, is_EAM_j
820 +  
821 +   integer :: me_i, me_j
822 +  
823 +   r = sqrt(rijsq)
824 +  
825 +
826 + #ifdef IS_MPI
827 +   if (tagRow(i) .eq. tagColumn(j)) then
828 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
829 +   endif
830 +  
831 +   me_i = atid_row(i)
832 +   me_j = atid_col(j)
833 +  
834 + #else
835 +  
836 +   me_i = atid(i)
837 +   me_j = atid(j)
838 +  
839 + #endif
840 +    
841 +   if (FF_uses_EAM .and. SimUsesEAM()) then
842 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
843 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
844 +      
845 +      if ( is_EAM_i .and. is_EAM_j ) &
846 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
847 +   endif
848  
849 + end subroutine do_prepair
850  
851 +
852 +
853 +
854 +  subroutine do_preforce(nlocal,pot)
855 +    integer :: nlocal
856 +    real( kind = dp ) :: pot
857 +
858 +    if (FF_uses_EAM .and. SimUsesEAM()) then
859 +       call calc_EAM_preforce_Frho(nlocal,pot)
860 +    endif
861 +
862 +
863 +  end subroutine do_preforce
864 +  
865 +  
866    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
867      
868      real (kind = dp), dimension(3) :: q_i
# Line 649 | Line 920 | contains
920      error = 0
921      ! Make sure we are properly initialized.
922      if (.not. do_forces_initialized) then
923 +       write(*,*) "Forces not initialized"
924         error = -1
925         return
926      endif
# Line 696 | Line 968 | contains
968  
969   #endif
970  
971 +
972 +    if (FF_uses_EAM .and. SimUsesEAM()) then
973 +       call clean_EAM()
974 +    endif
975 +
976 +
977 +
978 +
979 +
980      rf = 0.0_dp
981      tau_Temp = 0.0_dp
982      virial_Temp = 0.0_dp
# Line 808 | Line 1089 | end module do_Forces
1089      doesit = FF_uses_RF
1090    end function FF_RequiresPostpairCalc
1091    
1092 + #ifdef PROFILE
1093 +  function getforcetime() result(totalforcetime)
1094 +    real(kind=dp) :: totalforcetime
1095 +    totalforcetime = forcetime
1096 +  end function getforcetime
1097 + #endif
1098 +
1099 + !! This cleans componets of force arrays belonging only to fortran
1100 +
1101   end module do_Forces

Diff Legend

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