1!
2! Copyright (C) 2001-2007 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!-----------------------------------------------------------------------
10SUBROUTINE set_rhoc
11  !-----------------------------------------------------------------------
12  !
13  !    This routine computes the core charge on the real space 3D mesh
14  !    (rho_core) and in reciprocal space (rhog_core)
15  !
16  !
17  USE kinds,     ONLY : dp
18  USE io_global, ONLY : stdout
19  USE atom,      ONLY : msh, rgrid
20  USE uspp_param,ONLY : upf
21  USE ions_base, ONLY : ntyp => nsp
22  USE cell_base, ONLY : omega, tpiba2
23  USE fft_base,  ONLY : dfftp
24  USE fft_rho,   ONLY : rho_g2r
25  USE gvect,     ONLY : ngm, ngl, gl, igtongl
26  USE vlocal,    ONLY : strf
27  USE mp_bands,  ONLY : intra_bgrp_comm
28  USE mp,        ONLY : mp_sum
29  USE scf,       ONLY : rho_core, rhog_core
30  !
31  IMPLICIT NONE
32  !
33  REAL(DP) , ALLOCATABLE ::  rhocg(:)
34  ! the radial fourier transform
35  REAL(DP) ::  rhoneg
36  ! used to check the core charge
37  INTEGER :: ir, nt, ng
38  ! counter on mesh points
39  ! counter on atomic types
40  ! counter on g vectors
41
42  rhog_core(:) = 0.0_DP
43  rho_core(:)  = 0.0_DP
44
45  IF ( ANY( upf(1:ntyp)%nlcc ) ) THEN
46
47     ALLOCATE (rhocg( ngl))
48     !
49     !    the sum is on atom types
50     !
51     DO nt = 1, ntyp
52        IF ( upf(nt)%nlcc ) THEN
53           !
54           ! drhoc computes the radial fourier transform for each shell of g vec
55           !
56           CALL drhoc (ngl, gl, omega, tpiba2, msh (nt), rgrid(nt)%r, &
57             rgrid(nt)%rab, upf(nt)%rho_atc, rhocg)
58           !
59           !     multiply by the structure factor and sum
60           !
61           DO ng = 1, ngm
62              rhog_core(ng) = rhog_core(ng) + strf(ng,nt) * rhocg(igtongl(ng))
63           END DO
64       ENDIF
65     ENDDO
66     DEALLOCATE (rhocg)
67     !
68     CALL rho_g2r( dfftp, rhog_core, rho_core )
69     !
70     !    test on the charge and computation of the core energy
71     !
72     rhoneg = 0.d0
73     DO ir = 1, dfftp%nnr
74        rhoneg = rhoneg + min (0.d0, rho_core (ir) )
75        !
76        ! NOTE: Core charge is computed in reciprocal space and brought to real
77        ! space by FFT. For non smooth core charges (or insufficient cut-off)
78        ! this may result in negative values in some grid points.
79        ! Up to October 1999 the core charge was forced to be positive definite.
80        ! This induces an error in force, and probably stress, calculation if
81        ! the number of grid points where the core charge would be negative is
82        ! large. The error disappears for sufficiently high cut-off, but may be
83        ! rather large and it is better to leave the core charge as it is.
84        ! If you insist to have it positive definite (with the possible problems
85        ! mentioned above) uncomment the following lines.  SdG, Oct 15 1999
86        !
87        !         rho_core(ir) = MAX (0.0_dp, rho_core(ir))
88     ENDDO
89     !
90     rhoneg = rhoneg / (dfftp%nr1 * dfftp%nr2 * dfftp%nr3)
91     CALL mp_sum(  rhoneg, intra_bgrp_comm )
92     !
93     IF (rhoneg < -1.0d-6) &
94        WRITE(stdout,'(/5x,"Check: negative core charge=",2f12.6)') rhoneg
95     !
96     ! calculate core_only exch-corr energy etxcc=E_xc[rho_core] if required
97     ! The term was present in previous versions of the code but it shouldn't
98     !
99     !   call create_scf_type(dum)
100     !   dum%of_r(:,:) = 0.0_DP
101     !   dum%of_g(:,:) = (0.0_DP, 0.0_DP)
102     !
103     !   call v_xc( dum, rho_core, rhog_core, etxcc, vtxcc, aux )
104     !
105     !   call destroy_scf_type(dum)
106     !   WRITE( stdout, 9000) etxcc
107     ! 9000 format (5x,'core-only xc energy         = ',f15.8,' Ry')
108     !   WRITE( stdout,  * ) 'BEWARE it will be subtracted from total energy !'
109     !
110  END IF
111  !
112  RETURN
113
114END SUBROUTINE set_rhoc
115
116