MAIN.F59

tau_orig = tauin
    FtsvdW_exact = FtsvdW

    do s = 1, nat, 1
      do i = 1, 3, 1
        do is = 1, 2, 1
          grad = 0.0_DP
          plus = 0.0_DP
          minus = 0.0_DP
          tau_temp = tau_orig
          if(is.eq.1) tau_temp(i, s) = tau_orig(i, s) + dx
          if(is.eq.2) tau_temp(i, s) = tau_orig(i, s) - dx
          CALL tsvdw_para_init()
          CALL tsvdw_pbc(tau_temp)
          CALL tsvdw_unique_pair()
          CALL tsvdw_rhotot( rhor )
          CALL tsvdw_screen()
          CALL tsvdw_veff()
          CALL tsvdw_dveff()
          CALL tsvdw_effqnts()
          CALL tsvdw_energy()
          CALL tsvdw_cleanup()
          if(is.eq.1) plus = EtsvdW
          if(is.eq.2) minus = EtsvdW
          end do
        grad = (plus-minus)/(2.0_DP*dx)
        write(stdout, *), i, s, minus, grad, FtsvdW_exact(i, s)
        end do
      end do

    end subroutine tsvdw_check_dEdr

  subroutine tsvdw_check_dEdh(tauin, rhor)
    implicit none
    REAL(DP), INTENT(IN) :: rhor(:,:)
    REAL(DP) :: tauin(3,nat)
    real(dp) :: dx, ainv_orig(3,3), h_orig(3,3), HtsvdW_exact(3,3), ainv_temp(3,3), h_temp(3,3), plus, minus, grad
    integer :: i, s, is
    dx = 0.00001_DP

    ainv_orig = ainv_
    h_orig = h_
    HtsvdW_exact = HtsvdW

    do s = 1, 3, 1
      do i = 1, 3, 1
        do is = 1, 2, 1
          grad = 0.0_DP
          plus = 0.0_DP
          minus = 0.0_DP
          h_temp = h_orig
          if(is.eq.1) h_temp(s, i) = h_orig(s, i) + dx
          if(is.eq.2) h_temp(s, i) = h_orig(s, i) - dx
          call inv3x3_mat(h_temp, ainv_temp)
          h_ = h_temp
          ainv_ = ainv_temp
          CALL tsvdw_para_init()
          CALL tsvdw_pbc(tauin)
          CALL tsvdw_unique_pair()
          CALL tsvdw_rhotot( rhor )
          CALL tsvdw_screen()
          CALL tsvdw_veff()
          CALL tsvdw_dveff()
          CALL tsvdw_effqnts()
          CALL tsvdw_energy()
          CALL tsvdw_cleanup()
          if(is.eq.1) plus = EtsvdW
          if(is.eq.2) minus = EtsvdW
          end do
        grad = (plus-minus)/(2.0_DP*dx)
        write(stdout, *), i, s, minus, grad, HtsvdW_exact(i, s)
        end do
      end do

    end subroutine tsvdw_check_dEdh

  subroutine tsvdw_check_dEdn(tauin, rhor)
    implicit none
    REAL(DP), INTENT(IN) :: rhor(:,:)
    REAL(DP) :: tauin(3,nat)
    real(dp) :: plus, minus, grad, rhotot_temp(nr1*nr2*nr3), UtsvdW_exact(nr1*nr2*nr3), dx
    integer  :: is, ir

    dx = 0.00001_DP
    UtsvdW_exact = UtsvdW

    write(stdout, *), "ir          Numerical              Analytic"
    write(stdout, *), "-------------------------------------------"
    do ir = 1, nr1*nr2*nr3, 1
      grad = 0.0_DP
      plus = 0.0_DP
      minus = 0.0_DP
      do is = 1, 2, 1
        vdw_lstep0 = .true.
        vdw_lscreen = .true.
        CALL tsvdw_para_init()
        CALL tsvdw_pbc(tauin)
        CALL tsvdw_unique_pair()
        CALL tsvdw_rhotot( rhor )
        if(is.eq.1) rhotot(ir) = rhotot(ir) + dx
        if(is.eq.2) rhotot(ir) = rhotot(ir) - dx
        CALL tsvdw_screen()
        CALL tsvdw_veff()
        CALL tsvdw_dveff()
        CALL tsvdw_effqnts()
        CALL tsvdw_energy()
        CALL tsvdw_cleanup()
        if(is.eq.1) plus = EtsvdW
        if(is.eq.2) minus = EtsvdW
        end do
      grad = (plus-minus)/(2.0_DP*dx)*(nr1*nr2*nr3)/(omega*8.0_DP)
      write(stdout, *), ir, grad, UtsvdW_exact(ir)
      end do
  end subroutine tsvdw_check_dEdn

  subroutine inv3x3_mat(A,T)
    use kinds, ONLY : DP
    IMPLICIT NONE
    real(DP), intent(in) :: A(3,3)
    real(DP) :: D
    real(DP), intent(out) :: T(3,3)
        D = A(1,1)*A(2,2)*A(3,3)-A(1,1)*A(2,3)*A(3,2)- &
        A(1,2)*A(2,1)*A(3,3)+A(1,2)*A(2,3)*A(3,1)+ &
        A(1,3)*A(2,1)*A(3,2)-A(1,3)*A(2,2)*A(3,1)
        T(1,1) = (A(2,2)*A(3,3)-A(2,3)*A(3,2))/D
        T(1,2) = -(A(1,2)*A(3,3)-A(1,3)*A(3,2))/D
        T(1,3) = (A(1,2)*A(2,3)-A(1,3)*A(2,2))/D
        T(2,1) = -(A(2,1)*A(3,3)-A(2,3)*A(3,1))/D
        T(2,2) = (A(1,1)*A(3,3)-A(1,3)*A(3,1))/D
        T(2,3) = -(A(1,1)*A(2,3)-A(1,3)*A(2,1))/D
        T(3,1) = (A(2,1)*A(3,2)-A(2,2)*A(3,1))/D
        T(3,2) = -(A(1,1)*A(3,2)-A(1,2)*A(3,1))/D
        T(3,3) = (A(1,1)*A(2,2)-A(1,2)*A(2,1))/D
      return
         end subroutine inv3x3_mat