!
! Copyright (C) 2010 Davide Ceresoli
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!---------------------------------------------------------------------
SUBROUTINE impose_deviatoric_strain ( at_old, at )
  !---------------------------------------------------------------------
  !
  !     Impose a pure deviatoric (volume-conserving) deformation
  !     Needed to enforce volume conservation in variable-cell MD/optimization
  !
  USE kinds, ONLY: dp
! RICHARD ADDED LOGICAL VARIABLE FOR 2DSHAPE
  USE cell_base, ONLY: fix_area ! ADDED FOR 2DSHAPE logical .true. for cell_dofree='2Dshape'
! RICHARD
  IMPLICIT NONE
  REAL(dp), INTENT(in)    :: at_old(3,3)
  REAL(dp), INTENT(inout) :: at(3,3)
  REAL(dp) :: tr, omega, omega_old
  INTEGER :: i,j  

  tr = (at(1,1)+at(2,2)+at(3,3))/3.d0
  tr = tr - (at_old(1,1)+at_old(2,2)+at_old(3,3))/3.d0
  ! Commented out, while waiting for better idea:
  ! it breaks the symmetry of hexagonal lattices - PG
  ! at(1,1) = at(1,1) - tr
  ! at(2,2) = at(2,2) - tr
  ! at(3,3) = at(3,3) - tr
  ! print '("difference in trace: ",e12.4)', tr

  call volume (1.d0, at_old(1,1), at_old(1,2), at_old(1,3), omega_old)
  call volume (1.d0, at(1,1), at(1,2), at(1,3), omega)
! RICHARD REVISED CODE FOR cell_dofree='2Dshape'********************
  do i = 1,3
  do j = 1,3
  if (fix_area .and. j.eq.3) then
    at(i,j) = at(i,j) ! DON'T CHANGE IN z- DIRECTION IF 2DSHAPE
  else  
  at(i,j) = at(i,j) * (omega_old / omega)**(1.d0/3.d0)
  endif
  enddo
  enddo
! RICHARD
END SUBROUTINE impose_deviatoric_strain

!---------------------------------------------------------------------
SUBROUTINE impose_deviatoric_stress ( sigma )
  !---------------------------------------------------------------------
  !
  !     Impose a pure deviatoric stress
  !
  USE kinds, ONLY: dp
  USE cell_base, ONLY: fix_area ! ADDED FOR 2DSHAPE
  USE io_global, ONLY: stdout
  IMPLICIT NONE
  REAL(dp), INTENT(inout) :: sigma(3,3)
  REAL(dp) :: tr
! RICHARD MODIFIED CODE TO ADD 2DSHAPE
  if(fix_area) then ! IF cell_dofree='2Dshape' IMPOSE FOR 2D CASE FOR FIXED AREA IN xy PLANE
  tr = (sigma(1,1)+sigma(2,2))/2.d0
  sigma(1,1) = sigma(1,1) - tr
  sigma(2,2) = sigma(2,2) - tr
  WRITE (stdout,'(5x,"Area is kept fixed: isostatic in-plane pressure in xy set to zero")')
  else ! IMPOSE FOR cell_dofree='shape' 3D CASE AT FIXED VOLUME
  tr = (sigma(1,1)+sigma(2,2)+sigma(3,3))/3.d0
  sigma(1,1) = sigma(1,1) - tr
  sigma(2,2) = sigma(2,2) - tr
  sigma(3,3) = sigma(3,3) - tr
  WRITE (stdout,'(5x,"Volume is kept fixed: isostatic pressure set to zero")')
  endif
! RICHARD
END SUBROUTINE impose_deviatoric_stress