Branch data 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 : 612 : subroutine pawinit(effmass_free,gnt_option,gsqcut_eff,hyb_range_fock,lcutdens,lmix,mpsang,&
166 [ - + - + ]: 612 : & 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 : 612 : integer,allocatable :: indl(:,:),klm_diag(:),kmix_tmp(:)
193 : 612 : integer, contiguous, pointer :: indlmn(:,:)
194 : : real(dp) :: tsec(2)
195 : 612 : real(dp),allocatable :: der(:),ff(:),gg(:),hh(:),indklmn_(:,:),intvhatl(:)
196 : 612 : real(dp),allocatable :: rad(:),rgl(:,:),vhatijl(:,:),vhatl(:),work(:)
197 : 612 : real(dp),pointer :: eijkl(:,:)
198 : 612 : real(dp), allocatable :: old_phiphj(:,:),old_tphitphj(:,:),old_qijl(:,:),old_vhatijl(:,:)
199 : :
200 : : !************************************************************************
201 : :
202 : : DBG_ENTER("COLL")
203 : :
204 : 612 : call timab(553,1,tsec)
205 : :
206 : 612 : ntypat=size(pawtab)
207 [ - + ]: 612 : if (size(pawrad)/=ntypat) then
208 : 0 : ABI_BUG('pawrad and pawtab should have the same size!')
209 : : end if
210 : :
211 : 612 : rcpaw_update_=.false.
212 [ + + ]: 612 : if(present(rcpaw_update)) rcpaw_update_=rcpaw_update
213 : :
214 [ + + ]: 1503 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
215 [ + + ]: 612 : if(.not.rcpaw_update_) then
216 : : !Immediately set the value of usepotzero
217 : : !it will be used later on in this subroutine
218 [ + + ]: 1495 : 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 : 608 : usekden=pawxc_get_usekden(ixc)
227 [ + + ]: 608 : nabgnt_option=0;if (usekden>0) nabgnt_option=1 ! If kin. ene. density is used, need nabla Gaunt coeffs
228 [ + + ]: 608 : use_angular_grid=0;if (pawxcdev==0) use_angular_grid=1
229 [ + + ]: 608 : use_ylm=0;if (pawxcdev==0) use_ylm=1
230 [ + + ]: 608 : use_ls_ylm=0;if (pawspnorb>0) use_ls_ylm=1
231 [ + + ]: 608 : ngrad2_ylm=0;if (pawxc_get_xclevel(ixc)>=2) ngrad2_ylm=1
232 [ + + ]: 608 : if (pawxc_get_uselaplacian(ixc)>0) ngrad2_ylm=2
233 : 608 : call pawang_free(pawang)
234 : : call pawang_init(pawang,gnt_option,nabgnt_option,mpsang-1,nphi,ntheta,nsym,ngrad2_ylm,&
235 : 608 : & use_angular_grid,use_ylm,use_ls_ylm)
236 : : endif
237 : :
238 : : !*******************
239 : : !Loop on atom types
240 : : !*******************
241 [ + + ]: 1503 : do itypat=1,ntypat
242 : 891 : mesh_size=pawtab(itypat)%mesh_size
243 : 891 : l_size=pawtab(itypat)%l_size
244 : 891 : lmn_size=pawtab(itypat)%lmn_size
245 : 891 : lmn2_size=pawtab(itypat)%lmn2_size
246 : 891 : basis_size=pawtab(itypat)%basis_size
247 : 891 : ij_size=pawtab(itypat)%ij_size
248 : 891 : indlmn => pawtab(itypat)%indlmn(:,:)
249 [ - + + - : 2673 : ABI_MALLOC(indklmn_,(8,lmn2_size))
+ - + - -
+ - + ]
250 [ - + + - : 2673 : ABI_MALLOC(klm_diag,(lmn2_size))
+ - + - -
+ - + ]
251 [ - + + - : 2673 : ABI_MALLOC(ff,(mesh_size))
+ - + - -
+ - + -
+ ]
252 [ - + + - : 1782 : ABI_MALLOC(gg,(mesh_size))
- + - + -
+ ]
253 [ - + + - : 1782 : ABI_MALLOC(hh,(mesh_size))
- + - + -
+ ]
254 [ - + + - : 1782 : ABI_MALLOC(rad,(mesh_size))
- + - + ]
255 [ + + ]: 824894 : rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
256 : :
257 [ - + ]: 891 : 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 [ + + ]: 891 : 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 [ + + ]: 887 : if (pawtab(itypat)%shape_type/=-1) then
283 [ + + ]: 884 : if (allocated(pawtab(itypat)%shapefunc)) then
284 : 124 : ABI_FREE(pawtab(itypat)%shapefunc)
285 : : end if
286 [ - + + - : 3536 : 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 [ + + ]: 887 : if (allocated(pawtab(itypat)%gnorm)) then
292 : 148 : ABI_FREE(pawtab(itypat)%gnorm)
293 : : end if
294 [ - + + - : 2661 : ABI_MALLOC(pawtab(itypat)%gnorm,(l_size))
+ - + - -
+ - + -
+ ]
295 : :
296 : : ! Compute shape function
297 [ + + ]: 4034 : do il=1,l_size
298 : 3147 : ll=il-1
299 : 3147 : call atompaw_shpfun(ll,pawrad(itypat),gnrm,pawtab(itypat),ff)
300 [ + + ]: 2899184 : pawtab(itypat)%shapefunc(1:mesh_size,il)=ff(1:mesh_size)
301 : 7181 : pawtab(itypat)%gnorm(il)=gnrm
302 : : end do
303 : : ! In case of numerical shape function, compute some derivatives
304 [ + + ]: 887 : 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 [ + + ]: 887 : 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 [ - + ]: 887 : 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 : 887 : 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 [ + + ]: 887 : if (allocated(pawtab(itypat)%indklmn)) then
370 : 148 : ABI_FREE(pawtab(itypat)%indklmn)
371 : : end if
372 [ - + + - : 2661 : ABI_MALLOC(pawtab(itypat)%indklmn,(8,lmn2_size))
+ - + - -
+ - + -
+ ]
373 : :
374 [ + + ]: 72593 : 3 klm_diag=0
375 [ + + ]: 10534 : do jlmn=1,lmn_size
376 : 9643 : jl= indlmn(1,jlmn);jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
377 : 9643 : j0lmn=jlmn*(jlmn-1)/2
378 : 9643 : j0lm =jlm *(jlm -1)/2
379 : 9643 : j0ln =jln *(jln -1)/2
380 [ + + ]: 82236 : do ilmn=1,jlmn
381 : 71702 : il= indlmn(1,ilmn);ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
382 : 71702 : klmn=j0lmn+ilmn
383 [ + + ]: 71702 : if (ilm<=jlm) then
384 : 66261 : indklmn_(1,klmn)=j0lm+ilm
385 : : else
386 : 5441 : i0lm=ilm*(ilm-1)/2
387 : 5441 : indklmn_(1,klmn)=i0lm+jlm
388 : : end if
389 [ + - ]: 71702 : if (iln<=jln) then
390 : 71702 : 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 : 71702 : indklmn_(3,klmn)=min(abs(il-jl),lcutdens)
396 : 71702 : indklmn_(4,klmn)=min(il+jl,lcutdens)
397 : 71702 : indklmn_(5,klmn)=ilm
398 : 71702 : indklmn_(6,klmn)=jlm
399 : 71702 : indklmn_(7,klmn)=ilmn
400 : 71702 : indklmn_(8,klmn)=jlmn
401 [ + + ]: 645318 : pawtab(itypat)%indklmn(:,klmn)=indklmn_(:,klmn)
402 [ + + ]: 81345 : if (ilm==jlm) klm_diag(klmn)=1
403 : : end do
404 : : end do
405 : :
406 [ + + ]: 891 : if(rcpaw_update_) goto 5
407 : :
408 : : ! ==================================================
409 : : ! 4- COMPUTE various FACTORS/SIZES (depending on (l,m,n))
410 : :
411 : 887 : pawtab(itypat)%usespnorb=pawspnorb
412 : 887 : pawtab(itypat)%lcut_size=min(l_size,lcutdens+1)
413 : :
414 [ + + ]: 887 : if (allocated(pawtab(itypat)%dltij)) then
415 : 148 : ABI_FREE(pawtab(itypat)%dltij)
416 : : end if
417 [ - + + - : 2661 : ABI_MALLOC(pawtab(itypat)%dltij,(lmn2_size))
+ - + - -
+ - + -
+ ]
418 [ + + ]: 69895 : pawtab(itypat)%dltij(:)=two
419 [ + + ]: 10398 : do ilmn=1,lmn_size
420 : 10398 : pawtab(itypat)%dltij(ilmn*(ilmn+1)/2)=one
421 : : end do
422 : :
423 : 887 : lmnmix=zero
424 [ - + + - : 1774 : ABI_MALLOC(kmix_tmp,(lmn2_size))
- + - + ]
425 [ + + ]: 10398 : do jlmn=1,lmn_size
426 : 9511 : jl=indlmn(1,jlmn)
427 [ + + ]: 10398 : if (jl<=lmix) then
428 : 9501 : j0lmn=jlmn*(jlmn-1)/2
429 [ + + ]: 78374 : do ilmn=1,jlmn
430 : 68873 : il=indlmn(1,ilmn)
431 [ + - ]: 78374 : if (il<=lmix) then
432 : 68873 : lmnmix=lmnmix+1
433 : 68873 : kmix_tmp(lmnmix)=j0lmn+ilmn
434 : : end if
435 : : end do
436 : : end if
437 : : end do
438 [ + + ]: 887 : if (allocated(pawtab(itypat)%kmix)) then
439 : 148 : ABI_FREE(pawtab(itypat)%kmix)
440 : : end if
441 [ - + + - : 2661 : ABI_MALLOC(pawtab(itypat)%kmix,(lmnmix))
+ - + - -
+ - + -
+ ]
442 : 887 : pawtab(itypat)%lmnmix_sz=lmnmix
443 [ + + ]: 69760 : pawtab(itypat)%kmix(1:lmnmix)=kmix_tmp(1:lmnmix)
444 [ - + ]: 887 : ABI_FREE(kmix_tmp)
445 : :
446 : : ! ==================================================
447 : : ! 5- STORE SOME USEFUL QUANTITIES FROM PARTIAL WAVES
448 : :
449 [ + + ]: 887 : if (allocated(pawtab(itypat)%phiphj)) then
450 : 148 : ABI_FREE(pawtab(itypat)%phiphj)
451 : : end if
452 [ + + ]: 887 : if (allocated(pawtab(itypat)%tphitphj)) then
453 : 148 : ABI_FREE(pawtab(itypat)%tphitphj)
454 : : end if
455 [ - + + - : 3548 : ABI_MALLOC(pawtab(itypat)%phiphj,(mesh_size,ij_size))
- + + - +
- + - - +
- + - + ]
456 [ - + + - : 2665 : ABI_MALLOC(pawtab(itypat)%tphitphj,(mesh_size,ij_size))
- + + - -
+ - + -
+ ]
457 [ + + ]: 4944 : 5 do jln=1,basis_size
458 : 4053 : j0ln=jln*(jln-1)/2
459 [ + + ]: 16845 : do iln=1,jln
460 : 11901 : kln=j0ln+iln
461 [ + + ]: 11901 : 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 [ + + ]: 11217746 : & *pawtab(itypat)%phi(1:mesh_size,jln)
467 : : pawtab(itypat)%tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphi(1:mesh_size,iln)&
468 [ + + ]: 11221799 : & *pawtab(itypat)%tphi(1:mesh_size,jln)
469 : : end do
470 : : end do
471 : :
472 [ + + ]: 891 : if(rcpaw_update_) goto 6
473 : :
474 [ + + ]: 887 : 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 [ + + ]: 887 : if (allocated(pawtab(itypat)%qijl)) then
504 : 148 : ABI_FREE(pawtab(itypat)%qijl)
505 : : end if
506 [ + + ]: 887 : if (allocated(pawtab(itypat)%sij)) then
507 : 148 : ABI_FREE(pawtab(itypat)%sij)
508 : : end if
509 [ - + + - : 3548 : ABI_MALLOC(pawtab(itypat)%qijl,(l_size*l_size,lmn2_size))
- + + - +
- + - - +
- + - + ]
510 [ - + + - : 2661 : ABI_MALLOC(pawtab(itypat)%sij,(lmn2_size))
+ - + - -
+ - + -
+ ]
511 [ + + + + ]: 1680911 : pawtab(itypat)%qijl=zero
512 [ + + ]: 69899 : pawtab(itypat)%sij=zero
513 [ + + ]: 72593 : 6 do klmn=1,lmn2_size
514 : 71702 : klm=indklmn_(1,klmn);kln=indklmn_(2,klmn)
515 : 71702 : lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn)
516 [ - + + + ]: 216393 : do ll=lmin,lmax,2
517 : 144691 : 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 [ + + ]: 138798256 : & *rad(2:mesh_size)**ll
521 : 144691 : call simp_gen(intg,ff,pawrad(itypat))
522 [ + + ]: 880690 : do mm=-ll,ll
523 : 664297 : isel=pawang%gntselect(lm0+mm,klm)
524 [ + + ]: 664297 : if (isel>0.and.rcpaw_update_) old_qijl(lm0+mm,klmn)=pawtab(itypat)%qijl(lm0+mm,klmn)
525 [ + + ]: 808988 : if (isel>0) pawtab(itypat)%qijl(lm0+mm,klmn)=intg*pawang%realgnt(isel)
526 : : end do
527 : : end do
528 [ + + ]: 71702 : if (klm_diag(klmn)==1) pawtab(itypat)%sij(klmn)= &
529 : 15452 : & pawtab(itypat)%qijl(1,klmn)*sqrt(four_pi)
530 : : end do
531 : :
532 [ + + ]: 891 : 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 [ + + ]: 887 : if (allocated(pawtab(itypat)%eijkl)) then
539 : 148 : ABI_FREE(pawtab(itypat)%eijkl)
540 : : end if
541 [ - + + - : 3548 : ABI_MALLOC(pawtab(itypat)%eijkl,(lmn2_size,lmn2_size))
- + + - +
- + - - +
- + - + ]
542 [ + + ]: 887 : 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 [ + + ]: 2673 : 7 do iloop=1,2
552 [ + + + + ]: 1782 : if (iloop==2.and.abs(hyb_range_fock)<=tol8) cycle
553 [ + + ]: 895 : if (iloop==1) eijkl => pawtab(itypat)%eijkl
554 [ + + ]: 895 : 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 [ - + + - : 1790 : ABI_MALLOC(vhatl,(mesh_size))
- + - + -
+ ]
562 [ - + + - : 3580 : ABI_MALLOC(vhatijl,(lmn2_size,l_size))
- + + - +
- + - - +
- + ]
563 [ - + + - : 2685 : ABI_MALLOC(intvhatl,(l_size))
+ - + - -
+ - + ]
564 [ + + + + : 1172834 : intvhatl(:)=zero;vhatl(:)=zero;vhatijl(:,:)=zero
+ + + + ]
565 [ + + ]: 4076 : do il=1,l_size
566 : 3181 : vhatl(1)=zero;ff(1)=zero
567 [ + + ]: 2937864 : ff(2:mesh_size)=pawtab(itypat)%shapefunc(2:mesh_size,il)*rad(2:mesh_size)**2
568 [ + + ]: 3181 : if (iloop==1) call poisson(ff,il-1,pawrad(itypat),vhatl)
569 [ + + ]: 3181 : if (iloop==2) call poisson(ff,il-1,pawrad(itypat),vhatl,screened_sr_separation=hyb_range_fock)
570 [ + + ]: 2937864 : vhatl(2:mesh_size)=two*vhatl(2:mesh_size)/rad(2:mesh_size)
571 [ + + ]: 2941045 : gg(1:mesh_size)=vhatl(1:mesh_size)*ff(1:mesh_size)
572 : 3181 : call simp_gen(intvhatl(il),gg,pawrad(itypat))
573 [ + + ]: 342511 : do klmn=1,lmn2_size
574 : 338435 : kln=indklmn_(2,klmn)
575 [ + + ]: 326242119 : hh(1:mesh_size)=vhatl(1:mesh_size)*pawtab(itypat)%tphitphj(1:mesh_size,kln)
576 : 338435 : call simp_gen(vhatijl(klmn,il),hh,pawrad(itypat))
577 [ + + ]: 341616 : 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 [ - + ]: 895 : 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 [ + + + + ]: 17467985 : eijkl(:,:)=zero
597 [ + + + + ]: 4638 : meshsz=pawrad(itypat)%int_meshsz;if (mesh_size>meshsz) ff(meshsz+1:mesh_size)=zero
598 [ + + ]: 72708 : do klmn=1,lmn2_size
599 : 71813 : klm=indklmn_(1,klmn);kln=indklmn_(2,klmn)
600 : 71813 : lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn)
601 [ - + + + ]: 217573 : do ll=lmin,lmax,2
602 : 144865 : lm0=ll*ll+ll+1
603 [ + + ]: 138466562 : ff(1:meshsz)=pawtab(itypat)%phiphj (1:meshsz,kln)
604 [ + + ]: 144865 : if (iloop==1) call poisson(ff,ll,pawrad(itypat),gg)
605 [ + + ]: 144865 : if (iloop==2) call poisson(ff,ll,pawrad(itypat),gg,screened_sr_separation=hyb_range_fock)
606 [ + + ]: 138466562 : ff(1:meshsz)=pawtab(itypat)%tphitphj(1:meshsz,kln)
607 [ + + ]: 144865 : if (iloop==1) call poisson(ff,ll,pawrad(itypat),hh)
608 [ + + ]: 144865 : if (iloop==2) call poisson(ff,ll,pawrad(itypat),hh,screened_sr_separation=hyb_range_fock)
609 : 144865 : k1min=klmn
610 [ + + ]: 144865 : if(rcpaw_update_) k1min=1
611 [ + + ]: 23221226 : do klmn1=k1min,lmn2_size
612 : 23004548 : klm1=indklmn_(1,klmn1);kln1=indklmn_(2,klmn1)
613 : 23004548 : lmin1=indklmn_(3,klmn1);lmax1=indklmn_(4,klmn1)
614 : 23004548 : vh1=zero
615 [ + + + + ]: 23004548 : if ((ll.ge.lmin1).and.(ll.le.lmax1)) then
616 : 17193962 : ff(1)=zero
617 [ + + ]: 17193962 : 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 [ + + ]: 12392919933 : & *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 : 17193962 : call simp_gen(vh1,ff,pawrad(itypat))
627 : : end if
628 [ + + ]: 151017735 : do mm=-ll,ll
629 : 127868322 : isel =pawang%gntselect(lm0+mm,klm)
630 : 127868322 : isel1=pawang%gntselect(lm0+mm,klm1)
631 [ + + ]: 150872870 : if (isel>0.and.isel1>0) then
632 : 2159975 : rg =pawang%realgnt(isel)
633 : 2159975 : rg1=pawang%realgnt(isel1)
634 : 2159975 : ql =pawtab(itypat)%qijl(lm0+mm,klmn)
635 [ + + ]: 2159975 : if(.not.rcpaw_update_) then
636 : 1580383 : 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 : 1580383 : & )*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 [ - + ]: 895 : ABI_FREE(vhatijl)
658 [ - + ]: 2673 : ABI_FREE(intvhatl)
659 : : end do ! iloop
660 : :
661 [ + + ]: 891 : if(rcpaw_update_) goto 10
662 : :
663 : : ! ==================================================
664 : : ! 8- COMPUTE gamma_ij TERMS
665 : : ! Corrections to get the background right
666 : :
667 [ + + ]: 887 : 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 [ - + ]: 887 : 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 [ - + ]: 891 : 10 ABI_FREE(ff)
717 [ - + ]: 891 : ABI_FREE(gg)
718 [ - + ]: 891 : ABI_FREE(hh)
719 [ - + ]: 891 : ABI_FREE(indklmn_)
720 [ - + ]: 891 : ABI_FREE(klm_diag)
721 [ - + ]: 891 : ABI_FREE(rad)
722 [ + + ]: 1503 : 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 : 612 : call timab(553,2,tsec)
731 : :
732 : : DBG_EXIT("COLL")
733 : :
734 [ - + - + : 1224 : 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 : 2085 : 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 : 2085 : call_pawinit = .False.
791 : 1477 : 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 [ + + + + : 1477 : gencond(9)/=Dtset%usepotzero.or.gencond(10)/=Dtset%usekden) call_pawinit = .True.
+ - + - +
+ + + + +
+ + + + -
+ ]
799 : :
800 : : case ("save")
801 : : ! Update internal values
802 : 608 : gencond(1)=Dtset%pawlcutd ; gencond(2) =Dtset%pawlmix
803 : 608 : gencond(3)=Dtset%pawnphi ; gencond(4) =Dtset%pawntheta
804 : 608 : gencond(5)=Dtset%pawspnorb ; gencond(6) =Dtset%pawxcdev
805 : 608 : gencond(7)=Dtset%nsym ; gencond(8) =gnt_option
806 : 608 : gencond(9)=Dtset%usepotzero; gencond(10)=Dtset%usekden
807 : :
808 : : case ("reset")
809 [ # # ]: 0 : gencond = -1
810 : :
811 : : case default
812 [ + + - - : 2085 : ABI_BUG("Wrong value for mode: "//trim(mode))
- - ]
813 : : end select
814 : :
815 : 2085 : end subroutine paw_gencond
816 : : !!***
817 : :
818 : : !----------------------------------------------------------------------
819 : :
820 : : END MODULE m_paw_init
821 : : !!***
|
