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 634 by mmeineke, Thu Jul 17 19:38:23 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.23 2003-07-17 19:38:23 mmeineke Exp $, $Date: 2003-07-17 19:38:23 $, $Name: not supported by cvs2svn $, $Revision: 1.23 $
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)  
261 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
262      endif
263      
264   #ifdef IS_MPI
# Line 268 | Line 281 | contains
281              
282               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
283              
284 <             if (rijsq <  rlistsq) then            
284 >             if (rijsq < rlistsq) then            
285                  
286                  nlist = nlist + 1
287                  
# Line 284 | Line 297 | contains
297                  
298                  list(nlist) = j
299                                  
300 <                if (rijsq <  rcutsq) then
301 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
302 <                        u_l, A, f, t, pot_local)
290 <                endif
300 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
301 >                     u_l, A, f, t, pot_local)
302 >                
303               endif
304            enddo inner
305         enddo
# Line 337 | Line 349 | contains
349               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
350            
351  
352 <             if (rijsq <  rlistsq) then
352 >             if (rijsq < rlistsq) then
353                  
354                  nlist = nlist + 1
355                
# Line 353 | Line 365 | contains
365                  
366                  list(nlist) = j
367                  
368 <                if (rijsq <  rcutsq) then
357 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
368 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
369                          u_l, A, f, t, pot)
370 <                endif
370 >                
371               endif
372            enddo inner
373         enddo
# Line 525 | Line 536 | contains
536  
537      r = sqrt(rijsq)
538  
528
529
539   #ifdef IS_MPI
540      if (tagRow(i) .eq. tagColumn(j)) then
541         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 589 | Line 598 | contains
598         endif
599      endif
600      
601 +
602 +
603    end subroutine do_pair
604  
605  
606 +
607 +  subroutine do_preforce(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
608 +   real( kind = dp ) :: pot
609 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
610 +   real (kind=dp), dimension(9,getNlocal()) :: A
611 +   real (kind=dp), dimension(3,getNlocal()) :: f
612 +   real (kind=dp), dimension(3,getNlocal()) :: t
613 +  
614 +   logical, intent(inout) :: do_pot, do_stress
615 +   integer, intent(in) :: i, j
616 +   real ( kind = dp ), intent(inout)    :: rijsq
617 +   real ( kind = dp )                :: r
618 +   real ( kind = dp ), intent(inout) :: d(3)
619 +  
620 +   logical :: is_EAM_i, is_EAM_j
621 +  
622 +   integer :: me_i, me_j
623 +  
624 +   r = sqrt(rijsq)
625 +  
626 + #ifdef IS_MPI
627 +   if (tagRow(i) .eq. tagColumn(j)) then
628 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
629 +   endif
630 +  
631 +   me_i = atid_row(i)
632 +   me_j = atid_col(j)
633 +  
634 + #else
635 +  
636 +   me_i = atid(i)
637 +   me_j = atid(j)
638 +  
639 + #endif
640 +  
641 +   if (FF_uses_EAM .and. SimUsesEAM()) then
642 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
643 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
644 +      
645 + !!$      if ( is_EAM_i .and. is_EAM_j ) &
646 + !!$           call calc_EAM_prepair(i, j, d, r, rijsq )
647 +   endif
648 +  end subroutine do_preforce
649 +  
650 +  
651    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
652      
653      real (kind = dp), dimension(3) :: q_i

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines