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 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC vs.
Revision 669 by chuckv, Thu Aug 7 00:47:33 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.15 2003-04-11 15:16:59 gezelter Exp $, $Date: 2003-04-11 15:16:59 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $
7 > !! @version $Id: do_Forces.F90,v 1.27 2003-08-07 00:47:33 chuckv Exp $, $Date: 2003-08-07 00:47:33 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $
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   contains
49  
50 +  subroutine setRlistDF( this_rlist )
51 +    
52 +    real(kind=dp) :: this_rlist
53 +
54 +    rlist = this_rlist
55 +    rlistsq = rlist * rlist
56 +    
57 +    haveRlist = .true.
58 +    if( havePolicies ) do_forces_initialized = .true.
59 +
60 +  end subroutine setRlistDF    
61 +
62    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
63  
64      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 106 | contains
106      !! check to make sure the FF_uses_RF setting makes sense
107      
108      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
109         if (FF_uses_RF) then
110            dielect = getDielect()
111 <          call initialize_rf(rrf, rt, dielect)
111 >          call initialize_rf(dielect)
112         endif
113      else
114         if (FF_uses_RF) then          
# Line 100 | Line 116 | contains
116            thisStat = -1
117            return
118         endif
119 <    endif
119 >    endif
120  
121      if (FF_uses_LJ) then
122        
107 –       call getRcut(rcut)
108 –
123         select case (LJMIXPOLICY)
124         case (LB_MIXING_RULE)
125 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
125 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
126         case (EXPLICIT_MIXING_RULE)
127 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
127 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
128         case default
129            write(default_error,*) 'unknown LJ Mixing Policy!'
130            thisStat = -1
# Line 129 | Line 143 | contains
143            return
144         end if
145      endif
146 +
147 +
148 +    if (FF_uses_EAM) then
149 +       call init_EAM_FF(my_status)
150 +       if (my_status /= 0) then
151 +          thisStat = -1
152 +          return
153 +       end if
154 +    endif
155 +
156 +
157      
158      if (FF_uses_GB) then
159         call check_gb_pair_FF(my_status)
# Line 149 | Line 174 | contains
174            return
175         endif
176      endif
177 +    
178  
179 <    do_forces_initialized = .true.    
180 <
179 >    havePolicies = .true.
180 >    if( haveRlist ) do_forces_initialized = .true.
181 >
182    end subroutine init_FF
183    
184  
# Line 185 | Line 212 | contains
212      logical :: update_nlist  
213      integer :: i, j, jbeg, jend, jnab
214      integer :: nlist
215 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
215 >    real( kind = DP ) ::  rijsq
216      real(kind=dp),dimension(3) :: d
217      real(kind=dp) :: rfpot, mu_i, virial
218      integer :: me_i
# Line 194 | Line 221 | contains
221      integer :: listerror, error
222      integer :: localError
223  
224 +    real(kind=dp) :: listSkin = 1.0
225 +    
226 +
227      !! initialize local variables  
228  
229   #ifdef IS_MPI
# Line 205 | Line 235 | contains
235      nlocal = getNlocal()
236      natoms = nlocal
237   #endif
238 <  
209 <    call getRcut(rcut,rc2=rcutsq)
210 <    call getRlist(rlist,rlistsq)
211 <    
238 >
239      call check_initialization(localError)
240      if ( localError .ne. 0 ) then
241 +       call handleError("do_force_loop","Not Initialized")
242         error = -1
243         return
244      end if
# Line 219 | Line 247 | contains
247      do_pot = do_pot_c
248      do_stress = do_stress_c
249  
250 +
251      ! Gather all information needed by all force loops:
252      
253   #ifdef IS_MPI    
# Line 235 | Line 264 | contains
264      endif
265      
266   #endif
267 <    
267 >  
268      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
269         !! See if we need to update neighbor lists
270 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
270 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
271         !! if_mpi_gather_stuff_for_prepair
272         !! do_prepair_loop_if_needed
273         !! if_mpi_scatter_stuff_from_prepair
274         !! if_mpi_gather_stuff_from_prepair_to_main_loop
275 <    else
276 <       !! See if we need to update neighbor lists
277 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
275 >
276 > !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
277 > #ifdef IS_MPI
278 >    
279 >    if (update_nlist) then
280 >      
281 >       !! save current configuration, construct neighbor list,
282 >       !! and calculate forces
283 >       call saveNeighborList(nlocal, q)
284 >      
285 >       neighborListSize = size(list)
286 >       nlist = 0      
287 >      
288 >       do i = 1, nrow
289 >          point(i) = nlist + 1
290 >          
291 >          prepair_inner: do j = 1, ncol
292 >            
293 >             if (skipThisPair(i,j)) cycle prepair_inner
294 >            
295 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
296 >            
297 >             if (rijsq < rlistsq) then            
298 >                
299 >                nlist = nlist + 1
300 >                
301 >                if (nlist > neighborListSize) then
302 >                   call expandNeighborList(nlocal, listerror)
303 >                   if (listerror /= 0) then
304 >                      error = -1
305 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
306 >                      return
307 >                   end if
308 >                   neighborListSize = size(list)
309 >                endif
310 >                
311 >                list(nlist) = j
312 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
313 >             endif
314 >          enddo prepair_inner
315 >       enddo
316 >
317 >       point(nrow + 1) = nlist + 1
318 >      
319 >    else  !! (of update_check)
320 >
321 >       ! use the list to find the neighbors
322 >       do i = 1, nrow
323 >          JBEG = POINT(i)
324 >          JEND = POINT(i+1) - 1
325 >          ! check thiat molecule i has neighbors
326 >          if (jbeg .le. jend) then
327 >            
328 >             do jnab = jbeg, jend
329 >                j = list(jnab)
330 >
331 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
332 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
333 >                     u_l, A, f, t, pot_local)
334 >
335 >             enddo
336 >          endif
337 >       enddo
338      endif
339      
340 + #else
341 +    
342 +    if (update_nlist) then
343 +      
344 +       ! save current configuration, contruct neighbor list,
345 +       ! and calculate forces
346 +       call saveNeighborList(natoms, q)
347 +      
348 +       neighborListSize = size(list)
349 +  
350 +       nlist = 0
351 +      
352 +       do i = 1, natoms-1
353 +          point(i) = nlist + 1
354 +          
355 +          prepair_inner: do j = i+1, natoms
356 +            
357 +             if (skipThisPair(i,j))  cycle prepair_inner
358 +                          
359 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
360 +          
361 +
362 +             if (rijsq < rlistsq) then
363 +                
364 +                nlist = nlist + 1
365 +              
366 +                if (nlist > neighborListSize) then
367 +                   call expandNeighborList(natoms, listerror)
368 +                   if (listerror /= 0) then
369 +                      error = -1
370 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
371 +                      return
372 +                   end if
373 +                   neighborListSize = size(list)
374 +                endif
375 +                
376 +                list(nlist) = j
377 +                
378 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
379 +                        u_l, A, f, t, pot)
380 +                
381 +             endif
382 +          enddo prepair_inner
383 +       enddo
384 +      
385 +       point(natoms) = nlist + 1
386 +      
387 +    else !! (update)
388 +      
389 +       ! use the list to find the neighbors
390 +       do i = 1, natoms-1
391 +          JBEG = POINT(i)
392 +          JEND = POINT(i+1) - 1
393 +          ! check thiat molecule i has neighbors
394 +          if (jbeg .le. jend) then
395 +            
396 +             do jnab = jbeg, jend
397 +                j = list(jnab)
398 +
399 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
400 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
401 +                     u_l, A, f, t, pot)
402 +
403 +             enddo
404 +          endif
405 +       enddo
406 +    endif    
407 + #endif
408 +    !! Do rest of preforce calculations
409 +   call do_preforce(nlocal,pot)
410 +    else
411 +       !! See if we need to update neighbor lists for non pre-pair
412 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
413 +    endif
414 +
415 +
416 +
417 +
418 +
419 + !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
420 +
421 +
422 +
423 +
424 +  
425   #ifdef IS_MPI
426      
427      if (update_nlist) then
# Line 268 | Line 442 | contains
442              
443               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
444              
445 <             if (rijsq <  rlistsq) then            
445 >             if (rijsq < rlistsq) then            
446                  
447                  nlist = nlist + 1
448                  
# Line 284 | Line 458 | contains
458                  
459                  list(nlist) = j
460                                  
461 <                if (rijsq <  rcutsq) then
462 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
463 <                        u_l, A, f, t, pot_local)
290 <                endif
461 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
462 >                     u_l, A, f, t, pot_local)
463 >                
464               endif
465            enddo inner
466         enddo
# Line 337 | Line 510 | contains
510               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
511            
512  
513 <             if (rijsq <  rlistsq) then
513 >             if (rijsq < rlistsq) then
514                  
515                  nlist = nlist + 1
516                
# Line 353 | Line 526 | contains
526                  
527                  list(nlist) = j
528                  
529 <                if (rijsq <  rcutsq) then
357 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
529 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
530                          u_l, A, f, t, pot)
531 <                endif
531 >                
532               endif
533            enddo inner
534         enddo
# Line 520 | Line 692 | contains
692      logical :: is_LJ_i, is_LJ_j
693      logical :: is_DP_i, is_DP_j
694      logical :: is_GB_i, is_GB_j
695 +    logical :: is_EAM_i,is_EAM_j
696      logical :: is_Sticky_i, is_Sticky_j
697      integer :: me_i, me_j
698  
699      r = sqrt(rijsq)
700  
528 –
529 –
701   #ifdef IS_MPI
702      if (tagRow(i) .eq. tagColumn(j)) then
703         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 589 | Line 760 | contains
760         endif
761      endif
762      
763 +
764 +  
765 +   if (FF_uses_EAM .and. SimUsesEAM()) then
766 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
767 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
768 +      
769 +      if ( is_EAM_i .and. is_EAM_j ) &
770 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
771 +   endif
772 +
773 +
774 +
775 +
776    end subroutine do_pair
777  
778  
779 +
780 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
781 +   real( kind = dp ) :: pot
782 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
783 +   real (kind=dp), dimension(9,getNlocal()) :: A
784 +   real (kind=dp), dimension(3,getNlocal()) :: f
785 +   real (kind=dp), dimension(3,getNlocal()) :: t
786 +  
787 +   logical, intent(inout) :: do_pot, do_stress
788 +   integer, intent(in) :: i, j
789 +   real ( kind = dp ), intent(inout)    :: rijsq
790 +   real ( kind = dp )                :: r
791 +   real ( kind = dp ), intent(inout) :: d(3)
792 +  
793 +   logical :: is_EAM_i, is_EAM_j
794 +  
795 +   integer :: me_i, me_j
796 +  
797 +   r = sqrt(rijsq)
798 +  
799 +
800 + #ifdef IS_MPI
801 +   if (tagRow(i) .eq. tagColumn(j)) then
802 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
803 +   endif
804 +  
805 +   me_i = atid_row(i)
806 +   me_j = atid_col(j)
807 +  
808 + #else
809 +  
810 +   me_i = atid(i)
811 +   me_j = atid(j)
812 +  
813 + #endif
814 +  
815 +   if (FF_uses_EAM .and. SimUsesEAM()) then
816 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
817 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
818 +      
819 +      if ( is_EAM_i .and. is_EAM_j ) &
820 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
821 +   endif
822 +  end subroutine do_prepair
823 +
824 +
825 +
826 +
827 +  subroutine do_preforce(nlocal,pot)
828 +    integer :: nlocal
829 +    real( kind = dp ) :: pot
830 +
831 +    if (FF_uses_EAM .and. SimUsesEAM()) then
832 +       call calc_EAM_preforce_Frho(nlocal,pot)
833 +    endif
834 +
835 +
836 +  end subroutine do_preforce
837 +  
838 +  
839    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
840      
841      real (kind = dp), dimension(3) :: q_i
842      real (kind = dp), dimension(3) :: q_j
843      real ( kind = dp ), intent(out) :: r_sq
844 <    real( kind = dp ) :: d(3)
845 <    real( kind = dp ) :: d_old(3)
844 >    real( kind = dp ) :: d(3), scaled(3)
845 >    integer i
846 >
847      d(1:3) = q_j(1:3) - q_i(1:3)
848 <    d_old = d
848 >
849      ! Wrap back into periodic box if necessary
850      if ( SimUsesPBC() ) then
851        
852 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
853 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
852 >       if( .not.boxIsOrthorhombic ) then
853 >          ! calc the scaled coordinates.
854 >          
855 >          scaled = matmul(HmatInv, d)
856 >          
857 >          ! wrap the scaled coordinates
858 >
859 >          scaled = scaled  - anint(scaled)
860 >          
861 >
862 >          ! calc the wrapped real coordinates from the wrapped scaled
863 >          ! coordinates
864 >
865 >          d = matmul(Hmat,scaled)
866 >
867 >       else
868 >          ! calc the scaled coordinates.
869 >          
870 >          do i = 1, 3
871 >             scaled(i) = d(i) * HmatInv(i,i)
872 >            
873 >             ! wrap the scaled coordinates
874 >            
875 >             scaled(i) = scaled(i) - anint(scaled(i))
876 >            
877 >             ! calc the wrapped real coordinates from the wrapped scaled
878 >             ! coordinates
879 >
880 >             d(i) = scaled(i)*Hmat(i,i)
881 >          enddo
882 >       endif
883        
884      endif
885 +    
886      r_sq = dot_product(d,d)
887 <        
887 >    
888    end subroutine get_interatomic_vector
889 <
889 >  
890    subroutine check_initialization(error)
891      integer, intent(out) :: error
892      
893      error = 0
894      ! Make sure we are properly initialized.
895      if (.not. do_forces_initialized) then
896 +       write(*,*) "Forces not initialized"
897         error = -1
898         return
899      endif

Diff Legend

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