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 441 by chuckv, Tue Apr 1 16:50:14 2003 UTC vs.
Revision 572 by mmeineke, Wed Jul 2 21:26:55 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.7 2003-04-01 16:50:14 chuckv Exp $, $Date: 2003-04-01 16:50:14 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
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 $
8  
9   module do_Forces
10    use force_globals
# Line 140 | Line 140 | contains
140  
141      if (FF_uses_GB .and. FF_uses_LJ) then
142      endif
143 +    if (.not. do_forces_initialized) then
144 +       !! Create neighbor lists
145 +       call expandNeighborList(getNlocal(), my_status)
146 +       if (my_Status /= 0) then
147 +          write(default_error,*) "SimSetup: ExpandNeighborList returned error."
148 +          thisStat = -1
149 +          return
150 +       endif
151 +    endif
152  
144
153      do_forces_initialized = .true.    
154  
155    end subroutine init_FF
# Line 246 | Line 254 | contains
254        
255         !! save current configuration, construct neighbor list,
256         !! and calculate forces
257 <       call saveNeighborList(q)
257 >       call saveNeighborList(nlocal, q)
258        
259         neighborListSize = size(list)
260         nlist = 0      
# Line 313 | Line 321 | contains
321        
322         ! save current configuration, contruct neighbor list,
323         ! and calculate forces
324 <       call saveNeighborList(q)
324 >       call saveNeighborList(natoms, q)
325        
326         neighborListSize = size(list)
327    
# Line 479 | Line 487 | contains
487      endif
488  
489      if (do_stress) then
490 <       call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
490 >      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
491              mpi_comm_world,mpi_err)
492 <       call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
492 >       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
493              mpi_comm_world,mpi_err)
494      endif
495  
# Line 499 | Line 507 | contains
507    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
508  
509      real( kind = dp ) :: pot
510 <    real( kind = dp ), dimension(:,:) :: u_l
511 <    real (kind=dp), dimension(:,:) :: A
512 <    real (kind=dp), dimension(:,:) :: f
513 <    real (kind=dp), dimension(:,:) :: t
510 >    real( kind = dp ), dimension(3,getNlocal()) :: u_l
511 >    real (kind=dp), dimension(9,getNlocal()) :: A
512 >    real (kind=dp), dimension(3,getNlocal()) :: f
513 >    real (kind=dp), dimension(3,getNlocal()) :: t
514  
515      logical, intent(inout) :: do_pot, do_stress
516      integer, intent(in) :: i, j
# Line 517 | Line 525 | contains
525  
526      r = sqrt(rijsq)
527  
528 +
529 +
530   #ifdef IS_MPI
531 +    if (tagRow(i) .eq. tagColumn(j)) then
532 +       write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
533 +    endif
534  
535      me_i = atid_row(i)
536      me_j = atid_col(j)
# Line 584 | Line 597 | contains
597      real (kind = dp), dimension(3) :: q_i
598      real (kind = dp), dimension(3) :: q_j
599      real ( kind = dp ), intent(out) :: r_sq
600 <    real( kind = dp ) :: d(3)
601 <    real( kind = dp ) :: d_old(3)
602 <    d(1:3) = q_i(1:3) - q_j(1:3)
603 <    d_old = d
600 >    real( kind = dp ) :: d(3), scaled(3)
601 >    integer i
602 >
603 >    d(1:3) = q_j(1:3) - q_i(1:3)
604 >
605      ! Wrap back into periodic box if necessary
606      if ( SimUsesPBC() ) then
607        
608 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
609 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
608 >       if( .not.boxIsOrthorhombic ) then
609 >          ! calc the scaled coordinates.
610 >          
611 >          scaled = matmul(HmatInv, d)
612 >          
613 >          ! wrap the scaled coordinates
614 >
615 >          scaled = scaled  - anint(scaled)
616 >          
617 >
618 >          ! calc the wrapped real coordinates from the wrapped scaled
619 >          ! coordinates
620 >
621 >          d = matmul(Hmat,scaled)
622 >
623 >       else
624 >          ! calc the scaled coordinates.
625 >          
626 >          do i = 1, 3
627 >             scaled(i) = d(i) * HmatInv(i,i)
628 >            
629 >             ! wrap the scaled coordinates
630 >            
631 >             scaled(i) = scaled(i) - anint(scaled(i))
632 >            
633 >             ! calc the wrapped real coordinates from the wrapped scaled
634 >             ! coordinates
635 >
636 >             d(i) = scaled(i)*Hmat(i,i)
637 >          enddo
638 >       endif
639        
640      endif
641 +    
642      r_sq = dot_product(d,d)
643 <        
643 >    
644    end subroutine get_interatomic_vector
645 <
645 >  
646    subroutine check_initialization(error)
647      integer, intent(out) :: error
648      
# Line 655 | Line 699 | contains
699      rf = 0.0_dp
700      tau_Temp = 0.0_dp
701      virial_Temp = 0.0_dp
658    
702    end subroutine zero_work_arrays
703    
704    function skipThisPair(atom1, atom2) result(skip_it)

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines