Branch data Line data Source code
1 : : !!****m* m_xc_tb09/m_xc_tb09
2 : : !! NAME
3 : : !! m_xc_tb09
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains routine(s) related to exchange-correlation Tran-Blaha 2009
7 : : !! functional (modified Becke-Johnson)
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2023-2026 ABINIT group (MT)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : MODULE m_xc_tb09
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_xmpi
28 : : use m_errors
29 : : use m_time, only : timab
30 : : use defs_abitypes, only : MPI_type
31 : :
32 : : use m_xctk, only : xcden
33 : : use m_drivexc, only : mkdenpos
34 : : use m_geometry, only : metric
35 : :
36 : : use m_pawang, only : pawang_type
37 : : use m_pawrad, only : pawrad_type,pawrad_deducer0,nderiv_gen,simp_gen
38 : : use m_pawtab, only : pawtab_type
39 : : use m_pawrhoij, only : pawrhoij_type
40 : : use m_paw_denpot, only : pawdensities
41 : : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
42 : :
43 : : use libxc_functionals, only : libxc_functional_type,libxc_functionals_set_c_tb09, &
44 : : & libxc_functionals_is_tb09,libxc_functionals_ixc
45 : :
46 : : implicit none
47 : :
48 : : private
49 : :
50 : : !public procedure(s)
51 : : public :: xc_tb09_update_c ! Compute the C parameter for the TB09 functional (NCPP+PAW)
52 : :
53 : : !Private variables
54 : : real(dp),save :: grho_over_rho_pw ! Contribution to TB09 c from the plane-wave pseudo-density
55 : : real(dp),save :: grho_over_rho_paw ! Contribution to TB09 c from the on-site PAW densities
56 : :
57 : : CONTAINS !========================================================================================
58 : : !!***
59 : :
60 : : !----------------------------------------------------------------------
61 : :
62 : : !!****f* m_xc_tb09/xc_tb09_update_c
63 : : !! NAME
64 : : !! xc_tb09_update_c
65 : : !!
66 : : !! FUNCTION
67 : : !! This routine computes from the density the C parameter of the Tran-Blaha 2009
68 : : !! XC functional (modified Becke-Johnson)
69 : : !! C = 1/V Int[|Grad(Rho)|/Rho] + Sum_PAW_aug_regions[ Natom/V [|Grad(Rho)|/Rho] ]
70 : : !!
71 : : !! INPUTS
72 : : !! intxc = 1 if the XC functional has to be interpolated on a more refined mesh
73 : : !! ixc= choice of exchange-correlation scheme
74 : : !! mpi_enreg= information about MPI parallelization
75 : : !! natom= total number of atoms in cell
76 : : !! nfft= (effective) number of FFT grid points (for this processor)
77 : : !! ngfft(18)= contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
78 : : !! nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
79 : : !! nhatdim= -PAW only- 0 if nhat array is not used ; 1 otherwise
80 : : !! nhatgr(nfft,xcdata%nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
81 : : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
82 : : !! nspden= number of spin-density components
83 : : !! ntypat= number of types of atoms in unit cell.
84 : : !! n3xccc= dimension of the xccc3d array (0 or nfft or cplx*nfft).
85 : : !! pawang <type(pawang_type)> =paw angular mesh and related data
86 : : !! pawrad(ntypat) <type(pawrad_type)> =paw radial mesh and related data
87 : : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
88 : : !! pawtab(ntypat) <type(pawtab_type)>= paw tabulated starting data
89 : : !! pawxcdev= choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
90 : : !! rhor(nfft,nspden)= electron density in real space in electrons/bohr**3
91 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
92 : : !! usepaw = 1 for PAW computation, 0 else
93 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
94 : : !! xc_denpos= lowest allowed density
95 : : !!
96 : : !! ----- Optional arguments -----
97 : : !! [computation_type]= 'only_pw' : update only plane-wave contribution
98 : : !! 'only_paw': update only PAW on-site contribution
99 : : !! 'all' (default): compute all contributions
100 : : !! [mpi_atmtab(:)]= indexes of the atoms treated by current proc
101 : : !! [comm_atom]= MPI communicator over atoms
102 : : !! [xc_funcs(2)]= <type(libxc_functional_type)>=libxc XC functionals (if not the global one)
103 : : !!
104 : : !! OUTPUT
105 : : !! No output
106 : : !! The C parameter is directly set in the xc_functional datastructure
107 : : !!
108 : : !! SOURCE
109 : :
110 : 111422 : subroutine xc_tb09_update_c(intxc,ixc,mpi_enreg,natom,nfft,ngfft,nhat,nhatdim, &
111 : 111422 : & nhatgr,nhatgrdim,nspden,ntypat,n3xccc,pawang,pawrad, &
112 [ - + ]: 55711 : & pawrhoij,pawtab,pawxcdev,rhor,rprimd,usepaw,xccc3d,xc_denpos, &
113 [ + - + + : 55711 : & computation_type,mpi_atmtab,comm_atom,xc_funcs) ! Optional arguments
- + ]
114 : :
115 : : !Arguments ------------------------------------
116 : : !scalars
117 : : integer,intent(in) :: intxc,ixc,n3xccc,nfft,natom,nhatdim,nhatgrdim
118 : : integer,intent(in) :: nspden,ntypat,pawxcdev,usepaw
119 : : integer,intent(in),optional :: comm_atom
120 : : real(dp),intent(in),optional :: xc_denpos
121 : : character(len=*),intent(in),optional :: computation_type
122 : : type(pawang_type),intent(in) :: pawang
123 : : type(MPI_type),intent(in) :: mpi_enreg
124 : : !arrays
125 : : integer,intent(in) :: ngfft(18)
126 : : integer,intent(in),optional,target :: mpi_atmtab(:)
127 : : real(dp),intent(in) :: nhat(nfft,nspden*nhatdim),nhatgr(nfft,nspden,3*nhatgrdim)
128 : : real(dp),intent(in) :: rhor(nfft,nspden)
129 : : real(dp),intent(in) :: rprimd(3,3),xccc3d(n3xccc)
130 : : type(libxc_functional_type),intent(inout),optional :: xc_funcs(2)
131 : : type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw)
132 : : type(pawrhoij_type),intent(in) :: pawrhoij(:)
133 : : type(pawtab_type),intent(in),target :: pawtab(ntypat*usepaw)
134 : :
135 : : !Local variables-------------------------------
136 : : !scalars
137 : : integer :: cplex,iatom,iatom_tot,ierr,ii,ilm,iloop,ir,itypat,ixc_from_lib,ipts
138 : : integer :: ishift,iwarn,lm_size,lm_size_eff,mesh_size,my_comm_atom,my_natom,ngrad
139 : : integer :: nfftot,npts,nspden1,nzlmopt,option,opt_compch,opt_dens,pawprtvol,usecore
140 : : integer :: usenhat,usexcnhat
141 : : real(dp) :: compch_dum,factor,grho_over_rho,grho2,rho,sumg,ucvol,xc_c_tb09
142 : : logical :: my_atmtab_allocated,paral_atom,test_tb09,use_exact_nhat_gradient
143 : : logical :: compute_pw,compute_paw
144 : : character(len=500) :: msg
145 : : !arrays
146 : 55711 : integer,pointer :: my_atmtab(:)
147 : 55711 : logical,allocatable :: lmselect(:),lmselect_tmp(:)
148 : : real(dp) :: gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3)
149 : : !real(dp) :: tsec(2)
150 : 55711 : real(dp),allocatable :: dylmdr(:,:,:),ff(:)
151 : 55711 : real(dp),allocatable :: rhotot(:),grhotot(:,:),rhonow(:,:,:)
152 : 55711 : real(dp),allocatable :: rhoxc(:,:),drhoxc(:),dcorexc(:)
153 : 55711 : real(dp),allocatable :: rho1(:,:,:),trho1(:,:,:),nhat1(:,:,:)
154 : 55711 : real(dp), contiguous, pointer :: corexc(:)
155 : :
156 : : ! *************************************************************************
157 : :
158 : : !call timab(xx,1,tsec)
159 : :
160 : : !Type of calculation
161 : 55711 : compute_pw=.true. ; compute_paw=.true.
162 : 55711 : if (present(computation_type)) then
163 [ - + ]: 54990 : if (computation_type=='only_pw') compute_paw=.false.
164 [ - + ]: 54990 : if (computation_type=='only_paw') compute_pw=.false.
165 : : end if
166 [ # # ]: 0 : if ((.not.compute_pw).and.(.not.compute_paw)) return
167 : :
168 : : !Nothing to do if XC functional is not TB09
169 [ - + ]: 55711 : if (present(xc_funcs)) then
170 : 0 : test_tb09=libxc_functionals_is_tb09(xc_functionals=xc_funcs)
171 : : else
172 : 55711 : test_tb09=libxc_functionals_is_tb09()
173 : : end if
174 [ + + ]: 55711 : if (.not.test_tb09) return
175 : :
176 : : !Check options
177 [ + - ]: 71 : if(ixc<0) then
178 [ - + ]: 71 : if (present(xc_funcs)) then
179 : 0 : ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs)
180 : : else
181 : 71 : ixc_from_lib=libxc_functionals_ixc()
182 : : end if
183 [ - + ]: 71 : if (ixc /= ixc_from_lib) then
184 : : write(msg, '(a,i0,a,a,i0)')&
185 : 0 : & 'The value of ixc specified in input, ixc = ',ixc,ch10,&
186 : 0 : & 'differs from the one used to initialize the functional ',ixc_from_lib
187 : 0 : ABI_BUG(msg)
188 : : end if
189 : : end if
190 [ + + + - ]: 71 : if (usepaw==1.and.compute_paw) then
191 [ - + ]: 31 : if (pawxcdev/=0) then
192 : 0 : msg='TB09 (mBJ) XC functional only available with pawxcdev=0!'
193 : 0 : ABI_BUG(msg)
194 : : end if
195 [ - + ]: 31 : if(pawang%angl_size==0) then
196 : 0 : msg='pawang%angl_size=0!'
197 : 0 : ABI_BUG(msg)
198 : : end if
199 [ - + ]: 31 : if(.not.allocated(pawang%ylmr)) then
200 : 0 : msg='pawang%ylmr must be allocated!'
201 : 0 : ABI_BUG(msg)
202 : : end if
203 [ - + ]: 31 : if(.not.allocated(pawang%ylmrgr)) then
204 : 0 : msg='pawang%ylmrgr must be allocated!'
205 : 0 : ABI_BUG(msg)
206 : : end if
207 [ - + ]: 31 : if (size(pawrhoij)/=mpi_enreg%my_natom) then
208 : 0 : msg='Invalid pawrhoij size!'
209 : 0 : ABI_BUG(msg)
210 : : end if
211 : : end if
212 : :
213 : : !Initialize cell data
214 : 71 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
215 : :
216 : : !Other initializations
217 [ + + + + ]: 102 : usexcnhat=0; if (usepaw==1) usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
218 : :
219 : : !Initialization of Int[|Grad(Rho)|/Rho] contributions (plane-wave and on-site if PAW)
220 [ + - ]: 71 : if (compute_pw) grho_over_rho_pw=zero
221 [ + - ]: 71 : if (compute_paw) grho_over_rho_paw=zero
222 : :
223 : : ! ========================================================================
224 : : ! =========== Plane wave contribution to Int[|Grad(Rho)|/Rho] ============
225 : :
226 [ + - ]: 71 : if (compute_pw) then
227 : :
228 : : ! Total density used in XC (with/without core density and compensation density)
229 [ - + + - : 213 : ABI_MALLOC(rhotot,(nfft))
+ - + - -
+ - + ]
230 [ + + ]: 1394031 : rhotot(1:nfft)=rhor(1:nfft,1)
231 [ + - + + : 428615 : if (usexcnhat==0.and.nhatdim==1) rhotot(1:nfft)=rhotot(1:nfft)-nhat(1:nfft,1)
+ + ]
232 [ + + + + ]: 428615 : if (n3xccc>0) rhotot(1:nfft)=rhotot(1:nfft)+xccc3d(1:nfft)
233 : :
234 : : ! When possible, it is better to use the exact (analytical) gradient of nhat
235 : : ! In that case, we substract nhat from the density now
236 : : ! If not possible, the gradient will be computed in Fourier space (see call to xcden)
237 [ + + - + : 71 : use_exact_nhat_gradient=(nhatdim==1.and.nhatgrdim==1.and.usexcnhat==1.and.intxc==0)
- - - - ]
238 [ # # ]: 0 : if (use_exact_nhat_gradient) rhotot(1:nfft)=rhotot(1:nfft)-nhat(1:nfft,1)
239 : :
240 : : ! Loop on unshifted or shifted grids
241 [ + + ]: 142 : do ishift=0,intxc
242 : :
243 : : ! Set up density on unshifted or shifted grid (will be in rhonow(:,:,1)),
244 : : ! as well as the gradient of the density, also on the unshifted
245 : : ! or shifted grid (will be in rhonow(:,:,2:4)), if needed.
246 [ - + + - : 213 : ABI_MALLOC(rhonow,(nfft,1,4))
+ - + - -
+ - + ]
247 : 71 : nspden1=1 ; cplex=1 ; ngrad=2 ; qphon(:)=zero
248 : 71 : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden1,qphon,rhotot,rhonow)
249 : :
250 : : ! PAW: if they were previously substracted, re-add nhat and its "exact" gradient
251 [ - + ]: 71 : if (use_exact_nhat_gradient) then
252 [ # # ]: 0 : rhonow(1:nfft,1,1)=rhonow(1:nfft,1,1)+nhat(1:nfft,1)
253 [ # # ]: 0 : do ii=1,3
254 [ # # ]: 0 : rhonow(1:nfft,1,ii+1)=rhonow(1:nfft,1,ii+1)+nhatgr(1:nfft,1,ii)
255 : : end do
256 : : end if
257 : :
258 : : ! Make the density positive everywhere (but do not care about gradients)
259 : 71 : nspden1=1 ; iwarn=1 ; option=1
260 : 71 : call mkdenpos(iwarn,nfft,nspden1,option,rhonow(:,1,1),xc_denpos)
261 : :
262 : : ! Accumulate integral of |Grad_rho|/Rho (to be used for TB09 XC)
263 [ + + ]: 1394031 : do ipts=1,nfft
264 : 1393960 : rho=rhonow(ipts,1,1)
265 [ + - ]: 1394031 : if (abs(rho)>tol10) then
266 : 1393960 : grho2=rhonow(ipts,1,2)**2+rhonow(ipts,1,3)**2+rhonow(ipts,1,4)**2
267 : 1393960 : grho_over_rho_pw=grho_over_rho_pw+sqrt(grho2)/rho
268 : : end if
269 : : end do
270 : :
271 [ - + ]: 142 : ABI_FREE(rhonow)
272 : :
273 : : ! End loop on unshifted or shifted grids
274 : : end do
275 : :
276 : : ! Normalize the integral
277 : 71 : nfftot=ngfft(1)*ngfft(2)*ngfft(3)
278 : 71 : grho_over_rho_pw=grho_over_rho_pw*ucvol/dble(nfftot)/dble(intxc+1)
279 : :
280 : : ! Reduction in case of FFT distribution
281 [ - + ]: 71 : if (mpi_enreg%nproc_fft>1) then
282 : 0 : call xmpi_sum(grho_over_rho_pw,mpi_enreg%comm_fft,ierr)
283 : : end if
284 : :
285 [ - + ]: 71 : ABI_FREE(rhotot)
286 : :
287 : : end if ! compute_pw
288 : :
289 : : ! ========================================================================
290 : : ! ========== PAW on-site contributions to Int[|Grad(Rho)|/Rho] ===========
291 : :
292 [ + + + - ]: 71 : if (usepaw==1.and.compute_paw) then
293 : :
294 : : ! Set up parallelism over atoms
295 : 31 : my_natom=mpi_enreg%my_natom
296 [ + - + - ]: 31 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
297 [ - + ]: 31 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
298 [ + - ]: 31 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
299 : 31 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
300 : :
301 : 31 : usecore=1 ; nzlmopt=-1
302 : 31 : npts=pawang%angl_size
303 : :
304 : : ! Compute Ylm gradients
305 : : ! dYlm/dr_i = { dYlm/dr_i^hat - r_i^hat * Sum_j[dYlm/dr_j^hat r_j^hat] } * (1/r)
306 [ - + + - : 124 : ABI_MALLOC(dylmdr,(3,npts,pawang%ylm_size))
- + + - +
- + - - +
- + ]
307 [ + + ]: 527 : do ilm=1,pawang%ylm_size
308 [ + + ]: 77903 : do ipts=1,npts
309 [ + + ]: 309504 : factor=sum(pawang%ylmrgr(1:3,ilm,ipts)*pawang%anginit(1:3,ipts))
310 [ + + ]: 310000 : dylmdr(1:3,ipts,ilm)=pawang%ylmrgr(1:3,ilm,ipts)-factor*pawang%anginit(1:3,ipts)
311 : : end do
312 : : end do
313 : :
314 : : ! Loop on atom sites
315 [ + + ]: 93 : do iatom=1,my_natom
316 [ - + ]: 62 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
317 : :
318 : 62 : itypat=pawrhoij(iatom)%itypat
319 : 62 : mesh_size=pawtab(itypat)%mesh_size
320 : 62 : lm_size=pawtab(itypat)%lcut_size**2
321 : 62 : lm_size_eff=min(lm_size,pawang%ylm_size)
322 [ - + + - : 186 : ABI_MALLOC(ff,(mesh_size))
+ - + - -
+ - + -
+ ]
323 : :
324 : : ! Compute on-site densities
325 [ - + + - : 248 : ABI_MALLOC(rho1 ,(mesh_size,lm_size,1))
- + + - +
- + - - +
- + ]
326 [ - + + - : 186 : ABI_MALLOC(trho1,(mesh_size,lm_size,1))
- + + - -
+ - + ]
327 [ - + + - : 310 : ABI_MALLOC(nhat1,(mesh_size,lm_size,usexcnhat))
- + + - +
- - - + -
- + - + -
+ ]
328 [ - + + - : 186 : ABI_MALLOC(lmselect,(lm_size))
+ - + - -
+ - + ]
329 [ - + + - : 124 : ABI_MALLOC(lmselect_tmp,(lm_size))
- + - + ]
330 [ + + ]: 620 : lmselect_tmp(:)=.true.
331 : 62 : cplex=1 ; nspden1=1 ; opt_compch=0 ; opt_dens=1-usexcnhat ; pawprtvol=0
332 : : call pawdensities(compch_dum,cplex,iatom_tot,lmselect_tmp,lmselect,&
333 : : & lm_size,nhat1,nspden1,nzlmopt,opt_compch,opt_dens,-1,0,pawang,pawprtvol,&
334 : 62 : & pawrad(itypat),pawrhoij(iatom),pawtab(itypat),rho1,trho1)
335 [ - + ]: 62 : ABI_FREE(lmselect_tmp)
336 : :
337 : : ! Allocation of temporary densities
338 [ - + + - : 124 : ABI_MALLOC(rhotot,(mesh_size))
- + - + ]
339 [ - + + - : 186 : ABI_MALLOC(grhotot,(mesh_size,3))
+ - + - -
+ - + ]
340 [ - + + - : 186 : ABI_MALLOC(rhoxc,(mesh_size,lm_size))
- + + - -
+ - + ]
341 [ - + + - : 124 : ABI_MALLOC(drhoxc,(mesh_size))
- + - + -
+ ]
342 [ - + + - : 124 : ABI_MALLOC(dcorexc,(mesh_size))
- + - + -
+ ]
343 : :
344 : : ! Loop over AE and PS contributions
345 [ + + ]: 186 : do iloop=1,2
346 [ + + ]: 124 : if (iloop==1) then
347 [ + + + + ]: 776240 : rhoxc(:,:)=rho1(:,:,1)
348 : 62 : corexc => pawtab(itypat)%coredens(:)
349 : : usenhat=0
350 : : else
351 [ + + + + ]: 776240 : rhoxc(:,:)=trho1(:,:,1)
352 : 62 : corexc => pawtab(itypat)%tcoredens(:,1)
353 : : usenhat=usexcnhat
354 : : end if
355 : :
356 : : ! Add hat density if needed
357 [ - + - - : 62 : if (usenhat>0) rhoxc(:,:)=rhoxc(:,:)+nhat1(:,:,1)
- - ]
358 : :
359 : : ! Need derivative of core density
360 : : if (usecore==1) then
361 : 124 : call nderiv_gen(dcorexc,corexc,pawrad(itypat))
362 : : end if
363 : :
364 : : ! Do loop on the angular part (theta,phi)
365 [ + + ]: 19530 : do ipts=1,npts
366 : :
367 : : ! Compute the density and its gradient for this (theta,phi)
368 [ + + + + : 107630016 : rhotot(:)=zero ; grhotot(:,:)=zero
+ + ]
369 [ + + ]: 193440 : do ilm=1,lm_size_eff
370 [ + + ]: 193440 : if (lmselect(ilm)) then
371 : : !Density
372 [ + + ]: 65966784 : rhotot(:)=rhotot(:)+rhoxc(:,ilm)*pawang%ylmr(ilm,ipts)
373 : : !Gradient
374 [ + + ]: 65966784 : ff(1:mesh_size)=rhoxc(1:mesh_size,ilm)
375 : 47424 : call nderiv_gen(drhoxc,ff,pawrad(itypat))
376 [ + + ]: 65919360 : ff(2:mesh_size)=ff(2:mesh_size)/pawrad(itypat)%rad(2:mesh_size)
377 : 47424 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
378 [ + + ]: 189696 : do ii=1,3
379 : : grhotot(1:mesh_size,ii)=grhotot(1:mesh_size,ii) &
380 : : & +drhoxc(1:mesh_size)*pawang%ylmr(ilm,ipts)*pawang%anginit(ii,ipts) &
381 [ + + ]: 197947776 : & +ff(1:mesh_size)*dylmdr(ii,ipts,ilm)
382 : : end do
383 : : end if
384 : : end do
385 : : if (usecore==1) then
386 [ + + ]: 26907504 : rhotot(:)=rhotot(:)+corexc(:)
387 [ + + ]: 77376 : do ii=1,3
388 [ + + ]: 80741856 : grhotot(:,ii)=grhotot(:,ii)+dcorexc(:)*pawang%anginit(ii,ipts)
389 : : end do
390 : : end if
391 : :
392 : : ! Make the density positive everywhere (but do not care about gradients)
393 : 19344 : nspden1=1 ; iwarn=1
394 : 19344 : call mkdenpos(iwarn,mesh_size,nspden1,0,rhotot,xc_denpos)
395 : :
396 : : ! Accumulate integral of |Grad_rho|/Rho
397 [ + + ]: 26907504 : do ir=1,mesh_size
398 : 26888160 : rho=rhotot(ir)
399 [ + - ]: 26907504 : if (abs(rho)>tol10) then
400 : 26888160 : grho2=grhotot(ir,1)**2+grhotot(ir,2)**2+grhotot(ir,3)**2
401 : 26888160 : ff(ir)=sqrt(grho2)/rho
402 : : end if
403 : : end do
404 [ + + ]: 26907504 : ff(1:mesh_size)=ff(1:mesh_size)*pawrad(itypat)%rad(1:mesh_size)**2
405 : 19344 : call simp_gen(sumg,ff,pawrad(itypat))
406 : :
407 : : ! Add the contribution of the atom
408 [ + + ]: 19468 : if (iloop==1) then
409 : 9672 : grho_over_rho_paw=grho_over_rho_paw+sumg*four_pi*pawang%angwgth(ipts)
410 : : else
411 : 9672 : grho_over_rho_paw=grho_over_rho_paw-sumg*four_pi*pawang%angwgth(ipts)
412 : : end if
413 : :
414 : : end do ! npts
415 : : end do ! AE/PS loop
416 : :
417 : : ! Deallocate temporary memory
418 [ - + ]: 62 : ABI_FREE(ff)
419 [ - + ]: 62 : ABI_FREE(rho1)
420 [ - + ]: 62 : ABI_FREE(trho1)
421 [ - + ]: 62 : ABI_FREE(nhat1)
422 [ - + ]: 62 : ABI_FREE(lmselect)
423 [ - + ]: 62 : ABI_FREE(rhotot)
424 [ - + ]: 62 : ABI_FREE(grhotot)
425 [ - + ]: 62 : ABI_FREE(rhoxc)
426 [ - + ]: 62 : ABI_FREE(drhoxc)
427 [ - + ]: 93 : ABI_FREE(dcorexc)
428 : :
429 : : end do ! atom sites
430 : :
431 : : ! Deallocate temporary memory
432 [ - + ]: 31 : ABI_FREE(dylmdr)
433 : :
434 : : ! Reduction in case of parallelism
435 [ - + ]: 31 : if (paral_atom) then
436 : 0 : call xmpi_sum(grho_over_rho_paw,my_comm_atom,ierr)
437 : : end if
438 : :
439 : : ! Destroy atom table used for parallelism
440 : 31 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
441 : :
442 : : end if ! PAW and compute_paw
443 : :
444 : : ! ========================================================================
445 : : ! ============ Assemble contributions and update mBJ C value ============
446 : :
447 : 71 : grho_over_rho = (grho_over_rho_pw + grho_over_rho_paw)/ucvol
448 : 71 : xc_c_tb09 = -0.012_dp + 1.023_dp * sqrt(grho_over_rho)
449 [ + + ]: 71 : if (xc_c_tb09<one) xc_c_tb09=one
450 [ - + ]: 71 : if (xc_c_tb09>two) xc_c_tb09=two
451 : :
452 : 71 : write(msg,'(2a,f10.5,a)') ch10, &
453 : 142 : & ' In the mGGA functional TB09, c (computed value) = ',xc_c_tb09,ch10
454 : 71 : call wrtout(std_out,msg,'COLL')
455 : :
456 [ - + ]: 71 : if (present(xc_funcs)) then
457 : 0 : call libxc_functionals_set_c_tb09(xc_c_tb09,xc_functionals=xc_funcs)
458 : : else
459 : 71 : call libxc_functionals_set_c_tb09(xc_c_tb09)
460 : : end if
461 : :
462 : : !call timab(xx,2,tsec)
463 : :
464 [ + + - + : 111422 : end subroutine xc_tb09_update_c
- + - + ]
465 : : !!***
466 : :
467 : : !----------------------------------------------------------------------
468 : :
469 : : END MODULE m_xc_tb09
470 : : !!***
|
