Line data Source code
1 : !!****m* ABINIT/m_paw_init
2 : !! NAME
3 : !! m_paw_init
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines related tp PAW calculations initialization.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2018-2026 ABINIT group (FJ, MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_paw_init
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_splines
28 : use m_dtset
29 :
30 : use m_time, only : timab
31 : use m_pawpsp, only : pawpsp_nl
32 : use m_paw_atom,only : atompaw_shpfun
33 : use m_pawang, only : pawang_type, pawang_init, pawang_free
34 : use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen, poisson, pawrad_deducer0
35 : use m_pawtab, only : pawtab_type
36 : use m_pawxc, only : pawxc_get_usekden,pawxc_get_uselaplacian,pawxc_get_xclevel
37 : use m_paw_numeric, only : paw_derfc
38 :
39 : implicit none
40 :
41 : private
42 :
43 : !public procedures.
44 : public :: pawinit ! Initialize some tabulated data for PAW calculations
45 : public :: paw_gencond ! Test whether we have to call pawinit to regenerate tabulated data.
46 :
47 : CONTAINS !========================================================================================
48 : !!***
49 :
50 : !----------------------------------------------------------------------
51 :
52 : !!****f* m_paw_init/pawinit
53 : !! NAME
54 : !! pawinit
55 : !!
56 : !! FUNCTION
57 : !! Initialize some starting values of several arrays used in PAW calculations.
58 : !!
59 : !! 1-Initialize data related to angular mesh
60 : !! 2-Tabulate normalized shape function g(r)
61 : !! 3-Compute indklmn indexes giving some l,m,n,lp,mp,np info
62 : !! from klmn=[(l,m,n),(lp,mp,np)]
63 : !! 4-Compute various factors/sizes (depending on (l,m,n))
64 : !! 5-Compute $q_ijL=\displaystyle
65 : !! \int_{0}^{r_c}{(\phi_i\phi_j-\widetilde{\phi_i}\widetilde{\phi_j}) r^l\,dr}
66 : !! Gaunt(l_i m_i,l_j m_j,l m))$
67 : !! $S_ij=\displaystyle \sqrt{4 \pi} q_ij0$
68 : !! 6-Compute $e_ijkl= vh1_ijkl - Vhatijkl - Bijkl - Cijkl$
69 : !! With:
70 : !! $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$
71 : !! $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$
72 : !! $B_ijkl =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$
73 : !! $C_ijkl =\sum_{L,m} {intvhatL*q_ijL*q_klL}$
74 : !! and:
75 : !! vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]]
76 : !! 7-Compute Ex-correlation energy for the core density
77 : !!
78 : !! COPYRIGHT
79 : !! Copyright (C) 1998-2026 ABINIT group (FJ, MT)
80 : !! This file is distributed under the terms of the
81 : !! GNU General Public License, see ~abinit/COPYING
82 : !! or http://www.gnu.org/copyleft/gpl.txt .
83 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
84 : !!
85 : !! INPUTS
86 : !! effmass_free=effective mass for electrons (1. in common case)
87 : !! gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated
88 : !! also determine the size of these pointers
89 : !! gsqcut_shp=effective cut-off to determine shape functions in reciprocal space
90 : !! hyb_range_fock=range coefficient for screened hybrid XC functionals
91 : !! ixc=choice of exchange-correlation functional
92 : !! lcutdens=max. l for densities/potentials moments computations
93 : !! lmix=max. l for which spherical terms will be mixed durinf SCF cycle
94 : !! mpsang=1+maximum angular momentum
95 : !! nphi="phi" dimension of paw angular mesh
96 : !! nsym=Number of symmetry elements in space group
97 : !! ntheta="theta" dimension of paw angular mesh
98 : !! pawang <type(pawang_type)>=paw angular mesh and related data
99 : !! %lmax=Maximum value of angular momentum l+1
100 : !! %gntselect((2*l_max-1)**2,l_max**2,l_max**2)=
101 : !! selection rules for Gaunt coefficients
102 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
103 : !! %mesh_size=Dimension of radial mesh
104 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
105 : !! %radfact(mesh_size)=Factor used to compute radial integrals
106 : !! pawspnorb=flag: 1 if spin-orbit coupling is activated
107 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
108 : !! %basis_size=Number of elements for the PAW nl basis
109 : !! %l_size=Maximum value of l+1 leading to non zero Gaunt coeffs
110 : !! %lmn_size=Number of (l,m,n) elements for the PAW basis
111 : !! %lmn2_size=lmn_size*(lmn_size+1)/2
112 : !! %dltij(lmn2_size)=factors used to compute sums over (ilmn,jlmn)
113 : !! %phi(mesh_size,basis_size)=PAW all electron wavefunctions
114 : !! %rshp=shape function radius (radius for compensation charge)
115 : !! %shape_type=Radial shape function type
116 : !! %shape_alpha=Alpha parameters in Bessel shape function
117 : !! %shape_lambda=Lambda parameter in gaussian shape function
118 : !! %shape_q=Q parameters in Bessel shape function
119 : !! %shape_sigma=Sigma parameter in gaussian shape function
120 : !! %tphi(mesh_size,basis_size)=PAW atomic pseudowavefunctions
121 : !! pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1=dev. on moments)
122 : !! usekden= 1 is kinetic energy density has to be computed, 0 otherwise
123 : !!
124 : !! OUTPUT
125 : !! pawang
126 : !! %gntselect(l_size_max**2,l_max**2*(l_max**2+1)/2)=selection rules for Gaunt coefficients
127 : !! %l_max=maximum value of angular momentum l+1
128 : !! %l_size_max=maximum value of angular momentum l_size=2*l_max-1
129 : !! %nsym=number of symmetry elements in space group
130 : !! %ngnt=number of non-zero Gaunt coefficients
131 : !! %realgnt(pawang%ngnt)=non-zero real Gaunt coefficients
132 : !! === only if pawxcdev==1 ==
133 : !! %anginit(3,angl_size)=for each point of the angular mesh, gives the coordinates
134 : !! of the corresponding point on an unitary sphere
135 : !! %angl_size=dimension of paw angular mesh (angl_size=ntheta*nphi)
136 : !! %angwgth(angl_size)=for each point of the angular mesh, gives the weight
137 : !! of the corresponding point on an unitary sphere
138 : !! %ntheta, nphi=dimensions of paw angular mesh
139 : !! %ylmr(l_size_max**2,angl_size)=real Ylm calculated in real space
140 : !! %ylmrgr(1:3,l_size_max**2,angl_size)=first gradients of real Ylm calculated in real space
141 : !! === only if pawspnorb==1 ==
142 : !! %ls_ylm(2,l_max**2,l_max**2,4)=LS operator in the real spherical harmonics basis
143 : !! %use_ls_ylm=flag activated if ls_ylm is allocated
144 : !! %usespnorb=flag activated for spin-orbit coupling
145 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated data read at start:
146 : !! %lcut_size_=max. value of l+1 leading to non zero Gaunt coeffs modified by lcutdens
147 : !! %lmnmix_sz=number of (lmn,lmn_prime) verifying l<=lmix and l_prime<=lmix
148 : !! %mqgrid_shp=number of points in reciprocal space for shape function
149 : !! %indklmn(8,lmn2_size)=array giving klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn for each klmn=(ilmn,jlmn)
150 : !! %dshpfunc(mesh_size,l_size,4)=derivatives of shape function (used only for numerical shape functions)
151 : !! %eijkl(lmn2_size,lmn2_size)=part of the Dij that depends only from the projected occupation coeffs
152 : !! %exccore=Exchange-correlation energy for the core density
153 : !! %gnorm(l_size)=normalization factor of radial shape function
154 : !! %phiphj(:,:)=useful product Phi(:,i)*Phi(:,j)
155 : !! %qgrid_shp(mqgrid_shp)=points in reciprocal space for shape function
156 : !! %qijl(l_size**2,lmn2_size)=moments of the difference charge density between AE and PS partial wave
157 : !! %rad_for_spline(mesh_size)=radial grid used for spline (copy of pawrad%rad)
158 : !! %shapefunc(mesh_size,l_size)=normalized radial shape function
159 : !! %shapefncg(mqgrid_shp,l_size)=normalized radial shape function in reciprocal space
160 : !! %sij(lmn2_size)=nonlocal part of the overlap operator
161 : !! %tphitphj(:,:)=useful product tPhi(:,i)*tPhi(:,j)
162 : !!
163 : !! SOURCE
164 :
165 596 : subroutine pawinit(effmass_free,gnt_option,gsqcut_eff,hyb_range_fock,lcutdens,lmix,mpsang,&
166 596 : & nphi,nsym,ntheta,pawang,pawrad,pawspnorb,pawtab,pawxcdev,ixc,usepotzero,&
167 : & rcpaw_update)
168 :
169 : !Arguments ---------------------------------------------
170 : !scalars
171 : integer,intent(in) :: gnt_option,ixc,lcutdens,lmix,mpsang,nphi,nsym,ntheta
172 : integer,intent(in) :: pawspnorb,pawxcdev,usepotzero
173 : logical,intent(in),optional :: rcpaw_update
174 : real(dp),intent(in) :: effmass_free,gsqcut_eff,hyb_range_fock
175 : type(pawang_type),intent(inout) :: pawang
176 : !arrays
177 : type(pawrad_type),intent(in) :: pawrad(:)
178 : type(pawtab_type),target,intent(inout) :: pawtab(:)
179 :
180 : !Local variables ------------------------------
181 : !scalars
182 : integer,parameter :: mqgrid_shp_default=300
183 : integer :: basis_size,i0lm,i0ln,ij_size,il,ilm,ilmn,iln,iloop,iq,isel,isel1
184 : integer :: itypat,j0lm,j0lmn,j0ln,jl,jlm,jlmn,jln,klm,klm1
185 : integer :: klmn,klmn1,kln,kln1,l_size,ll,lm0,lmax,lmax1,lmin,lmin1,lmn2_size
186 : integer :: lmn_size,lmnmix,mesh_size,meshsz,mm,nabgnt_option,ngrad2_ylm,ntypat,pw_mesh_size
187 : integer :: usexcnhat,use_angular_grid,use_ls_ylm,use_ylm,usekden,k1min
188 : logical :: rcpaw_update_
189 : real(dp) :: dq,gnrm,intg,ql,ql1,rg,rg1,vh1,yp1,ypn
190 : character(len=500) :: message
191 : !arrays
192 596 : integer,allocatable :: indl(:,:),klm_diag(:),kmix_tmp(:)
193 596 : integer, ABI_CONTIGUOUS pointer :: indlmn(:,:)
194 : real(dp) :: tsec(2)
195 596 : real(dp),allocatable :: der(:),ff(:),gg(:),hh(:),indklmn_(:,:),intvhatl(:)
196 596 : real(dp),allocatable :: rad(:),rgl(:,:),vhatijl(:,:),vhatl(:),work(:)
197 596 : real(dp),pointer :: eijkl(:,:)
198 596 : real(dp), allocatable :: old_phiphj(:,:),old_tphitphj(:,:),old_qijl(:,:),old_vhatijl(:,:)
199 :
200 : !************************************************************************
201 :
202 : DBG_ENTER("COLL")
203 :
204 596 : call timab(553,1,tsec)
205 :
206 596 : ntypat=size(pawtab)
207 596 : if (size(pawrad)/=ntypat) then
208 0 : ABI_BUG('pawrad and pawtab should have the same size!')
209 : end if
210 :
211 596 : rcpaw_update_=.false.
212 596 : if(present(rcpaw_update)) rcpaw_update_=rcpaw_update
213 :
214 1464 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
215 596 : if(.not.rcpaw_update_) then
216 : !Immediately set the value of usepotzero
217 : !it will be used later on in this subroutine
218 1456 : pawtab%usepotzero=usepotzero
219 :
220 : !==================================================
221 : !1- INITIALIZE DATA RELATED TO ANGULAR MESH
222 : !* ANGULAR GRID
223 : !* REAL SPHERICAL HARMONICS
224 : !* REAL GAUNT COEFFICIENTS
225 :
226 592 : usekden=pawxc_get_usekden(ixc)
227 592 : nabgnt_option=0;if (usekden>0) nabgnt_option=1 ! If kin. ene. density is used, need nabla Gaunt coeffs
228 592 : use_angular_grid=0;if (pawxcdev==0) use_angular_grid=1
229 592 : use_ylm=0;if (pawxcdev==0) use_ylm=1
230 592 : use_ls_ylm=0;if (pawspnorb>0) use_ls_ylm=1
231 592 : ngrad2_ylm=0;if (pawxc_get_xclevel(ixc)>=2) ngrad2_ylm=1
232 592 : if (pawxc_get_uselaplacian(ixc)>0) ngrad2_ylm=2
233 592 : call pawang_free(pawang)
234 : call pawang_init(pawang,gnt_option,nabgnt_option,mpsang-1,nphi,ntheta,nsym,ngrad2_ylm,&
235 592 : & use_angular_grid,use_ylm,use_ls_ylm)
236 : endif
237 :
238 : !*******************
239 : !Loop on atom types
240 : !*******************
241 1464 : do itypat=1,ntypat
242 868 : mesh_size=pawtab(itypat)%mesh_size
243 868 : l_size=pawtab(itypat)%l_size
244 868 : lmn_size=pawtab(itypat)%lmn_size
245 868 : lmn2_size=pawtab(itypat)%lmn2_size
246 868 : basis_size=pawtab(itypat)%basis_size
247 868 : ij_size=pawtab(itypat)%ij_size
248 868 : indlmn => pawtab(itypat)%indlmn(:,:)
249 2604 : ABI_MALLOC(indklmn_,(8,lmn2_size))
250 2604 : ABI_MALLOC(klm_diag,(lmn2_size))
251 2604 : ABI_MALLOC(ff,(mesh_size))
252 1736 : ABI_MALLOC(gg,(mesh_size))
253 1736 : ABI_MALLOC(hh,(mesh_size))
254 1736 : ABI_MALLOC(rad,(mesh_size))
255 797299 : rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
256 :
257 868 : if (pawtab(itypat)%usexcnhat/=usexcnhat) then
258 : write(message, '(7a)' )&
259 0 : & 'You cannot simultaneously use atomic data with different',ch10,&
260 0 : & 'formulation of XC [using compensation charge in XC or not] !',ch10,&
261 0 : & 'Action: change at least one of your atomic data (psp) file',ch10,&
262 0 : & ' or use usexcnhat keyword in input file.'
263 0 : ABI_ERROR(message)
264 : end if
265 :
266 868 : if(rcpaw_update_) then
267 16 : ABI_MALLOC(old_phiphj,(mesh_size,ij_size))
268 12 : ABI_MALLOC(old_tphitphj,(mesh_size,ij_size))
269 16 : ABI_MALLOC(old_qijl,(l_size**2,lmn2_size))
270 16 : ABI_MALLOC(old_vhatijl,(lmn2_size,l_size))
271 279754 : old_phiphj=zero
272 279754 : old_tphitphj=zero
273 126496 : old_qijl=zero
274 18202 : old_vhatijl=zero
275 : goto 3
276 : endif
277 :
278 : ! ==================================================
279 : ! 2- TABULATE SHAPE FUNCTION
280 :
281 : ! Allocated shape function
282 864 : if (pawtab(itypat)%shape_type/=-1) then
283 861 : if (allocated(pawtab(itypat)%shapefunc)) then
284 122 : ABI_FREE(pawtab(itypat)%shapefunc)
285 : end if
286 3444 : ABI_MALLOC(pawtab(itypat)%shapefunc,(mesh_size,l_size))
287 3 : else if (.not.allocated(pawtab(itypat)%shapefunc)) then
288 0 : message='shapefunc should be allocated with shape_type=-1'
289 0 : ABI_ERROR(message)
290 : end if
291 864 : if (allocated(pawtab(itypat)%gnorm)) then
292 143 : ABI_FREE(pawtab(itypat)%gnorm)
293 : end if
294 2592 : ABI_MALLOC(pawtab(itypat)%gnorm,(l_size))
295 :
296 : ! Compute shape function
297 3936 : do il=1,l_size
298 3072 : ll=il-1
299 3072 : call atompaw_shpfun(ll,pawrad(itypat),gnrm,pawtab(itypat),ff)
300 2811939 : pawtab(itypat)%shapefunc(1:mesh_size,il)=ff(1:mesh_size)
301 7008 : pawtab(itypat)%gnorm(il)=gnrm
302 : end do
303 : ! In case of numerical shape function, compute some derivatives
304 864 : if (pawtab(itypat)%shape_type==-1) then
305 3 : if (allocated(pawtab(itypat)%dshpfunc)) then
306 0 : ABI_FREE(pawtab(itypat)%dshpfunc)
307 : end if
308 15 : ABI_MALLOC(pawtab(itypat)%dshpfunc,(mesh_size,l_size,4))
309 6 : ABI_MALLOC(work,(mesh_size))
310 18 : do il=1,l_size
311 15 : call nderiv_gen(pawtab(itypat)%dshpfunc(:,il,1),pawtab(itypat)%shapefunc(:,il),pawrad(itypat))
312 15 : yp1=pawtab(itypat)%dshpfunc(1,il,1);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,1)
313 15 : call spline(rad,pawtab(itypat)%shapefunc(:,il),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,2))
314 15 : yp1=pawtab(itypat)%dshpfunc(1,il,2);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,2)
315 15 : call spline(rad,pawtab(itypat)%dshpfunc(:,il,1),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,3))
316 15 : yp1=pawtab(itypat)%dshpfunc(1,il,3);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,3)
317 18 : call spline(rad,pawtab(itypat)%dshpfunc(:,il,2),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,4))
318 : end do
319 3 : ABI_FREE(work)
320 : end if
321 :
322 : ! In some cases, has to store radial mesh for shape function in pawtab variable
323 864 : if (pawtab(itypat)%shape_type==-1) then
324 3 : if (allocated(pawtab(itypat)%rad_for_spline)) then
325 0 : ABI_FREE(pawtab(itypat)%rad_for_spline)
326 : end if
327 6 : ABI_MALLOC(pawtab(itypat)%rad_for_spline,(mesh_size))
328 1785 : pawtab(itypat)%rad_for_spline(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
329 : end if
330 :
331 : ! In some cases, has to store shape function in reciprocal space
332 864 : if (pawtab(itypat)%has_shapefncg>0) then
333 0 : if (gsqcut_eff<tol8) then
334 0 : message='Computation of shapefncg only possible when gsqcut>0!'
335 0 : ABI_BUG(message)
336 : end if
337 0 : pawtab(itypat)%mqgrid_shp=mqgrid_shp_default
338 0 : if (allocated(pawtab(itypat)%shapefncg)) then
339 0 : ABI_FREE(pawtab(itypat)%shapefncg)
340 : end if
341 0 : if (allocated(pawtab(itypat)%qgrid_shp)) then
342 0 : ABI_FREE(pawtab(itypat)%qgrid_shp)
343 : end if
344 0 : ABI_MALLOC(pawtab(itypat)%shapefncg,(pawtab(itypat)%mqgrid_shp,2,l_size))
345 0 : ABI_MALLOC(pawtab(itypat)%qgrid_shp,(pawtab(itypat)%mqgrid_shp))
346 0 : dq=1.1_dp*sqrt(gsqcut_eff)/dble(pawtab(itypat)%mqgrid_shp-1)
347 0 : do iq=1,pawtab(itypat)%mqgrid_shp
348 0 : pawtab(itypat)%qgrid_shp(iq)=dble(iq-1)*dq
349 : end do
350 0 : ABI_MALLOC(indl,(6,l_size))
351 0 : ABI_MALLOC(rgl,(mesh_size,il))
352 0 : do il=1,l_size
353 0 : indl(:,il)=0;indl(1,il)=il-1;indl(5,il)=il
354 0 : rgl(1:mesh_size,il)=rad(1:mesh_size)*pawtab(itypat)%shapefunc(1:mesh_size,il)
355 : end do
356 : call pawpsp_nl(pawtab(itypat)%shapefncg,indl,l_size,l_size,&
357 0 : & pawtab(itypat)%mqgrid_shp,pawtab(itypat)%qgrid_shp,pawrad(itypat),rgl)
358 0 : pawtab(itypat)%shapefncg=four_pi*pawtab(itypat)%shapefncg
359 0 : ABI_FREE(indl)
360 0 : ABI_FREE(rgl)
361 : else
362 864 : pawtab(itypat)%mqgrid_shp=0
363 : end if
364 :
365 : ! ==================================================
366 : ! 3- COMPUTE indklmn INDEXES GIVING klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn
367 : ! for each klmn=(ilmn,jlmn)
368 :
369 864 : if (allocated(pawtab(itypat)%indklmn)) then
370 143 : ABI_FREE(pawtab(itypat)%indklmn)
371 : end if
372 2592 : ABI_MALLOC(pawtab(itypat)%indklmn,(8,lmn2_size))
373 :
374 71497 : 3 klm_diag=0
375 10307 : do jlmn=1,lmn_size
376 9439 : jl= indlmn(1,jlmn);jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
377 9439 : j0lmn=jlmn*(jlmn-1)/2
378 9439 : j0lm =jlm *(jlm -1)/2
379 9439 : j0ln =jln *(jln -1)/2
380 80936 : do ilmn=1,jlmn
381 70629 : il= indlmn(1,ilmn);ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
382 70629 : klmn=j0lmn+ilmn
383 70629 : if (ilm<=jlm) then
384 65267 : indklmn_(1,klmn)=j0lm+ilm
385 : else
386 5362 : i0lm=ilm*(ilm-1)/2
387 5362 : indklmn_(1,klmn)=i0lm+jlm
388 : end if
389 70629 : if (iln<=jln) then
390 70629 : indklmn_(2,klmn)=j0ln+iln
391 : else
392 0 : i0ln=iln*(iln-1)/2
393 0 : indklmn_(2,klmn)=i0ln+jln
394 : end if
395 70629 : indklmn_(3,klmn)=min(abs(il-jl),lcutdens)
396 70629 : indklmn_(4,klmn)=min(il+jl,lcutdens)
397 70629 : indklmn_(5,klmn)=ilm
398 70629 : indklmn_(6,klmn)=jlm
399 70629 : indklmn_(7,klmn)=ilmn
400 70629 : indklmn_(8,klmn)=jlmn
401 635661 : pawtab(itypat)%indklmn(:,klmn)=indklmn_(:,klmn)
402 80068 : if (ilm==jlm) klm_diag(klmn)=1
403 : end do
404 : end do
405 :
406 868 : if(rcpaw_update_) goto 5
407 :
408 : ! ==================================================
409 : ! 4- COMPUTE various FACTORS/SIZES (depending on (l,m,n))
410 :
411 864 : pawtab(itypat)%usespnorb=pawspnorb
412 864 : pawtab(itypat)%lcut_size=min(l_size,lcutdens+1)
413 :
414 864 : if (allocated(pawtab(itypat)%dltij)) then
415 143 : ABI_FREE(pawtab(itypat)%dltij)
416 : end if
417 2592 : ABI_MALLOC(pawtab(itypat)%dltij,(lmn2_size))
418 68799 : pawtab(itypat)%dltij(:)=two
419 10171 : do ilmn=1,lmn_size
420 10171 : pawtab(itypat)%dltij(ilmn*(ilmn+1)/2)=one
421 : end do
422 :
423 864 : lmnmix=zero
424 1728 : ABI_MALLOC(kmix_tmp,(lmn2_size))
425 10171 : do jlmn=1,lmn_size
426 9307 : jl=indlmn(1,jlmn)
427 10171 : if (jl<=lmix) then
428 9297 : j0lmn=jlmn*(jlmn-1)/2
429 77097 : do ilmn=1,jlmn
430 67800 : il=indlmn(1,ilmn)
431 77097 : if (il<=lmix) then
432 67800 : lmnmix=lmnmix+1
433 67800 : kmix_tmp(lmnmix)=j0lmn+ilmn
434 : end if
435 : end do
436 : end if
437 : end do
438 864 : if (allocated(pawtab(itypat)%kmix)) then
439 143 : ABI_FREE(pawtab(itypat)%kmix)
440 : end if
441 2592 : ABI_MALLOC(pawtab(itypat)%kmix,(lmnmix))
442 864 : pawtab(itypat)%lmnmix_sz=lmnmix
443 68664 : pawtab(itypat)%kmix(1:lmnmix)=kmix_tmp(1:lmnmix)
444 864 : ABI_FREE(kmix_tmp)
445 :
446 : ! ==================================================
447 : ! 5- STORE SOME USEFUL QUANTITIES FROM PARTIAL WAVES
448 :
449 864 : if (allocated(pawtab(itypat)%phiphj)) then
450 143 : ABI_FREE(pawtab(itypat)%phiphj)
451 : end if
452 864 : if (allocated(pawtab(itypat)%tphitphj)) then
453 143 : ABI_FREE(pawtab(itypat)%tphitphj)
454 : end if
455 3456 : ABI_MALLOC(pawtab(itypat)%phiphj,(mesh_size,ij_size))
456 2596 : ABI_MALLOC(pawtab(itypat)%tphitphj,(mesh_size,ij_size))
457 4825 : 5 do jln=1,basis_size
458 3957 : j0ln=jln*(jln-1)/2
459 16475 : do iln=1,jln
460 11650 : kln=j0ln+iln
461 11650 : if(rcpaw_update_) then
462 279750 : old_phiphj(1:mesh_size,kln)=pawtab(itypat)%phiphj(1:mesh_size,kln)
463 279750 : old_tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphitphj(1:mesh_size,kln)
464 : endif
465 : pawtab(itypat)%phiphj(1:mesh_size,kln)=pawtab(itypat)%phi(1:mesh_size,iln)&
466 10921478 : & *pawtab(itypat)%phi(1:mesh_size,jln)
467 : pawtab(itypat)%tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphi(1:mesh_size,iln)&
468 10925435 : & *pawtab(itypat)%tphi(1:mesh_size,jln)
469 : end do
470 : end do
471 :
472 868 : if(rcpaw_update_) goto 6
473 :
474 864 : if (usekden==1) then
475 20 : pw_mesh_size=pawtab(itypat)%partialwave_mesh_size
476 20 : if (allocated(pawtab(itypat)%nablaphi)) then
477 3 : ABI_FREE(pawtab(itypat)%nablaphi)
478 : end if
479 80 : ABI_MALLOC(pawtab(itypat)%nablaphi,(pw_mesh_size,basis_size))
480 20 : if (allocated(pawtab(itypat)%tnablaphi)) then
481 3 : ABI_FREE(pawtab(itypat)%tnablaphi)
482 : end if
483 60 : ABI_MALLOC(pawtab(itypat)%tnablaphi,(pw_mesh_size,basis_size))
484 60 : ABI_MALLOC(der,(pw_mesh_size))
485 120 : do iln=1,basis_size
486 100 : call nderiv_gen(der,pawtab(itypat)%phi(1:pw_mesh_size,iln),pawrad(itypat))
487 : pawtab(itypat)%nablaphi(2:pw_mesh_size,iln)=der(2:pw_mesh_size) &
488 200040 : & -pawtab(itypat)%phi(2:pw_mesh_size,iln)/pawrad(itypat)%rad(2:pw_mesh_size)
489 100 : call nderiv_gen(der,pawtab(itypat)%tphi(1:pw_mesh_size,iln),pawrad(itypat))
490 : pawtab(itypat)%tnablaphi(2:pw_mesh_size,iln)=der(2:pw_mesh_size) &
491 200040 : & -pawtab(itypat)%tphi(2:pw_mesh_size,iln)/pawrad(itypat)%rad(2:pw_mesh_size)
492 100 : call pawrad_deducer0(pawtab(itypat)%nablaphi(1:pw_mesh_size,iln),pw_mesh_size,pawrad(itypat))
493 120 : call pawrad_deducer0(pawtab(itypat)%tnablaphi(1:pw_mesh_size,iln),pw_mesh_size,pawrad(itypat))
494 : end do
495 20 : ABI_FREE(der)
496 20 : pawtab(itypat)%has_nablaphi=2
497 : end if
498 :
499 : ! ==================================================
500 : ! 6- COMPUTE Qijl TERMS AND Sij MATRIX
501 :
502 : ! Compute q_ijL and S_ij=q_ij0
503 864 : if (allocated(pawtab(itypat)%qijl)) then
504 143 : ABI_FREE(pawtab(itypat)%qijl)
505 : end if
506 864 : if (allocated(pawtab(itypat)%sij)) then
507 143 : ABI_FREE(pawtab(itypat)%sij)
508 : end if
509 3456 : ABI_MALLOC(pawtab(itypat)%qijl,(l_size*l_size,lmn2_size))
510 2592 : ABI_MALLOC(pawtab(itypat)%sij,(lmn2_size))
511 1664510 : pawtab(itypat)%qijl=zero
512 68803 : pawtab(itypat)%sij=zero
513 71497 : 6 do klmn=1,lmn2_size
514 70629 : klm=indklmn_(1,klmn);kln=indklmn_(2,klmn)
515 70629 : lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn)
516 213474 : do ll=lmin,lmax,2
517 142845 : lm0=ll*ll+ll+1;ff(1)=zero
518 : ff(2:mesh_size)=(pawtab(itypat)%phiphj (2:mesh_size,kln)&
519 : & -pawtab(itypat)%tphitphj(2:mesh_size,kln))&
520 136680417 : & *rad(2:mesh_size)**ll
521 142845 : call simp_gen(intg,ff,pawrad(itypat))
522 871301 : do mm=-ll,ll
523 657827 : isel=pawang%gntselect(lm0+mm,klm)
524 657827 : if (isel>0.and.rcpaw_update_) old_qijl(lm0+mm,klmn)=pawtab(itypat)%qijl(lm0+mm,klmn)
525 800672 : if (isel>0) pawtab(itypat)%qijl(lm0+mm,klmn)=intg*pawang%realgnt(isel)
526 : end do
527 : end do
528 70629 : if (klm_diag(klmn)==1) pawtab(itypat)%sij(klmn)= &
529 15128 : & pawtab(itypat)%qijl(1,klmn)*sqrt(four_pi)
530 : end do
531 :
532 868 : if(rcpaw_update_) goto 7
533 :
534 : ! ==================================================
535 : ! 7- COMPUTE Eijkl TERMS (Hartree)
536 : ! Compute eventually short-range screened version of Eijkl (Fock)
537 :
538 864 : if (allocated(pawtab(itypat)%eijkl)) then
539 143 : ABI_FREE(pawtab(itypat)%eijkl)
540 : end if
541 3456 : ABI_MALLOC(pawtab(itypat)%eijkl,(lmn2_size,lmn2_size))
542 864 : if (abs(hyb_range_fock)>tol8) then
543 4 : if (allocated(pawtab(itypat)%eijkl_sr)) then
544 0 : ABI_FREE(pawtab(itypat)%eijkl_sr)
545 : end if
546 12 : ABI_MALLOC(pawtab(itypat)%eijkl_sr,(lmn2_size,lmn2_size))
547 : end if
548 :
549 : ! First loop is for eijkl (Hartree)
550 : ! 2nd loop is for eijkl_sr (short-range screened Fock exchange)
551 2604 : 7 do iloop=1,2
552 1736 : if (iloop==2.and.abs(hyb_range_fock)<=tol8) cycle
553 872 : if (iloop==1) eijkl => pawtab(itypat)%eijkl
554 872 : if (iloop==2) eijkl => pawtab(itypat)%eijkl_sr
555 :
556 : ! Compute:
557 : ! vhatL(r) according to eq. (A14) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]]
558 : ! intvhatL=$\int_{0}^{r_c}{vhatL(r) shapefunc_L(r) r^2\,dr}$
559 : ! vhatijL =$\int_{0}^{r_c}{vhatL(r) \tilde{\phi}_i \tilde{\phi}_j \,dr}$
560 : ! -----------------------------------------------------------------
561 1744 : ABI_MALLOC(vhatl,(mesh_size))
562 3488 : ABI_MALLOC(vhatijl,(lmn2_size,l_size))
563 2616 : ABI_MALLOC(intvhatl,(l_size))
564 1141164 : intvhatl(:)=zero;vhatl(:)=zero;vhatijl(:,:)=zero
565 3978 : do il=1,l_size
566 3106 : vhatl(1)=zero;ff(1)=zero
567 2850694 : ff(2:mesh_size)=pawtab(itypat)%shapefunc(2:mesh_size,il)*rad(2:mesh_size)**2
568 3106 : if (iloop==1) call poisson(ff,il-1,pawrad(itypat),vhatl)
569 3106 : if (iloop==2) call poisson(ff,il-1,pawrad(itypat),vhatl,screened_sr_separation=hyb_range_fock)
570 2850694 : vhatl(2:mesh_size)=two*vhatl(2:mesh_size)/rad(2:mesh_size)
571 2853800 : gg(1:mesh_size)=vhatl(1:mesh_size)*ff(1:mesh_size)
572 3106 : call simp_gen(intvhatl(il),gg,pawrad(itypat))
573 338488 : do klmn=1,lmn2_size
574 334510 : kln=indklmn_(2,klmn)
575 321876849 : hh(1:mesh_size)=vhatl(1:mesh_size)*pawtab(itypat)%tphitphj(1:mesh_size,kln)
576 334510 : call simp_gen(vhatijl(klmn,il),hh,pawrad(itypat))
577 337616 : if(rcpaw_update_) then
578 25628322 : hh(1:mesh_size)=vhatl(1:mesh_size)*old_tphitphj(1:mesh_size,kln)
579 18174 : call simp_gen(old_vhatijl(klmn,il),hh,pawrad(itypat))
580 : endif
581 : end do
582 : end do
583 872 : ABI_FREE(vhatl)
584 :
585 : ! Compute:
586 : ! eijkl=$ vh1_ijkl - Vhatijkl - Bijkl - Cijkl$
587 : ! With:
588 : ! $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$
589 : ! $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$
590 : ! $B_ijkl =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$
591 : ! $C_ijkl =\sum_{L,m} {intvhatL*q_ijL*q_klL}$
592 : ! and:
593 : ! vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]]
594 : ! Warning: compute only eijkl for (i,j)<=(k,l)
595 : ! -----------------------------------------------------------------
596 17395166 : eijkl(:,:)=zero
597 4520 : meshsz=pawrad(itypat)%int_meshsz;if (mesh_size>meshsz) ff(meshsz+1:mesh_size)=zero
598 71612 : do klmn=1,lmn2_size
599 70740 : klm=indklmn_(1,klmn);kln=indklmn_(2,klmn)
600 70740 : lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn)
601 214631 : do ll=lmin,lmax,2
602 143019 : lm0=ll*ll+ll+1
603 136353817 : ff(1:meshsz)=pawtab(itypat)%phiphj (1:meshsz,kln)
604 143019 : if (iloop==1) call poisson(ff,ll,pawrad(itypat),gg)
605 143019 : if (iloop==2) call poisson(ff,ll,pawrad(itypat),gg,screened_sr_separation=hyb_range_fock)
606 136353817 : ff(1:meshsz)=pawtab(itypat)%tphitphj(1:meshsz,kln)
607 143019 : if (iloop==1) call poisson(ff,ll,pawrad(itypat),hh)
608 143019 : if (iloop==2) call poisson(ff,ll,pawrad(itypat),hh,screened_sr_separation=hyb_range_fock)
609 143019 : k1min=klmn
610 143019 : if(rcpaw_update_) k1min=1
611 23156491 : do klmn1=k1min,lmn2_size
612 22942732 : klm1=indklmn_(1,klmn1);kln1=indklmn_(2,klmn1)
613 22942732 : lmin1=indklmn_(3,klmn1);lmax1=indklmn_(4,klmn1)
614 22942732 : vh1=zero
615 22942732 : if ((ll.ge.lmin1).and.(ll.le.lmax1)) then
616 17145354 : ff(1)=zero
617 17145354 : if(.not.rcpaw_update_) then
618 : ff(2:meshsz)=(pawtab(itypat)%phiphj (2:meshsz,kln1)*gg(2:meshsz)&
619 : & -pawtab(itypat)%tphitphj(2:meshsz,kln1)*hh(2:meshsz))&
620 12335465436 : & *two/rad(2:meshsz)
621 : else
622 : ff(2:meshsz)=(old_phiphj(2:meshsz,kln1)*gg(2:meshsz)&
623 : & -old_tphitphj(2:meshsz,kln1)*hh(2:meshsz))&
624 7088928960 : & *two/rad(2:meshsz)
625 : endif
626 17145354 : call simp_gen(vh1,ff,pawrad(itypat))
627 : end if
628 150723705 : do mm=-ll,ll
629 127637954 : isel =pawang%gntselect(lm0+mm,klm)
630 127637954 : isel1=pawang%gntselect(lm0+mm,klm1)
631 150580686 : if (isel>0.and.isel1>0) then
632 2150955 : rg =pawang%realgnt(isel)
633 2150955 : rg1=pawang%realgnt(isel1)
634 2150955 : ql =pawtab(itypat)%qijl(lm0+mm,klmn)
635 2150955 : if(.not.rcpaw_update_) then
636 1571363 : ql1=pawtab(itypat)%qijl(lm0+mm,klmn1)
637 : eijkl(klmn,klmn1)=eijkl(klmn,klmn1)&
638 : & +( vh1 *rg *rg1& ! vh1_ijkl
639 : & - vhatijl(klmn ,ll+1)*rg *ql1& ! Vhat_ijkl
640 : & - vhatijl(klmn1,ll+1)*rg1*ql & ! B_ijkl
641 : & - intvhatl(ll+1) *ql *ql1& ! C_ijkl
642 1571363 : & )*two_pi
643 : else
644 579592 : ql1=old_qijl(lm0+mm,klmn1)
645 : eijkl(klmn,klmn1)=eijkl(klmn,klmn1)&
646 : & +( vh1 *rg *rg1& ! vh1_ijkl
647 : & - vhatijl(klmn ,ll+1)*rg *ql1& ! Vhat_ijkl
648 : & - old_vhatijl(klmn1,ll+1)*rg1*ql & ! B_ijkl
649 : & - intvhatl(ll+1) *ql *ql1& ! C_ijkl
650 579592 : & )*two_pi
651 : endif
652 : end if
653 : end do
654 : end do
655 : end do
656 : end do
657 872 : ABI_FREE(vhatijl)
658 2604 : ABI_FREE(intvhatl)
659 : end do ! iloop
660 :
661 868 : if(rcpaw_update_) goto 10
662 :
663 : ! ==================================================
664 : ! 8- COMPUTE gamma_ij TERMS
665 : ! Corrections to get the background right
666 :
667 864 : if (pawtab(itypat)%usepotzero==1) then
668 2 : if (allocated(pawtab(itypat)%gammaij)) then
669 0 : ABI_FREE(pawtab(itypat)%gammaij)
670 : end if
671 6 : ABI_MALLOC(pawtab(itypat)%gammaij,(lmn2_size))
672 4 : ABI_MALLOC(work,(mesh_size))
673 39 : do klmn=1,lmn2_size
674 39 : if (klm_diag(klmn)==1) then
675 13 : kln=indklmn_(2,klmn)
676 13 : ff(1)=zero
677 15690 : ff(2:mesh_size)=pawtab(itypat)%phiphj(2:mesh_size,kln)-pawtab(itypat)%tphitphj(2:mesh_size,kln)
678 : ! First, compute q_ij^00
679 13 : call simp_gen(intg,ff,pawrad(itypat))
680 : ! Second, compute phi^2 - tphi^2 - 4pi*r^2 q_ij^00 g_0(r)
681 15690 : ff(2:mesh_size)= ff(2:mesh_size) - intg*pawtab(itypat)%shapefunc(2:mesh_size,1)*rad(2:mesh_size)**2
682 13 : call poisson(ff,0,pawrad(itypat),work)
683 : ! work is r*vh; should be then multiplied by r to prepare the integration in the sphere
684 15690 : work(1)=zero ; work(2:mesh_size)=work(2:mesh_size)*rad(2:mesh_size)
685 : ! Third, average it over the sphere
686 13 : call simp_gen(intg,work,pawrad(itypat))
687 : ! Finally, store it in pawtab%gammaij
688 13 : pawtab(itypat)%gammaij(klmn)=intg*four_pi
689 : else
690 24 : pawtab(itypat)%gammaij(klmn)=zero
691 : end if
692 : end do
693 2 : ABI_FREE(work)
694 : end if
695 :
696 : ! ==================================================
697 : ! 9- TAKE into account a modified effective mass for the electrons
698 :
699 864 : if (abs(effmass_free-one)>tol8) then
700 0 : if (pawtab(itypat)%has_kij/=2) then
701 0 : message='we need kij and has_kij/=2!'
702 0 : ABI_BUG(message)
703 : end if
704 0 : if (allocated(pawtab(itypat)%dij0)) then
705 0 : pawtab(itypat)%dij0(1:lmn2_size)=pawtab(itypat)%dij0(1:lmn2_size)-pawtab(itypat)%kij(1:lmn2_size)
706 : end if
707 0 : pawtab(itypat)%kij(1:lmn2_size)=pawtab(itypat)%kij(1:lmn2_size)/effmass_free
708 0 : if (allocated(pawtab(itypat)%dij0)) then
709 0 : pawtab(itypat)%dij0(1:lmn2_size)=pawtab(itypat)%dij0(1:lmn2_size)+pawtab(itypat)%kij(1:lmn2_size)
710 : end if
711 : end if
712 :
713 : ! ***********************
714 : ! End Loop on atom types
715 : ! ***********************
716 868 : 10 ABI_FREE(ff)
717 868 : ABI_FREE(gg)
718 868 : ABI_FREE(hh)
719 868 : ABI_FREE(indklmn_)
720 868 : ABI_FREE(klm_diag)
721 868 : ABI_FREE(rad)
722 1464 : if(rcpaw_update_) then
723 4 : ABI_FREE(old_phiphj)
724 4 : ABI_FREE(old_tphitphj)
725 4 : ABI_FREE(old_qijl)
726 4 : ABI_FREE(old_vhatijl)
727 : endif
728 : end do
729 :
730 596 : call timab(553,2,tsec)
731 :
732 : DBG_EXIT("COLL")
733 :
734 1192 : end subroutine pawinit
735 : !!***
736 :
737 : !----------------------------------------------------------------------
738 :
739 : !!****f* m_paw_init/paw_gencond
740 : !! NAME
741 : !! paw_gencond
742 : !!
743 : !! FUNCTION
744 : !! This routine tests whether we have to call pawinit in one of the optdriver
745 : !! routines since important values have changed wrt to the previous dataset.
746 : !! The function uses an internal array to store of the previous values
747 : !!
748 : !! Usage example:
749 : !!
750 : !! call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
751 : !!
752 : !! if (psp_gencond == 1 .or. call_pawinit) then
753 : !! call pawinit(...)
754 : !! call paw_gencond(Dtset, gnt_option, "save", call_pawinit)
755 : !! end if
756 : !!
757 : !! where psp_gencond is the value returned by pspini.
758 : !!
759 : !! INPUT
760 : !! Dtset<type(dataset_type)>=all input variables for this dataset
761 : !! gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated
762 : !! also determine the size of these pointers
763 : !! mode= "test" to test if pawinit must be called
764 : !! "save" to update the internal variables.
765 : !! "reset" to reset the internal variables
766 : !!
767 : !! OUTPUT
768 : !! call_pawinit=True if pawinit must be called. Meaninfull only if mode=="test"
769 : !!
770 : !! SIDE EFFECTS
771 : !! mode=="save" updates the internal variables.
772 : !! "reset" reset the internal variables to -1
773 : !!
774 : !! SOURCE
775 :
776 2032 : subroutine paw_gencond(Dtset,gnt_option,mode,call_pawinit)
777 :
778 : !Arguments ------------------------------------
779 : integer,intent(in) :: gnt_option
780 : logical,intent(out) :: call_pawinit
781 : character(len=*),intent(in) :: mode
782 : type(dataset_type),intent(in) :: Dtset
783 :
784 : !Local variables-------------------------------
785 : !scalars
786 : integer,save :: gencond(10)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/)
787 :
788 : ! *********************************************************************
789 :
790 2032 : call_pawinit = .False.
791 1440 : select case (mode)
792 : case ("test")
793 :
794 : if (gencond(1)/=Dtset%pawlcutd .or.gencond(2) /=Dtset%pawlmix .or.&
795 : gencond(3)/=Dtset%pawnphi .or.gencond(4) /=Dtset%pawntheta.or.&
796 : gencond(5)/=Dtset%pawspnorb .or.gencond(6) /=Dtset%pawxcdev .or.&
797 : gencond(7)/=Dtset%nsym .or.gencond(8) /=gnt_option .or.&
798 1440 : gencond(9)/=Dtset%usepotzero.or.gencond(10)/=Dtset%usekden) call_pawinit = .True.
799 :
800 : case ("save")
801 : ! Update internal values
802 592 : gencond(1)=Dtset%pawlcutd ; gencond(2) =Dtset%pawlmix
803 592 : gencond(3)=Dtset%pawnphi ; gencond(4) =Dtset%pawntheta
804 592 : gencond(5)=Dtset%pawspnorb ; gencond(6) =Dtset%pawxcdev
805 592 : gencond(7)=Dtset%nsym ; gencond(8) =gnt_option
806 592 : gencond(9)=Dtset%usepotzero; gencond(10)=Dtset%usekden
807 :
808 : case ("reset")
809 0 : gencond = -1
810 :
811 : case default
812 2032 : ABI_BUG("Wrong value for mode: "//trim(mode))
813 : end select
814 :
815 2032 : end subroutine paw_gencond
816 : !!***
817 :
818 : !----------------------------------------------------------------------
819 :
820 : END MODULE m_paw_init
821 : !!***
|
