Branch data Line data Source code
1 : : !!****m* ABINIT/m_psp6
2 : : !! NAME
3 : : !! m_psp6
4 : : !!
5 : : !! FUNCTION
6 : : !! Initialize pspcod=6 (Pseudopotentials from the fhi98pp code):
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (XG, AF, GJ,FJ,MT, DRH)
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_psp6
23 : :
24 : : use defs_basis
25 : : use m_splines
26 : : use m_errors
27 : : use m_abicore
28 : :
29 : : use m_numeric_tools, only : smooth, ctrap
30 : : use m_psptk, only : psp5lo, psp5nl, cc_derivatives
31 : :
32 : : implicit none
33 : :
34 : : private
35 : : !!***
36 : :
37 : : public :: psp6in
38 : : !!***
39 : :
40 : : contains
41 : : !!***
42 : :
43 : : !!****f* m_psp6/psp6in
44 : : !! NAME
45 : : !! psp6in
46 : : !!
47 : : !! FUNCTION
48 : : !! Initialize pspcod=6 (Pseudopotentials from the fhi98pp code):
49 : : !! continue to read the corresponding file, then compute the
50 : : !! local and non-local potentials.
51 : : !!
52 : : !! INPUTS
53 : : !! lloc=angular momentum choice of local pseudopotential
54 : : !! lmax=value of lmax mentioned at the second line of the psp file
55 : : !! lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
56 : : !! =if useylm=0, max number of (l,n) comp. over all type of psps
57 : : !! lnmax=max. number of (l,n) components over all type of psps
58 : : !! mmax=maximum number of points in real space grid in the psp file
59 : : !! angular momentum of nonlocal pseudopotential
60 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
61 : : !! mqgrid=dimension of q (or G) grid for arrays.
62 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
63 : : !! optnlxccc=option for nl XC core correction (input variable)
64 : : !! positron=0 if electron GS calculation
65 : : !! 1 if positron GS calculation
66 : : !! 2 if electron GS calculation in presence of the positron
67 : : !! qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
68 : : !! useylm=governs the way the nonlocal operator is to be applied:
69 : : !! 1=using Ylm, 0=using Legendre polynomials
70 : : !! zion=nominal valence of atom as specified in psp file
71 : : !! znucl=nuclear number of atom as specified in psp file
72 : : !!
73 : : !! OUTPUT
74 : : !! ekb(lnmax)=Kleinman-Bylander energy,
75 : : !! {{\ \begin{equation}
76 : : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
77 : : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
78 : : !! \end{equation} }}
79 : : !! for each (l,n)
80 : : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree)
81 : : !! ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
82 : : !! second derivative from spline fit for each angular momentum and
83 : : !! each projector
84 : : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
85 : : !! or i=lmn (if useylm=1)
86 : : !! nproj(mpsang)=number of projection functions for each angular momentum
87 : : !! qchrg is not used, and could be suppressed later
88 : : !! vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit
89 : : !! xcccrc=XC core correction cutoff radius (bohr)
90 : : !! xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
91 : : !!
92 : : !! SOURCE
93 : :
94 : 179 : subroutine psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
95 : 179 : & mmax,mpsang,mqgrid,nproj,n1xccc,optnlxccc,positron,qchrg,qgrid,&
96 : 179 : & useylm,vlspl,xcccrc,xccc1d,zion,znucl)
97 : :
98 : : !Arguments ------------------------------------
99 : : !scalars
100 : : integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc
101 : : integer,intent(in) :: optnlxccc,positron,useylm
102 : : real(dp),intent(in) :: zion,znucl
103 : : real(dp),intent(out) :: epsatm,qchrg,xcccrc
104 : : !arrays
105 : : integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang)
106 : : real(dp),intent(in) :: qgrid(mqgrid)
107 : : real(dp),intent(out) :: ekb(lnmax),vlspl(mqgrid,2) !vz_i
108 : : real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
109 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
110 : :
111 : : !Local variables-------------------------------
112 : : !scalars
113 : : integer :: ii,index,ipsang,irad,jj,jpsang,mm,mmax2
114 : : real(dp) :: al,al_announced,amesh,fchrg,ratio,rchrg,yp1,ypn
115 : : character(len=3) :: testxc
116 : : character(len=500) :: message,errmsg
117 : : !arrays
118 : 179 : real(dp),allocatable :: ekb_tmp(:),ffspl_tmp(:,:,:),rad(:),vloc(:)
119 : : !real(dp),allocatable :: radbis
120 : 179 : real(dp),allocatable :: vpspll(:,:),wfll(:,:),work_space(:),work_spl(:)
121 : :
122 : : ! ***************************************************************************
123 : :
124 : : !File format of formatted fhi psp input, as adapted for use
125 : : !by the ABINIT code (the 3 first lines
126 : : !have already been read in calling -pspatm- routine) :
127 : :
128 : : !(1) title (character) line
129 : : !(2) znucl,zion,pspdat
130 : : !(3) pspcod,pspxc,lmax,lloc,mmax,r2well
131 : : !(4) rchrg,fchrg,qchrg
132 : : !Note : prior to version 2.2, this 4th line started with 4-- ,
133 : : !and no core-correction was available.
134 : : !(5)-(18) -empty-
135 : : !(19) mmax, amesh ( mesh increment r(m+1)/r(m) )
136 : : !Then, for ll=0,lmax :
137 : : !for irad=1,mmax : irad, r(irad), upsp(irad,ll), vpsp(irad,ll)
138 : :
139 [ + - ]: 179 : read (tmp_unit, '(a3)', err=10, iomsg=errmsg) testxc
140 [ + + ]: 179 : if(testxc/='4--')then
141 : 102 : backspace(tmp_unit)
142 [ + - ]: 102 : read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
143 : 102 : write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
144 : 102 : call wrtout(ab_out,message,'COLL')
145 : 102 : call wrtout(std_out, message,'COLL')
146 : : else
147 : 77 : write(message, '(a)' ) ' No XC core correction.'
148 : 77 : call wrtout(ab_out,message,'COLL')
149 : 77 : call wrtout(std_out, message,'COLL')
150 : 77 : rchrg=zero ; fchrg=zero ; qchrg=zero
151 : : end if
152 [ + + ]: 2685 : do ii=5,18
153 [ + - ]: 2685 : read(tmp_unit,*, err=10, iomsg=errmsg)
154 : : end do
155 : :
156 [ + + - + ]: 179 : if (positron==1.and.abs(fchrg)<=tol14) then
157 : : write(message,'(5a)')&
158 : 0 : & 'You can only perform positronic ground-state calculations (positron=1)',ch10,&
159 : 0 : & 'using fhi pseudopotentials with a core density (fchrg>0)',ch10,&
160 : 0 : & 'Action: change your psp file (add fchrg>0).'
161 : 0 : ABI_ERROR(message)
162 : : end if
163 : : !--------------------------------------------------------------------
164 : : !Will now proceed at the reading of pots and wfs
165 : :
166 : : !rad(:)=radial grid r(i)
167 : : !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials
168 : : !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions
169 : :
170 [ - + + - : 537 : ABI_MALLOC(rad,(mmax))
+ - + - -
+ - + ]
171 [ - + + - : 716 : ABI_MALLOC(vpspll,(mmax,mpsang))
- + + - +
- + - - +
- + ]
172 [ - + + - : 537 : ABI_MALLOC(wfll,(mmax,mpsang))
- + + - -
+ - + ]
173 : :
174 : : !Read atomic pseudopotential for each l, filling up arrays vpspll
175 : : !and wfll. Also set up rad array (actually read more than once)
176 : : !Note: put each l into vpspll(:,l+1)
177 [ + + ]: 690 : do ipsang=1,lmax+1
178 : 511 : nproj(ipsang)=1
179 [ + - ]: 511 : read(tmp_unit,*, err=10, iomsg=errmsg)mmax2,amesh
180 [ + + ]: 511 : if(ipsang==1)then
181 : 179 : write(message, '(f10.6,t20,a)' ) amesh,' amesh (Hamman grid)'
182 : : al_announced=log(amesh)
183 : 179 : call wrtout(ab_out,message,'COLL')
184 : 179 : call wrtout(std_out, message,'COLL')
185 : : end if
186 [ + + ]: 275990 : do irad=1,mmax
187 [ + - ]: 275811 : read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),wfll(irad,ipsang),vpspll(irad,ipsang)
188 : : end do
189 : : end do
190 : :
191 : :
192 : : !Generate core charge function and derivatives, if needed
193 [ + + ]: 179 : if(fchrg>tol14)then
194 : :
195 [ + + ]: 26 : if (positron==1) then
196 : 1 : call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,vh_tnzc=vpspll(:,lloc+1))
197 [ + + ]: 25 : else if(optnlxccc==1)then
198 : 24 : call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl)
199 [ + - ]: 1 : else if(optnlxccc==2)then
200 : 1 : call psp6cc_drh(mmax,n1xccc,rchrg,xccc1d)
201 : : end if
202 : :
203 : : ! The core charge function for pspcod=6 becomes zero beyond rchrg.
204 : : ! Thus xcccrc must be set equal to rchrg.
205 : 26 : xcccrc=rchrg
206 : : else
207 [ + + + + ]: 676341 : xccc1d(:,:)=zero
208 : 153 : xcccrc=zero
209 : : end if
210 : :
211 : : !Compute in real(dp) al : the announced amesh is inaccurate.
212 : 179 : ratio=rad(mmax)/rad(1)
213 : 179 : al=log(ratio)/dble(mmax-1)
214 : :
215 : : !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg.
216 [ - + + - : 358 : ABI_MALLOC(vloc,(mmax))
- + - + ]
217 : : !Copy appropriate nonlocal psp for use as local one
218 [ + + ]: 96577 : vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 )
219 : :
220 : : !--------------------------------------------------------------------
221 : : !Carry out calculations for local (lloc) pseudopotential.
222 : : !Obtain Fourier transform (1-d sine transform) to get q^2 V(q).
223 : :
224 : : call psp5lo(al,epsatm,mmax,mqgrid,qgrid,&
225 : 179 : & vlspl(:,1),rad,vloc,yp1,ypn,zion)
226 : :
227 : : !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
228 [ - + + - : 537 : ABI_MALLOC(work_space,(mqgrid))
+ - + - -
+ - + ]
229 [ - + + - : 358 : ABI_MALLOC(work_spl,(mqgrid))
- + - + ]
230 : 179 : call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
231 [ + + ]: 535106 : vlspl(:,2)=work_spl(:)
232 [ - + ]: 179 : ABI_FREE(work_space)
233 [ - + ]: 179 : ABI_FREE(work_spl)
234 : :
235 : : !--------------------------------------------------------------------
236 : : !Take care of non-local part
237 : :
238 [ - + + - : 537 : ABI_MALLOC(ekb_tmp,(mpsang))
+ - + - -
+ - + ]
239 [ - + + - : 716 : ABI_MALLOC(ffspl_tmp,(mqgrid,2,mpsang))
- + + - +
- + - - +
- + ]
240 : :
241 : : !Zero out all Kleinman-Bylander energies to initialize
242 [ + + ]: 696 : ekb_tmp(:)=zero
243 [ + + ]: 714 : ekb(:)=zero
244 : :
245 : : !Allow for option of no nonlocal corrections (lloc=lmax=0)
246 [ + + + + ]: 179 : if (lloc==0.and.lmax==0) then
247 : 9 : write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
248 : 9 : call wrtout(ab_out,message,'COLL')
249 : 9 : call wrtout(std_out, message,'COLL')
250 : : else
251 : :
252 : : ! ----------------------------------------------------------------------
253 : : ! Compute KB form factors and fit splines
254 : :
255 : 170 : call psp5nl(al,ekb_tmp,ffspl_tmp,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc, vpspll,wfll)
256 : :
257 : : end if
258 : :
259 [ + + + + ]: 4540 : jj=0;index=0;indlmn(:,:)=0
260 [ + + ]: 690 : do ipsang=1,lmax+1
261 : : ! nproj had been set at 1, by default
262 [ + + ]: 511 : if(abs(ekb_tmp(ipsang))<tol10)then
263 : 179 : nproj(ipsang)=0
264 : : end if
265 : : ! Possible values for nproj in this routine : 0 or 1.
266 [ + + ]: 690 : if(nproj(ipsang)==1)then
267 [ + + ]: 332 : if (useylm==1) then
268 : 24 : jj=jj+1
269 [ + + ]: 84 : do mm=1,2*ipsang-1
270 : 60 : index=index+1
271 : 60 : indlmn(1,index)=ipsang-1
272 : 60 : indlmn(2,index)=mm-ipsang
273 : 60 : indlmn(3,index)=1
274 : 60 : indlmn(4,index)=mm+(ipsang-1)*(ipsang-1)
275 : 60 : indlmn(5,index)=jj
276 : 84 : indlmn(6,index)=1
277 : : end do
278 : : else
279 : 308 : jj=jj+1
280 : 308 : index=index+1
281 : 308 : indlmn(1,index)=ipsang-1
282 : 308 : indlmn(2,index)=0
283 : 308 : indlmn(3,index)=1
284 : 308 : indlmn(4,index)=ipsang+(ipsang-1)*(ipsang-1)
285 : 308 : indlmn(5,index)=jj
286 : 308 : indlmn(6,index)=1
287 : : end if
288 : : end if
289 : : end do
290 : : !Transfer ekb and ffspl to their definitive location
291 : 179 : jpsang=1
292 [ + + ]: 690 : do ipsang=1,lmax+1
293 [ + + ]: 690 : if(nproj(ipsang)/=0)then
294 : 332 : ekb(jpsang)=ekb_tmp(ipsang)
295 [ + + + + ]: 1980156 : ffspl(:,:,jpsang)=ffspl_tmp(:,:,ipsang)
296 : 332 : jpsang=jpsang+1
297 [ - + ]: 332 : if(jpsang>lnmax)then
298 : : write(message,'(3a,2i6)')&
299 : 0 : & 'Problem with the dimension of the ekb and ffspl arrays.',ch10,&
300 : 0 : & 'ipsang,lnmax=',ipsang,lnmax
301 : : end if
302 : : end if
303 : : end do
304 : :
305 [ - + ]: 179 : ABI_FREE(ekb_tmp)
306 [ - + ]: 179 : ABI_FREE(ffspl_tmp)
307 [ - + ]: 179 : ABI_FREE(vpspll)
308 : 179 : ABI_FREE(rad)
309 [ - + ]: 179 : ABI_FREE(vloc)
310 [ - + ]: 179 : ABI_FREE(wfll)
311 : :
312 : 179 : return
313 : :
314 : : ! Handle IO error
315 : : 10 continue
316 : 0 : ABI_ERROR(errmsg)
317 : :
318 [ - + - + : 179 : end subroutine psp6in
- + ]
319 : : !!***
320 : :
321 : : !!****f* m_psp6/psp6cc
322 : : !! NAME
323 : : !! psp6cc
324 : : !!
325 : : !! FUNCTION
326 : : !! Compute the core charge density, for use in the XC core
327 : : !! correction, following the function definition valid
328 : : !! for the format 6 of pseudopotentials.
329 : : !!
330 : : !! INPUTS
331 : : !! mmax=maximum number of points in real space grid in the psp file
332 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
333 : : !! rchrg=cut-off radius for the core density
334 : : !! znucl=nuclear number of atom as specified in psp file
335 : : !!
336 : : !! OUTPUT
337 : : !! xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
338 : : !! Optional output:
339 : : !! vh_tnzc(mmax) = Hartree potential induced by density tild_[n_Z+n_core]
340 : : !! (pseudized [n_Z+n_core], where n_Z=ions, n_core=core electrons)
341 : : !! using a simple pseudization scheme
342 : : !!
343 : : !! SOURCE
344 : :
345 : 25 : subroutine psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,&
346 : : & vh_tnzc) ! optional argument
347 : :
348 : : !Arguments ------------------------------------
349 : : !scalars
350 : : integer,intent(in) :: mmax,n1xccc
351 : : real(dp),intent(in) :: rchrg,znucl
352 : : !arrays
353 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
354 : : real(dp),intent(out),optional :: vh_tnzc(mmax)
355 : :
356 : : !Local variables-------------------------------
357 : : !scalars
358 : : integer :: i1xccc,irad
359 : : real(dp) :: der1,dern
360 : : character(len=500) :: errmsg
361 : : !arrays
362 : 25 : real(dp),allocatable :: ff(:),ff1(:),ff2(:),ff3(:),gg(:),gg1(:),gg2(:),gg3(:)
363 : 25 : real(dp),allocatable :: gg4(:),nc(:),rad(:),work(:),xx(:)
364 : :
365 : : !**********************************************************************
366 : :
367 [ - + + - : 75 : ABI_MALLOC(ff,(mmax))
+ - + - -
+ - + ]
368 [ - + + - : 50 : ABI_MALLOC(ff1,(mmax))
- + - + ]
369 [ - + + - : 50 : ABI_MALLOC(ff2,(mmax))
- + - + ]
370 [ - + + - : 50 : ABI_MALLOC(ff3,(mmax))
- + - + ]
371 [ - + + - : 50 : ABI_MALLOC(rad,(mmax))
- + - + ]
372 [ - + + - : 75 : ABI_MALLOC(gg,(n1xccc))
+ - + - -
+ - + ]
373 [ - + + - : 50 : ABI_MALLOC(gg1,(n1xccc))
- + - + ]
374 [ - + + - : 50 : ABI_MALLOC(gg2,(n1xccc))
- + - + ]
375 [ - + + - : 50 : ABI_MALLOC(gg3,(n1xccc))
- + - + ]
376 [ - + + - : 50 : ABI_MALLOC(gg4,(n1xccc))
- + - + ]
377 [ - + + - : 50 : ABI_MALLOC(work,(n1xccc))
- + - + ]
378 [ - + + - : 50 : ABI_MALLOC(xx,(n1xccc))
- + + - ]
379 : :
380 : : !read from pp file the model core charge (ff) and first (ff1) and
381 : : !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid
382 : : !the input functions contain the 4pi factor, it must be rescaled.
383 : :
384 [ + + ]: 16179 : do irad=1,mmax
385 [ + - ]: 16154 : read(tmp_unit,*, err=10, iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad)
386 : 16154 : ff(irad)=ff(irad)/four_pi
387 : 16154 : ff1(irad)=ff1(irad)/four_pi
388 : 16179 : ff2(irad)=ff2(irad)/four_pi
389 : : end do
390 : :
391 : : !Optional output: VHartree(tild_[n_Z+n_core])
392 [ + + ]: 25 : if (present(vh_tnzc)) then
393 [ - + + - : 2 : ABI_MALLOC(nc,(mmax))
- + - + ]
394 [ - + + + ]: 495 : nc=ff ! n_core
395 : 1 : call psden(1,ff,mmax,nc,rchrg,rad,ff1=ff1,ff2=ff2)
396 : 1 : call vhtnzc(ff,rchrg,vh_tnzc,mmax,rad,znucl)
397 [ - + ]: 1 : ABI_FREE(nc)
398 : : end if
399 : :
400 : 25 : rad(1)=zero
401 : :
402 : : !calculate third derivative ff3 on logarithmic grid
403 : 25 : der1=ff2(1)
404 : 25 : dern=ff2(mmax)
405 : 25 : call spline(rad,ff1,mmax,der1,dern,ff3)
406 : :
407 : : !generate uniform mesh xx in the box cut by rchrg:
408 : :
409 [ + + ]: 62550 : do i1xccc=1,n1xccc
410 : 62550 : xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1)
411 : : end do
412 : :
413 : : !now interpolate core charge and derivatives on the uniform grid
414 : : !core charge, input=ff, output=gg
415 : 25 : call splint(mmax,rad,ff,ff2,n1xccc,xx,gg)
416 : :
417 : : !first derivative input=ff1, output=gg1
418 : 25 : call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1)
419 : :
420 : : !normalize gg1
421 [ + + ]: 62550 : gg1(:)=gg1(:)*rchrg
422 : :
423 : : !now calculate second to fourth derivative by forward differences
424 : : !to avoid numerical noise uses a smoothing function
425 : :
426 : 25 : call smooth(gg1,n1xccc,10)
427 : :
428 : 25 : gg2(n1xccc)=zero
429 [ + + ]: 62525 : do i1xccc=1,n1xccc-1
430 : 62525 : gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1)
431 : : end do
432 : :
433 : 25 : call smooth(gg2,n1xccc,10)
434 : :
435 : 25 : gg3(n1xccc)=zero
436 [ + + ]: 62525 : do i1xccc=1,n1xccc-1
437 : 62525 : gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1)
438 : : end do
439 : :
440 : 25 : call smooth(gg3,n1xccc,10)
441 : :
442 : 25 : gg4(n1xccc)=zero
443 [ + + ]: 62525 : do i1xccc=1,n1xccc-1
444 : 62525 : gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1)
445 : : end do
446 : :
447 : 25 : call smooth(gg4,n1xccc,10)
448 : :
449 : : !write on xcc1d
450 [ + + ]: 62550 : xccc1d(:,1)=gg(:)
451 [ + + ]: 62550 : xccc1d(:,2)=gg1(:)
452 [ + + ]: 62550 : xccc1d(:,3)=gg2(:)
453 [ + + ]: 62550 : xccc1d(:,4)=gg3(:)
454 [ + + ]: 62550 : xccc1d(:,5)=gg4(:)
455 : :
456 : : !WARNING : fifth derivative not yet computed
457 [ + + ]: 62550 : xccc1d(:,6)=zero
458 : :
459 : : !note: the normalization condition is the following:
460 : : !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg
461 : : !
462 : : !norm=zero
463 : : !do i1xccc=1,n1xccc
464 : : !norm = norm + four_pi*rchrg/dble(n1xccc-1)*&
465 : : !& xx(i1xccc)**2*xccc1d(i1xccc,1)
466 : : !end do
467 : : !write(std_out,*) ' norm=',norm
468 : : !
469 : : !write(std_out,*)' psp1cc : output of core charge density and derivatives '
470 : : !write(std_out,*)' xx gg gg1 '
471 : : !do i1xccc=1,n1xccc
472 : : !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
473 : : !end do
474 : : !write(std_out,*)' xx gg2 gg3 '
475 : : !do i1xccc=1,n1xccc
476 : : !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
477 : : !end do
478 : : !write(std_out,*)' xx gg4 gg5 '
479 : : !do i1xccc=1,n1xccc
480 : : !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
481 : : !end do
482 : : !write(std_out,*)' psp1cc : debug done, stop '
483 : : !stop
484 : : !ENDDEBUG
485 : :
486 [ - + ]: 25 : ABI_FREE(ff)
487 [ - + ]: 25 : ABI_FREE(ff1)
488 : 25 : ABI_FREE(ff2)
489 [ - + ]: 25 : ABI_FREE(ff3)
490 [ - + ]: 25 : ABI_FREE(gg)
491 [ - + ]: 25 : ABI_FREE(gg1)
492 : 25 : ABI_FREE(gg2)
493 : 25 : ABI_FREE(gg3)
494 : 25 : ABI_FREE(gg4)
495 : 25 : ABI_FREE(rad)
496 [ - + ]: 25 : ABI_FREE(work)
497 [ - + ]: 25 : ABI_FREE(xx)
498 : :
499 : 25 : return
500 : :
501 : : ! Handle IO error
502 : : 10 continue
503 : 0 : ABI_ERROR(errmsg)
504 : :
505 [ - + - + : 25 : end subroutine psp6cc
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
506 : : !!***
507 : :
508 : : !!****f* m_psp6/psden
509 : : !! NAME
510 : : !! psden
511 : : !!
512 : : !! FUNCTION
513 : : !! Calculate a pseudo-density from an original density on a radial grid (regular or logarithmic)
514 : : !!
515 : : !! INPUTS
516 : : !! ilog=1 if grid is logarithmic, else 0
517 : : !! mesh= dimension of nc
518 : : !! nc(mesh)= density to be pseudized
519 : : !! rc= cut-off radius
520 : : !! rad(mesh) = radial mesh
521 : : !!
522 : : !! OUTPUT
523 : : !! ff(mesh)= pseudized density
524 : : !!
525 : : !!SIDE EFFECTS
526 : : !! Optional:
527 : : !! ff1(mesh)= 1st derivative of pseudo density (only r<rc modified)
528 : : !! ff2(mesh)= 2nd derivative of pseudo density (only r<rc modified)
529 : : !!
530 : : !! NOTES
531 : : !! ff=exp(-(a+b.r^2+c.r^4))
532 : : !!
533 : : !! SOURCE
534 : :
535 : 1 : subroutine psden(ilog,ff,mesh,nc,rc,rad,ff1,ff2)
536 : :
537 : : !Arguments ------------------------------------
538 : : !scalars
539 : : integer,intent(in) :: ilog,mesh
540 : : real(dp),intent(in) :: rc
541 : : !arrays
542 : : real(dp),intent(in) :: nc(mesh),rad(mesh)
543 : : real(dp),intent(out) :: ff(mesh)
544 : : real(dp),intent(inout),optional :: ff1(mesh),ff2(mesh)
545 : :
546 : : !Local variables-------------------------------
547 : : !scalars
548 : : integer :: ii,nc1
549 : : real(dp) :: aa,aa1,aa2,bb,cc,c1,c3,f0,f0p,norm1,norm2,rc1,step
550 : : !arrays
551 : 1 : real(dp),allocatable :: fpir(:),gg(:)
552 : :
553 : : ! *************************************************************************
554 : :
555 : 1 : rc1=rc/four
556 : :
557 [ - + + - : 3 : ABI_MALLOC(fpir,(mesh))
+ - + - -
+ - + ]
558 [ + + ]: 494 : fpir(1:mesh)=four_pi*rad(1:mesh)**2
559 [ + - + + ]: 494 : if (ilog==1) fpir(1:mesh)=fpir(1:mesh)*rad(1:mesh)
560 : :
561 [ - + ]: 1 : if (ilog==0) then
562 : 0 : step=rad(2)-rad(1)
563 : 0 : nc1=int(rc1/step)+1
564 : 0 : rc1=(nc1-1)*step
565 [ + - ]: 1 : else if (ilog==1) then
566 : 1 : step=log(rad(2)/rad(1))
567 : 1 : nc1=int(log(rc1/rad(1))/step)+1
568 : 1 : rc1=rad(nc1)
569 : : end if
570 [ + + ]: 319 : ff(1:nc1)=nc(1:nc1)*fpir(1:nc1)
571 : 1 : call ctrap(nc1,ff(1:nc1),step,c3)
572 [ + - ]: 1 : if (ilog==1) c3=c3+half*ff(1)
573 : 1 : f0=nc(nc1);c1=-log(f0)
574 : 1 : f0p=half*(nc(nc1+1)-nc(nc1-1))/step
575 : :
576 : 1 : ii=0;aa1=zero;norm1=c3+one
577 [ + + - + ]: 2 : do while (norm1>c3.and.ii<100)
578 : 1 : ii=ii+1;aa1=aa1+one
579 : 1 : aa=c1-aa1*rc1**4+rc1*(f0p/f0+four*aa1*rc1**3)*half
580 : 1 : bb=-half*(f0p/f0+four*aa1*rc1**3)/rc1
581 [ + + ]: 319 : ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa1*rad(1:nc1)**4)
582 : 1 : call ctrap(nc1,ff(1:nc1),step,norm1)
583 [ - + ]: 2 : if (ilog==1) norm1=norm1+half*ff(1)
584 : : end do
585 [ - + ]: 1 : if (ii==100) then
586 : 0 : ABI_ERROR('Big pb 1 in psden !')
587 : : end if
588 : :
589 : 1 : ii=0;aa2=zero;norm2=c3-one
590 [ + + - + ]: 2 : do while (norm2<c3.and.ii<100)
591 : 1 : ii=ii+1;aa2=aa2-one
592 : 1 : aa=c1-aa2*rc1**4+rc1*(f0p/f0+four*aa2*rc1**3)*half
593 : 1 : bb=-half*(f0p/f0+four*aa2*rc1**3)/rc1
594 [ + + ]: 319 : ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa2*rad(1:nc1)**4)
595 : 1 : call ctrap(nc1,ff(1:nc1),step,norm2)
596 [ - + ]: 2 : if (ilog==1) norm2=norm2+half*ff(1)
597 : : end do
598 [ - + ]: 1 : if (ii==100) then
599 : 0 : ABI_ERROR('Big pb 2 in psden !')
600 : : end if
601 : :
602 [ + + ]: 36 : do while (abs(norm2-c3)>tol10)
603 : :
604 : 35 : cc=(aa1+aa2)*half
605 : 35 : aa=c1-cc*rc1**4+rc1*(f0p/f0+four*cc*rc1**3)*half
606 : 35 : bb=-half*(f0p/f0+four*cc*rc1**3)/rc1
607 [ + + ]: 11165 : ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-cc*rad(1:nc1)**4)
608 : 35 : call ctrap (nc1,ff(1:nc1),step,norm2)
609 [ + - ]: 35 : if (ilog==1) norm2=norm2+half*ff(1)
610 [ + + ]: 36 : if ((norm1-c3)*(norm2-c3)>zero) then
611 : 19 : aa1=cc
612 : 19 : norm1=norm2
613 : : else
614 : : aa2=cc
615 : : end if
616 : :
617 : : end do ! while
618 : :
619 [ + - ]: 1 : ff(1)=exp(-aa);if (ilog==1) ff(1)=ff(1)*exp(-bb*rad(1)**2-cc*rad(1)**4)
620 [ + + ]: 318 : ff(2:nc1)=ff(2:nc1)/fpir(2:nc1)
621 [ + - + + ]: 176 : if (nc1<mesh) ff(nc1+1:mesh)=nc(nc1+1:mesh)
622 [ + - + + ]: 319 : if (present(ff1)) ff1(1:nc1)=-(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)*ff(1:nc1)
623 [ + - ]: 1 : if (present(ff2)) ff2(1:nc1)=-(two*bb+12.0_dp*cc*rad(1:nc1)**2)*ff(1:nc1) &
624 [ + + ]: 319 : & +(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)**2*ff(1:nc1)
625 : :
626 [ - + + - : 2 : ABI_MALLOC(gg,(mesh))
- + - + ]
627 [ + + ]: 494 : gg(1:mesh)=fpir(1:mesh)*ff(1:mesh)
628 : 1 : call ctrap(mesh,gg(1:mesh),step,norm1)
629 : : if (ilog==1) norm1=norm1+half*gg(1)
630 : : !write(std_out,*) 'psden: tild_nc integral= ',norm1
631 [ - + ]: 1 : ABI_FREE(gg)
632 : :
633 [ - + ]: 1 : ABI_FREE(fpir)
634 : :
635 : 1 : end subroutine psden
636 : : !!***
637 : :
638 : : !!****f* m_psp6/vhtnzc
639 : : !! NAME
640 : : !! vhtnzc
641 : : !!
642 : : !! FUNCTION
643 : : !! Compute VHartree(tild[n_Z+n_core]) from input ncore
644 : : !!
645 : : !! INPUTS
646 : : !! mesh=dimension of radial mesh
647 : : !! nc= core density (to be pseudized)
648 : : !! rad(mesh)=radial mesh
649 : : !! rc=cut-off radius
650 : : !! znucl=nuclear number of atom as specified in psp file
651 : : !!
652 : : !! OUTPUT
653 : : !! vhtnzc(mesh) = hartree potential induced by density tild[n_Z+n_core] (pseudo core density + nucleus)
654 : : !!
655 : : !! SOURCE
656 : :
657 : 1 : subroutine vhtnzc(nc,rc,vh_tnzc,mesh,rad,znucl)
658 : :
659 : : !Arguments ------------------------------------
660 : : !scalars
661 : : integer,intent(in) :: mesh
662 : : real(dp),intent(in) :: znucl
663 : : real(dp),intent(in) :: rc
664 : : !arrays
665 : : real(dp),intent(in) :: nc(mesh),rad(mesh)
666 : : real(dp),intent(out) :: vh_tnzc(mesh)
667 : :
668 : : !Local variables-------------------------------
669 : : !scalars
670 : : integer :: ir,nc1
671 : : real(dp) :: gnorm,rc1,step,yp1,yp2,yp3
672 : : !arrays
673 : 1 : real(dp),allocatable :: den1(:),den2(:),den3(:),den4(:),nzc(:),rvhn(:),shapefunc(:)
674 : :
675 : : ! *************************************************************************
676 : :
677 : 1 : rc1=rc/four
678 : :
679 : 1 : step=log(rad(2)/rad(1))
680 : 1 : nc1=int(log(rc1/rad(1))/step)+1
681 : 1 : rc1=rad(nc1)
682 : :
683 [ - + + - : 3 : ABI_MALLOC(shapefunc,(mesh))
+ - + - -
+ - + ]
684 : 1 : shapefunc(1)=one
685 [ + + ]: 318 : shapefunc(2:nc1)=(sin(pi*rad(2:nc1)/rc1)/(pi*rad(2:nc1)/rc1))**2
686 [ + - + + ]: 176 : if (nc1<mesh) shapefunc(nc1+1:mesh)=zero
687 : :
688 [ - + + - : 2 : ABI_MALLOC(den1,(mesh))
- + - + ]
689 [ + + ]: 494 : den1(1:mesh)=four_pi*shapefunc(1:mesh)*rad(1:mesh)**3
690 : 1 : call ctrap(mesh,den1,step,gnorm)
691 : 1 : gnorm =one/gnorm
692 [ - + ]: 1 : ABI_FREE(den1)
693 : :
694 [ - + + - : 2 : ABI_MALLOC(nzc,(mesh))
- + - + ]
695 [ + + ]: 494 : nzc(1:mesh)=four*pi*nc(1:mesh)*rad(1:mesh)**2-four_pi*shapefunc(1:mesh)*rad(1:mesh)**2*znucl*gnorm
696 : 1 : ABI_FREE(shapefunc)
697 : :
698 [ - + + - : 2 : ABI_MALLOC(rvhn,(mesh))
- + - + ]
699 : 1 : rvhn(1)=zero
700 : :
701 [ - + + - : 2 : ABI_MALLOC(den1,(mesh))
- + - + ]
702 [ - + + - : 2 : ABI_MALLOC(den2,(mesh))
- + - + ]
703 [ - + + - : 2 : ABI_MALLOC(den3,(mesh))
- + - + ]
704 [ - + + - : 2 : ABI_MALLOC(den4,(mesh))
- + - + ]
705 : :
706 : 1 : den1(1)=zero;den2(1)=zero
707 [ + + ]: 493 : do ir=2,mesh
708 : 492 : den1(ir)= rad(ir)*nzc(ir)
709 : 493 : den2(ir)= den1(ir)/rad(ir)
710 : : end do
711 : :
712 : : !For first few points do stupid integral
713 : 1 : den3(1)=zero;den4(1)=zero
714 [ + + ]: 493 : do ir=2,mesh
715 : 492 : call ctrap(ir,den1(1:ir),step,den3(ir))
716 : 493 : call ctrap(ir,den2(1:ir),step,den4(ir))
717 : : end do
718 : :
719 [ + + ]: 494 : do ir=1,mesh
720 : 494 : rvhn(ir)=den3(ir)+rad(ir)*(den4(mesh)-den4(ir))
721 : : end do
722 : :
723 : 1 : ABI_FREE(den1)
724 : 1 : ABI_FREE(den2)
725 : 1 : ABI_FREE(den3)
726 : 1 : ABI_FREE(den4)
727 : :
728 [ + + ]: 493 : vh_tnzc(2:mesh)=rvhn(2:mesh)/rad(2:mesh)
729 : 1 : yp2=(vh_tnzc(3)-vh_tnzc(2))/(rad(3)-rad(2))
730 : 1 : yp3=(vh_tnzc(4)-vh_tnzc(3))/(rad(4)-rad(3))
731 : 1 : yp1=yp2+(yp2-yp3)*rad(2)/(rad(3)-rad(2))
732 : 1 : vh_tnzc(1)=vh_tnzc(2)-(yp1+yp2)*rad(2)
733 : :
734 [ - + ]: 1 : ABI_FREE(nzc)
735 : 1 : ABI_FREE(rvhn)
736 : :
737 : 1 : end subroutine vhtnzc
738 : : !!***
739 : :
740 : : !!****f* m_psp6/psp6cc_drh
741 : : !! NAME
742 : : !! psp6cc_drh
743 : : !!
744 : : !! FUNCTION
745 : : !! Compute the core charge density, for use in the XC core
746 : : !! correction, following the function definition valid
747 : : !! for the format 6 of pseudopotentials.
748 : : !! Version modified by DHamann, with consistent treatment
749 : : !! of the derivatives in this routine and the remaining of the code.
750 : : !!
751 : : !! INPUTS
752 : : !! mmax=maximum number of points in real space grid in the psp file
753 : : !! n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
754 : : !! rchrg=cut-off radius for the core density
755 : : !!
756 : : !! OUTPUT
757 : : !! xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
758 : : !!
759 : : !! NOTES
760 : : !! Test version by DRH - requires very smooth model core charge
761 : : !!
762 : : !! SOURCE
763 : :
764 : 1 : subroutine psp6cc_drh(mmax,n1xccc,rchrg,xccc1d)
765 : :
766 : : !Arguments ------------------------------------
767 : : !scalars
768 : : integer,intent(in) :: mmax,n1xccc
769 : : real(dp),intent(in) :: rchrg
770 : : !arrays
771 : : real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
772 : :
773 : : !Local variables-------------------------------
774 : : !scalars
775 : : integer :: irad
776 : : character(len=500) :: errmsg
777 : : !arrays
778 : 1 : real(dp),allocatable :: ff(:),ff1(:),ff2(:),rad(:)
779 : :
780 : : !**********************************************************************
781 : :
782 [ - + + - : 3 : ABI_MALLOC(ff,(mmax))
+ - + - -
+ - + ]
783 [ - + + - : 2 : ABI_MALLOC(ff1,(mmax))
- + - + ]
784 [ - + + - : 2 : ABI_MALLOC(ff2,(mmax))
- + - + ]
785 [ - + + - : 2 : ABI_MALLOC(rad,(mmax))
- + + - ]
786 : :
787 : : !
788 : : !read from pp file the model core charge (ff) and first (ff1) and
789 : : !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid
790 : : !the input functions contain the 4pi factor, it must be rescaled.
791 : :
792 : : !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc
793 [ + + ]: 601 : do irad=1,mmax
794 [ + - ]: 600 : read(tmp_unit,*,err=10,iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad)
795 : 600 : ff(irad)=ff(irad)/4.d0/pi
796 : 600 : ff1(irad)=ff1(irad)/4.d0/pi
797 : 601 : ff2(irad)=ff2(irad)/4.d0/pi
798 : : end do
799 : 1 : rad(1)=0.d0
800 : :
801 : 1 : call cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d)
802 : :
803 [ - + ]: 1 : ABI_FREE(ff)
804 [ - + ]: 1 : ABI_FREE(ff1)
805 [ - + ]: 1 : ABI_FREE(ff2)
806 : 1 : ABI_FREE(rad)
807 : :
808 : : return
809 : :
810 : : ! Handle IO error
811 : : 10 continue
812 : 0 : ABI_ERROR(errmsg)
813 : :
814 [ # # # # : 0 : end subroutine psp6cc_drh
# # # # ]
815 : : !!***
816 : :
817 : : end module m_psp6
818 : : !!***
|
