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

Diff Legend

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