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 470 by chuckv, Mon Apr 7 20:50:46 2003 UTC vs.
Revision 626 by mmeineke, Wed Jul 16 21:30:56 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.11 2003-04-07 20:50:46 chuckv Exp $, $Date: 2003-04-07 20:50:46 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
7 > !! @version $Id: do_Forces.F90,v 1.21 2003-07-16 21:30:55 mmeineke Exp $, $Date: 2003-07-16 21:30:55 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $
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 140 | Line 152 | contains
152  
153      if (FF_uses_GB .and. FF_uses_LJ) then
154      endif
155 +    if (.not. do_forces_initialized) then
156 +       !! Create neighbor lists
157 +       call expandNeighborList(getNlocal(), my_status)
158 +       if (my_Status /= 0) then
159 +          write(default_error,*) "SimSetup: ExpandNeighborList returned error."
160 +          thisStat = -1
161 +          return
162 +       endif
163 +    endif
164  
165 <
166 <    do_forces_initialized = .true.    
167 <
165 >    havePolicies = .true.
166 >    if( haveRlist ) do_forces_initialized = .true.
167 >    
168    end subroutine init_FF
169    
170  
# Line 177 | 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 186 | 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 198 | Line 222 | contains
222      natoms = nlocal
223   #endif
224    
201 –    call getRcut(rcut,rc2=rcutsq)
202 –    call getRlist(rlist,rlistsq)
203 –    
225      call check_initialization(localError)
226      if ( localError .ne. 0 ) then
227         error = -1
# Line 210 | Line 231 | contains
231  
232      do_pot = do_pot_c
233      do_stress = do_stress_c
213 –    
234  
235      ! Gather all information needed by all force loops:
236      
# Line 231 | 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 261 | 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 277 | 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)
283 <                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 330 | 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 346 | Line 365 | contains
365                  
366                  list(nlist) = j
367                  
368 <                if (rijsq <  rcutsq) then
350 <                   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 480 | Line 498 | contains
498      endif
499  
500      if (do_stress) then
501 <       call mpi_allreduce(tau_Temp, tau,9,mpi_double_precision,mpi_sum, &
501 >      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
502              mpi_comm_world,mpi_err)
503         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
504              mpi_comm_world,mpi_err)
# Line 519 | Line 537 | contains
537      r = sqrt(rijsq)
538  
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)
542 +    endif
543  
544      me_i = atid_row(i)
545      me_j = atid_col(j)
# Line 577 | Line 598 | contains
598         endif
599      endif
600      
601 +
602 +
603    end subroutine do_pair
604  
605  
# Line 585 | Line 608 | contains
608      real (kind = dp), dimension(3) :: q_i
609      real (kind = dp), dimension(3) :: q_j
610      real ( kind = dp ), intent(out) :: r_sq
611 <    real( kind = dp ) :: d(3)
612 <    real( kind = dp ) :: d_old(3)
613 <    d(1:3) = q_i(1:3) - q_j(1:3)
614 <    d_old = d
611 >    real( kind = dp ) :: d(3), scaled(3)
612 >    integer i
613 >
614 >    d(1:3) = q_j(1:3) - q_i(1:3)
615 >
616      ! Wrap back into periodic box if necessary
617      if ( SimUsesPBC() ) then
618        
619 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
620 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
619 >       if( .not.boxIsOrthorhombic ) then
620 >          ! calc the scaled coordinates.
621 >          
622 >          scaled = matmul(HmatInv, d)
623 >          
624 >          ! wrap the scaled coordinates
625 >
626 >          scaled = scaled  - anint(scaled)
627 >          
628 >
629 >          ! calc the wrapped real coordinates from the wrapped scaled
630 >          ! coordinates
631 >
632 >          d = matmul(Hmat,scaled)
633 >
634 >       else
635 >          ! calc the scaled coordinates.
636 >          
637 >          do i = 1, 3
638 >             scaled(i) = d(i) * HmatInv(i,i)
639 >            
640 >             ! wrap the scaled coordinates
641 >            
642 >             scaled(i) = scaled(i) - anint(scaled(i))
643 >            
644 >             ! calc the wrapped real coordinates from the wrapped scaled
645 >             ! coordinates
646 >
647 >             d(i) = scaled(i)*Hmat(i,i)
648 >          enddo
649 >       endif
650        
651      endif
652 +    
653      r_sq = dot_product(d,d)
654 <        
654 >    
655    end subroutine get_interatomic_vector
656 <
656 >  
657    subroutine check_initialization(error)
658      integer, intent(out) :: error
659      

Diff Legend

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