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