*** debug.F90	2013-12-05 18:00:37.549843156 -0500
--- debug.modified	2013-12-05 17:58:53.438841500 -0500
***************
*** 194,202 ****
--- 194,220 ----
     if ( master )then
        write(6,*)'DEBUG FORCE!; calling force routine'
     end if
+    do k=1,3
+       do j=1,3
+          subvir = 0.d0
+       end do
+    end do
     call get_analfrc(xx,ix,ih,ipairs,x,f, &
           vir,ener,qsetup)
     call get_intvir(natom,x,f,subvir,intvir)
+ 
+    WRITE(6,'(A,100ES12.3)')"rec_vir ",rec_vir
+    WRITE(6,'(A,100ES12.3)')"dir_vir ",dir_vir
+    WRITE(6,'(A,100ES12.3)')"adj_vir ",adj_vir
+    WRITE(6,'(A,100ES12.3)')"rec_vird",rec_vird
+    WRITE(6,'(A,100ES12.3)')"self_vir",self_vir
+    WRITE(6,'(A,100ES12.3)')"atvir   ",atvir
+    WRITE(6,'(A,100ES12.3)')"molvir  ",molvir
+    WRITE(6,'(A,100ES12.3)')"subvir  ",subvir
+    WRITE(6,'(A,100ES12.3)')"intvir  ",intvir
+ 
+ 
+ 
     do k = 1,3
        do j = 1,3
           atvir(j,k) = atvir(j,k) + intvir(j,k)
***************
*** 326,339 ****
              ene,del,dudv,type, qsetup)
        apmvir = 3*volume*dudv
        call check_vtens(xx,ix,ih,ipairs,x,f, &
!             ene,del,mduda,molvir,mapduda,type, qsetup)
        ! Next the atvir, type = 2
        type = 2
        call check_virial(xx,ix,ih,ipairs,x,f, &
              ene,del,dudv,type, qsetup)
        apavir = 3*volume*dudv
        call check_vtens(xx,ix,ih,ipairs,x,f, &
!             ene,del,duda,atvir,apduda,type, qsetup)
        if ( master )then
           write(6,*)'Checking analytic virial trace versus'
           write(6,*)'Numerical calculation of 3V dU/dV'
--- 344,357 ----
              ene,del,dudv,type, qsetup)
        apmvir = 3*volume*dudv
        call check_vtens(xx,ix,ih,ipairs,x,f, &
!             ene,del,mduda,mvten,mapduda,type, qsetup)
        ! Next the atvir, type = 2
        type = 2
        call check_virial(xx,ix,ih,ipairs,x,f, &
              ene,del,dudv,type, qsetup)
        apavir = 3*volume*dudv
        call check_vtens(xx,ix,ih,ipairs,x,f, &
!             ene,del,duda,avten,apduda,type, qsetup)
        if ( master )then
           write(6,*)'Checking analytic virial trace versus'
           write(6,*)'Numerical calculation of 3V dU/dV'
***************
*** 1004,1010 ****
     den = max(abs(a),abs(b))+small
     err = abs(a-b)/den
     write(6,60)string,a,b,err
!    60 format(1x,a,e12.5,' vs ',e12.5,' relative error: ',e9.3)
     return
  end subroutine compare 
  !--------------------------------------------------
--- 1022,1028 ----
     den = max(abs(a),abs(b))+small
     err = abs(a-b)/den
     write(6,60)string,a,b,err
!    60 format(1x,a,e12.5,' vs ',e12.5,' relative error: ',e12.5)
     return
  end subroutine compare 
  !--------------------------------------------------
***************
*** 1438,1449 ****
  subroutine check_virial(xx,ix,ih,ipairs,x,f, &
        ene,del,dudv,type, qsetup)
     use nblist,only: volume
     implicit none
  #  include "memory.h"
  #  include "md.h"
  
     _REAL_ xx(*),x(*),f(*), &
!          ene(*),del,dudv
     integer ipairs(*),ix(*),type
     character(len=4) ih(*)
     _REAL_ factor(3),term,vir(4)
--- 1456,1469 ----
  subroutine check_virial(xx,ix,ih,ipairs,x,f, &
        ene,del,dudv,type, qsetup)
     use nblist,only: volume
+    use state
     implicit none
  #  include "memory.h"
  #  include "md.h"
  
+    TYPE(state_rec) :: ene
     _REAL_ xx(*),x(*),f(*), &
!          del,dudv
     integer ipairs(*),ix(*),type
     character(len=4) ih(*)
     _REAL_ factor(3),term,vir(4)
***************
*** 1458,1466 ****
     factor(3) = term
     call redo_ucell(factor)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene(23),vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update,0 )
!    enep = ene(23)
     volp = volume
     ! put cell back
     term = 1.d0/term
--- 1478,1486 ----
     factor(3) = term
     call redo_ucell(factor)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene,vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update,0 )
!    enep = ene%pot%tot
     volp = volume
     ! put cell back
     term = 1.d0/term
***************
*** 1476,1484 ****
     factor(3) = term
     call redo_ucell(factor)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene(23),vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0 )
!    enem = ene(23)
     volm = volume
     dudv = (enep - enem)/(volp-volm)
     ! put cell back
--- 1496,1504 ----
     factor(3) = term
     call redo_ucell(factor)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene,vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0 )
!    enem = ene%pot%tot
     volm = volume
     dudv = (enep - enem)/(volp-volm)
     ! put cell back
***************
*** 1488,1494 ****
     factor(3) = term
     call redo_ucell(factor)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene(23),vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0 )
     return
  end subroutine check_virial 
--- 1508,1514 ----
     factor(3) = term
     call redo_ucell(factor)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene,vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0 )
     return
  end subroutine check_virial 
***************
*** 1499,1510 ****
  subroutine check_vtens(xx,ix,ih,ipairs,x,f, &
        ene,del,duda,vten,apduda,type, qsetup)
     use nblist, only: ucell, recip
     implicit none
  #  include "memory.h"
  #  include "md.h"
  
     _REAL_ xx(*),x(*),f(*), &
!          ene(*),del,duda(3,3),vir(4),vten(3,3),apduda(3,3)
     integer ipairs(*),ix(*),type
     character(len=4) ih(*)
     _REAL_ enep,enem,save(3,3),new(3,3)
--- 1519,1533 ----
  subroutine check_vtens(xx,ix,ih,ipairs,x,f, &
        ene,del,duda,vten,apduda,type, qsetup)
     use nblist, only: ucell, recip
+    use nblist, only : volume
+    use state
     implicit none
  #  include "memory.h"
  #  include "md.h"
  
+    TYPE(state_rec) :: ene
     _REAL_ xx(*),x(*),f(*), &
!          del,duda(3,3),vir(4),vten(3,3),apduda(3,3)
     integer ipairs(*),ix(*),type
     character(len=4) ih(*)
     _REAL_ enep,enem,save(3,3),new(3,3)
***************
*** 1529,1543 ****
           new(i,j) = save(i,j) + del
           call new_ucell_gen(new)
           call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!          call force(xx,ix,ih,ipairs,x,f,ene(23),vir, &
                 xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0)
!          enep = ene(23)
           new(i,j) = save(i,j) - del
           call new_ucell_gen(new)
           call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!          call force(xx,ix,ih,ipairs,x,f,ene(23),vir, &
                 xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0)
!          enem = ene(23)
           apduda(i,j) = (enep-enem)/(2.d0*del)
        end do
     end do
--- 1552,1566 ----
           new(i,j) = save(i,j) + del
           call new_ucell_gen(new)
           call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!          call force(xx,ix,ih,ipairs,x,f,ene,vir, &
                 xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0)
!          enep = ene%pot%tot
           new(i,j) = save(i,j) - del
           call new_ucell_gen(new)
           call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!          call force(xx,ix,ih,ipairs,x,f,ene,vir, &
                 xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0)
!          enem = ene%pot%tot
           apduda(i,j) = (enep-enem)/(2.d0*del)
        end do
     end do
***************
*** 1549,1557 ****
     end do
     call new_ucell_gen(new)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene(23),vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0)
     ! calc approx duda
     do i = 1,3
        do j = 1,3
           duda(i,j) = 0.d0
--- 1572,1592 ----
     end do
     call new_ucell_gen(new)
     call ew_pscale(natom,x,xx(lmass),nspm,ix(i70),type)
!    call force(xx,ix,ih,ipairs,x,f,ene,vir, &
           xx(l96), xx(l97), xx(l98), xx(l99), qsetup, do_list_update, 0)
+ 
+    do i = 1,3
+       do j = 1,3
+          new(i,j) = 0.d0
+          do k = 1,3
+             new(i,j) = new(i,j) + apduda(i,k)*ucell(k,j) !/volume
+          end do
+       end do
+    end do
+    WRITE(6,'(A,100ES12.3)')"new=",new
+ 
     ! calc approx duda
+    WRITE(6,'(A,100ES12.3)')"vten=",vten
     do i = 1,3
        do j = 1,3
           duda(i,j) = 0.d0