MODULE esm_stres_mod

  USE kinds,    ONLY : DP
  USE esm_common_mod, ONLY : esm_w, esm_bc, &
                             mill_2d, imill_2d, ngm_2d, exp_erfc
  IMPLICIT NONE

  ! Workaround for Cray bug - note that exp, cosh, sinh with complex argument
  ! are in the F2008 standard so qe_exp, qe_cosh, qe_sinh are no longer needed

#if defined(__CRAY)
#define QE_EXP  exp
#define QE_COSH cosh
#define QE_SINH sinh
#endif  
CONTAINS

  !-----------------------------------------------------------------------
  !--------------ESM STRESS SUBROUTINE------------------------------------
  !-----------------------------------------------------------------------
  SUBROUTINE esm_stres_har(sigmahar, rhog)
    USE kinds,    ONLY : DP
    USE gvect,    ONLY : ngm
    IMPLICIT NONE
    REAL(DP), INTENT(out)   :: sigmahar(3, 3)
    COMPLEX(DP), INTENT(in) :: rhog(ngm)   !  n(G)

    SELECT CASE (esm_bc)
    CASE ('pbc')
      STOP 'esm_stres_har must not be called for esm_bc = pbc'
    CASE ('bc1')
      CALL esm_stres_har_bc1(sigmahar, rhog)
    CASE ('bc2')
      CALL esm_stres_har_bc2(sigmahar, rhog)
    CASE ('bc3')
      CALL esm_stres_har_bc3(sigmahar, rhog)
    CASE ('bc4')
      STOP 'esm_stres_har has not yet implemented for esm_bc = bc4'
    END SELECT

    RETURN
  END SUBROUTINE esm_stres_har

  SUBROUTINE esm_stres_ewa(sigmaewa)
    !-----------------------------------------------------------------------
    !
    ! Calculates Ewald stresswith both G- and R-space terms.
    ! Determines optimal alpha. Should hopefully work for any structure.
    !
    USE kinds,     ONLY : DP
    USE constants, ONLY : tpi
    USE cell_base, ONLY : tpiba2
    USE ions_base, ONLY : zv, nat, ityp
    USE gvect,     ONLY : gcutm
    IMPLICIT NONE
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    ! output: the ewald stress
    !
    !    here the local variables
    !
    INTEGER :: ia
    ! counter on atoms

    REAL(DP) :: charge, alpha, upperbound
    ! total ionic charge in the cell
    ! alpha term in ewald sum
    ! the maximum radius to consider real space sum
    REAL(DP) :: sigmaewg(3, 3), sigmaewr(3, 3)
    ! ewald stress computed in reciprocal space
    ! ewald stress computed in real space

    charge = sum(zv(ityp(:)))

    ! choose alpha in order to have convergence in the sum over G
    ! upperbound is a safe upper bound for the error in the sum over G
    alpha = 2.9d0
    DO
      alpha = alpha - 0.1d0
      IF (alpha .le. 0.d0) CALL errore('esm_stres_ewa', 'optimal alpha not found', 1)
      upperbound = 2.d0*charge**2*sqrt(2.d0*alpha/tpi)* &
                   erfc(sqrt(tpiba2*gcutm/4.d0/alpha))
      IF (upperbound < 1.0d-7) EXIT
    END DO

    ! G-space sum here.
    ! Determine if this processor contains G=0 and set the constant term
    CALL esm_stres_ewg(alpha, sigmaewg)

    ! R-space sum here (only for the processor that contains G=0)
    CALL esm_stres_ewr(alpha, sigmaewr)

    sigmaewa(:, :) = sigmaewg(:, :) + sigmaewr(:, :)

    RETURN
  END SUBROUTINE esm_stres_ewa

  SUBROUTINE esm_stres_loclong(sigmaloclong, rhog)
    USE kinds,    ONLY : DP
    USE gvect,    ONLY : ngm
    IMPLICIT NONE
    REAL(DP), INTENT(out)   :: sigmaloclong(3, 3)
    COMPLEX(DP), INTENT(in) :: rhog(ngm)   !  n(G)

    SELECT CASE (esm_bc)
    CASE ('pbc')
      STOP 'esm_stres_loclong must not be called for esm_bc = pbc'
    CASE ('bc1')
      CALL esm_stres_loclong_bc1(sigmaloclong, rhog)
    CASE ('bc2')
      CALL esm_stres_loclong_bc2(sigmaloclong, rhog)
    CASE ('bc3')
      CALL esm_stres_loclong_bc3(sigmaloclong, rhog)
    CASE ('bc4')
      STOP 'esm_stres_loclong has not yet implemented for esm_bc = bc4'
    END SELECT

    RETURN
  END SUBROUTINE esm_stres_loclong

  SUBROUTINE esm_stres_har_bc1(sigmahar, rhog)
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : ngm, mill
    USE constants,     ONLY : tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE control_flags, ONLY : gamma_only
    USE fft_base,      ONLY : dfftp
    USE fft_scalar,    ONLY : cft_1z
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(out)   :: sigmahar(3, 3)
    COMPLEX(DP), INTENT(in) :: rhog(ngm)   !  n(G)

    INTEGER :: ig, iga, igb, igz, igp, la, mu, iz, jz
    REAL(DP) :: L, S, z0, z
    REAL(DP) :: g(2), gp, gz
    COMPLEX(DP), PARAMETER :: ci = dcmplx(0.0d0, 1.0d0)
    COMPLEX(DP) :: rg3
    COMPLEX(DP) :: sum1p, sum1m, sum1c, sum2c, sum2p, sum2m
    REAL(DP)    :: z_l, z_r
    COMPLEX(DP) :: f1, f2, f3, f4, a0, a1, a2, a3
    COMPLEX(DP) :: poly_fr, poly_fl, poly_dfr, poly_dfl
    COMPLEX(DP) :: poly_a, poly_b, poly_c, poly_d
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dgp2_deps(2, 2)  !! dgp^2/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps

    COMPLEX(DP), ALLOCATABLE :: rhog3(:, :)
    COMPLEX(DP), ALLOCATABLE :: dVr_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: dVg_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: Vr(:)
    COMPLEX(DP), ALLOCATABLE :: Vg(:)

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0

    ! initialize
    sigmahar(:, :) = 0.0d0

    ALLOCATE (rhog3(dfftp%nr3, ngm_2d))
    ALLOCATE (dVr_deps(dfftp%nr3, 2, 2))
    ALLOCATE (dVg_deps(dfftp%nr3, 2, 2))
    ALLOCATE (Vr(dfftp%nr3))
    ALLOCATE (Vg(dfftp%nr3))

    ! reconstruct rho(gz,gp)
    rhog3(:, :) = (0.d0, 0.d0)

    DO ig = 1, ngm
      iga = mill(1, ig)
      igb = mill(2, ig)
      igz = mill(3, ig) + 1
      igp = imill_2d(iga, igb)
      IF (igz < 1) THEN
        igz = igz + dfftp%nr3
      END IF

      rg3 = rhog(ig)
      rhog3(igz, igp) = rg3

      ! expand function symmetrically to gz<0
      IF (gamma_only .and. iga == 0 .and. igb == 0) THEN
        igz = 1 - mill(3, ig)
        IF (igz < 1) THEN
          igz = igz + dfftp%nr3
        END IF
        rhog3(igz, igp) = CONJG(rg3)
      END IF
    END DO ! ig

    !****For gp!=0 case ********************
    DO igp = 1, ngm_2d
      iga = mill_2d(1, igp)
      igb = mill_2d(2, igp)
      g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
      gp = sqrt(g(1)*g(1) + g(2)*g(2))

      IF (gp == 0.0d0) CYCLE ! skip gp=0

      ! derivatives by strain tensor
      DO la = 1, 2
        DO mu = 1, 2
          dgp_deps(la, mu) = -g(la)*g(mu)/gp
          dgp2_deps(la, mu) = -g(la)*g(mu)*2.0d0
          dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
        END DO
      END DO

      ! summations over gz
      sum1p = (0.d0, 0.d0)
      sum1m = (0.d0, 0.d0)
      sum2p = (0.d0, 0.d0)
      sum2m = (0.d0, 0.d0)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, igp)
        sum1p = sum1p + rg3*QE_EXP(+ci*gz*z0)/(gp - ci*gz)
        sum1m = sum1m + rg3*QE_EXP(-ci*gz*z0)/(gp + ci*gz)
        sum2p = sum2p + rg3*QE_EXP(+ci*gz*z0)/(gp - ci*gz)**2
        sum2m = sum2m + rg3*QE_EXP(-ci*gz*z0)/(gp + ci*gz)**2
      END DO ! igz

      ! calculate dV(z)/deps
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        dVr_deps(iz, :, :) = &
          -(dgp_deps(:, :)*tpi/gp**2*(gp*(z - z0) - 1.0d0) &
          - delta(:, :)*tpi/gp) &
          *EXP(+gp*(z - z0))*sum1p &
          + dgp_deps(:, :)*tpi/gp*EXP(+gp*(z - z0))*sum2p &
          + (dgp_deps(:, :)*tpi/gp**2*(gp*(z + z0) + 1.0d0) &
          + delta(:, :)*tpi/gp) &
          *EXP(-gp*(z + z0))*sum1m &
          + dgp_deps(:, :)*tpi/gp*EXP(-gp*(z + z0))*sum2m
      END DO ! iz

      ! convert dV(z)/deps to dV(gz)/deps
      DO la = 1, 2
        DO mu = 1, 2
          CALL cft_1z(dVr_deps(:, la, mu), 1, dfftp%nr3, dfftp%nr3, -1, dVg_deps(:, la, mu))
        END DO
      END DO

      ! add bare coulomn terms to dV(gz)/deps
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, igp)

        dVg_deps(iz, :, :) = dVg_deps(iz, :, :) &
                              - delta(:, :)*fpi*rg3/(gp**2 + gz**2) &
                              - dgp2_deps(:, :)*fpi*rg3/(gp**2 + gz**2)**2
      END DO ! igz

      ! modifications
      IF (gamma_only) THEN
        dVg_deps(:, :, :) = dVg_deps(:, :, :)*2.0d0
      END IF

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, igp)
        sigmahar(1:2, 1:2) = sigmahar(1:2, 1:2) &
                             + REAL(CONJG(rg3)*dVg_deps(igz, :, :))
      END DO ! igz
    END DO ! igp

    !****For gp=0 case ********************
    IF (imill_2d(0, 0) > 0) THEN
      ! summations over gz
      sum1c = (0.d0, 0.d0)
      sum2c = (0.d0, 0.d0)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, imill_2d(0, 0))
        sum1c = sum1c + rg3*ci*cos(gz*z0)/gz
        sum2c = sum2c + rg3*cos(gz*z0)/gz**2
      END DO ! igz

      ! calculate V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        rg3 = rhog3(1, imill_2d(0, 0))
        Vr(iz) = &
          -tpi*z**2*rg3 &
          - tpi*z0**2*rg3 &
          - fpi*z*sum1c &
          - fpi*sum2c
      END DO ! iz

      ! separation by polynomial
      z_l = -z0
      z_r = +z0
      f1 = -tpi*z_r**2*rg3 &
           - tpi*z0**2*rg3 &
           - fpi*z_r*sum1c &
           - fpi*sum2c
      f2 = -tpi*z_l**2*rg3 &
           - tpi*z0**2*rg3 &
           - fpi*z_l*sum1c &
           - fpi*sum2c
      f3 = -fpi*z_r*rg3 &
           - fpi*sum1c
      f4 = -fpi*z_l*rg3 &
           - fpi*sum1c
      a0 = (f1*z_l**2*(z_l - 3.d0*z_r) + z_r*(f3*z_l**2*(-z_l + z_r) &
           + z_r*(f2*(3.d0*z_l - z_r) + f4*z_l*(-z_l + z_r))))/(z_l - z_r)**3
      a1 = (f3*z_l**3 + z_l*(6.d0*f1 - 6.d0*f2 + (f3 + 2.d0*f4)*z_l)*z_r &
            - (2*f3 + f4)*z_l*z_r**2 - f4*z_r**3)/(z_l - z_r)**3
      a2 = (-3*f1*(z_l + z_r) + 3.d0*f2*(z_l + z_r) - (z_l - z_r)*(2*f3*z_l &
           + f4*z_l + f3*z_r + 2*f4*z_r))/(z_l - z_r)**3
      a3 = (2.d0*f1 - 2.d0*f2 + (f3 + f4)*(z_l - z_r))/(z_l - z_r)**3

      ! remove polynomial from V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L
        Vr(iz) = Vr(iz) - (a0 + a1*z + a2*z**2 + a3*z**3)
      ENDDO

      ! convert V(z) to V(gz) without polynomial
      CALL cft_1z(Vr, 1, dfftp%nr3, dfftp%nr3, -1, Vg)

      ! add polynomial to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        Vg(iz) = Vg(iz) &
                  + a1*ci*cos(gz*z0)/gz &
                  + a2*2.0d0*cos(gz*z0)/gz**2 &
                  + a3*ci*z0**2*cos(gz*z0)/gz &
                  - a3*ci*6.0d0*cos(gz*z0)/gz**3
      END DO
      Vg(1) = Vg(1) + a0*1.0d0 + a2*z0**2/3.0d0

      ! add bare coulomn terms to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, imill_2d(0, 0))
        Vg(iz) = Vg(iz) + fpi*rg3/gz**2
      END DO ! igz

      ! calculate dV/deps(gz)
      DO igz = 1, dfftp%nr3
        dVg_deps(igz, :, :) = -delta(:, :)*Vg(igz)
      END DO ! igz

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, imill_2d(0, 0))
        sigmahar(1:2, 1:2) = sigmahar(1:2, 1:2) &
                             + REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz

    END IF ! imill_2d(0,0) > 0

    ! half means removing duplications.
    ! e2 means hartree -> Ry.
    sigmahar(:, :) = sigmahar(:, :)*(-0.5d0*e2)

    CALL mp_sum(sigmahar, intra_bgrp_comm)

    DEALLOCATE (rhog3)
    DEALLOCATE (dVr_deps)
    DEALLOCATE (dVg_deps)
    DEALLOCATE (Vr)
    DEALLOCATE (Vg)

    RETURN
  END SUBROUTINE esm_stres_har_bc1

  SUBROUTINE esm_stres_har_bc2(sigmahar, rhog)
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : ngm, mill
    USE constants,     ONLY : tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE control_flags, ONLY : gamma_only
    USE fft_base,      ONLY : dfftp
    USE fft_scalar,    ONLY : cft_1z
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(out)   :: sigmahar(3, 3)
    COMPLEX(DP), INTENT(in) :: rhog(ngm)   !  n(G)

    INTEGER :: ig, iga, igb, igz, igp, la, mu, iz, jz
    REAL(DP) :: L, S, z0, z1, z
    REAL(DP) :: g(2), gp, gz
    COMPLEX(DP), PARAMETER :: ci = dcmplx(0.0d0, 1.0d0)
    COMPLEX(DP) :: rg3
    COMPLEX(DP) :: sum1p, sum1m, sum1c, sum2c, sum2p, sum2m
    COMPLEX(DP) :: sum1sp, sum1sm, sum1cp, sum1cm, sum2sp, sum2sm
    REAL(DP)    :: z_l, z_r
    COMPLEX(DP) :: f1, f2, f3, f4, a0, a1, a2, a3
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dgp2_deps(2, 2)  !! dgp^2/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps

    COMPLEX(DP), ALLOCATABLE :: rhog3(:, :)
    COMPLEX(DP), ALLOCATABLE :: dVr_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: dVg_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: Vr(:)
    COMPLEX(DP), ALLOCATABLE :: Vg(:)

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0
    z1 = z0 + esm_w

    ! initialize
    sigmahar(:, :) = 0.0d0

    ALLOCATE (rhog3(dfftp%nr3, ngm_2d))
    ALLOCATE (dVr_deps(dfftp%nr3, 2, 2))
    ALLOCATE (dVg_deps(dfftp%nr3, 2, 2))
    ALLOCATE (Vr(dfftp%nr3))
    ALLOCATE (Vg(dfftp%nr3))

    ! reconstruct rho(gz,gp)
    rhog3(:, :) = (0.d0, 0.d0)

    DO ig = 1, ngm
      iga = mill(1, ig)
      igb = mill(2, ig)
      igz = mill(3, ig) + 1
      igp = imill_2d(iga, igb)
      IF (igz < 1) THEN
        igz = igz + dfftp%nr3
      END IF

      rg3 = rhog(ig)
      rhog3(igz, igp) = rg3

      ! expand function symmetrically to gz<0
      IF (gamma_only .and. iga == 0 .and. igb == 0) THEN
        igz = 1 - mill(3, ig)
        IF (igz < 1) THEN
          igz = igz + dfftp%nr3
        END IF
        rhog3(igz, igp) = CONJG(rg3)
      END IF
    END DO ! ig

    !****For gp!=0 case ********************
    DO igp = 1, ngm_2d
      iga = mill_2d(1, igp)
      igb = mill_2d(2, igp)
      g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
      gp = sqrt(g(1)*g(1) + g(2)*g(2))

      IF (gp == 0.0d0) CYCLE ! skip gp=0

      ! derivatives by strain tensor
      DO la = 1, 2
        DO mu = 1, 2
          dgp_deps(la, mu) = -g(la)*g(mu)/gp
          dgp2_deps(la, mu) = -g(la)*g(mu)*2.0d0
          dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
        END DO
      END DO

      ! summations over gz
      sum1p = (0.d0, 0.d0)
      sum1m = (0.d0, 0.d0)
      sum2p = (0.d0, 0.d0)
      sum2m = (0.d0, 0.d0)
      sum1sp = (0.d0, 0.d0)
      sum1sm = (0.d0, 0.d0)
      sum1cp = (0.d0, 0.d0)
      sum1cm = (0.d0, 0.d0)
      sum2sp = (0.d0, 0.d0)
      sum2sm = (0.d0, 0.d0)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, igp)

        sum1p = sum1p + rg3*QE_EXP(+ci*gz*z0)/(gp - ci*gz)
        sum1m = sum1m + rg3*QE_EXP(-ci*gz*z0)/(gp + ci*gz)
        sum2p = sum2p + rg3*QE_EXP(+ci*gz*z0)/(gp - ci*gz)**2
        sum2m = sum2m + rg3*QE_EXP(-ci*gz*z0)/(gp + ci*gz)**2

        sum1sp = sum1sp + rg3*QE_SINH(gp*z0 + ci*gz*z0)/(gp + ci*gz)
        sum1sm = sum1sm + rg3*QE_SINH(gp*z0 - ci*gz*z0)/(gp - ci*gz)

        sum1cp = sum1cp + rg3*QE_COSH(gp*z0 + ci*gz*z0)/(gp + ci*gz)*z0
        sum1cm = sum1cm + rg3*QE_COSH(gp*z0 - ci*gz*z0)/(gp - ci*gz)*z0

        sum2sp = sum2sp + rg3*QE_SINH(gp*z0 + ci*gz*z0)/(gp + ci*gz)**2
        sum2sm = sum2sm + rg3*QE_SINH(gp*z0 - ci*gz*z0)/(gp - ci*gz)**2
      END DO ! igz

      ! calculate dV(z)/deps
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        !! BC1 terms
        dVr_deps(iz, :, :) = &
          -(dgp_deps(:, :)*tpi/gp**2*(gp*(z - z0) - 1.0d0) &
          - delta(:, :)*tpi/gp) &
          *EXP(+gp*(z - z0))*sum1p &
          + dgp_deps(:, :)*tpi/gp*EXP(+gp*(z - z0))*sum2p &
          + (dgp_deps(:, :)*tpi/gp**2*(gp*(z + z0) + 1.0d0) &
          + delta(:, :)*tpi/gp) &
          *EXP(-gp*(z + z0))*sum1m &
          + dgp_deps(:, :)*tpi/gp*EXP(-gp*(z + z0))*sum2m

        !! BC2 terms
        dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                             + dgp_deps(:, :)*( &
                             - tpi/gp**2*(EXP(-gp*(z + 2*z1)) - EXP(+gp*z))/sinh(2*gp*z1) &
                             + tpi/gp*((-z - 2*z1)*EXP(-gp*(z + 2*z1)) - z*EXP(+gp*z))/sinh(2*gp*z1) &
                             - tpi/gp*(EXP(-gp*(z + 2*z1)) - EXP(+gp*z))/sinh(2*gp*z1)**2*2*z1*cosh(2*gp*z1)) &
                             * sum1sp &
                             + tpi/gp*(EXP(-gp*(z + 2*z1)) - EXP(+gp*z))/sinh(2*gp*z1) &
                             * (-delta(:, :)*sum1sp + dgp_deps(:, :)*(sum1cp - sum2sp)) &
                             + dgp_deps(:, :)*( &
                             - tpi/gp**2*(EXP(+gp*(z - 2*z1)) - EXP(-gp*z))/sinh(2*gp*z1) &
                             + tpi/gp*((+z - 2*z1)*EXP(+gp*(z - 2*z1)) + z*EXP(-gp*z))/sinh(2*gp*z1) &
                             - tpi/gp*(EXP(+gp*(z - 2*z1)) - EXP(-gp*z))/sinh(2*gp*z1)**2*2*z1*cosh(2*gp*z1)) &
                             * sum1sm &
                             + tpi/gp*(EXP(+gp*(z - 2*z1)) - EXP(-gp*z))/sinh(2*gp*z1) &
                             * (-delta(:, :)*sum1sm + dgp_deps(:, :)*(sum1cm - sum2sm))

      END DO ! iz

      ! convert dV(z)/deps to dV(gz)/deps
      DO la = 1, 2
        DO mu = 1, 2
          CALL cft_1z(dVr_deps(:, la, mu), 1, dfftp%nr3, dfftp%nr3, -1, dVg_deps(:, la, mu))
        END DO
      END DO

      ! add bare couloum terms to dV(gz)/deps
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, igp)

        dVg_deps(iz, :, :) = dVg_deps(iz, :, :) &
                              - delta(:, :)*fpi*rg3/(gp**2 + gz**2) &
                              - dgp2_deps(:, :)*fpi*rg3/(gp**2 + gz**2)**2
      END DO ! igz

      ! modifications
      IF (gamma_only) THEN
        dVg_deps(:, :, :) = dVg_deps(:, :, :)*2.0d0
      END IF

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, igp)
        sigmahar(1:2, 1:2) = sigmahar(1:2, 1:2) &
                             + REAL(CONJG(rg3)*dVg_deps(igz, :, :))
      END DO ! igz
    END DO ! igp

    !****For gp=0 case ********************
    IF (imill_2d(0, 0) > 0) THEN
      ! summations over gz
      sum1c = (0.d0, 0.d0)
      sum2c = (0.d0, 0.d0)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, imill_2d(0, 0))
        sum1c = sum1c + rg3*ci*cos(gz*z0)/gz
        sum2c = sum2c + rg3*cos(gz*z0)/gz**2
      END DO ! igz

      ! calculate V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        rg3 = rhog3(1, imill_2d(0, 0))

        !! BC1 terms
        Vr(iz) = &
          - tpi*z**2*rg3 &
          - tpi*z0**2*rg3 &
          - fpi*z*sum1c &
          - fpi*sum2c

        !! BC2 terms
        Vr(iz) = Vr(iz) &
                 + tpi*z1*2*z0*rg3 - tpi*(-z/z1)*2*z0*sum1c
      END DO ! iz

      ! separation by polynomial
      z_l = -z0
      z_r = +z0
      f1 = -tpi*z_r**2*rg3 &
           - tpi*z0**2*rg3 &
           - fpi*z_r*sum1c &
           - fpi*sum2c
      f1 = f1 &
           + tpi*z1*2*z0*rg3 - tpi*(-z_r/z1)*2*z0*sum1c

      f2 = -tpi*z_l**2*rg3 &
           - tpi*z0**2*rg3 &
           - fpi*z_l*sum1c &
           - fpi*sum2c
      f2 = f2 &
           + tpi*z1*2*z0*rg3 - tpi*(-z_l/z1)*2*z0*sum1c

      f3 = -fpi*z_r*rg3 &
           - fpi*sum1c
      f3 = f3 &
           - tpi*(-1.0d0/z1)*2*z0*sum1c

      f4 = -fpi*z_l*rg3 &
           - fpi*sum1c
      f4 = f4 &
           - tpi*(-1.0d0/z1)*2*z0*sum1c

      a0 = (f1*z_l**2*(z_l - 3.d0*z_r) + z_r*(f3*z_l**2*(-z_l + z_r) &
                                              + z_r*(f2*(3.d0*z_l - z_r) + f4*z_l*(-z_l + z_r))))/(z_l - z_r)**3
      a1 = (f3*z_l**3 + z_l*(6.d0*f1 - 6.d0*f2 + (f3 + 2.d0*f4)*z_l)*z_r &
            - (2*f3 + f4)*z_l*z_r**2 - f4*z_r**3)/(z_l - z_r)**3
      a2 = (-3*f1*(z_l + z_r) + 3.d0*f2*(z_l + z_r) - (z_l - z_r)*(2*f3*z_l &
                                                                   + f4*z_l + f3*z_r + 2*f4*z_r))/(z_l - z_r)**3
      a3 = (2.d0*f1 - 2.d0*f2 + (f3 + f4)*(z_l - z_r))/(z_l - z_r)**3

      ! remove polynomial from V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L
        Vr(iz) = Vr(iz) - (a0 + a1*z + a2*z**2 + a3*z**3)
      ENDDO

      ! convert V(z) to V(gz) without polynomial
      CALL cft_1z(Vr, 1, dfftp%nr3, dfftp%nr3, -1, Vg)

      ! add polynomial to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        Vg(iz) = Vg(iz) &
                  + a1*ci*cos(gz*z0)/gz &
                  + a2*2.0d0*cos(gz*z0)/gz**2 &
                  + a3*ci*z0**2*cos(gz*z0)/gz &
                  - a3*ci*6.0d0*cos(gz*z0)/gz**3
      END DO
      Vg(1) = Vg(1) + a0*1.0d0 + a2*z0**2/3.0d0

      ! add bare coulomn terms to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, imill_2d(0, 0))
        Vg(iz) = Vg(iz) + fpi*rg3/gz**2
      END DO ! igz

      ! calculate dV/deps(gz)
      DO igz = 1, dfftp%nr3
        dVg_deps(igz, :, :) = -delta(:, :)*Vg(igz)
      END DO ! igz

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, imill_2d(0, 0))
        sigmahar(1:2, 1:2) = sigmahar(1:2, 1:2) &
                             + REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    END IF ! imill_2d(0,0) > 0

    ! half means removing duplications.
    ! e2 means hartree -> Ry.
    sigmahar(:, :) = sigmahar(:, :)*(-0.5d0*e2)

    CALL mp_sum(sigmahar, intra_bgrp_comm)

    DEALLOCATE (rhog3)
    DEALLOCATE (dVr_deps)
    DEALLOCATE (dVg_deps)
    DEALLOCATE (Vr)
    DEALLOCATE (Vg)

    RETURN
  END SUBROUTINE esm_stres_har_bc2

  SUBROUTINE esm_stres_har_bc3(sigmahar, rhog)
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : ngm, mill
    USE constants,     ONLY : tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE control_flags, ONLY : gamma_only
    USE fft_base,      ONLY : dfftp
    USE fft_scalar,    ONLY : cft_1z
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(out)   :: sigmahar(3, 3)
    COMPLEX(DP), INTENT(in) :: rhog(ngm)   !  n(G)

    INTEGER :: ig, iga, igb, igz, igp, la, mu, iz, jz
    REAL(DP) :: L, S, z0, z1, z
    REAL(DP) :: g(2), gp, gz
    COMPLEX(DP), PARAMETER :: ci = dcmplx(0.0d0, 1.0d0)
    COMPLEX(DP) :: rg3
    COMPLEX(DP) :: sum1p, sum1m, sum2p, sum2m, sum1c, sum2c
    COMPLEX(DP) :: sum1sh, sum1ch, sum2sh
    REAL(DP)    :: z_l, z_r
    COMPLEX(DP) :: f1, f2, f3, f4, a0, a1, a2, a3
    COMPLEX(DP) :: poly_fr, poly_fl, poly_dfr, poly_dfl
    COMPLEX(DP) :: poly_a, poly_b, poly_c, poly_d
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dgp2_deps(2, 2)  !! dgp^2/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps

    COMPLEX(DP), ALLOCATABLE :: rhog3(:, :)
    COMPLEX(DP), ALLOCATABLE :: dVr_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: dVg_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: Vr(:)
    COMPLEX(DP), ALLOCATABLE :: Vg(:)

    REAL(DP) :: sigmahar_bc1(3, 3)

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0
    z1 = z0 + esm_w

    ! initialize
    sigmahar(:, :) = 0.0d0

    ALLOCATE (rhog3(dfftp%nr3, ngm_2d))
    ALLOCATE (dVr_deps(dfftp%nr3, 2, 2))
    ALLOCATE (dVg_deps(dfftp%nr3, 2, 2))
    ALLOCATE (Vr(dfftp%nr3))
    ALLOCATE (Vg(dfftp%nr3))

    ! reconstruct rho(gz,gp)
    rhog3(:, :) = (0.d0, 0.d0)

    DO ig = 1, ngm
      iga = mill(1, ig)
      igb = mill(2, ig)
      igz = mill(3, ig) + 1
      igp = imill_2d(iga, igb)
      IF (igz < 1) THEN
        igz = igz + dfftp%nr3
      END IF

      rg3 = rhog(ig)
      rhog3(igz, igp) = rg3

      ! expand function symmetrically to gz<0
      IF (gamma_only .and. iga == 0 .and. igb == 0) THEN
        igz = 1 - mill(3, ig)
        IF (igz < 1) THEN
          igz = igz + dfftp%nr3
        END IF
        rhog3(igz, igp) = CONJG(rg3)
      END IF
    END DO ! ig

    !****For gp!=0 case ********************
    DO igp = 1, ngm_2d
      iga = mill_2d(1, igp)
      igb = mill_2d(2, igp)
      g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
      gp = sqrt(g(1)*g(1) + g(2)*g(2))

      IF (gp == 0.0d0) CYCLE ! skip gp=0

      ! derivatives by strain tensor
      DO la = 1, 2
        DO mu = 1, 2
          dgp_deps(la, mu) = -g(la)*g(mu)/gp
          dgp2_deps(la, mu) = -g(la)*g(mu)*2.0d0
          dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
        END DO
      END DO

      ! summations over gz
      sum1p = (0.d0, 0.d0)
      sum1m = (0.d0, 0.d0)
      sum2p = (0.d0, 0.d0)
      sum2m = (0.d0, 0.d0)
      sum1sh = (0.d0, 0.d0)
      sum1ch = (0.d0, 0.d0)
      sum2sh = (0.d0, 0.d0)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, igp)
        sum1p = sum1p + rg3*QE_EXP(+ci*gz*z0)/(gp - ci*gz)
        sum1m = sum1m + rg3*QE_EXP(-ci*gz*z0)/(gp + ci*gz)
        sum2p = sum2p + rg3*QE_EXP(+ci*gz*z0)/(gp - ci*gz)**2
        sum2m = sum2m + rg3*QE_EXP(-ci*gz*z0)/(gp + ci*gz)**2
        sum1sh = sum1sh + rg3*QE_SINH(gp*z0 + ci*gz*z0)/(gp + ci*gz)
        sum1ch = sum1ch + rg3*QE_COSH(gp*z0 + ci*gz*z0)/(gp + ci*gz)*z0
        sum2sh = sum2sh + rg3*QE_SINH(gp*z0 + ci*gz*z0)/(gp + ci*gz)**2
      END DO ! igz

      ! calculate dV(z)/deps
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        !! BC1 terms
        dVr_deps(iz, :, :) = &
          -(dgp_deps(:, :)*tpi/gp**2*(gp*(z - z0) - 1.0d0) &
            - delta(:, :)*tpi/gp) &
          *EXP(+gp*(z - z0))*sum1p &
          + dgp_deps(:, :)*tpi/gp*EXP(+gp*(z - z0))*sum2p &
          + (dgp_deps(:, :)*tpi/gp**2*(gp*(z + z0) + 1.0d0) &
             + delta(:, :)*tpi/gp) &
          *EXP(-gp*(z + z0))*sum1m &
          + dgp_deps(:, :)*tpi/gp*EXP(-gp*(z + z0))*sum2m

        !! BC3 termn
        dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                             - dgp_deps(:, :)*( &
                             -fpi/gp**2*EXP(-gp*(-z + 2*z1)) &
                             - fpi/gp*(-z + 2*z1)*EXP(-gp*(-z + 2*z1)) &
                             )*sum1sh &
                             - fpi/gp*EXP(-gp*(-z + 2*z1))*( &
                             -delta(:, :)*sum1sh &
                             + dgp_deps(:, :)*(sum1ch - sum2sh))
      END DO ! iz

      ! convert dV(z)/deps to dV(gz)/deps
      DO la = 1, 2
        DO mu = 1, 2
          CALL cft_1z(dVr_deps(:, la, mu), 1, dfftp%nr3, dfftp%nr3, -1, dVg_deps(:, la, mu))
        END DO
      END DO

      ! add bare coulomn terms to dV(gz)/deps
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, igp)

        dVg_deps(iz, :, :) = dVg_deps(iz, :, :) &
                              - delta(:, :)*fpi*rg3/(gp**2 + gz**2) &
                              - dgp2_deps(:, :)*fpi*rg3/(gp**2 + gz**2)**2
      END DO ! igz

      ! modifications
      IF (gamma_only) THEN
        dVg_deps(:, :, :) = dVg_deps(:, :, :)*2.0d0
      END IF

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, igp)
        sigmahar(1:2, 1:2) = sigmahar(1:2, 1:2) &
                             + REAL(CONJG(rg3)*dVg_deps(igz, :, :))
      END DO ! igz
    END DO ! igp

    !****For gp=0 case ********************
    IF (imill_2d(0, 0) > 0) THEN
      ! summations over gz
      sum1c = (0.d0, 0.d0)
      sum2c = (0.d0, 0.d0)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, imill_2d(0, 0))
        sum1c = sum1c + rg3*ci*cos(gz*z0)/gz
        sum2c = sum2c + rg3*cos(gz*z0)/gz**2
      END DO ! igz

      ! calculate V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        rg3 = rhog3(1, imill_2d(0, 0))
        !! BC1 terms
        Vr(iz) = &
          - tpi*z**2*rg3 &
          - tpi*z0**2*rg3 &
          - fpi*z*sum1c &
          - fpi*sum2c

        !! BC3 terms
        Vr(iz) = Vr(iz) - tpi*(z - 2*z1)*2*z0*rg3 + fpi*z0*sum1c
      END DO ! iz

      ! separation by polynomial
      z_l = -z0
      z_r = +z0
      f1 = -tpi*z_r**2*rg3 &
           - tpi*z0**2*rg3 &
           - fpi*z_r*sum1c &
           - fpi*sum2c
      f1 = f1 &
           - tpi*(z_r - 2*z1)*2*z0*rg3 + fpi*z0*sum1c

      f2 = -tpi*z_l**2*rg3 &
           - tpi*z0**2*rg3 &
           - fpi*z_l*sum1c &
           - fpi*sum2c
      f2 = f2 &
           - tpi*(z_l - 2*z1)*2*z0*rg3 + fpi*z0*sum1c

      f3 = -fpi*z_r*rg3 &
           - fpi*sum1c
      f3 = f3 &
           - tpi*(1.0d0)*2*z0*rg3

      f4 = -fpi*z_l*rg3 &
           - fpi*sum1c
      f4 = f4 &
           - tpi*(1.0d0)*2*z0*rg3

      a0 = (f1*z_l**2*(z_l - 3.d0*z_r) + z_r*(f3*z_l**2*(-z_l + z_r) &
                                              + z_r*(f2*(3.d0*z_l - z_r) + f4*z_l*(-z_l + z_r))))/(z_l - z_r)**3
      a1 = (f3*z_l**3 + z_l*(6.d0*f1 - 6.d0*f2 + (f3 + 2.d0*f4)*z_l)*z_r &
            - (2*f3 + f4)*z_l*z_r**2 - f4*z_r**3)/(z_l - z_r)**3
      a2 = (-3*f1*(z_l + z_r) + 3.d0*f2*(z_l + z_r) - (z_l - z_r)*(2*f3*z_l &
                                                                   + f4*z_l + f3*z_r + 2*f4*z_r))/(z_l - z_r)**3
      a3 = (2.d0*f1 - 2.d0*f2 + (f3 + f4)*(z_l - z_r))/(z_l - z_r)**3

      ! remove polynomial from V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L
        Vr(iz) = Vr(iz) - (a0 + a1*z + a2*z**2 + a3*z**3)
      ENDDO

      ! convert V(z) to V(gz) without polynomial
      CALL cft_1z(Vr, 1, dfftp%nr3, dfftp%nr3, -1, Vg)

      ! add polynomial to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        Vg(iz) = Vg(iz) &
                  + a1*ci*cos(gz*z0)/gz &
                  + a2*2.0d0*cos(gz*z0)/gz**2 &
                  + a3*ci*z0**2*cos(gz*z0)/gz &
                  - a3*ci*6.0d0*cos(gz*z0)/gz**3
      END DO
      Vg(1) = Vg(1) + a0*1.0d0 + a2*z0**2/3.0d0

      ! add bare coulomn terms to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        rg3 = rhog3(iz, imill_2d(0, 0))
        Vg(iz) = Vg(iz) + fpi*rg3/gz**2
      END DO ! igz

      ! calculate dV/deps(gz)
      DO igz = 1, dfftp%nr3
        dVg_deps(igz, :, :) = -delta(:, :)*Vg(igz)
      END DO ! igz

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, imill_2d(0, 0))

        sigmahar(1:2, 1:2) = sigmahar(1:2, 1:2) &
                             + REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    END IF ! imill_2d(0,0) > 0

    ! half means removing duplications.
    ! e2 means hartree -> Ry.
    sigmahar(:, :) = sigmahar(:, :)*(-0.5d0*e2)

    CALL mp_sum(sigmahar, intra_bgrp_comm)

    DEALLOCATE (rhog3)
    DEALLOCATE (dVr_deps)
    DEALLOCATE (dVg_deps)
    DEALLOCATE (Vr)
    DEALLOCATE (Vg)

    RETURN
  END SUBROUTINE esm_stres_har_bc3

  SUBROUTINE esm_stres_ewr(alpha, sigmaewa)
    USE kinds, ONLY : DP
    IMPLICIT NONE

    REAL(DP), INTENT(in)  :: alpha
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    SELECT CASE (esm_bc)
    CASE ('pbc')
      STOP 'esm_stres_ewa must not be called for esm_bc = pbc'
    CASE ('bc1')
      CALL esm_stres_ewr_pbc(alpha, sigmaewa)
    CASE ('bc2')
      CALL esm_stres_ewr_pbc(alpha, sigmaewa)
    CASE ('bc3')
      CALL esm_stres_ewr_pbc(alpha, sigmaewa)
    CASE ('bc4')
      STOP 'esm_stres_ewa has not yet implemented for esm_bc = bc4'
    END SELECT

    RETURN
  END SUBROUTINE esm_stres_ewr

  SUBROUTINE esm_stres_ewr_pbc(alpha, sigmaewa)
    USE kinds,     ONLY : DP
    USE constants, ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base, ONLY : omega, alat, at, tpiba, bg
    USE ions_base, ONLY : zv, nat, tau, ityp
    USE gvect,     ONLY : gstart
    USE mp_bands,  ONLY : intra_bgrp_comm
    USE mp,        ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(in)  :: alpha
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    INTEGER, PARAMETER :: mxr = 50
    ! the maximum number of R vectors included in r sum
    INTEGER  :: ia, ib, nr, nrm, la, mu
    REAL(DP) :: Qa, Qb, dtau(3), rmax
    REAL(DP) :: salp, r(3, mxr), r2(mxr), rr, fac

    salp = sqrt(alpha)

    ! initialize
    sigmaewa(:, :) = 0.d0

    !
    ! R-space sum here (only for the processor that contains G=0)
    !
    IF (gstart == 2) THEN
      rmax = 4.0d0/salp/alat
      !
      ! with this choice terms up to ZiZj*erfc(5) are counted (erfc(5)=2x10^-1
      !
      DO ib = 1, nat
        Qb = (-1.0d0)*zv(ityp(ib))
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          !
          !     generates nearest-neighbors shells r(i)=R(i)-dtau(i)
          !
          dtau(:) = tau(:, ib) - tau(:, ia)
          CALL rgen(dtau, rmax, mxr, at, bg, r, r2, nrm)

          DO nr = 1, nrm
            rr = sqrt(r2(nr))*alat
            r(:, nr) = r(:, nr)*alat

            fac = Qb*Qa/rr**3 &
                  *(erfc(salp*rr) &
                  + rr*2.0d0*salp*sqrtpm1*EXP(-alpha*rr**2))
            DO la = 1, 3
              DO mu = 1, 3
                sigmaewa(la, mu) = sigmaewa(la, mu) + fac*r(la, nr)*r(mu, nr)
              END DO ! mu
            END DO ! la
          END DO ! nr
        END DO ! ia
      END DO ! ib
    END IF

    sigmaewa(:, :) = sigmaewa(:, :)*(e2/2.0d0/omega)

    CALL mp_sum(sigmaewa, intra_bgrp_comm)

    RETURN
  END SUBROUTINE esm_stres_ewr_pbc

  SUBROUTINE esm_stres_ewg(alpha, sigmaewa)
    USE kinds, ONLY : DP
    IMPLICIT NONE

    REAL(DP), INTENT(in)  :: alpha
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    SELECT CASE (esm_bc)
    CASE ('pbc')
      STOP 'esm_stres_ewa must not be called for esm_bc = pbc'
    CASE ('bc1')
      CALL esm_stres_ewg_bc1(alpha, sigmaewa)
    CASE ('bc2')
      CALL esm_stres_ewg_bc2(alpha, sigmaewa)
    CASE ('bc3')
      CALL esm_stres_ewg_bc3(alpha, sigmaewa)
    CASE ('bc4')
      STOP 'esm_stres_ewa must not be called for esm_bc = bc4'
    END SELECT

    RETURN
  END SUBROUTINE esm_stres_ewg

  SUBROUTINE esm_stres_ewg_bc1(alpha, sigmaewa)
    USE kinds,         ONLY : DP
    USE constants,     ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE ions_base,     ONLY : zv, nat, tau, ityp
    USE control_flags, ONLY : gamma_only
    USE gvect,         ONLY : gstart
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(in)  :: alpha
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    INTEGER  :: ia, ib, igp, iga, igb, la, mu, iz
    REAL(DP) :: L, S, salp
    REAL(DP) :: Qa, Qb, ra(2), rb(2), za, zb
    REAL(DP) :: g(2), gp, Vr
    REAL(DP) :: cosgpr, experfcm, experfcp, dexperfcm_dgp, dexperfcp_dgp
    REAL(DP) :: zbza, isalp, gpzbza, gp2a, mgazz, pgazz, fact

    REAL(DP) :: dE_deps(2, 2)
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    salp  = sqrt(alpha)
    isalp = 1.0_DP/salp
    fact  = 0.5_DP * isalp

    ! initialize
    sigmaewa(:, :) = 0.0d0

    !****For gp!=0 case ********************
    DO ib = 1, nat
      Qb = (-1.0d0)*zv(ityp(ib))
      rb(1:2) = tau(1:2, ib)*alat
      zb = tau(3, ib)*alat
      IF (zb > L*0.5d0) THEN
        zb = zb - L
      END IF

      DO ia = 1, nat
        Qa = (-1.0d0)*zv(ityp(ia))
        ra(1:2) = tau(1:2, ia)*alat
        za = tau(3, ia)*alat
        IF (za > L*0.5d0) THEN
          za = za - L
        END IF

        ! distance between atoms is defined
        zbza = zb - za

        ! summations over gp
        dE_deps(:, :) = 0.0d0
        DO igp = 1, ngm_2d
          iga = mill_2d(1, igp)
          igb = mill_2d(2, igp)
          g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
          gp = sqrt(g(1)*g(1) + g(2)*g(2))

          IF (gp == 0.0d0) CYCLE ! skip gp=0

          ! define exp phases
          gpzbza = gp * zbza
          gp2a   = gp * 0.5_DP * isalp
          mgazz  = gp2a - salp * zbza
          pgazz  = gp2a + salp * zbza

          ! derivatives by strain tensor
          DO la = 1, 2
            DO mu = 1, 2
              dgp_deps(la, mu) = -g(la)*g(mu)/gp
              dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
            END DO
          END DO

          ! coefficients
          cosgpr = cos(g(1)*(rb(1) - ra(1)) + g(2)*(rb(2) - ra(2)))
          experfcm = exp_erfc(-gpzbza, mgazz)
          experfcp = exp_erfc(+gpzbza, pgazz)
          dexperfcm_dgp = -zbza*experfcm &
                          -exp_gauss( -gpzbza, mgazz ) * fact
          dexperfcp_dgp = +zbza*experfcp &
                          -exp_gauss( +gpzbza, pgazz ) * fact
          !
          ! Old code is not safe, because diverged terms are included.
          ! However, this code is a faithful for original formula.
          ! For this reason, we leave following old codes as comment.
          !
          ! experfcm = exp_erfc(-gp*(zb - za), gp/2.d0/salp - salp*(zb - za))
          ! experfcp = exp_erfc(+gp*(zb - za), gp/2.d0/salp + salp*(zb - za))
          ! dexperfcm_dgp = -(zb - za)*exp_erfc(-gp*(zb - za), gp/2.d0/salp - salp*(zb - za)) &
          !                 - EXP(-gp*(zb - za))*qe_gauss(gp/2.d0/salp - salp*(zb - za))/2.d0/salp
          ! dexperfcp_dgp = +(zb - za)*exp_erfc(+gp*(zb - za), gp/2.d0/salp + salp*(zb - za)) &
          !                 - EXP(+gp*(zb - za))*qe_gauss(gp/2.d0/salp + salp*(zb - za))/2.d0/salp
          !
          dE_deps(:, :) = dE_deps(:, :) &
                          + gp*dinvgp_deps(:, :)*pi/gp*Qb*Qa/S*cosgpr*experfcm &
                          - pi/gp*delta(:, :)*Qb*Qa/S*cosgpr*experfcm &
                          + pi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*dexperfcm_dgp &
                          + gp*dinvgp_deps(:, :)*pi/gp*Qb*Qa/S*cosgpr*experfcp &
                          - pi/gp*delta(:, :)*Qb*Qa/S*cosgpr*experfcp &
                          + pi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*dexperfcp_dgp
        END DO ! igp

        ! modifications
        IF (gamma_only) THEN
          dE_deps(:, :) = dE_deps(:, :)*2.0d0
        END IF

        ! calculate stress tensor
        sigmaewa(1:2, 1:2) = sigmaewa(1:2, 1:2) - dE_deps(1:2, 1:2)/omega

      END DO ! ia
    END DO ! ib

    !****For gp=0 case ********************
    IF (gstart == 2) THEN
      DO ib = 1, nat
        Qb = (-1.0d0)*zv(ityp(ib))
        rb(1:2) = tau(1:2, ib)*alat
        zb = tau(3, ib)*alat
        IF (zb > L*0.5d0) THEN
          zb = zb - L
        END IF

        Vr = 0.0d0
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          ra(1:2) = tau(1:2, ia)*alat
          za = tau(3, ia)*alat
          IF (za > L*0.5d0) THEN
            za = za - L
          END IF

          Vr = Vr - tpi*Qa/S &
               *((zb - za)*erf(salp*(zb - za)) &
                 + EXP(-alpha*(zb - za)**2)*sqrtpm1/salp)
        END DO ! ia

        dE_deps(1:2, 1:2) = -delta(1:2, 1:2)*Vr*Qb

        ! calculate stress tensor
        sigmaewa(1:2, 1:2) = sigmaewa(1:2, 1:2) - dE_deps(1:2, 1:2)/omega
      END DO ! ib
    END IF

    ! half means removing duplications.
    ! e2 means hartree -> Ry.
    sigmaewa(:, :) = sigmaewa(:, :)*(0.5d0*e2)

    CALL mp_sum(sigmaewa, intra_bgrp_comm)

    RETURN
  END SUBROUTINE esm_stres_ewg_bc1

  SUBROUTINE esm_stres_ewg_bc2(alpha, sigmaewa)
    USE kinds,         ONLY : DP
    USE constants,     ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE ions_base,     ONLY : zv, nat, tau, ityp
    USE control_flags, ONLY : gamma_only
    USE gvect,         ONLY : gstart
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(in)  :: alpha
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    INTEGER  :: ia, ib, igp, iga, igb, la, mu
    REAL(DP) :: L, S, salp, z0, z1
    REAL(DP) :: Qa, Qb, ra(2), rb(2), za, zb
    REAL(DP) :: g(2), gp, Vr
    REAL(DP) :: cosgpr, experfcm, experfcp, dexperfcm_dgp, dexperfcp_dgp
    REAL(DP) :: exph1, exph2, exph3
    REAL(DP) :: zbza, isalp, gpzbza, gp2a, mgazz, pgazz, fact

    REAL(DP) :: dE_deps(2, 2)
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0
    z1 = z0 + esm_w
    salp = sqrt(alpha)
    isalp = 1.0_DP/salp
    fact  = 0.5_DP * isalp

    ! initialize
    sigmaewa(:, :) = 0.0d0

    !****For gp!=0 case ********************
    DO ib = 1, nat
      Qb = (-1.0d0)*zv(ityp(ib))
      rb(1:2) = tau(1:2, ib)*alat
      zb = tau(3, ib)*alat
      IF (zb > L*0.5d0) THEN
        zb = zb - L
      END IF

      DO ia = 1, nat
        Qa = (-1.0d0)*zv(ityp(ia))
        ra(1:2) = tau(1:2, ia)*alat
        za = tau(3, ia)*alat
        IF (za > L*0.5d0) THEN
          za = za - L
        END IF

        ! distance between atoms is defined
        zbza = zb - za

        ! summations over gp
        dE_deps(:, :) = 0.0d0
        DO igp = 1, ngm_2d
          iga = mill_2d(1, igp)
          igb = mill_2d(2, igp)
          g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
          gp = sqrt(g(1)*g(1) + g(2)*g(2))

          IF (gp == 0.0d0) CYCLE ! skip gp=0

          ! define exp phases
          gpzbza = gp * zbza
          gp2a   = gp * 0.5_DP * isalp
          mgazz  = gp2a - salp * zbza
          pgazz  = gp2a + salp * zbza

          ! derivatives by strain tensor
          DO la = 1, 2
            DO mu = 1, 2
              dgp_deps(la, mu) = -g(la)*g(mu)/gp
              dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
            END DO
          END DO

          ! coefficients
          cosgpr = cos(g(1)*(rb(1) - ra(1)) + g(2)*(rb(2) - ra(2)))
          experfcm = exp_erfc(-gpzbza, mgazz)
          experfcp = exp_erfc(+gpzbza, pgazz)
          dexperfcm_dgp = -zbza*experfcm &
                          -exp_gauss( -gpzbza, mgazz ) * fact
          dexperfcp_dgp = +zbza*experfcp &
                          -exp_gauss( +gpzbza, pgazz ) * fact
          !
          ! Old code is not safe, because diverged terms are included.
          ! However, this code is a faithful for original formula.
          ! For this reason, we leave following old codes as comment.
          !
          ! experfcm = exp_erfc(-gp*(zb - za), gp/2.d0/salp - salp*(zb - za))
          ! experfcp = exp_erfc(+gp*(zb - za), gp/2.d0/salp + salp*(zb - za))
          ! dexperfcm_dgp = -(zb - za)*exp_erfc(-gp*(zb - za), gp/2.d0/salp - salp*(zb - za)) &
          !                 - EXP(-gp*(zb - za))*qe_gauss(gp/2.d0/salp - salp*(zb - za))/2.d0/salp
          ! dexperfcp_dgp = +(zb - za)*exp_erfc(+gp*(zb - za), gp/2.d0/salp + salp*(zb - za)) &
          !                 - EXP(+gp*(zb - za))*qe_gauss(gp/2.d0/salp + salp*(zb - za))/2.d0/salp
          !
          exph1 = (cosh(gp*(zb - za))*EXP(-2*gp*z1) - cosh(gp*(zb + za)))/sinh(2*gp*z1)
          exph2 = ((zb - za)*sinh(gp*(zb - za))*EXP(-2*gp*z1) &
                  - 2*z1*cosh(gp*(zb - za))*EXP(-2*gp*z1) &
                  - (zb + za)*sinh(gp*(zb + za)))/sinh(2*gp*z1)
          exph3 = -(cosh(gp*(zb - za))*EXP(-2*gp*z1) - cosh(gp*(zb + za)))/sinh(2*gp*z1)**2*2*z1*cosh(2*gp*z1)

          !! BC1 terms
          dE_deps(:, :) = dE_deps(:, :) &
                          + gp*dinvgp_deps(:, :)*pi/gp*Qb*Qa/S*cosgpr*experfcm &
                          - pi/gp*delta(:, :)*Qb*Qa/S*cosgpr*experfcm &
                          + pi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*dexperfcm_dgp &
                          + gp*dinvgp_deps(:, :)*pi/gp*Qb*Qa/S*cosgpr*experfcp &
                          - pi/gp*delta(:, :)*Qb*Qa/S*cosgpr*experfcp &
                          + pi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*dexperfcp_dgp

          !! BC2 terms
          dE_deps(:, :) = dE_deps(:, :) &
                          + gp*dinvgp_deps(:, :)*tpi/gp*Qb*Qa/S*cosgpr*exph1 &
                          - tpi/gp*delta(:, :)*Qb*Qa/S*cosgpr*exph1 &
                          + tpi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*(exph2 + exph3)
        END DO ! igp

        ! modifications
        IF (gamma_only) THEN
          dE_deps(:, :) = dE_deps(:, :)*2.0d0
        END IF

        ! calculate stress tensor
        sigmaewa(1:2, 1:2) = sigmaewa(1:2, 1:2) - dE_deps(1:2, 1:2)/omega

      END DO ! ia
    END DO ! ib

    !****For gp=0 case ********************
    IF (gstart == 2) THEN
      DO ib = 1, nat
        Qb = (-1.0d0)*zv(ityp(ib))
        rb(1:2) = tau(1:2, ib)*alat
        zb = tau(3, ib)*alat
        IF (zb > L*0.5d0) THEN
          zb = zb - L
        END IF

        ! [note] this Vr does not contain a term due to efield z*efield
        ! because it vanishes in the differentiation with respect to strain.
        Vr = 0.0d0
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          ra(1:2) = tau(1:2, ia)*alat
          za = tau(3, ia)*alat
          IF (za > L*0.5d0) THEN
            za = za - L
          END IF

          !! BC1 terms
          Vr = Vr - tpi*Qa/S &
               *((zb - za)*erf(salp*(zb - za)) &
                 + EXP(-alpha*(zb - za)**2)*sqrtpm1/salp)

          !! BC2 terms
          Vr = Vr + tpi*Qa/S*(-zb*za + z1*z1)/z1
        END DO ! ia

        dE_deps(1:2, 1:2) = -delta(1:2, 1:2)*Vr*Qb

        ! calculate stress tensor
        sigmaewa(1:2, 1:2) = sigmaewa(1:2, 1:2) - dE_deps(1:2, 1:2)/omega
      END DO ! ib
    END IF

    ! half means removing duplications.
    ! e2 means hartree -> Ry.
    sigmaewa(:, :) = sigmaewa(:, :)*(0.5d0*e2)

    CALL mp_sum(sigmaewa, intra_bgrp_comm)

    RETURN
  END SUBROUTINE esm_stres_ewg_bc2

  SUBROUTINE esm_stres_ewg_bc3(alpha, sigmaewa)
    USE kinds,         ONLY : DP
    USE constants,     ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE ions_base,     ONLY : zv, nat, tau, ityp
    USE control_flags, ONLY : gamma_only
    USE gvect,         ONLY : gstart
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(in)  :: alpha
    REAL(DP), INTENT(out) :: sigmaewa(3, 3)

    INTEGER  :: ia, ib, igp, iga, igb, la, mu
    REAL(DP) :: L, S, salp, z0, z1
    REAL(DP) :: Qa, Qb, ra(2), rb(2), za, zb
    REAL(DP) :: g(2), gp, Vr
    REAL(DP) :: cosgpr, experfcm, experfcp, dexperfcm_dgp, dexperfcp_dgp, expm
    REAL(DP) :: zbza, isalp, gpzbza, gp2a, mgazz, pgazz, fact

    REAL(DP) :: dE_deps(2, 2)
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0
    z1 = z0 + esm_w
    salp = sqrt(alpha)
    isalp = 1.0_DP/salp
    fact  = 0.5_DP * isalp

    ! initialize
    sigmaewa(:, :) = 0.0d0

    !****For gp!=0 case ********************
    DO ib = 1, nat
      Qb = (-1.0d0)*zv(ityp(ib))
      rb(1:2) = tau(1:2, ib)*alat
      zb = tau(3, ib)*alat
      IF (zb > L*0.5d0) THEN
        zb = zb - L
      END IF

      DO ia = 1, nat
        Qa = (-1.0d0)*zv(ityp(ia))
        ra(1:2) = tau(1:2, ia)*alat
        za = tau(3, ia)*alat
        IF (za > L*0.5d0) THEN
          za = za - L
        END IF

        ! distance between atoms is defined
        zbza = zb - za

        ! summations over gp
        dE_deps(:, :) = 0.0d0
        DO igp = 1, ngm_2d
          iga = mill_2d(1, igp)
          igb = mill_2d(2, igp)
          g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
          gp = sqrt(g(1)*g(1) + g(2)*g(2))

          IF (gp == 0.0d0) CYCLE ! skip gp=0

          ! define exp phases
          gpzbza = gp * zbza
          gp2a   = gp * 0.5_DP * isalp
          mgazz  = gp2a - salp * zbza
          pgazz  = gp2a + salp * zbza

          ! derivatives by strain tensor
          DO la = 1, 2
            DO mu = 1, 2
              dgp_deps(la, mu) = -g(la)*g(mu)/gp
              dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
            END DO
          END DO

          ! coefficients
          cosgpr = cos(g(1)*(rb(1) - ra(1)) + g(2)*(rb(2) - ra(2)))
          experfcm = exp_erfc(-gpzbza, mgazz)
          experfcp = exp_erfc(+gpzbza, pgazz)
          dexperfcm_dgp = -zbza*experfcm &
                          -exp_gauss( -gpzbza, mgazz ) * fact
          dexperfcp_dgp = +zbza*experfcp &
                          -exp_gauss( +gpzbza, pgazz ) * fact
          expm = exp(-gp*(-zb + 2*z1 - za))
          !
          ! Old code is not safe, because diverged terms are included.
          ! However, this code is a faithful for original formula.
          ! For this reason, we leave following old codes as comment.
          !
          ! experfcm = exp_erfc(-gp*(zb - za), gp/2.d0/salp - salp*(zb - za))
          ! experfcp = exp_erfc(+gp*(zb - za), gp/2.d0/salp + salp*(zb - za))
          ! dexperfcm_dgp = -(zb - za)*exp_erfc(-gp*(zb - za), gp/2.d0/salp - salp*(zb - za)) &
          !                 - EXP(-gp*(zb - za))*qe_gauss(gp/2.d0/salp - salp*(zb - za))/2.d0/salp
          ! dexperfcp_dgp = +(zb - za)*exp_erfc(+gp*(zb - za), gp/2.d0/salp + salp*(zb - za)) &
          !                 - EXP(+gp*(zb - za))*qe_gauss(gp/2.d0/salp + salp*(zb - za))/2.d0/salp
          ! expm = EXP(-gp*(-zb + 2*z1 - za))
          !
          !! BC1 terms
          dE_deps(:, :) = dE_deps(:, :) &
                          + gp*dinvgp_deps(:, :)*pi/gp*Qb*Qa/S*cosgpr*experfcm &
                          - pi/gp*delta(:, :)*Qb*Qa/S*cosgpr*experfcm &
                          + pi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*dexperfcm_dgp &
                          + gp*dinvgp_deps(:, :)*pi/gp*Qb*Qa/S*cosgpr*experfcp &
                          - pi/gp*delta(:, :)*Qb*Qa/S*cosgpr*experfcp &
                          + pi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*dexperfcp_dgp

          !! BC3 terms
          !call errore('esm_stress_ewg_bc2','FIXME: expm used but not defined',1)
          dE_deps(:, :) = dE_deps(:, :) &
                          - gp*dinvgp_deps(:, :)*tpi/gp*Qb*Qa/S*cosgpr*expm &
                          + tpi/gp*delta(:, :)*Qb*Qa/S*cosgpr*expm &
                          + tpi/gp*Qb*Qa/S*cosgpr*dgp_deps(:, :)*(-zb + 2*z1 - za)*expm
        END DO ! igp

        ! modifications
        IF (gamma_only) THEN
          dE_deps(:, :) = dE_deps(:, :)*2.0d0
        END IF

        ! calculate stress tensor
        sigmaewa(1:2, 1:2) = sigmaewa(1:2, 1:2) - dE_deps(1:2, 1:2)/omega

      END DO ! ia
    END DO ! ib

    !****For gp=0 case ********************
    IF (gstart == 2) THEN
      DO ib = 1, nat
        Qb = (-1.0d0)*zv(ityp(ib))
        rb(1:2) = tau(1:2, ib)*alat
        zb = tau(3, ib)*alat
        IF (zb > L*0.5d0) THEN
          zb = zb - L
        END IF

        Vr = 0.0d0
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          ra(1:2) = tau(1:2, ia)*alat
          za = tau(3, ia)*alat
          IF (za > L*0.5d0) THEN
            za = za - L
          END IF

          !! BC1 terms
          Vr = Vr - tpi*Qa/S &
               *((zb - za)*erf(salp*(zb - za)) &
               + EXP(-alpha*(zb - za)**2)*sqrtpm1/salp)

          !! BC3 terms
          Vr = Vr + tpi*Qa/S*(-zb + 2*z1 - za)
        END DO ! ia

        dE_deps(1:2, 1:2) = -delta(1:2, 1:2)*Vr*Qb

        ! calculate stress tensor
        sigmaewa(1:2, 1:2) = sigmaewa(1:2, 1:2) - dE_deps(1:2, 1:2)/omega
      END DO ! ib
    END IF

    ! half means removing duplications.
    ! e2 means hartree -> Ry.
    sigmaewa(:, :) = sigmaewa(:, :)*(0.5d0*e2)

    CALL mp_sum(sigmaewa, intra_bgrp_comm)

    RETURN
  END SUBROUTINE esm_stres_ewg_bc3

  SUBROUTINE esm_stres_loclong_bc1(sigmaloclong, rhog)
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : ngm, mill
    USE constants,     ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE ions_base,     ONLY : zv, nat, tau, ityp
    USE control_flags, ONLY : gamma_only
    USE fft_base,      ONLY : dfftp
    USE fft_scalar,    ONLY : cft_1z
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(out) :: sigmaloclong(3, 3)
    COMPLEX(DP) :: rhog(ngm)   !  n(G)

    INTEGER  :: ig, iga, igb, igz, igp, la, mu, iz, jz, ia
    REAL(DP) :: L, S, z0, alpha, salp, z
    REAL(DP) :: Qa, ra(2), za
    REAL(DP) :: g(2), gp, gz
    COMPLEX(DP), PARAMETER :: ci = dcmplx(0.0d0, 1.0d0)
    COMPLEX(DP) :: rg3
    COMPLEX(DP) :: expimgpr, experfcm, experfcp, dexperfcm_dgp, dexperfcp_dgp
    REAL(DP)    :: z_r, z_l
    COMPLEX(DP) :: a0, a1, a2, a3, f1, f2, f3, f4
    COMPLEX(DP) :: poly_fr, poly_fl, poly_dfr, poly_dfl
    COMPLEX(DP) :: poly_a, poly_b, poly_c, poly_d
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps
    REAL(DP) :: isalp, fact
    REAL(DP) :: zza, mgza, pgza, g2a_maza, g2a_paza

    COMPLEX(DP), ALLOCATABLE :: rhog3(:, :)
    COMPLEX(DP), ALLOCATABLE :: dVr_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: dVg_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: Vr(:)
    COMPLEX(DP), ALLOCATABLE :: Vg(:)

    ALLOCATE (rhog3(dfftp%nr3, ngm_2d))
    ALLOCATE (dVr_deps(dfftp%nr3, 2, 2))
    ALLOCATE (dVg_deps(dfftp%nr3, 2, 2))
    ALLOCATE (Vr(dfftp%nr3))
    ALLOCATE (Vg(dfftp%nr3))

    ! reconstruct rho(gz,gp)
    rhog3(:, :) = (0.d0, 0.d0)

    DO ig = 1, ngm
      iga = mill(1, ig)
      igb = mill(2, ig)
      igz = mill(3, ig) + 1
      igp = imill_2d(iga, igb)
      IF (igz < 1) THEN
        igz = igz + dfftp%nr3
      END IF

      rg3 = rhog(ig)
      rhog3(igz, igp) = rg3

      IF (gamma_only .and. iga == 0 .and. igb == 0) THEN
        igz = 1 - mill(3, ig)
        IF (igz < 1) THEN
          igz = igz + dfftp%nr3
        END IF
        rhog3(igz, igp) = CONJG(rg3)
      ENDIF
    END DO ! ig

    ! cell settings
    L = at(3, 3)*alat
    S     = omega/L
    z0    = L/2.d0
    alpha = 1.0d0
    salp  = sqrt(alpha)
    ! useful values are setted
    isalp = 1.d0/sqrt(alpha)
    fact  = 1.d0/2.d0/salp
    ! initialize
    sigmaloclong(:, :) = 0.0d0

    !****For gp!=0 case ********************
    DO igp = 1, ngm_2d
      iga = mill_2d(1, igp)
      igb = mill_2d(2, igp)
      g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
      gp = sqrt(g(1)*g(1) + g(2)*g(2))

      IF (gp == 0.0d0) CYCLE ! skip gp=0

      ! derivatives by strain tensor
      DO la = 1, 2
        DO mu = 1, 2
          dgp_deps(la, mu) = -g(la)*g(mu)/gp
          dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
        END DO
      END DO

      ! calculate dV(z)/deps
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        ! summations over all atoms
        dVr_deps(iz, :, :) = (0.0d0, 0.0d0)
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          ra(1:2) = tau(1:2, ia)*alat
          za = tau(3, ia)*alat
          IF (za > L*0.5d0) THEN
            za = za - L
          END IF
! --------------------------------------------------------------------------------------------------
!         Following code is old version. However, this code explicitly shows the formulation of 
!         stress tensor within ESM scheme. For this reason, we left the code as comment.
!          expimgpr = qe_exp(-ci*(g(1)*ra(1) + g(2)*ra(2)))
!          experfcm = exp_erfc(-gp*(z - za), gp/2.d0/salp - salp*(z - za))
!          experfcp = exp_erfc(+gp*(z - za), gp/2.d0/salp + salp*(z - za))
!          dexperfcm_dgp = -(z - za)*exp_erfc(-gp*(z - za), gp/2.d0/salp - salp*(z - za)) &
!                          - EXP(-gp*(z - za))*qe_gauss(gp/2.d0/salp - salp*(z - za))/2.d0/salp
!          dexperfcp_dgp = +(z - za)*exp_erfc(+gp*(z - za), gp/2.d0/salp + salp*(z - za)) &
!                          - EXP(+gp*(z - za))*qe_gauss(gp/2.d0/salp + salp*(z - za))/2.d0/salp
!---------------------------------------------------------------------------------------------------
          !
          ! ... Set useful values
          zza = z - za
          mgza = -gp*zza
          pgza = +gp*zza
          g2a_maza = gp*0.5d0*isalp - salp * zza
          g2a_paza = gp*0.5d0*isalp + salp * zza
          !
          expimgpr = QE_EXP(-ci*(g(1)*ra(1) + g(2)*ra(2)))
          experfcm = exp_erfc(mgza, g2a_maza)
          experfcp = exp_erfc(pgza, g2a_paza)
          dexperfcm_dgp = -zza * experfcm &
                          -exp_gauss( mgza, g2a_maza ) * fact
          dexperfcp_dgp = +zza * experfcp &
                          -exp_gauss( pgza, g2a_paza ) * fact
          !
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*pi/gp*Qa/S*expimgpr*experfcm &
                               - pi/gp*delta(:, :)*Qa/S*expimgpr*experfcm &
                               + pi/gp*Qa/S*expimgpr*dgp_deps(:, :)*dexperfcm_dgp

          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*pi/gp*Qa/S*expimgpr*experfcp &
                               - pi/gp*delta(:, :)*Qa/S*expimgpr*experfcp &
                               + pi/gp*Qa/S*expimgpr*dgp_deps(:, :)*dexperfcp_dgp
        END DO ! ia
      END DO ! iz

      ! convert dV(z)/deps to dV(gz)/deps
      DO la = 1, 2
        DO mu = 1, 2
          CALL cft_1z(dVr_deps(:, la, mu), 1, dfftp%nr3, dfftp%nr3, -1, dVg_deps(:, la, mu))
        END DO
      END DO

      ! modifications
      IF (gamma_only) THEN
        dVg_deps(:, :, :) = dVg_deps(:, :, :)*2.0d0
      END IF

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, igp)
        sigmaloclong(1:2, 1:2) = sigmaloclong(1:2, 1:2) &
                                 - REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    END DO ! igp

    !****For gp=0 case ********************
    IF (imill_2d(0, 0) > 0) THEN
      ! calculate V(z)
      Vr(:) = 0.0d0
      ! separation by polynomial
      f1 = (0.d0, 0.d0); f2 = (0.d0, 0.d0); f3 = (0.d0, 0.d0); f4 = (0.d0, 0.d0)
      z_l = -z0
      z_r = +z0
      DO ia = 1, nat
        Qa = (-1.0d0)*zv(ityp(ia))
        ra(1:2) = tau(1:2, ia)*alat
        za = tau(3, ia)*alat
        IF (za > L*0.5d0) THEN
          za = za - L
        END IF

        DO iz = 1, dfftp%nr3
          jz = iz - 1
          IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
            jz = jz - dfftp%nr3
          END IF
          z = DBLE(jz) / DBLE(dfftp%nr3) * L

          Vr(iz) = Vr(iz) - tpi*Qa/S &
                   *((z - za)*erf(salp*(z - za)) &
                     + EXP(-alpha*(z - za)**2)*sqrtpm1/salp)
        END DO ! iz

        f1 = f1 - tpi*Qa/S &
             *((z_r - za)*erf(salp*(z_r - za)) &
               + EXP(-alpha*(z_r - za)**2)*sqrtpm1/salp)
        f2 = f2 - tpi*Qa/S &
             *((z_l - za)*erf(salp*(z_l - za)) &
               + EXP(-alpha*(z_l - za)**2)*sqrtpm1/salp)
        f3 = f3 - tpi*Qa/S &
             *erf(salp*(z_r - za))
        f4 = f4 - tpi*Qa/S &
             *erf(salp*(z_l - za))
      END DO ! ia

      a0 = (f1*z_l**2*(z_l - 3.d0*z_r) + z_r*(f3*z_l**2*(-z_l + z_r) &
                                              + z_r*(f2*(3.d0*z_l - z_r) + f4*z_l*(-z_l + z_r))))/(z_l - z_r)**3
      a1 = (f3*z_l**3 + z_l*(6.d0*f1 - 6.d0*f2 + (f3 + 2.d0*f4)*z_l)*z_r &
            - (2*f3 + f4)*z_l*z_r**2 - f4*z_r**3)/(z_l - z_r)**3
      a2 = (-3*f1*(z_l + z_r) + 3.d0*f2*(z_l + z_r) - (z_l - z_r)*(2*f3*z_l &
                                                                   + f4*z_l + f3*z_r + 2*f4*z_r))/(z_l - z_r)**3
      a3 = (2.d0*f1 - 2.d0*f2 + (f3 + f4)*(z_l - z_r))/(z_l - z_r)**3

      ! remove polynomial from V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L
        Vr(iz) = Vr(iz) - (a0 + a1*z + a2*z**2 + a3*z**3)
      ENDDO

      ! convert V(z) to V(gz) without polynomial
      CALL cft_1z(Vr, 1, dfftp%nr3, dfftp%nr3, -1, Vg)

      ! add polynomial to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        Vg(iz) = Vg(iz) &
                  + a1*ci*cos(gz*z0)/gz &
                  + a2*2.0d0*cos(gz*z0)/gz**2 &
                  + a3*ci*z0**2*cos(gz*z0)/gz &
                  - a3*ci*6.0d0*cos(gz*z0)/gz**3
      END DO
      Vg(1) = Vg(1) + a0*1.0d0 + a2*z0**2/3.0d0

      ! calculate dV/deps(gz)
      DO igz = 1, dfftp%nr3
        dVg_deps(igz, :, :) = -delta(:, :)*Vg(igz)
      END DO ! igz

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, imill_2d(0, 0))
        sigmaloclong(1:2, 1:2) = sigmaloclong(1:2, 1:2) &
                                 - REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    ENDIF ! imill_2d(0,0) > 0

    ! e2 means hartree -> Ry.
    sigmaloclong(:, :) = sigmaloclong(:, :)*(e2)

    CALL mp_sum(sigmaloclong, intra_bgrp_comm)

    DEALLOCATE (rhog3)
    DEALLOCATE (dVr_deps)
    DEALLOCATE (dVg_deps)
    DEALLOCATE (Vr)
    DEALLOCATE (Vg)

    RETURN
  END SUBROUTINE esm_stres_loclong_bc1

  SUBROUTINE esm_stres_loclong_bc2(sigmaloclong, rhog)
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : ngm, mill
    USE constants,     ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE ions_base,     ONLY : zv, nat, tau, ityp
    USE control_flags, ONLY : gamma_only
    USE fft_base,      ONLY : dfftp
    USE fft_scalar,    ONLY : cft_1z
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(out) :: sigmaloclong(3, 3)
    COMPLEX(DP) :: rhog(ngm)   !  n(G)

    INTEGER  :: ig, iga, igb, igz, igp, la, mu, iz, jz, ia
    REAL(DP) :: L, S, z0, z1, alpha, salp, z
    REAL(DP) :: Qa, ra(2), za
    REAL(DP) :: g(2), gp, gz
    COMPLEX(DP), PARAMETER :: ci = dcmplx(0.0d0, 1.0d0)
    COMPLEX(DP) :: rg3
    COMPLEX(DP) :: expimgpr, experfcm, experfcp, dexperfcm_dgp, dexperfcp_dgp
    COMPLEX(DP) :: exph1, exph2, exph3
    REAL(DP)    :: z_r, z_l
    COMPLEX(DP) :: a0, a1, a2, a3, f1, f2, f3, f4
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps
    REAL(DP) :: isalp, fact
    REAL(DP) :: zza, mgza, pgza, g2a_maza, g2a_paza

    COMPLEX(DP), ALLOCATABLE :: rhog3(:, :)
    COMPLEX(DP), ALLOCATABLE :: dVr_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: dVg_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: Vr(:)
    COMPLEX(DP), ALLOCATABLE :: Vg(:)

    ALLOCATE (rhog3(dfftp%nr3, ngm_2d))
    ALLOCATE (dVr_deps(dfftp%nr3, 2, 2))
    ALLOCATE (dVg_deps(dfftp%nr3, 2, 2))
    ALLOCATE (Vr(dfftp%nr3))
    ALLOCATE (Vg(dfftp%nr3))

    ! reconstruct rho(gz,gp)
    rhog3(:, :) = (0.d0, 0.d0)

    DO ig = 1, ngm
      iga = mill(1, ig)
      igb = mill(2, ig)
      igz = mill(3, ig) + 1
      igp = imill_2d(iga, igb)
      IF (igz < 1) THEN
        igz = igz + dfftp%nr3
      END IF

      rg3 = rhog(ig)
      rhog3(igz, igp) = rg3

      IF (gamma_only .and. iga == 0 .and. igb == 0) THEN
        igz = 1 - mill(3, ig)
        IF (igz < 1) THEN
          igz = igz + dfftp%nr3
        END IF
        rhog3(igz, igp) = CONJG(rg3)
      ENDIF
    END DO ! ig

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0
    z1 = z0 + esm_w
    alpha = 1.0d0
    salp = sqrt(alpha)
    ! useful values are setted
    isalp = 1.d0/sqrt(alpha)
    fact  = 1.d0/2.d0/salp
    ! initialize
    sigmaloclong(:, :) = 0.0d0

    !****For gp!=0 case ********************
    DO igp = 1, ngm_2d
      iga = mill_2d(1, igp)
      igb = mill_2d(2, igp)
      g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
      gp = sqrt(g(1)*g(1) + g(2)*g(2))

      IF (gp == 0.0d0) CYCLE ! skip gp=0

      ! derivatives by strain tensor
      DO la = 1, 2
        DO mu = 1, 2
          dgp_deps(la, mu) = -g(la)*g(mu)/gp
          dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
        END DO
      END DO

      ! calculate dV(z)/deps
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        ! summations over all atoms
        dVr_deps(iz, :, :) = (0.0d0, 0.0d0)
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          ra(1:2) = tau(1:2, ia)*alat
          za = tau(3, ia)*alat
          IF (za > L*0.5d0) THEN
            za = za - L
          END IF
! --------------------------------------------------------------------------------------------------
!         Following code is old version. However, this code explicitly shows the formulation of 
!         stress tensor within ESM scheme. For this reason, we left the code as comment.
!          expimgpr = qe_exp(-ci*(g(1)*ra(1) + g(2)*ra(2)))
!          experfcm = exp_erfc(-gp*(z - za), gp/2.d0/salp - salp*(z - za))
!          experfcp = exp_erfc(+gp*(z - za), gp/2.d0/salp + salp*(z - za))
!          dexperfcm_dgp = -(z - za)*exp_erfc(-gp*(z - za), gp/2.d0/salp - salp*(z - za)) &
!                          - EXP(-gp*(z - za))*qe_gauss(gp/2.d0/salp - salp*(z - za))/2.d0/salp
!          dexperfcp_dgp = +(z - za)*exp_erfc(+gp*(z - za), gp/2.d0/salp + salp*(z - za)) &
!                          - EXP(+gp*(z - za))*qe_gauss(gp/2.d0/salp + salp*(z - za))/2.d0/salp
!---------------------------------------------------------------------------------------------------
          !
          ! ... Set useful values
          zza = z - za
          mgza = -gp*zza
          pgza = +gp*zza
          g2a_maza = gp*0.5d0*isalp - salp * zza
          g2a_paza = gp*0.5d0*isalp + salp * zza
          !
          expimgpr = QE_EXP(-ci*(g(1)*ra(1) + g(2)*ra(2)))
          experfcm = exp_erfc(mgza, g2a_maza)
          experfcp = exp_erfc(pgza, g2a_paza)
          dexperfcm_dgp = -zza * experfcm &
                          -exp_gauss( mgza, g2a_maza ) * fact
          dexperfcp_dgp = +zza * experfcp &
                          -exp_gauss( pgza, g2a_paza ) * fact
          !
          exph1 = (cosh(gp*(z - za))*EXP(-2*gp*z1) - cosh(gp*(z + za)))/sinh(2*gp*z1)
          exph2 = ((z - za)*sinh(gp*(z - za))*EXP(-2*gp*z1) &
                   - 2*z1*cosh(gp*(z - za))*EXP(-2*gp*z1) &
                   - (z + za)*sinh(gp*(z + za)))/sinh(2*gp*z1)
          exph3 = -(cosh(gp*(z - za))*EXP(-2*gp*z1) - cosh(gp*(z + za)))/sinh(2*gp*z1)**2*2*z1*cosh(2*gp*z1)

          !! BC1 terms
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*pi/gp*Qa/S*expimgpr*experfcm &
                               - pi/gp*delta(:, :)*Qa/S*expimgpr*experfcm &
                               + pi/gp*Qa/S*expimgpr*dgp_deps(:, :)*dexperfcm_dgp

          !! BC1 terms
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*pi/gp*Qa/S*expimgpr*experfcp &
                               - pi/gp*delta(:, :)*Qa/S*expimgpr*experfcp &
                               + pi/gp*Qa/S*expimgpr*dgp_deps(:, :)*dexperfcp_dgp

          !! BC2 terms
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*tpi/gp*Qa/S*expimgpr*exph1 &
                               - tpi/gp*delta(:, :)*Qa/S*expimgpr*exph1 &
                               + tpi/gp*Qa/S*expimgpr*dgp_deps(:, :)*(exph2 + exph3)
        END DO ! ia
      END DO ! iz

      ! convert dV(z)/deps to dV(gz)/deps
      DO la = 1, 2
        DO mu = 1, 2
          CALL cft_1z(dVr_deps(:, la, mu), 1, dfftp%nr3, dfftp%nr3, -1, dVg_deps(:, la, mu))
        END DO
      END DO

      ! modifications
      IF (gamma_only) THEN
        dVg_deps(:, :, :) = dVg_deps(:, :, :)*2.0d0
      END IF

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, igp)
        sigmaloclong(1:2, 1:2) = sigmaloclong(1:2, 1:2) &
                                 - REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    END DO ! igp

    !****For gp=0 case ********************
    IF (imill_2d(0, 0) > 0) THEN
      ! calculate V(z)
      ! [note] this Vr does not contain a term due to efield z*efield
      ! because it vanishes in the differentiation with respect to strain.
      Vr(:) = 0.0d0
      ! separation by polynomial
      f1 = (0.d0, 0.d0); f2 = (0.d0, 0.d0); f3 = (0.d0, 0.d0); f4 = (0.d0, 0.d0)
      z_l = -z0
      z_r = +z0
      DO ia = 1, nat
        Qa = (-1.0d0)*zv(ityp(ia))
        ra(1:2) = tau(1:2, ia)*alat
        za = tau(3, ia)*alat
        IF (za > L*0.5d0) THEN
          za = za - L
        END IF

        DO iz = 1, dfftp%nr3
          jz = iz - 1
          IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
            jz = jz - dfftp%nr3
          END IF
          z = DBLE(jz) / DBLE(dfftp%nr3) * L

          !! BC1 terms
          Vr(iz) = Vr(iz) - tpi*Qa/S &
                   *((z - za)*erf(salp*(z - za)) &
                     + EXP(-alpha*(z - za)**2)*sqrtpm1/salp)

          !! BC2 terms
          Vr(iz) = Vr(iz) + tpi*Qa/S*(-z*za + z1*z1)/z1
        END DO ! iz

        f1 = f1 - tpi*Qa/S &
             *((z_r - za)*erf(salp*(z_r - za)) &
               + EXP(-alpha*(z_r - za)**2)*sqrtpm1/salp)
        f1 = f1 + tpi*Qa/S*(-z_r*za + z1*z1)/z1

        f2 = f2 - tpi*Qa/S &
             *((z_l - za)*erf(salp*(z_l - za)) &
               + EXP(-alpha*(z_l - za)**2)*sqrtpm1/salp)
        f2 = f2 + tpi*Qa/S*(-z_l*za + z1*z1)/z1

        f3 = f3 - tpi*Qa/S &
             *erf(salp*(z_r - za))
        f3 = f3 + tpi*Qa/S*(-za)/z1

        f4 = f4 - tpi*Qa/S &
             *erf(salp*(z_l - za))
        f4 = f4 + tpi*Qa/S*(-za)/z1

      END DO ! ia

      a0 = (f1*z_l**2*(z_l - 3.d0*z_r) + z_r*(f3*z_l**2*(-z_l + z_r) &
           + z_r*(f2*(3.d0*z_l - z_r) + f4*z_l*(-z_l + z_r))))/(z_l - z_r)**3
      a1 = (f3*z_l**3 + z_l*(6.d0*f1 - 6.d0*f2 + (f3 + 2.d0*f4)*z_l)*z_r &
           - (2*f3 + f4)*z_l*z_r**2 - f4*z_r**3)/(z_l - z_r)**3
      a2 = (-3*f1*(z_l + z_r) + 3.d0*f2*(z_l + z_r) - (z_l - z_r)*(2*f3*z_l &
           + f4*z_l + f3*z_r + 2*f4*z_r))/(z_l - z_r)**3
      a3 = (2.d0*f1 - 2.d0*f2 + (f3 + f4)*(z_l - z_r))/(z_l - z_r)**3

      ! remove polynomial from V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L
        Vr(iz) = Vr(iz) - (a0 + a1*z + a2*z**2 + a3*z**3)
      ENDDO

      ! convert V(z) to V(gz) without polynomial
      CALL cft_1z(Vr, 1, dfftp%nr3, dfftp%nr3, -1, Vg)

      ! add polynomial to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        Vg(iz) = Vg(iz) &
                  + a1*ci*cos(gz*z0)/gz &
                  + a2*2.0d0*cos(gz*z0)/gz**2 &
                  + a3*ci*z0**2*cos(gz*z0)/gz &
                  - a3*ci*6.0d0*cos(gz*z0)/gz**3
      END DO
      Vg(1) = Vg(1) + a0*1.0d0 + a2*z0**2/3.0d0

      ! calculate dV/deps(gz)
      DO igz = 1, dfftp%nr3
        dVg_deps(igz, :, :) = -delta(:, :)*Vg(igz)
      END DO ! igz

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, imill_2d(0, 0))
        sigmaloclong(1:2, 1:2) = sigmaloclong(1:2, 1:2) &
                                 - REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    ENDIF ! imill_2d(0,0) > 0

    ! e2 means hartree -> Ry.
    sigmaloclong(:, :) = sigmaloclong(:, :)*(e2)

    CALL mp_sum(sigmaloclong, intra_bgrp_comm)

    DEALLOCATE (rhog3)
    DEALLOCATE (dVr_deps)
    DEALLOCATE (dVg_deps)
    DEALLOCATE (Vr)
    DEALLOCATE (Vg)

    RETURN
  END SUBROUTINE esm_stres_loclong_bc2

  SUBROUTINE esm_stres_loclong_bc3(sigmaloclong, rhog)
    USE kinds,         ONLY : DP
    USE gvect,         ONLY : ngm, mill
    USE constants,     ONLY : pi, sqrtpm1, tpi, fpi, e2
    USE cell_base,     ONLY : omega, alat, at, tpiba, bg
    USE ions_base,     ONLY : zv, nat, tau, ityp
    USE control_flags, ONLY : gamma_only
    USE fft_base,      ONLY : dfftp
    USE fft_scalar,    ONLY : cft_1z
    USE mp_bands,      ONLY : intra_bgrp_comm
    USE mp,            ONLY : mp_sum
    IMPLICIT NONE

    REAL(DP), INTENT(out) :: sigmaloclong(3, 3)
    COMPLEX(DP) :: rhog(ngm)   !  n(G)

    INTEGER  :: ig, iga, igb, igz, igp, la, mu, iz, jz, ia
    REAL(DP) :: L, S, z0, z1, alpha, salp, z
    REAL(DP) :: Qa, ra(2), za
    REAL(DP) :: g(2), gp, gz
    COMPLEX(DP), PARAMETER :: ci = dcmplx(0.0d0, 1.0d0)
    COMPLEX(DP) :: rg3
    COMPLEX(DP) :: expimgpr, experfcm, experfcp, dexperfcm_dgp, dexperfcp_dgp
    COMPLEX(DP) :: expm
    REAL(DP)    :: z_r, z_l
    COMPLEX(DP) :: a0, a1, a2, a3, f1, f2, f3, f4
    REAL(DP), PARAMETER :: delta(2, 2) = reshape((/1.0d0, 0.0d0, 0.0d0, 1.0d0/), (/2, 2/))
    REAL(DP) :: dgp_deps(2, 2)  !! dgp/deps
    REAL(DP) :: dinvgp_deps(2, 2)  !! dgp^-1/deps
    REAL(DP) :: isalp, fact
    REAL(DP) :: zza, mgza, pgza, g2a_maza, g2a_paza, tz1

    COMPLEX(DP), ALLOCATABLE :: rhog3(:, :)
    COMPLEX(DP), ALLOCATABLE :: dVr_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: dVg_deps(:, :, :)
    COMPLEX(DP), ALLOCATABLE :: Vr(:)
    COMPLEX(DP), ALLOCATABLE :: Vg(:)

    REAL(DP) :: sigmaloclong_bc1(3, 3)

    ALLOCATE (rhog3(dfftp%nr3, ngm_2d))
    ALLOCATE (dVr_deps(dfftp%nr3, 2, 2))
    ALLOCATE (dVg_deps(dfftp%nr3, 2, 2))
    ALLOCATE (Vr(dfftp%nr3))
    ALLOCATE (Vg(dfftp%nr3))

    ! reconstruct rho(gz,gp)
    rhog3(:, :) = (0.d0, 0.d0)

    DO ig = 1, ngm
      iga = mill(1, ig)
      igb = mill(2, ig)
      igz = mill(3, ig) + 1
      igp = imill_2d(iga, igb)
      IF (igz < 1) THEN
        igz = igz + dfftp%nr3
      END IF

      rg3 = rhog(ig)
      rhog3(igz, igp) = rg3

      IF (gamma_only .and. iga == 0 .and. igb == 0) THEN
        igz = 1 - mill(3, ig)
        IF (igz < 1) THEN
          igz = igz + dfftp%nr3
        END IF
        rhog3(igz, igp) = CONJG(rg3)
      ENDIF
    END DO ! ig

    ! cell settings
    L = at(3, 3)*alat
    S = omega/L
    z0 = L/2.d0
    z1 = z0 + esm_w
    alpha = 1.0d0
    salp = sqrt(alpha)
    ! useful values are setted
    isalp = 1.d0/sqrt(alpha)
    fact  = 1.d0/2.d0/salp
    tz1 = 2.d0 * z1
    !
    ! initialize
    sigmaloclong(:, :) = 0.0d0

    !****For gp!=0 case ********************
    DO igp = 1, ngm_2d
      iga = mill_2d(1, igp)
      igb = mill_2d(2, igp)
      g(1:2) = (iga*bg(1:2, 1) + igb*bg(1:2, 2))*tpiba
      gp = sqrt(g(1)*g(1) + g(2)*g(2))

      IF (gp == 0.0d0) CYCLE ! skip gp=0

      ! derivatives by strain tensor
      DO la = 1, 2
        DO mu = 1, 2
          dgp_deps(la, mu) = -g(la)*g(mu)/gp
          dinvgp_deps(la, mu) = +g(la)*g(mu)/gp**3
        END DO
      END DO

      ! calculate dV(z)/deps
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L

        ! summations over all atoms
        dVr_deps(iz, :, :) = (0.0d0, 0.0d0)
        DO ia = 1, nat
          Qa = (-1.0d0)*zv(ityp(ia))
          ra(1:2) = tau(1:2, ia)*alat
          za = tau(3, ia)*alat
          IF (za > L*0.5d0) THEN
            za = za - L
          END IF
! --------------------------------------------------------------------------------------------------
!         Following code is old version. However, this code explicitly shows the formulation of 
!         stress tensor within ESM scheme. For this reason, we left the code as comment.
!          expimgpr = qe_exp(-ci*(g(1)*ra(1) + g(2)*ra(2)))
!          experfcm = exp_erfc(-gp*(z - za), gp/2.d0/salp - salp*(z - za))
!          experfcp = exp_erfc(+gp*(z - za), gp/2.d0/salp + salp*(z - za))
!          dexperfcm_dgp = -(z - za)*exp_erfc(-gp*(z - za), gp/2.d0/salp - salp*(z - za)) &
!                          - EXP(-gp*(z - za))*qe_gauss(gp/2.d0/salp - salp*(z - za))/2.d0/salp
!          dexperfcp_dgp = +(z - za)*exp_erfc(+gp*(z - za), gp/2.d0/salp + salp*(z - za)) &
!                          - EXP(+gp*(z - za))*qe_gauss(gp/2.d0/salp + salp*(z - za))/2.d0/salp
!---------------------------------------------------------------------------------------------------
          !
          ! ... Set useful values
          zza = z - za
          mgza = -gp*zza
          pgza = +gp*zza
          g2a_maza = gp*0.5d0*isalp - salp * zza
          g2a_paza = gp*0.5d0*isalp + salp * zza
          !
          expimgpr = QE_EXP(-ci*(g(1)*ra(1) + g(2)*ra(2)))
          experfcm = exp_erfc(mgza, g2a_maza)
          experfcp = exp_erfc(pgza, g2a_paza)
          dexperfcm_dgp = -zza * experfcm &
                          -exp_gauss( mgza, g2a_maza ) * fact
          dexperfcp_dgp = +zza * experfcp &
                          -exp_gauss( pgza, g2a_paza ) * fact
          !
          expm = EXP(-gp*(-z + tz1 - za))

          !! BC1 terms
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*pi/gp*Qa/S*expimgpr*experfcm &
                               - pi/gp*delta(:, :)*Qa/S*expimgpr*experfcm &
                               + pi/gp*Qa/S*expimgpr*dgp_deps(:, :)*dexperfcm_dgp

          !! BC1 terms
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               + gp*dinvgp_deps(:, :)*pi/gp*Qa/S*expimgpr*experfcp &
                               - pi/gp*delta(:, :)*Qa/S*expimgpr*experfcp &
                               + pi/gp*Qa/S*expimgpr*dgp_deps(:, :)*dexperfcp_dgp

          !! BC3 terms
          dVr_deps(iz, :, :) = dVr_deps(iz, :, :) &
                               - gp*dinvgp_deps(:, :)*tpi/gp*Qa/S*expimgpr*expm &
                               + tpi/gp*delta(:, :)*Qa/S*expimgpr*expm &
                               + tpi/gp*Qa/S*expimgpr*dgp_deps(:, :)*(-z + 2*z1 - za)*expm
        END DO ! ia
      END DO ! iz

      ! convert dV(z)/deps to dV(gz)/deps
      DO la = 1, 2
        DO mu = 1, 2
          CALL cft_1z(dVr_deps(:, la, mu), 1, dfftp%nr3, dfftp%nr3, -1, dVg_deps(:, la, mu))
        END DO
      END DO

      ! modifications
      IF (gamma_only) THEN
        dVg_deps(:, :, :) = dVg_deps(:, :, :)*2.0d0
      END IF

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, igp)
        sigmaloclong(1:2, 1:2) = sigmaloclong(1:2, 1:2) &
                                 - REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    END DO ! igp

    !****For gp=0 case ********************
    IF (imill_2d(0, 0) > 0) THEN
      ! calculate V(z)
      Vr(:) = 0.0d0
      ! separation by polynomial
      f1 = (0.d0, 0.d0); f2 = (0.d0, 0.d0); f3 = (0.d0, 0.d0); f4 = (0.d0, 0.d0)
      z_l = -z0
      z_r = +z0
      DO ia = 1, nat
        Qa = (-1.0d0)*zv(ityp(ia))
        ra(1:2) = tau(1:2, ia)*alat
        za = tau(3, ia)*alat
        IF (za > L*0.5d0) THEN
          za = za - L
        END IF

        DO iz = 1, dfftp%nr3
          jz = iz - 1
          IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
            jz = jz - dfftp%nr3
          END IF
          z = DBLE(jz) / DBLE(dfftp%nr3) * L

          !! BC1 terms
          Vr(iz) = Vr(iz) - tpi*Qa/S &
                   *((z - za)*erf(salp*(z - za)) &
                   + EXP(-alpha*(z - za)**2)*sqrtpm1/salp)

          !! BC3 terms
          Vr(iz) = Vr(iz) + tpi*Qa/S*(-z + 2*z1 - za)
        END DO ! iz

        f1 = f1 - tpi*Qa/S &
             *((z_r - za)*erf(salp*(z_r - za)) &
             + EXP(-alpha*(z_r - za)**2)*sqrtpm1/salp)
        f1 = f1 + tpi*Qa/S*(-z_r + 2*z1 - za)

        f2 = f2 - tpi*Qa/S &
             *((z_l - za)*erf(salp*(z_l - za)) &
             + EXP(-alpha*(z_l - za)**2)*sqrtpm1/salp)
        f2 = f2 + tpi*Qa/S*(-z_l + 2*z1 - za)

        f3 = f3 - tpi*Qa/S &
             *erf(salp*(z_r - za))
        f3 = f3 + tpi*Qa/S*(-1.0d0)

        f4 = f4 - tpi*Qa/S &
             *erf(salp*(z_l - za))
        f4 = f4 + tpi*Qa/S*(-1.0d0)
      END DO ! ia

      a0 = (f1*z_l**2*(z_l - 3.d0*z_r) + z_r*(f3*z_l**2*(-z_l + z_r) &
           + z_r*(f2*(3.d0*z_l - z_r) + f4*z_l*(-z_l + z_r))))/(z_l - z_r)**3
      a1 = (f3*z_l**3 + z_l*(6.d0*f1 - 6.d0*f2 + (f3 + 2.d0*f4)*z_l)*z_r &
           - (2*f3 + f4)*z_l*z_r**2 - f4*z_r**3)/(z_l - z_r)**3
      a2 = (-3*f1*(z_l + z_r) + 3.d0*f2*(z_l + z_r) - (z_l - z_r)*(2*f3*z_l &
           + f4*z_l + f3*z_r + 2*f4*z_r))/(z_l - z_r)**3
      a3 = (2.d0*f1 - 2.d0*f2 + (f3 + f4)*(z_l - z_r))/(z_l - z_r)**3

      ! remove polynomial from V(z)
      DO iz = 1, dfftp%nr3
        jz = iz - 1
        IF (jz >= (dfftp%nr3 - dfftp%nr3/2)) THEN
          jz = jz - dfftp%nr3
        END IF
        z = DBLE(jz) / DBLE(dfftp%nr3) * L
        Vr(iz) = Vr(iz) - (a0 + a1*z + a2*z**2 + a3*z**3)
      ENDDO

      ! convert V(z) to V(gz) without polynomial
      CALL cft_1z(Vr, 1, dfftp%nr3, dfftp%nr3, -1, Vg)

      ! add polynomial to V(gz)
      DO igz = -(dfftp%nr3 - 1)/2, (dfftp%nr3 - 1)/2
        IF (igz == 0) CYCLE
        iz = igz + 1
        IF (iz < 1) iz = iz + dfftp%nr3
        gz = dble(igz)*tpi/L

        Vg(iz) = Vg(iz) &
                  + a1*ci*cos(gz*z0)/gz &
                  + a2*2.0d0*cos(gz*z0)/gz**2 &
                  + a3*ci*z0**2*cos(gz*z0)/gz &
                  - a3*ci*6.0d0*cos(gz*z0)/gz**3
      END DO
      Vg(1) = Vg(1) + a0*1.0d0 + a2*z0**2/3.0d0

      ! calculate dV/deps(gz)
      DO igz = 1, dfftp%nr3
        dVg_deps(igz, :, :) = -delta(:, :)*Vg(igz)
      END DO ! igz

      ! calculate stress tensor
      DO igz = 1, dfftp%nr3
        rg3 = rhog3(igz, imill_2d(0, 0))
        sigmaloclong(1:2, 1:2) = sigmaloclong(1:2, 1:2) &
                                 - REAL(CONJG(rg3)*dVg_deps(igz, 1:2, 1:2))
      END DO ! igz
    ENDIF ! imill_2d(0,0) > 0

    ! e2 means hartree -> Ry.
    sigmaloclong(:, :) = sigmaloclong(:, :)*(e2)

    CALL mp_sum(sigmaloclong, intra_bgrp_comm)

    DEALLOCATE (rhog3)
    DEALLOCATE (dVr_deps)
    DEALLOCATE (dVg_deps)
    DEALLOCATE (Vr)
    DEALLOCATE (Vg)

    RETURN
  END SUBROUTINE esm_stres_loclong_bc3

  COMPLEX(DP) FUNCTION qe_exp(x)
    COMPLEX(DP), INTENT(in) :: x
    REAL(DP) :: r, i, c, s

    r = dreal(x)
    i = dimag(x)
    c = cos(i)
    s = sin(i)

    qe_exp = EXP(r)*cmplx(c, s, kind=DP)

  END FUNCTION qe_exp

  COMPLEX(DP) FUNCTION qe_sinh(x)
    COMPLEX(DP), INTENT(in) :: x
    REAL(DP) :: r, i, c, s

    r = dreal(x)
    i = dimag(x)
    c = cos(i)
    s = sin(i)

    qe_sinh = 0.5d0*(EXP(r)*cmplx(c, s, kind=DP) - EXP(-r)*cmplx(c, -s, kind=DP))

  END FUNCTION qe_sinh

  COMPLEX(DP) FUNCTION qe_cosh(x)
    COMPLEX(DP), INTENT(in) :: x
    REAL(DP) :: r, i, c, s

    r = dreal(x)
    i = dimag(x)
    c = cos(i)
    s = sin(i)

    qe_cosh = 0.5d0*(EXP(r)*cmplx(c, s, kind=DP) + EXP(-r)*cmplx(c, -s, kind=DP))

  END FUNCTION qe_cosh

  FUNCTION qe_gauss(x) result(gauss)
    USE kinds,     ONLY : DP
    USE constants, ONLY : sqrtpm1  ! 1/sqrt(pi)
    IMPLICIT NONE

    REAL(DP), INTENT(in) :: x
    REAL(DP) :: gauss

    gauss = 2.0d0*sqrtpm1*EXP(-x*x)

  END FUNCTION qe_gauss

  FUNCTION exp_gauss( x, y )
    USE kinds,     ONLY : DP
    USE constants, ONLY : sqrtpm1 !1/sqrt(pi)

    REAL(DP), INTENT(IN) :: x, y
    REAL(DP) :: exp_gauss

    exp_gauss = 2._DP*sqrtpm1*EXP( x - y*y )

  END FUNCTION exp_gauss
END MODULE esm_stres_mod
