1
2!
3! Copyright (C) 2005-2014 Quantum ESPRESSO group
4! This file is distributed under the terms of the
5! GNU General Public License. See the file `License'
6! in the root directory of the present distribution,
7! or http://www.gnu.org/copyleft/gpl.txt .
8!
9!
10SUBROUTINE tpssmeta(nnr, nspin,grho,rho,kedtau,etxc)
11  !     ===================
12  !--------------------------------------------------------------------
13  use kinds,   only: dp
14  use funct,   only: get_meta
15  use xc_mgga, only: xc_metagcx, change_threshold_mgga !, &
16                     !tau_xc_array, tau_xc_array_spin
17  !
18  IMPLICIT NONE
19  !
20  ! input
21  integer :: nspin, nnr
22  real(dp) :: grho(3,nnr,nspin), rho(nnr,nspin), kedtau(nnr,nspin)
23  ! output: excrho: exc * rho ;  E_xc = \int excrho(r) d_r
24  ! output: rhor:   contains the exchange-correlation potential
25  real(dp) :: etxc
26  REAL(dp) :: zeta, rh, grh2
27  INTEGER :: k, ipol, is
28  REAL(dp), PARAMETER :: epsr = 1.0d-6, epsg = 1.0d-10
29  INTEGER :: imeta
30  !
31  etxc = 0.d0
32  ! calculate the gradient of rho+rho_core in real space
33  imeta = get_meta()
34  !
35  call exch_corr_meta() !HK/MCA
36  !
37  RETURN
38  !
39  !
40 contains
41  !
42  !
43  subroutine  exch_corr_meta()
44    !
45    implicit none
46    !
47    INTEGER :: np
48    REAL(DP), ALLOCATABLE :: sx(:), v1x(:,:), v2x(:,:), v3x(:,:), &
49                             sc(:), v1c(:,:), v2c(:,:,:), v3c(:,:)
50    !
51    np=1
52    if (nspin==2) np=3
53    !
54    allocate( sx(nnr), v1x(nnr,nspin), v2x(nnr,nspin),    v3x(nnr,nspin) )
55    allocate( sc(nnr), v1c(nnr,nspin), v2c(np,nnr,nspin), v3c(nnr,nspin) )
56    !
57    if (nspin==1) then
58      !
59      call change_threshold_mgga( epsr, epsg, epsr )
60      !
61      call xc_metagcx( nnr, 1, np, rho, grho, kedtau, sx, sc, &
62                       v1x, v2x, v3x, v1c, v2c, v3c )
63      !
64      rho(:,1) = v1x(:,1) + v1c(:,1)
65      kedtau(:,1) = (v3x(:,1) + v3c(:,1)) * 0.5d0
66      ! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
67      do ipol = 1, 3
68         grho(ipol,:,1) = (v2x(:,1) + v2c(1,:,1))*grho(ipol,:,1)
69      enddo
70      etxc = SUM( (sx(:) + sc(:)) * SIGN(1.d0,rho(:,1)) )
71      !
72    else
73      !
74      call change_threshold_mgga( epsr )
75      !
76      call xc_metagcx( nnr, 2, np, rho, grho, kedtau, sx, sc, &
77                       v1x, v2x, v3x, v1c, v2c, v3c )
78      !
79      rho(:,1) = v1x(:,1) + v1c(:,1)
80      rho(:,2) = v1x(:,2) + v1c(:,2)
81      !
82      ! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
83      !
84      do ipol = 1, 3
85        grho(ipol,:,1) = (v2x(:,1)*grho(ipol,:,1) + v2c(ipol,:,1))
86        grho(ipol,:,2) = (v2x(:,2)*grho(ipol,:,2) + v2c(ipol,:,2))
87      enddo
88      !
89      kedtau(:,1) = (v3x(:,1) + v3c(:,1)) * 0.5d0
90      kedtau(:,2) = (v3x(:,2) + v3c(:,2)) * 0.5d0
91      etxc = etxc + SUM( sx(:) + sc(:) )
92      !
93    endif
94    !
95    deallocate( sx, v1x, v2x, v3x )
96    deallocate( sc, v1c, v2c, v3c )
97    !
98    !
99    return
100    !
101  end subroutine exch_corr_meta
102  !
103END SUBROUTINE tpssmeta
104
105!-----------------------------------------------------------------------
106