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 619 by mmeineke, Tue Jul 15 22:22:41 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.20 2003-07-15 22:22:41 mmeineke Exp $, $Date: 2003-07-15 22:22:41 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
8  
9   module do_Forces
10    use force_globals
# Line 185 | Line 185 | contains
185      logical :: update_nlist  
186      integer :: i, j, jbeg, jend, jnab
187      integer :: nlist
188 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
188 >    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut, rrf, rt, dielect
189      real(kind=dp),dimension(3) :: d
190      real(kind=dp) :: rfpot, mu_i, virial
191      integer :: me_i
# Line 193 | Line 193 | contains
193      integer :: neighborListSize
194      integer :: listerror, error
195      integer :: localError
196 +    
197  
198      !! initialize local variables  
199  
# Line 208 | Line 209 | contains
209    
210      call getRcut(rcut,rc2=rcutsq)
211      call getRlist(rlist,rlistsq)
212 +    rt = getRt()
213 +    rrf = getRrf()
214 +    dielect = getDielect()
215      
216 +    if( FF_uses_LJ) then      
217 +       call lj_new_rcut( rcut, localError )
218 +       if ( localError .ne. 0 ) then
219 +          error = -1
220 +          return
221 +       end if
222 +    end if
223 +    
224 +    
225 +    if( FF_uses_dipoles ) then
226 +      
227 +       if( rcut .lt. rrf ) then
228 +          rcut = rrf
229 +          rlist = rcut + 1.0_dp
230 +          rcutsq = rcut * rcut
231 +          rlistsq = rlist * rlist
232 +       end if
233 +      
234 +       call initialize_dipole( rrf, rt )
235 +    end if
236 +    
237 +    if( FF_uses_RF )call initialize_rf( rrf, rt, dielect )
238 +    
239 +    
240      call check_initialization(localError)
241      if ( localError .ne. 0 ) then
242         error = -1
# Line 524 | Line 552 | contains
552      integer :: me_i, me_j
553  
554      r = sqrt(rijsq)
527
528
555  
556   #ifdef IS_MPI
557      if (tagRow(i) .eq. tagColumn(j)) then
# Line 589 | Line 615 | contains
615         endif
616      endif
617      
618 +
619 +
620    end subroutine do_pair
621  
622  
# Line 597 | Line 625 | contains
625      real (kind = dp), dimension(3) :: q_i
626      real (kind = dp), dimension(3) :: q_j
627      real ( kind = dp ), intent(out) :: r_sq
628 <    real( kind = dp ) :: d(3)
629 <    real( kind = dp ) :: d_old(3)
628 >    real( kind = dp ) :: d(3), scaled(3)
629 >    integer i
630 >
631      d(1:3) = q_j(1:3) - q_i(1:3)
632 <    d_old = d
632 >
633      ! Wrap back into periodic box if necessary
634      if ( SimUsesPBC() ) then
635        
636 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
637 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
636 >       if( .not.boxIsOrthorhombic ) then
637 >          ! calc the scaled coordinates.
638 >          
639 >          scaled = matmul(HmatInv, d)
640 >          
641 >          ! wrap the scaled coordinates
642 >
643 >          scaled = scaled  - anint(scaled)
644 >          
645 >
646 >          ! calc the wrapped real coordinates from the wrapped scaled
647 >          ! coordinates
648 >
649 >          d = matmul(Hmat,scaled)
650 >
651 >       else
652 >          ! calc the scaled coordinates.
653 >          
654 >          do i = 1, 3
655 >             scaled(i) = d(i) * HmatInv(i,i)
656 >            
657 >             ! wrap the scaled coordinates
658 >            
659 >             scaled(i) = scaled(i) - anint(scaled(i))
660 >            
661 >             ! calc the wrapped real coordinates from the wrapped scaled
662 >             ! coordinates
663 >
664 >             d(i) = scaled(i)*Hmat(i,i)
665 >          enddo
666 >       endif
667        
668      endif
669 +    
670      r_sq = dot_product(d,d)
671 <        
671 >    
672    end subroutine get_interatomic_vector
673 <
673 >  
674    subroutine check_initialization(error)
675      integer, intent(out) :: error
676      

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines