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 : !!***
|
