--- trunk/OOPSE/libmdtools/do_Forces.F90 2003/04/11 15:16:59 490 +++ trunk/OOPSE/libmdtools/do_Forces.F90 2003/07/14 21:28:54 597 @@ -4,7 +4,7 @@ !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @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 $ +!! @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 $ module do_Forces use force_globals @@ -502,6 +502,14 @@ contains #endif + write(*,*) 'T(1) = ' + write(*,'(3es12.3)') t(1,1), t(1,2), t(1,3) + write(*,*) + + write(*,*) 'T(2) = ' + write(*,'(3es12.3)') t(2,1), t(2,2), t(2,3) + write(*,*) + end subroutine do_force_loop subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot) @@ -525,8 +533,27 @@ contains r = sqrt(rijsq) + write(*,*) 'ul(1) = ' + write(*,'(3es12.3)') u_l(1,1), u_l(1,2), u_l(1,3) + write(*,*) + write(*,*) 'ul(2) = ' + write(*,'(3es12.3)') u_l(2,1), u_l(2,2), u_l(2,3) + write(*,*) + + write(*,*) 'A(1) = ' + write(*,'(3es12.3)') A(1,1), A(2,1), A(3,1) + write(*,'(3es12.3)') A(4,1), A(5,1), A(6,1) + write(*,'(3es12.3)') A(7,1), A(8,1), A(9,1) + write(*,*) + write(*,*) 'A(2) = ' + write(*,'(3es12.3)') A(1,2), A(2,2), A(3,2) + write(*,'(3es12.3)') A(4,2), A(5,2), A(6,2) + write(*,'(3es12.3)') A(7,2), A(8,2), A(9,2) + write(*,*) + + #ifdef IS_MPI if (tagRow(i) .eq. tagColumn(j)) then write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j) @@ -589,6 +616,8 @@ contains endif endif + + end subroutine do_pair @@ -597,21 +626,52 @@ contains real (kind = dp), dimension(3) :: q_i real (kind = dp), dimension(3) :: q_j real ( kind = dp ), intent(out) :: r_sq - real( kind = dp ) :: d(3) - real( kind = dp ) :: d_old(3) + real( kind = dp ) :: d(3), scaled(3) + integer i + d(1:3) = q_j(1:3) - q_i(1:3) - d_old = d + ! Wrap back into periodic box if necessary if ( SimUsesPBC() ) then - d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * & - int(abs(d(1:3)/box(1:3)) + 0.5_dp) + if( .not.boxIsOrthorhombic ) then + ! calc the scaled coordinates. + + scaled = matmul(HmatInv, d) + + ! wrap the scaled coordinates + + scaled = scaled - anint(scaled) + + + ! calc the wrapped real coordinates from the wrapped scaled + ! coordinates + + d = matmul(Hmat,scaled) + + else + ! calc the scaled coordinates. + + do i = 1, 3 + scaled(i) = d(i) * HmatInv(i,i) + + ! wrap the scaled coordinates + + scaled(i) = scaled(i) - anint(scaled(i)) + + ! calc the wrapped real coordinates from the wrapped scaled + ! coordinates + + d(i) = scaled(i)*Hmat(i,i) + enddo + endif endif + r_sq = dot_product(d,d) - + end subroutine get_interatomic_vector - + subroutine check_initialization(error) integer, intent(out) :: error