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 648 by chuckv, Wed Jul 23 22:13:59 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.24 2003-07-23 22:13:59 chuckv Exp $, $Date: 2003-07-23 22:13:59 $, $Name: not supported by cvs2svn $, $Revision: 1.24 $
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   #ifdef IS_MPI
22    use mpiSimulation
23   #endif
# Line 27 | Line 28 | module do_Forces
28   #define __FORTRAN90
29   #include "fForceField.h"
30  
31 <  logical, save :: do_forces_initialized = .false.
31 >  logical, save :: do_forces_initialized = .false., haveRlist = .false.
32 >  logical, save :: havePolicies = .false.
33    logical, save :: FF_uses_LJ
34    logical, save :: FF_uses_sticky
35    logical, save :: FF_uses_dipoles
# Line 35 | Line 37 | module do_Forces
37    logical, save :: FF_uses_GB
38    logical, save :: FF_uses_EAM
39  
40 +  real(kind=dp), save :: rlist, rlistsq
41 +
42    public :: init_FF
43    public :: do_force_loop
44 +  public :: setRlistDF
45  
46   contains
47  
48 +  subroutine setRlistDF( this_rlist )
49 +    
50 +    real(kind=dp) :: this_rlist
51 +
52 +    rlist = this_rlist
53 +    rlistsq = rlist * rlist
54 +    
55 +    haveRlist = .true.
56 +    if( havePolicies ) do_forces_initialized = .true.
57 +
58 +  end subroutine setRlistDF    
59 +
60    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
61  
62      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 104 | contains
104      !! check to make sure the FF_uses_RF setting makes sense
105      
106      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
107         if (FF_uses_RF) then
108            dielect = getDielect()
109 <          call initialize_rf(rrf, rt, dielect)
109 >          call initialize_rf(dielect)
110         endif
111      else
112         if (FF_uses_RF) then          
# Line 100 | Line 114 | contains
114            thisStat = -1
115            return
116         endif
117 <    endif
117 >    endif
118  
119      if (FF_uses_LJ) then
120        
107 –       call getRcut(rcut)
108 –
121         select case (LJMIXPOLICY)
122         case (LB_MIXING_RULE)
123 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
123 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
124         case (EXPLICIT_MIXING_RULE)
125 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
125 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
126         case default
127            write(default_error,*) 'unknown LJ Mixing Policy!'
128            thisStat = -1
# Line 150 | Line 162 | contains
162         endif
163      endif
164  
165 <    do_forces_initialized = .true.    
166 <
165 >    havePolicies = .true.
166 >    if( haveRlist ) do_forces_initialized = .true.
167 >    
168    end subroutine init_FF
169    
170  
# Line 185 | Line 198 | contains
198      logical :: update_nlist  
199      integer :: i, j, jbeg, jend, jnab
200      integer :: nlist
201 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
201 >    real( kind = DP ) ::  rijsq
202      real(kind=dp),dimension(3) :: d
203      real(kind=dp) :: rfpot, mu_i, virial
204      integer :: me_i
# Line 194 | Line 207 | contains
207      integer :: listerror, error
208      integer :: localError
209  
210 +    real(kind=dp) :: listSkin = 1.0
211 +    
212 +
213      !! initialize local variables  
214  
215   #ifdef IS_MPI
# Line 206 | Line 222 | contains
222      natoms = nlocal
223   #endif
224    
209 –    call getRcut(rcut,rc2=rcutsq)
210 –    call getRlist(rlist,rlistsq)
211 –    
225      call check_initialization(localError)
226      if ( localError .ne. 0 ) then
227         error = -1
# Line 238 | Line 251 | contains
251      
252      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
253         !! See if we need to update neighbor lists
254 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
254 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
255         !! if_mpi_gather_stuff_for_prepair
256         !! do_prepair_loop_if_needed
257         !! if_mpi_scatter_stuff_from_prepair
258         !! if_mpi_gather_stuff_from_prepair_to_main_loop
259 <    else
260 <       !! See if we need to update neighbor lists
261 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
259 >
260 > !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
261 > #ifdef IS_MPI
262 >    
263 >    if (update_nlist) then
264 >      
265 >       !! save current configuration, construct neighbor list,
266 >       !! and calculate forces
267 >       call saveNeighborList(nlocal, q)
268 >      
269 >       neighborListSize = size(list)
270 >       nlist = 0      
271 >      
272 >       do i = 1, nrow
273 >          point(i) = nlist + 1
274 >          
275 >          prepair_inner: do j = 1, ncol
276 >            
277 >             if (skipThisPair(i,j)) cycle prepair_inner
278 >            
279 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
280 >            
281 >             if (rijsq < rlistsq) then            
282 >                
283 >                nlist = nlist + 1
284 >                
285 >                if (nlist > neighborListSize) then
286 >                   call expandNeighborList(nlocal, listerror)
287 >                   if (listerror /= 0) then
288 >                      error = -1
289 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
290 >                      return
291 >                   end if
292 >                   neighborListSize = size(list)
293 >                endif
294 >                
295 >                list(nlist) = j
296 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
297 >             endif
298 >          enddo prepair_inner
299 >       enddo
300 >
301 >       point(nrow + 1) = nlist + 1
302 >      
303 >    else  !! (of update_check)
304 >
305 >       ! use the list to find the neighbors
306 >       do i = 1, nrow
307 >          JBEG = POINT(i)
308 >          JEND = POINT(i+1) - 1
309 >          ! check thiat molecule i has neighbors
310 >          if (jbeg .le. jend) then
311 >            
312 >             do jnab = jbeg, jend
313 >                j = list(jnab)
314 >
315 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
316 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
317 >                     u_l, A, f, t, pot_local)
318 >
319 >             enddo
320 >          endif
321 >       enddo
322      endif
323      
324 + #else
325 +    
326 +    if (update_nlist) then
327 +      
328 +       ! save current configuration, contruct neighbor list,
329 +       ! and calculate forces
330 +       call saveNeighborList(natoms, q)
331 +      
332 +       neighborListSize = size(list)
333 +  
334 +       nlist = 0
335 +      
336 +       do i = 1, natoms-1
337 +          point(i) = nlist + 1
338 +          
339 +          prepair_inner: do j = i+1, natoms
340 +            
341 +             if (skipThisPair(i,j))  cycle prepair_inner
342 +                          
343 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
344 +          
345 +
346 +             if (rijsq < rlistsq) then
347 +                
348 +                nlist = nlist + 1
349 +              
350 +                if (nlist > neighborListSize) then
351 +                   call expandNeighborList(natoms, listerror)
352 +                   if (listerror /= 0) then
353 +                      error = -1
354 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
355 +                      return
356 +                   end if
357 +                   neighborListSize = size(list)
358 +                endif
359 +                
360 +                list(nlist) = j
361 +                
362 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
363 +                        u_l, A, f, t, pot)
364 +                
365 +             endif
366 +          enddo prepair_inner
367 +       enddo
368 +      
369 +       point(natoms) = nlist + 1
370 +      
371 +    else !! (update)
372 +      
373 +       ! use the list to find the neighbors
374 +       do i = 1, natoms-1
375 +          JBEG = POINT(i)
376 +          JEND = POINT(i+1) - 1
377 +          ! check thiat molecule i has neighbors
378 +          if (jbeg .le. jend) then
379 +            
380 +             do jnab = jbeg, jend
381 +                j = list(jnab)
382 +
383 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
384 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
385 +                     u_l, A, f, t, pot)
386 +
387 +             enddo
388 +          endif
389 +       enddo
390 +    endif    
391 + #endif
392 +    !! Do rest of preforce calculations
393 +   call do_preforce(nlocal,pot)
394 +    else
395 +       !! See if we need to update neighbor lists for non pre-pair
396 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
397 +    endif
398 +
399 +
400 +
401 +
402 +
403 + !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
404 +
405 +
406 +
407 +
408 +  
409   #ifdef IS_MPI
410      
411      if (update_nlist) then
# Line 268 | Line 426 | contains
426              
427               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
428              
429 <             if (rijsq <  rlistsq) then            
429 >             if (rijsq < rlistsq) then            
430                  
431                  nlist = nlist + 1
432                  
# Line 284 | Line 442 | contains
442                  
443                  list(nlist) = j
444                                  
445 <                if (rijsq <  rcutsq) then
446 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
447 <                        u_l, A, f, t, pot_local)
290 <                endif
445 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
446 >                     u_l, A, f, t, pot_local)
447 >                
448               endif
449            enddo inner
450         enddo
# Line 337 | Line 494 | contains
494               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
495            
496  
497 <             if (rijsq <  rlistsq) then
497 >             if (rijsq < rlistsq) then
498                  
499                  nlist = nlist + 1
500                
# Line 353 | Line 510 | contains
510                  
511                  list(nlist) = j
512                  
513 <                if (rijsq <  rcutsq) then
357 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
513 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
514                          u_l, A, f, t, pot)
515 <                endif
515 >                
516               endif
517            enddo inner
518         enddo
# Line 502 | Line 658 | contains
658  
659   #endif
660      
505 –    write(*,*) 'T(1) = '
506 –    write(*,'(3es12.3)') t(1,1), t(1,2), t(1,3)
507 –    write(*,*)
508 –
509 –    write(*,*) 'T(2) = '
510 –    write(*,'(3es12.3)') t(2,1), t(2,2), t(2,3)
511 –    write(*,*)
512 –
661    end subroutine do_force_loop
662  
663    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
# Line 528 | Line 676 | contains
676      logical :: is_LJ_i, is_LJ_j
677      logical :: is_DP_i, is_DP_j
678      logical :: is_GB_i, is_GB_j
679 +    logical :: is_EAM_i,is_EAM_j
680      logical :: is_Sticky_i, is_Sticky_j
681      integer :: me_i, me_j
682  
683      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(*,*)
543 –
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(*,*)
684  
556 –
685   #ifdef IS_MPI
686      if (tagRow(i) .eq. tagColumn(j)) then
687         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 617 | Line 745 | contains
745      endif
746      
747  
748 +  
749 +   if (FF_uses_EAM .and. SimUsesEAM()) then
750 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
751 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
752 +      
753 +      if ( is_EAM_i .and. is_EAM_j ) &
754 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
755 +   endif
756 +
757  
758 +
759 +
760    end subroutine do_pair
761  
762  
763 +
764 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
765 +   real( kind = dp ) :: pot
766 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
767 +   real (kind=dp), dimension(9,getNlocal()) :: A
768 +   real (kind=dp), dimension(3,getNlocal()) :: f
769 +   real (kind=dp), dimension(3,getNlocal()) :: t
770 +  
771 +   logical, intent(inout) :: do_pot, do_stress
772 +   integer, intent(in) :: i, j
773 +   real ( kind = dp ), intent(inout)    :: rijsq
774 +   real ( kind = dp )                :: r
775 +   real ( kind = dp ), intent(inout) :: d(3)
776 +  
777 +   logical :: is_EAM_i, is_EAM_j
778 +  
779 +   integer :: me_i, me_j
780 +  
781 +   r = sqrt(rijsq)
782 +  
783 + #ifdef IS_MPI
784 +   if (tagRow(i) .eq. tagColumn(j)) then
785 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
786 +   endif
787 +  
788 +   me_i = atid_row(i)
789 +   me_j = atid_col(j)
790 +  
791 + #else
792 +  
793 +   me_i = atid(i)
794 +   me_j = atid(j)
795 +  
796 + #endif
797 +  
798 +   if (FF_uses_EAM .and. SimUsesEAM()) then
799 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
800 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
801 +      
802 +      if ( is_EAM_i .and. is_EAM_j ) &
803 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
804 +   endif
805 +  end subroutine do_prepair
806 +
807 +
808 +
809 +
810 +  subroutine do_preforce(nlocal,pot)
811 +    integer :: nlocal
812 +    real( kind = dp ) :: pot
813 +
814 +   if (FF_uses_EAM .and. SimUsesEAM()) then
815 +      call calc_EAM_preforce_Frho(nlocal,pot)
816 +   endif
817 +
818 +
819 +  end subroutine do_preforce
820 +  
821 +  
822    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
823      
824      real (kind = dp), dimension(3) :: q_i

Diff Legend

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