Line data Source code
1 : !!****m* ABINIT/m_pspini
2 : !! NAME
3 : !! m_pspini
4 : !!
5 : !! FUNCTION
6 : !! Initialize pseudopotential datastructures from files.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MT, FrD, AF, DRH, YP)
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_pspini
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_xmpi
28 : use m_psxml2ab
29 : use m_dtset
30 : use m_dtfil
31 : !use m_psps
32 :
33 : use defs_datatypes, only : pseudopotential_type, nctab_t, pspheader_type
34 : use m_fstrings, only : ftoa, sjoin
35 : use m_time, only : timab
36 : use m_io_tools, only : open_file
37 : use m_pawrad, only : pawrad_type
38 : use m_pawtab, only : pawtab_type, pawtab_set_flags
39 : use m_psps, only : psps_print, psps_ncwrite_path, nctab_init, nctab_free, nctab_mixalch, test_xml_xmlpaw_upf, &
40 : nctab_eval_tcorespl
41 : use m_pawpsp, only : pawpsp_bcast, pawpsp_read_pawheader, pawpsp_read_header_xml,&
42 : pawpsp_header_type, pawpsp_wvl, pawpsp_7in, pawpsp_17in
43 : use m_pawxmlps, only : paw_setup_free,paw_setuploc
44 : use m_pspheads, only : pawpsxml2ab
45 : #if defined HAVE_BIGDFT
46 : use BigDFT_API, only : dictionary, atomic_info, dict_init, dict_free, UNINITIALIZED
47 : #endif
48 :
49 : use m_psp1, only : psp1in
50 : use m_psp5, only : psp5in
51 : use m_psp6, only : psp6in
52 : use m_psp8, only : psp8in
53 : use m_psp9, only : psp9in
54 : use m_upf2abinit, only : upf1_to_abinit, upf2_to_abinit
55 : use m_psp_hgh, only : psp2in, psp3in, psp10in
56 : use m_wvl_descr_psp, only : wvl_descr_psp_fill
57 :
58 : implicit none
59 :
60 : private
61 : !!***
62 :
63 : public :: pspini
64 : public :: pspcor
65 : !!***
66 :
67 : contains
68 : !!***
69 :
70 : !!****f* ABINIT/pspini
71 : !! NAME
72 : !! pspini
73 : !!
74 : !! FUNCTION
75 : !! Looping over atom types 1 ... ntypat,
76 : !! read pseudopotential data filename, then call pspatm for each psp.
77 : !! Might combine the psps to generate pseudoatoms, thanks to alchemy.
78 : !! Also compute ecore=[Sum(i) zion(i)] * [Sum(i) epsatm(i)] by calling pspcor.
79 : !!
80 : !! INPUTS
81 : !! dtset <type(dataset_type)>=all input variables in this dataset
82 : !! | iscf=parameter controlling scf or non-scf calculations
83 : !! | ixc=exchange-correlation choice as input to main routine
84 : !! | natom=number of atoms in unit cell
85 : !! | pawxcdev=choice of XC development in PAW formalism
86 : !! | prtvol= control output volume
87 : !! | typat(natom)=type (integer) for each atom
88 : !! | main routine, for each type of atom
89 : !! gsqcut=cutoff for G^2 based on ecut for basis sphere (bohr^-2)
90 : !! gsqcutdg=PAW only - cutoff for G^2 based on ecutdg (fine grid) for basis sphere (bohr^-2)
91 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
92 : !! used to estimate real space mesh (if necessary)
93 : !!
94 : !! OUTPUT
95 : !! ecore=total psp core correction energy*ucvol (hartree*bohr^3)
96 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
97 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
98 : !! gencond=general condition for new computation of pseudopotentials
99 : !! (if gencond=1, new psps have been re-computed)
100 : !!
101 : !! SIDE EFFECTS
102 : !! psps <type(pseudopotential_type)>=at output, psps is completely initialized
103 : !! At the input, it is already partially or completely initialized.
104 : !!
105 : !! NOTES
106 : !! The interplay with the multi-dataset mode is interesting:
107 : !! the pseudopotentials are independent of the dataset, but the largest q vector, the
108 : !! spin-orbit characteristics, the use of Ylm as well as ixc
109 : !! play a role in the set up of pseudopotentials (ixc plays a very minor
110 : !! role, however). So, the pseudopotential data ought not be recomputed
111 : !! when gsqcut, gsqcutdg, mqgrid_ff, mqgrid_vl, npspso, ixc, dimekb and useylm do not change.
112 : !! In many cases, this routine is also called just to write the psp line
113 : !! of the header, without reading again the psp. This psp line
114 : !! is constant throughout run.
115 : !!
116 : !! SOURCE
117 :
118 6935 : subroutine pspini(dtset,dtfil,ecore,gencond,gsqcut,gsqcutdg,pawrad,pawtab,psps,rprimd,comm_mpi)
119 :
120 : !Arguments ------------------------------------
121 : !scalars
122 : integer, optional,intent(in) :: comm_mpi
123 : integer,intent(out) :: gencond
124 : real(dp),intent(in) :: gsqcut,gsqcutdg
125 : real(dp),intent(out) :: ecore
126 : type(dataset_type),intent(in) :: dtset
127 : type(datafiles_type),intent(in) :: dtfil
128 : !arrays
129 : real(dp),intent(in) :: rprimd(3,3)
130 : type(pseudopotential_type), target,intent(inout) :: psps
131 : type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
132 : type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
133 :
134 : !Local variables-------------------------------
135 : !scalars
136 : integer,parameter :: npspmax=50
137 : integer,save :: dimekb_old=0,ifirst=1,ixc_old=-1,lmnmax_old=0,lnmax_old=0,use_rcpaw_old=0
138 : integer,save :: mpssoang_old=0,mqgridff_old=0,mqgridvl_old=0,optnlxccc_old=-1
139 : integer,save :: paw_size_old=-1,pawxcdev_old=-1,positron_old=-2,usekden_old=-1,usepaw_old=-1
140 : integer,save :: usexcnhat_old=-1,usewvl_old=-1,useylm_old=-1
141 : integer :: comm_mpi_,ierr,ii,ilang,ilmn,ilmn0,iln,iproj,ipsp,ipspalch
142 : integer :: ispin,itypalch,itypat,mtypalch,npsp,npspalch,ntypalch
143 : integer :: ntypat,ntyppure,paw_size
144 : logical :: has_coretau,has_kij,has_tproj,has_tvale,has_nabla
145 : logical :: has_shapefncg,has_vminushalf,has_wvl,paw_add_core
146 : real(dp),save :: ecore_old=zero,gsqcut_old=zero,gsqcutdg_old=zero, spnorbscl_old=-one,hyb_mixing_old=-999.0_dp
147 : real(dp) :: dq,epsatm_psp,qmax,rmax,xcccrc
148 : character(len=500) :: msg
149 6935 : type(pawrad_type) :: pawrad_dum
150 6935 : type(pawtab_type) :: pawtab_dum
151 6935 : type(nctab_t) :: nctab_dum
152 : type(nctab_t),pointer :: nctab_ptr
153 : !arrays
154 : integer :: paw_options(11)
155 : integer,save :: paw_options_old(11)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/)
156 : integer,save :: pspso_old(npspmax),pspso_zero(npspmax)
157 6935 : integer,allocatable :: indlmn_alch(:,:,:),new_pspso(:)
158 6935 : integer,pointer :: indlmn(:,:)
159 : real(dp),save :: epsatm(npspmax)
160 : real(dp) :: tsec(2)
161 6935 : real(dp),allocatable :: dvlspl(:,:),dvlspl_alch(:,:,:),ekb(:),ekb_alch(:,:)
162 6935 : real(dp),allocatable :: epsatm_alch(:),ffspl(:,:,:),ffspl_alch(:,:,:,:)
163 6935 : real(dp),allocatable :: vlspl(:,:),vlspl_alch(:,:,:),xccc1d(:,:)
164 6935 : real(dp),allocatable :: xcctau1d(:,:)
165 6935 : real(dp),allocatable :: xccc1d_alch(:,:,:),xcccrc_alch(:)
166 6935 : real(dp),allocatable :: xcctau1d_alch(:,:,:)
167 6935 : type(nctab_t),target,allocatable :: nctab_alch(:)
168 : ! *************************************************************************
169 :
170 : DBG_ENTER("COLL")
171 :
172 : ! Keep track of time spent in this subroutine
173 6935 : call timab(15,1,tsec)
174 :
175 : !-------------------------------------------------------------
176 : ! Some initializations
177 : !-------------------------------------------------------------
178 :
179 : !Useful sizes
180 6935 : ntypat=psps%ntypat
181 6935 : mtypalch=psps%mtypalch
182 6935 : npsp=psps%npsp
183 6935 : if (npsp>npspmax) then
184 0 : ABI_BUG("npsp>npspmax in pspini !")
185 : end if
186 :
187 : ! Set up q grids for atomic data represented in reciprocal space
188 : ! make qmax 20% larger than largest expected:
189 6935 : qmax=1.2d0 * sqrt(gsqcut)
190 : !ffnl is always computed in reciprocal space
191 6935 : dq=qmax/(one*(psps%mqgrid_ff-1))
192 20963353 : do ii=1,psps%mqgrid_ff
193 20963353 : psps%qgrid_ff(ii)=(ii-1)*dq
194 : end do
195 6935 : if (psps%usepaw==1) qmax=1.2d0 * sqrt(gsqcutdg)
196 : !If vlspl is computed in real space, qgrid contains a real space mesh
197 : !the max is taken as the biggest distance in the box.
198 6935 : if (psps%vlspl_recipSpace) then
199 6935 : dq=qmax/(one*(psps%mqgrid_vl-1))
200 : else
201 0 : rmax = (rprimd(1, 1) + rprimd(1, 2) + rprimd(1, 3)) ** 2
202 0 : rmax = rmax + (rprimd(2, 1) + rprimd(2, 2) + rprimd(2, 3)) ** 2
203 0 : rmax = rmax + (rprimd(3, 1) + rprimd(3, 2) + rprimd(3, 3)) ** 2
204 0 : rmax = sqrt(rmax)
205 0 : dq = rmax /(one*(psps%mqgrid_vl-1))
206 : end if
207 21000524 : do ii=1,psps%mqgrid_vl
208 21000524 : psps%qgrid_vl(ii)=(ii-1)*dq
209 : end do
210 :
211 : !Determine whether new optional data requests have changed
212 6935 : paw_options=0;paw_size=0
213 6935 : if (psps%usepaw==1) then
214 1440 : paw_size=size(pawtab)
215 1440 : has_kij=(dtset%positron/=0.or.abs(dtset%effmass_free-one)>tol8.or.dtset%orbmag>0.or.dtset%use_rcpaw==1)
216 1440 : has_tvale=.true. ! Will be modified later (depending on PAW dataset format)
217 1440 : has_nabla=.false.
218 1440 : has_shapefncg=(dtset%optdriver==RUNL_GSTATE.and.((dtset%iprcel>=20.and.dtset%iprcel<70).or.dtset%iprcel>=80))
219 1440 : has_wvl=(dtset%usewvl==1.or.dtset%icoulomb/=0)
220 1440 : has_tproj=(dtset%usewvl==1.or.dtset%use_rcpaw==1) ! projectors will be free at the end of the psp reading
221 3502 : has_vminushalf=(maxval(dtset%ldaminushalf)==1)
222 1440 : has_coretau=(dtset%usekden>=1)
223 1440 : paw_add_core=(dtset%paw_add_core==1)
224 1440 : if (has_kij) paw_options(1)=1
225 1440 : if (has_tvale) paw_options(2)=1
226 : if (has_nabla) paw_options(5)=1
227 1440 : if (has_shapefncg) paw_options(6)=1
228 1440 : if (has_wvl) paw_options(7)=1
229 1440 : if (has_tproj) paw_options(8)=1
230 1440 : if (has_vminushalf)paw_options(9)=1
231 1440 : if (has_coretau) paw_options(10)=1
232 1440 : if (paw_add_core) paw_options(11)=1
233 : !if (dtset%prtvclmb /= 0) then
234 1440 : paw_options(3) = 1
235 1440 : paw_options(4) = 1
236 : !end if
237 : end if
238 :
239 : !Determine whether the spin-orbit characteristic has changed
240 : !Do not forget that the SO is not consistent with alchemy presently
241 20805 : ABI_MALLOC(new_pspso,(npsp))
242 76485 : if (ifirst==1) pspso_old(:)=-1
243 76485 : if (ifirst==1) pspso_zero(:)=-1
244 16764 : do ipsp=1,npsp
245 9829 : new_pspso(ipsp)=1
246 : ! No new characteristics if it is equal to the old one,
247 : ! or, if it is one, the old one is equal to the intrinsic characteristic one.
248 9829 : if (psps%pspso(ipsp)==pspso_old(ipsp).or. &
249 : & (psps%pspso(ipsp)==1.and.pspso_old(ipsp)==pspso_zero(ipsp))) then
250 7898 : new_pspso(ipsp)=0
251 : end if
252 : ! No new characteristics if PAW
253 9829 : if (psps%usepaw==1) new_pspso(ipsp)=0
254 : ! Prepare the saving of the intrinsic pseudopotential characteristics
255 16764 : if(psps%pspso(ipsp)==1) pspso_zero(ipsp)=0
256 : end do
257 :
258 : !Compute the general condition for new computation of pseudopotentials
259 6935 : gencond=0
260 : if( ixc_old /= dtset%ixc &
261 : & .or. ABS(hyb_mixing_old - dtset%hyb_mixing)>tol8 &
262 : & .or. mqgridff_old /= psps%mqgrid_ff &
263 : & .or. mqgridvl_old /= psps%mqgrid_vl &
264 : & .or. mpssoang_old /= psps%mpssoang &
265 : & .or. abs(gsqcut_old-gsqcut)>1.0d-10 &
266 : & .or. (psps%usepaw==1.and.abs(gsqcutdg_old-gsqcutdg)>1.0d-10) &
267 : & .or. dimekb_old /= psps%dimekb &
268 : & .or. lmnmax_old /= psps%lmnmax &
269 : & .or. lnmax_old /= psps%lnmax &
270 : & .or. optnlxccc_old /= psps%optnlxccc &
271 : & .or. usepaw_old /= psps%usepaw &
272 : & .or. useylm_old /= psps%useylm &
273 : & .or. pawxcdev_old /= dtset%pawxcdev &
274 : & .or. (psps%usepaw == 0 .and. spnorbscl_old /= dtset%spnorbscl) &
275 : & .or. positron_old /= dtset%positron &
276 : & .or. usewvl_old /= dtset%usewvl &
277 : & .or. paw_size_old /= paw_size &
278 : & .or. usekden_old/=dtset%usekden &
279 : & .or. usexcnhat_old/=dtset%usexcnhat_orig &
280 : & .or. any(paw_options_old(:)/=paw_options(:)) &
281 : & .or. sum(new_pspso(:))/=0 &
282 : & .or. mtypalch>0 &
283 : & .or. (dtset%usewvl==1.and.psps%usepaw==1)&
284 84427 : & .or. (use_rcpaw_old==1) &
285 2261 : & ) gencond=1
286 :
287 6935 : if (present(comm_mpi).and.psps%usepaw==1) then
288 1440 : if(xmpi_comm_size(comm_mpi)>1) call xmpi_sum(gencond,comm_mpi,ierr)
289 1440 : if (gencond/=0) gencond=1
290 : end if
291 6935 : ABI_FREE(new_pspso)
292 :
293 : !-------------------------------------------------------------
294 : ! Following section is only reached when new computation
295 : ! of pseudopotentials is needed
296 : !-------------------------------------------------------------
297 :
298 6935 : if (gencond==1) then
299 :
300 2261 : write(msg, '(a,a)' ) ch10,'--- Pseudopotential description ------------------------------------------------'
301 2261 : call wrtout(ab_out,msg)
302 :
303 6783 : ABI_MALLOC(ekb,(psps%dimekb*(1-psps%usepaw)))
304 6783 : ABI_MALLOC(xccc1d,(psps%n1xccc*(1-psps%usepaw),6))
305 4522 : ABI_MALLOC(xcctau1d,(psps%n1xccc*(1-psps%usepaw),6))
306 9044 : ABI_MALLOC(ffspl,(psps%mqgrid_ff,2,psps%lnmax))
307 6783 : ABI_MALLOC(vlspl,(psps%mqgrid_vl,2))
308 2261 : if (.not.psps%vlspl_recipSpace) then
309 0 : ABI_MALLOC(dvlspl,(psps%mqgrid_vl,2))
310 : else
311 2261 : ABI_MALLOC(dvlspl,(0,0))
312 : end if
313 :
314 : ! PAW: reset flags for optional data
315 2261 : if (psps%usepaw==1) then
316 : call pawtab_set_flags(pawtab,has_kij=paw_options(1),has_tvale=paw_options(2),&
317 : & has_vhnzc=paw_options(3),has_vhtnzc=paw_options(4),&
318 : & has_nabla=paw_options(5),has_shapefncg=paw_options(6),&
319 : & has_wvl=paw_options(7),has_tproj=paw_options(8),&
320 : & has_vminushalf=paw_options(9),has_coretau=paw_options(10),&
321 522 : & add_core_energy=paw_options(11))
322 : end if
323 :
324 : ! Read atomic pseudopotential data and get transforms
325 : ! for each atom type: two cases, alchemy or not.
326 :
327 2261 : if(mtypalch==0)then
328 : ! No alchemical pseudoatom, in all datasets, npsp=ntypat
329 :
330 5090 : do ipsp=1,npsp
331 :
332 2855 : xcccrc=zero
333 92496844 : ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero
334 2855 : if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero
335 2855 : if (psps%usepaw==0) then
336 16986577 : xccc1d(:,:)=zero
337 16986577 : xcctau1d(:,:)=zero
338 : end if
339 2855 : indlmn=>psps%indlmn(:,:,ipsp)
340 134392 : indlmn(:,:)=0
341 :
342 2855 : write(msg, '(a,i4,a,t38,a)' )'- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp))
343 8565 : call wrtout([std_out, ab_out], msg)
344 :
345 : ! Read atomic psp V(r) and wf(r) to get local and nonlocal psp:
346 : ! Cannot use the same call in case of bound checking, because of pawrad/pawtab
347 2855 : if(psps%usepaw==0)then
348 : call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,&
349 2113 : pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,psps%nctab(ipsp))
350 11019 : psps%ekb(:,ipsp)=ekb(:)
351 16986577 : psps%xccc1d(:,:,ipsp)=xccc1d(:,:)
352 16986577 : psps%xcctau1d(:,:,ipsp)=xcctau1d(:,:)
353 : else
354 742 : comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi
355 : call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,&
356 742 : pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,nctab_dum,comm_mpi=comm_mpi_)
357 742 : if (dtset%usefock==1.and.pawtab(ipsp)%has_fock==0) then
358 0 : ABI_BUG('The PAW data file does not contain Fock information. Change the PAW data file!')
359 : end if
360 : end if
361 :
362 : ! Copy data to psps datastructure.
363 2855 : psps%xcccrc(ipsp)=xcccrc
364 2855 : psps%znucltypat(ipsp)=psps%znuclpsp(ipsp)
365 75043638 : psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:)
366 17447155 : psps%vlspl(:,:,ipsp)=vlspl(:,:)
367 5090 : if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp) = dvlspl(:, :)
368 : end do ! ipsp
369 :
370 : else
371 : ! mtypalch/=0
372 :
373 26 : npspalch=psps%npspalch
374 26 : ntyppure=npsp-npspalch
375 26 : ntypalch=psps%ntypalch
376 78 : ABI_MALLOC(epsatm_alch,(npspalch))
377 104 : ABI_MALLOC(ekb_alch,(psps%dimekb,npspalch*(1-psps%usepaw)))
378 130 : ABI_MALLOC(ffspl_alch,(psps%mqgrid_ff,2,psps%lnmax,npspalch))
379 104 : ABI_MALLOC(xccc1d_alch,(psps%n1xccc*(1-psps%usepaw),6,npspalch))
380 78 : ABI_MALLOC(xcctau1d_alch,(psps%n1xccc*(1-psps%usepaw),6,npspalch))
381 52 : ABI_MALLOC(xcccrc_alch,(npspalch))
382 104 : ABI_MALLOC(vlspl_alch,(psps%mqgrid_vl,2,npspalch))
383 26 : if (.not.psps%vlspl_recipSpace) then
384 0 : ABI_MALLOC(dvlspl_alch,(psps%mqgrid_vl,2,npspalch))
385 : end if
386 78 : ABI_MALLOC(indlmn,(6,psps%lmnmax))
387 104 : ABI_MALLOC(indlmn_alch,(6,psps%lmnmax,npspalch))
388 :
389 : ! Allocate NC tables used for mixing.
390 26 : if (psps%usepaw == 0) then
391 130 : ABI_MALLOC(nctab_alch, (npspalch))
392 78 : do ipspalch=1,npspalch
393 78 : call nctab_init(nctab_alch(ipspalch), psps%mqgrid_vl, .False., .False.)
394 : end do
395 : end if
396 :
397 94 : do ipsp=1,npsp
398 68 : write(msg, '(a,i4,a,t38,a)' ) '- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp))
399 68 : call wrtout(ab_out,msg)
400 :
401 68 : xcccrc=zero
402 2768698 : ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero
403 68 : if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero
404 68 : if (psps%usepaw==0) then
405 180548 : xccc1d(:,:)=zero
406 180548 : xcctau1d(:,:)=zero
407 : end if
408 2819 : indlmn(:,:)=0
409 :
410 : ! Read atomic psp V(r) and wf(r) to get local and nonlocal psp:
411 68 : if (psps%usepaw==0) then
412 68 : if (ipsp <= ntyppure) then
413 : ! Store data in nctab if pure atom.
414 16 : nctab_ptr => psps%nctab(ipsp)
415 : else
416 : ! Store data in nctab_alch (to be mixed afterwards).
417 52 : nctab_ptr => nctab_alch(ipsp-ntyppure)
418 : end if
419 :
420 : call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,&
421 68 : & pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,nctab_ptr)
422 :
423 0 : else if (psps%usepaw==1) then
424 0 : comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi
425 : call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,&
426 : & pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,nctab_dum,&
427 0 : & comm_mpi=comm_mpi_)
428 : end if
429 :
430 94 : if (ipsp<=ntyppure) then
431 : ! Pure pseudopotentials, leading to pure pseudoatoms
432 16 : epsatm(ipsp)=epsatm_psp
433 16 : psps%znucltypat(ipsp)=psps%znuclpsp(ipsp)
434 129 : if (psps%usepaw==0) psps%ekb(:,ipsp)=ekb(:)
435 678581 : psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:)
436 96080 : psps%vlspl(:,:,ipsp)=vlspl(:,:)
437 16 : if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp)=dvlspl(:, :)
438 16 : if (psps%usepaw==0) then
439 60136 : psps%xccc1d(:,:,ipsp)=xccc1d(:,:)
440 60136 : psps%xcctau1d(:,:,ipsp)=xcctau1d(:,:)
441 : end if
442 16 : psps%xcccrc(ipsp)=xcccrc
443 1598 : psps%indlmn(:,:,ipsp)=indlmn(:,:)
444 :
445 : else
446 : ! Pseudopotentials for alchemical generation
447 52 : ipspalch=ipsp-ntyppure
448 52 : epsatm_alch(ipspalch)=epsatm_psp
449 1681452 : ffspl_alch(:,:,:,ipspalch)=ffspl(:,:,:)
450 312260 : vlspl_alch(:,:,ipspalch)=vlspl(:,:)
451 52 : if (.not.psps%vlspl_recipSpace) dvlspl_alch(:,:,ipspalch)=dvlspl(:,:)
452 52 : if (psps%usepaw==0) then
453 332 : ekb_alch(:,ipspalch)=ekb(:)
454 120412 : xccc1d_alch(:,:,ipspalch)=xccc1d(:,:)
455 120412 : xcctau1d_alch(:,:,ipspalch)=xcctau1d(:,:)
456 : end if
457 52 : xcccrc_alch(ipspalch)=xcccrc
458 2012 : indlmn_alch(:,:,ipspalch)=indlmn(:,:)
459 : ! write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,1,ipspalch)=',indlmn_alch(:,1,ipspalch)
460 : ! write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,2,ipspalch)=',indlmn_alch(:,2,ipspalch)
461 : end if
462 :
463 : end do ! ipsp
464 :
465 : ! Generate data for alchemical pseudos.
466 52 : do itypalch=1,ntypalch
467 26 : itypat=itypalch+ntyppure
468 26 : psps%znucltypat(itypat)=200.0+itypalch ! Convention for alchemical pseudoatoms
469 156130 : vlspl(:,:)=zero
470 26 : if (.not.psps%vlspl_recipSpace) dvlspl(:, :) = zero
471 26 : epsatm(itypat)=zero
472 26 : xcccrc=zero
473 26 : if (psps%usepaw==0) then
474 60206 : xccc1d(:,:)=zero
475 60206 : xcctau1d(:,:)=zero
476 : end if
477 :
478 : ! Here, linear combination of the quantities
479 : ! MG: FIXME I think that the mixing of xcccrc is wrong when the xxccrc are different!
480 : ! but this is minor bug since alchemical pseudos should not have XCCC (?)
481 78 : do ipspalch=1,npspalch
482 52 : epsatm(itypat) = epsatm(itypat) + epsatm_alch(ipspalch) * psps%mixalch(ipspalch,itypalch)
483 312260 : vlspl(:,:) = vlspl(:,:) + vlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
484 52 : if (.not.psps%vlspl_recipSpace) then
485 0 : dvlspl(:,:) = dvlspl(:,:) + dvlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
486 : end if
487 52 : xcccrc = xcccrc + xcccrc_alch(ipspalch) * psps%mixalch(ipspalch,itypalch)
488 78 : if (psps%usepaw==0) then
489 120412 : xccc1d(:,:) = xccc1d(:,:) + xccc1d_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
490 120412 : xcctau1d(:,:) = xcctau1d(:,:) + xcctau1d_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
491 : end if
492 : end do ! ipspalch
493 :
494 156130 : psps%vlspl(:,:,itypat)=vlspl(:,:)
495 26 : if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, itypat) = dvlspl(:, :)
496 26 : if (psps%usepaw==0) then
497 60206 : psps%xccc1d(:,:,itypat)=xccc1d(:,:)
498 60206 : psps%xcctau1d(:,:,itypat)=xcctau1d(:,:)
499 : end if
500 26 : psps%xcccrc(itypat)=xcccrc
501 :
502 26 : if (abs(xcccrc) > tol6) then
503 4 : write(std_out, *)"xcccrc", xcccrc
504 4 : ABI_WARNING("Alchemical pseudopotential with nlcc!")
505 : end if
506 :
507 : ! Combine the different non-local projectors : for the scalar part then
508 : ! the spin-orbit part, treat the different angular momenta
509 : ! WARNING : this coding does not work for PAW
510 26 : ilmn=0; iln=0
511 1006 : psps%indlmn(:,:,itypat)=0
512 104 : do ispin=1,2
513 286 : do ilang=0,3
514 208 : if(ispin==2 .and. ilang==0)cycle
515 : iproj=0
516 598 : do ipspalch=1,npspalch
517 572 : if(abs(psps%mixalch(ipspalch,itypalch))>tol10)then
518 2002 : do ilmn0=1,psps%lmnmax
519 2002 : if(indlmn_alch(5,ilmn0,ipspalch)/=0)then
520 805 : if(indlmn_alch(6,ilmn0,ipspalch)==ispin)then
521 448 : if(indlmn_alch(1,ilmn0,ipspalch)==ilang)then
522 115 : ilmn=ilmn+1 ! increment the counter
523 115 : if (indlmn_alch(2,ilmn0,ipspalch)==-ilang*psps%useylm)then
524 115 : iln = iln+1
525 115 : iproj = iproj+1
526 : end if
527 115 : if(ilmn>psps%lmnmax)then
528 0 : ABI_BUG('Problem with the alchemical pseudopotentials : ilmn>lmnmax.')
529 : end if
530 115 : psps%indlmn(1,ilmn,itypat)=ilang
531 115 : psps%indlmn(2,ilmn,itypat)=indlmn_alch(2,ilmn0,ipspalch)
532 115 : psps%indlmn(3,ilmn,itypat)=iproj
533 115 : psps%indlmn(4,ilmn,itypat)=ilmn
534 115 : psps%indlmn(5,ilmn,itypat)=iln
535 115 : psps%indlmn(6,ilmn,itypat)=ispin
536 : ! The two lines below do not work for PAW
537 115 : if (psps%usepaw==0) then
538 115 : psps%ekb(iln,itypat)=psps%mixalch(ipspalch,itypalch) *ekb_alch(indlmn_alch(5,ilmn0,ipspalch),ipspalch)
539 : end if
540 690575 : psps%ffspl(:,:,iln,itypat)=ffspl_alch(:,:,indlmn_alch(5,ilmn0,ipspalch),ipspalch)
541 :
542 115 : psps%indlmn(1,ilmn,itypat)=ilang
543 115 : psps%indlmn(2,ilmn,itypat)=indlmn_alch(2,ilmn0,ipspalch)
544 115 : psps%indlmn(3,ilmn,itypat)=iproj ! This does not work for PAW
545 115 : psps%indlmn(4,ilmn,itypat)=indlmn_alch(4,ilmn0,ipspalch) ! This does not work for PAW
546 115 : psps%indlmn(5,ilmn,itypat)=iln
547 115 : psps%indlmn(6,ilmn,itypat)=ispin
548 : end if ! ilang is OK
549 : end if ! ispin is OK
550 : end if ! ilmn0 exist
551 : end do ! ilmn0
552 : end if ! mixalch>tol10
553 : end do ! ipspalch
554 : end do ! ilang
555 : end do ! ispin
556 :
557 : end do ! itypalch
558 :
559 26 : ABI_FREE(epsatm_alch)
560 26 : ABI_FREE(ekb_alch)
561 26 : ABI_FREE(ffspl_alch)
562 26 : ABI_FREE(xccc1d_alch)
563 26 : ABI_FREE(xcctau1d_alch)
564 26 : ABI_FREE(xcccrc_alch)
565 26 : ABI_FREE(vlspl_alch)
566 26 : if (.not.psps%vlspl_recipSpace) then
567 0 : ABI_FREE(dvlspl_alch)
568 : end if
569 26 : ABI_FREE(indlmn_alch)
570 26 : ABI_FREE(indlmn)
571 :
572 : ! Mix NC tables.
573 26 : if (psps%usepaw == 0) then
574 26 : call nctab_mixalch(nctab_alch, npspalch, ntypalch, psps%algalch, psps%mixalch, psps%nctab(ntyppure+1:))
575 78 : do ipspalch=1,npspalch
576 78 : call nctab_free(nctab_alch(ipspalch))
577 : end do
578 78 : ABI_FREE(nctab_alch)
579 : end if
580 : end if ! mtypalch
581 :
582 2261 : ABI_FREE(ekb)
583 2261 : ABI_FREE(ffspl)
584 2261 : ABI_FREE(vlspl)
585 2261 : ABI_FREE(xccc1d)
586 2261 : ABI_FREE(xcctau1d)
587 2261 : ABI_FREE(dvlspl)
588 : end if ! End condition of new computation needed
589 :
590 : !-------------------------------------------------------------
591 : ! Following section is always executed
592 : !-------------------------------------------------------------
593 : !One should move this section of code outside of pspini,
594 : !but epsatm is needed, so should be in the psp datastructure.
595 : !Compute pseudo correction energy. Will differ from an already
596 : !computed one if the number of atom differ ...
597 16764 : do ipsp=1,npsp
598 16764 : psps%epsatm(ipsp)=epsatm(ipsp)
599 : enddo
600 6935 : call pspcor(ecore,epsatm,dtset%natom,ntypat,dtset%typat,psps%ziontypat)
601 6935 : if(abs(ecore_old-ecore)>tol8*abs(ecore_old+ecore))then
602 1466 : write(msg, '(2x,es15.8,t50,a)' ) ecore,'ecore*ucvol(ha*bohr**3)'
603 : ! ecore is useless if iscf<=0, but at least it has been initialized
604 1466 : if(dtset%iscf>=0) call wrtout(ab_out,msg)
605 1466 : call wrtout(std_out,msg)
606 : end if
607 :
608 : !End of pseudopotential output section
609 6935 : write(msg, '(2a)' )'--------------------------------------------------------------------------------',ch10
610 6935 : call wrtout(ab_out,msg)
611 :
612 : !-------------------------------------------------------------
613 : ! Keep track of this call to the routine
614 : !-------------------------------------------------------------
615 :
616 6935 : if (ifirst==1) ifirst=0
617 :
618 6935 : mqgridff_old=psps%mqgrid_ff
619 6935 : mqgridvl_old=psps%mqgrid_vl
620 6935 : mpssoang_old=psps%mpssoang
621 6935 : spnorbscl_old = dtset%spnorbscl
622 6935 : ixc_old=dtset%ixc
623 6935 : hyb_mixing_old=dtset%hyb_mixing
624 6935 : gsqcut_old=gsqcut;if (psps%usepaw==1) gsqcutdg_old=gsqcutdg
625 6935 : lmnmax_old=psps%lmnmax
626 6935 : lnmax_old=psps%lnmax
627 6935 : optnlxccc_old=psps%optnlxccc
628 6935 : usepaw_old=psps%usepaw
629 6935 : dimekb_old=psps%dimekb
630 6935 : useylm_old=psps%useylm
631 6935 : pawxcdev_old=dtset%pawxcdev
632 6935 : positron_old=dtset%positron
633 6935 : usewvl_old = dtset%usewvl
634 6935 : usekden_old = dtset%usekden
635 6935 : usexcnhat_old=dtset%usexcnhat_orig
636 6935 : paw_size_old=paw_size
637 6935 : use_rcpaw_old=dtset%use_rcpaw
638 6935 : ecore_old=ecore
639 6935 : paw_options_old(:)=paw_options(:)
640 :
641 16764 : do ipsp=1,npsp
642 9829 : pspso_old(ipsp)=psps%pspso(ipsp)
643 16764 : if(pspso_zero(ipsp)==0)pspso_zero(ipsp)=psps%pspso(ipsp)
644 : end do
645 62535 : psps%mproj = maxval(psps%indlmn(3,:,:))
646 :
647 9196 : if (gencond == 1) call psps_print(psps,[std_out], prtvol=dtset%prtvol)
648 :
649 : ! Write the PSPS.nc file and exit here if requested by the user.
650 6935 : if (abs(dtset%prtpsps) == 1) then
651 1 : if (xmpi_comm_rank(xmpi_world) == 0) call psps_ncwrite_path(psps, trim(dtfil%filnam_ds(4))//"_PSPS.nc")
652 1 : if (dtset%prtpsps == -1) then
653 0 : ABI_ERROR_NODUMP("prtpsps == -1 ==> aborting now")
654 : end if
655 : end if
656 :
657 6935 : call timab(15,2,tsec)
658 :
659 : DBG_EXIT("COLL")
660 :
661 20805 : end subroutine pspini
662 : !!***
663 :
664 : !!****f* ABINIT/pspcor
665 : !! NAME
666 : !! pspcor
667 : !!
668 : !! FUNCTION
669 : !! Compute ecore pseudoion-pseudoion correction energy from epsatm for
670 : !! different types of atoms in unit cell.
671 : !!
672 : !! INPUTS
673 : !! natom=number of atoms in cell
674 : !! ntypat=number of types of atoms
675 : !! typat(natom)=integer label of 'typat' for each atom in cell
676 : !! epsatm(ntypat)=pseudoatom energy for each type of atom
677 : !! zion(ntypat)=valence charge on each type of atom in cell
678 : !!
679 : !! OUTPUT
680 : !! ecore=resulting psion-psion energy in Hartrees
681 : !!
682 : !! SOURCE
683 :
684 6953 : subroutine pspcor(ecore,epsatm,natom,ntypat,typat,zion)
685 :
686 : !Arguments ------------------------------------
687 : !scalars
688 : integer,intent(in) :: natom,ntypat
689 : real(dp),intent(out) :: ecore
690 : !arrays
691 : integer,intent(in) :: typat(natom)
692 : real(dp),intent(in) :: epsatm(ntypat),zion(ntypat)
693 :
694 : !Local variables-------------------------------
695 : integer :: ia
696 : real(dp) :: charge,esum
697 : ! *************************************************************************
698 :
699 6953 : charge = 0.d0
700 6953 : esum = 0.d0
701 23101 : do ia=1,natom
702 : ! compute pseudocharge:
703 16148 : charge=charge+zion(typat(ia))
704 : ! add pseudocore energies together:
705 23101 : esum = esum + epsatm(typat(ia))
706 : end do
707 :
708 6953 : ecore=charge*esum
709 :
710 6953 : end subroutine pspcor
711 : !!***
712 :
713 : !!****f* ABINIT/pspatm
714 : !! NAME
715 : !! pspatm
716 : !!
717 : !! FUNCTION
718 : !! Open atomic pseudopotential data file for a given atom,
719 : !! read the three first lines, make some checks, then
720 : !! call appropriate subroutine for the reading of
721 : !! V(r) and wf R(r) data for each angular momentum, and subsequent
722 : !! Fourier and Bessel function transforms for local and nonlocal potentials.
723 : !! Close psp file at end.
724 : !!
725 : !! Handles pseudopotential files produced by (pspcod=1 or 4) Teter code,
726 : !! or from the Goedecker-Teter-Hutter paper (pspcod=2),
727 : !! or from the Hartwigsen-Goedecker-Hutter paper (pspcod=3 or 10)
728 : !! or "Phoney pseudopotentials" (Hamman grid in real space) (pspcod=5)
729 : !! or "Troullier-Martins pseudopotentials" from the FHI (pspcod=6)
730 : !! or "XML format" (pspcod=9)
731 : !! or "UPF1 PWSCF format" (pspcod=11)
732 : !! or "UPF2 PWSCF format" (pspcod=12)
733 : !!
734 : !! INPUTS
735 : !! dq= spacing of the q-grid
736 : !! dtset <type(dataset_type)>=all input variables in this dataset
737 : !! | ixc=exchange-correlation choice from main routine data file
738 : !! | pawxcdev=choice of XC development in PAW formalism
739 : !! | usexcnhat_orig=choice for use of nhat in Vxc in PAW formalism
740 : !! | xclevel= XC functional level
741 : !! ipsp=id in the array of the currently read pseudo.
742 : !!
743 : !! OUTPUT
744 : !! ekb(dimekb)=
745 : !! ->NORM-CONSERVING PSPS ONLY (pspcod/=7):
746 : !! (Real) Kleinman-Bylander energies (hartree)
747 : !! {{\ \begin{equation}
748 : !! \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]}
749 : !! {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]}
750 : !! \end{equation} }}
751 : !! for number of basis functions (l,n) (dimekb=lnmax)
752 : !! If any, spin-orbit components begin at l=mpsang+1
753 : !! epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
754 : !! indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
755 : !! or i=lmn (if useylm=1)
756 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data
757 : !! pawtab <type(pawtab_type)>=paw tabulated starting data
758 : !! vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
759 : !! ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
760 : !! second derivative from spline fit for each angular momentum and
761 : !! each projector; if any, spin-orbit components begin at l=mpsang+1
762 : !! xcccrc=XC core correction cutoff radius (bohr) from psp file
763 : !! xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives, from psp file (used in NC only)
764 : !! xcctau1d(n1xccc*(1-usepaw),6)=1D core charge kinetic energy, and five derivatives, from psp file (used in NC only)
765 : !! nctab=<nctab_t>
766 : !! has_tvale=True if the pseudo provides the valence density (used in NC only)
767 : !! tvalespl(mqgrid_vl(1-usepaw),2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
768 : !! (used in NC only)
769 : !!
770 : !! SIDE EFFECTS
771 : !! Input/Output :
772 : !! psps <type(pseudopotential_type)>=at output, values depending on the read
773 : !! pseudo are set.
774 : !! | dimekb(IN)=dimension of ekb (see module defs_psp.f)
775 : !! | filpsp(IN)=name of formatted external file containing atomic psp data.
776 : !! | lmnmax(IN)=if useylm=1, max number of (l,m,n) comp. over all type of psps
777 : !! | =if useylm=0, max number of (l,n) comp. over all type of psps
778 : !! | lnmax(IN)=max. number of (l,n) components over all type of psps
779 : !! | angular momentum of nonlocal pseudopotential
780 : !! | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
781 : !! | mpssoang(IN)= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
782 : !! | mqgrid_ff(IN)=dimension of q (or G) grid for nl form factors (array ffspl)
783 : !! | mqgrid_vl(IN)=dimension of q (or G) grid for Vloc (array vlspl)
784 : !! | n1xccc(IN)=dimension of xccc1d ; 0 if no XC core correction is used
785 : !! | optnlxccc(IN)=option for nl XC core correction
786 : !! | positron(IN)=0 if electron GS calculation
787 : !! | 1 if positron GS calculation
788 : !! | 2 if electron GS calculation in presence of the positron
789 : !! | pspso(INOUT)=spin-orbit characteristics, govern the content of ffspl and ekb
790 : !! | if =0 : this input requires NO spin-orbit characteristics of the psp
791 : !! | if =2 : this input requires HGH characteristics of the psp
792 : !! | if =3 : this input requires HFN characteristics of the psp
793 : !! | if =1 : this input will be changed at output to 1, 2, 3, according
794 : !! | to the intrinsic characteristics of the psp file
795 : !! | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
796 : !! | qgrid_vl(mqgrid_vl)(IN)=values of q on grid from 0 to qmax (bohr^-1) for Vloc
797 : !! | usepaw(IN)= 0 for non paw calculation; =1 for paw calculation
798 : !! | useylm(IN)=governs the way the nonlocal operator is to be applied:
799 : !! | 1=using Ylm, 0=using Legendre polynomials
800 : !! | vlspl_recipSpace(IN)=.true. if pseudo are expressed in reciprocal space.
801 : !! | znuclpsp(IN)=atomic number of atom as specified in input file to main routine
802 : !!
803 : !! NOTES
804 : !! Format expected for the three first lines of pseudopotentials
805 : !! (1) title (character) line
806 : !! (2) znucl,zion,pspdat
807 : !! (3) pspcod,pspxc,lmax,lloc,mmax,r2well
808 : !!
809 : !! Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case
810 : !!
811 : !! SOURCE
812 :
813 2923 : subroutine pspatm(dq,dtset,dtfil,ekb,epsatm,ffspl,indlmn,ipsp,pawrad,pawtab,&
814 2923 : psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,nctab,comm_mpi)
815 :
816 : !Arguments ---------------------------------------------
817 : !scalars
818 : integer,intent(in) :: ipsp
819 : integer, optional,intent(in) :: comm_mpi
820 : real(dp),intent(in) :: dq
821 : real(dp),intent(out) :: epsatm,xcccrc
822 : type(dataset_type),intent(in) :: dtset
823 : type(datafiles_type),intent(in) :: dtfil
824 : type(pawrad_type),intent(inout) :: pawrad
825 : type(pawtab_type),intent(inout) :: pawtab
826 : type(nctab_t),intent(inout) :: nctab
827 : type(pseudopotential_type),intent(inout) :: psps
828 : !arrays
829 : integer,intent(out) :: indlmn(6,psps%lmnmax)
830 : real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2)
831 : real(dp),intent(inout) :: ekb(psps%dimekb*(1-psps%usepaw))
832 : real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax)
833 : real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2)
834 : real(dp),intent(inout) :: xccc1d(psps%n1xccc*(1-psps%usepaw),6)
835 : real(dp),intent(inout) :: xcctau1d(psps%n1xccc*(1-psps%usepaw),6)
836 :
837 : !Local variables ---------------------------------------
838 : !scalars
839 : integer :: ii,il,ilmn,iln,iln0,lloc,lmax,me,mmax
840 : integer :: paral_mode,pspcod,pspdat,pspxc,useupf,usexml,xmlpaw,unt
841 : real(dp) :: maxrad,qchrg,r2well,zion,znucl,el_temp
842 : logical,parameter :: nc_debug = .False.
843 : !logical,parameter :: nc_debug = .True.
844 : character(len=500) :: msg,errmsg
845 : character(len=fnlen) :: title, filnam
846 : type(pawpsp_header_type):: pawpsp_header
847 2923 : type(pspheader_type) :: pspheads_tmp
848 : !arrays
849 2923 : integer,allocatable :: nproj(:)
850 : real(dp) :: tsec(2),ecut_tmp(3,2)
851 2923 : real(dp),allocatable :: e990(:),e999(:),ekb1(:),ekb2(:),epspsp(:),rcpsp(:)
852 2923 : real(dp),allocatable :: rms(:)
853 : #if defined HAVE_LIBPSML
854 : !! usexml= 0 for non xml ps format ; =1 for xml ps format
855 : character(len=3) :: atmsymb
856 : character(len=30) :: creator
857 2923 : type(pspheader_type) :: psphead
858 : #endif
859 : ! ******************************************************************************
860 :
861 : !paral_mode defines how we access to the psp file
862 : ! paral_mode=0: all processes access to the file (sequentially)
863 : ! paral_mode=1: only proc. 0 access to the file and then broadcast
864 2923 : paral_mode=0
865 2923 : if (present(comm_mpi)) then
866 742 : if (psps%usepaw==1.and.xmpi_comm_size(comm_mpi)>1) paral_mode=1
867 : end if
868 276 : me=0;if (paral_mode==1) me=xmpi_comm_rank(comm_mpi)
869 :
870 2457 : if (paral_mode == 1) then
871 276 : ABI_CHECK(psps%usepaw==1, "paral_mode==1 is only compatible with PAW, see call to pawpsp_bcast below")
872 : end if
873 :
874 2923 : nctab%has_tvale = .False.; nctab%has_tcore = .False.
875 2923 : pspcod = -1
876 : !Get electronic temperature from dtset
877 2923 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
878 :
879 2923 : if (me==0) then
880 : ! Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case
881 2717 : if (psps%usepaw==0 .and. psps%mqgrid_ff/=psps%mqgrid_vl) then
882 : write(msg, '(a,a,a,a,a)' )&
883 0 : 'Dimension of q-grid for nl form factors (mqgrid_ff)',ch10,&
884 0 : 'is different from dimension of q-grid for Vloc (mqgrid_vl) !',ch10,&
885 0 : 'This is not allowed for norm-conserving psp.'
886 0 : ABI_ERROR(msg)
887 : end if
888 :
889 2717 : write(msg, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(psps%filpsp(ipsp))
890 8151 : call wrtout([std_out, ab_out], msg)
891 :
892 : ! Check if the file pseudopotential file is written in (XML| XML-PAW | UPF)
893 2717 : call test_xml_xmlpaw_upf(psps%filpsp(ipsp), usexml, xmlpaw, useupf)
894 :
895 : ! ----------------------------------------------------------------------------
896 : ! allocate nproj here: can be read in now for UPF
897 8151 : ABI_MALLOC(nproj,(psps%mpssoang))
898 10651 : nproj(:)=0
899 :
900 2717 : if (usexml /= 1 .and. useupf == 0) then
901 :
902 : ! Open the atomic data file, and read the three first lines
903 : ! These three first lines have a similar format in all allowed psp files
904 : ! Open atomic data file (note: formatted input file)
905 2386 : if (open_file(psps%filpsp(ipsp), msg, unit=tmp_unit, form='formatted', status='old') /= 0) then
906 0 : ABI_ERROR(msg)
907 : end if
908 2386 : rewind (unit=tmp_unit,err=10,iomsg=errmsg)
909 :
910 : ! Read and write some description of file from first line (character data)
911 2386 : read (tmp_unit,'(a)',err=10,iomsg=errmsg) title
912 2386 : write(msg, '(a,a)' ) '- ',trim(title)
913 7158 : call wrtout([std_out, ab_out], msg)
914 :
915 : ! Read and write more data describing psp parameters
916 2386 : read (tmp_unit,*,err=10,iomsg=errmsg) znucl,zion,pspdat
917 2386 : write(msg,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat'
918 7158 : call wrtout([std_out, ab_out], msg)
919 :
920 2386 : read (tmp_unit,*,err=10,iomsg=errmsg) pspcod,pspxc,lmax,lloc,mmax,r2well
921 2386 : if(pspxc<0) then
922 : write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) &
923 193 : pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
924 : else
925 : write(msg, '(4i5,i10,f10.5,t47,a)' ) &
926 2193 : pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
927 : end if
928 7158 : call wrtout([std_out, ab_out], msg)
929 :
930 331 : else if (usexml == 1 .and. xmlpaw == 0) then
931 :
932 : ! the following is probably useless - already read in everything in inpspheads
933 : #if defined HAVE_LIBPSML
934 15 : call psxml2abheader( psps%filpsp(ipsp), psphead, atmsymb, creator, 0 )
935 15 : ABI_FREE(psphead%nproj)
936 15 : ABI_FREE(psphead%nprojso)
937 15 : znucl = psphead%znuclpsp
938 15 : zion = psphead%zionpsp
939 15 : pspdat = psphead%pspdat
940 15 : pspcod = psphead%pspcod
941 15 : pspxc = psphead%pspxc
942 15 : lmax = psphead%lmax
943 : !lloc = 0 ! does this mean s? in psml case the local potential can be different from any l channel
944 15 : lloc = -1
945 15 : mmax = -1
946 15 : r2well = 0
947 :
948 15 : write(msg,'(a,1x,a3,3x,a)') "-",atmsymb,trim(creator)
949 45 : call wrtout([std_out, ab_out], msg)
950 15 : write(msg,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat'
951 45 : call wrtout([std_out, ab_out], msg)
952 15 : if(pspxc<0) then
953 : write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) &
954 15 : pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
955 : else
956 : write(msg, '(4i5,i10,f10.5,t47,a)' ) &
957 0 : pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
958 : end if
959 45 : call wrtout([std_out, ab_out], msg)
960 : #else
961 : write(msg,'(a,a)') &
962 : 'ABINIT is not compiled with XML support for reading this type of pseudopotential ', trim(psps%filpsp(ipsp))
963 : ABI_BUG(msg)
964 : #endif
965 : ! END useless
966 316 : else if (usexml == 1 .and. xmlpaw == 1) then
967 : write(msg,'(a,a)') &
968 305 : '- pspatm : Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp))
969 915 : call wrtout([std_out, ab_out], msg)
970 :
971 : ! Return header information
972 305 : call pawpsxml2ab(psps%filpsp(ipsp),ecut_tmp, pspheads_tmp,0)
973 305 : ABI_FREE(pspheads_tmp%nproj)
974 305 : ABI_FREE(pspheads_tmp%nprojso)
975 305 : lmax=pspheads_tmp%lmax
976 305 : pspxc=pspheads_tmp%pspxc
977 305 : znucl=pspheads_tmp%znuclpsp
978 305 : pawpsp_header%basis_size=pspheads_tmp%pawheader%basis_size
979 305 : pawpsp_header%l_size=pspheads_tmp%pawheader%l_size
980 305 : pawpsp_header%lmn_size=pspheads_tmp%pawheader%lmn_size
981 305 : pawpsp_header%mesh_size=pspheads_tmp%pawheader%mesh_size
982 305 : pawpsp_header%pawver=pspheads_tmp%pawheader%pawver
983 305 : pawpsp_header%shape_type=pspheads_tmp%pawheader%shape_type
984 305 : pawpsp_header%rpaw=pspheads_tmp%pawheader%rpaw
985 305 : pawpsp_header%rshp=pspheads_tmp%pawheader%rshp
986 305 : lloc=0; pspcod=17
987 :
988 11 : else if (useupf /= 0) then
989 11 : if (psps%usepaw /= 0) then
990 0 : ABI_ERROR("UPF format not allowed with PAW (USPP part not read yet)")
991 : end if
992 :
993 11 : r2well = 0; qchrg=zero
994 :
995 : ! should initialize znucl,zion,pspxc,lmax,lloc,mmax
996 11 : if (useupf == 1) then
997 1 : pspcod = 11
998 : call upf1_to_abinit(psps%filpsp(ipsp), znucl, zion, pspxc, lmax, lloc, mmax, &
999 1 : psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj, vlspl, xccc1d)
1000 : else
1001 10 : pspcod = 12
1002 : call upf2_to_abinit(ipsp, psps%filpsp(ipsp), dtset%vloc_rcut, znucl, zion, pspxc, lmax, lloc, mmax, &
1003 10 : psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj, vlspl, xccc1d, xcctau1d, nctab, maxrad)
1004 :
1005 : if (nc_debug) then
1006 : call psp_dump_outputs("UPF2", pspcod, psps%lmnmax, psps%lnmax, psps%mpssoang, &
1007 : psps%mqgrid_ff, psps%n1xccc, mmax, maxrad, epsatm, qchrg, xcccrc, nctab, &
1008 : indlmn, nproj, ekb, ffspl, vlspl, xccc1d, xcctau1d)
1009 : end if
1010 : end if
1011 :
1012 : else
1013 0 : ABI_ERROR("You should not be here! erroneous type or pseudopotential input")
1014 : end if
1015 :
1016 : ! ------------------------------------------------------------------------------
1017 : ! Check data for consistency against main routine input
1018 :
1019 : ! Does required spin-orbit characteristics agree with format
1020 : ! TODO: in case of pspcod 5 (phoney) and 8 (oncvpsp) this is not specific enough.
1021 : ! they can be non-SOC as well.
1022 : ! HGH is ok - can always turn SOC on or off.
1023 : ! PAW is ok - can be used with or without SOC
1024 : ! write(std_out,*) pspso
1025 : if((pspcod/=3).and.(pspcod/=5).and.(pspcod/=8).and.(pspcod/=10).and. pspcod /= 12 .and. &
1026 2717 : (pspcod/=7).and.(pspcod/=17))then
1027 : ! If pspso requires internal characteristics, set it to 1 for non-HGH psps
1028 1217 : if(psps%pspso(ipsp)==1) psps%pspso(ipsp)=0
1029 1217 : if(psps%pspso(ipsp) /= 0)then
1030 : write(msg, '(3a,i0,3a)' )&
1031 0 : 'Pseudopotential file cannot give spin-orbit characteristics,',ch10,&
1032 0 : 'while pspso(itypat)= ',psps%pspso(ipsp),'.',ch10,&
1033 0 : 'Action: check your pseudopotential and input files for consistency.'
1034 0 : ABI_ERROR(msg)
1035 : end if
1036 : end if
1037 :
1038 : ! Does nuclear charge znuclpsp agree with psp input znucl
1039 2717 : if (abs(psps%znuclpsp(ipsp)-znucl) > tol8) then
1040 : write(msg, '(a,f10.5,2a,f10.5,5a)' )&
1041 0 : 'Pseudopotential file znucl: ',znucl,ch10,&
1042 0 : 'does not equal input znuclpsp: ',psps%znuclpsp(ipsp),' better than 1e-08 .',ch10,&
1043 0 : 'znucl is read from the psp file in pspatm, while',ch10,&
1044 0 : 'znuclpsp is read in iofn2.'
1045 0 : ABI_BUG(msg)
1046 : end if
1047 :
1048 : ! Is the highest angular momentum within limits?
1049 : ! Recall mpsang is 1+highest l for nonlocal correction.
1050 : ! Nonlocal corrections for s, p, d, and f are supported.
1051 2717 : if (lmax+1>psps%mpsang) then
1052 : write(msg, '(a,i0,a,i0,a,a)' )&
1053 0 : 'input lmax+1: ',lmax+1,' exceeds mpsang: ',psps%mpsang,ch10,&
1054 0 : 'indicates input lmax too large for dimensions.'
1055 0 : ABI_BUG(msg)
1056 : end if
1057 :
1058 : ! Check several choices for ixc against pspxc
1059 : ! ixc is from ABINIT code; pspxc is from atomic psp file
1060 2717 : if (dtset%ixc==0) then
1061 15 : ABI_WARNING('Note that input ixc=0 => no xc is being used.')
1062 2702 : else if(dtset%ixc/=pspxc) then
1063 : write(msg, '(a,i0,a,i0,8a)' )&
1064 763 : 'Pseudopotential file pspxc: ',pspxc,', not equal to input ixc: ',dtset%ixc,'.',ch10,&
1065 763 : 'These parameters must agree to get the same xc in ABINIT code as in psp construction.',ch10,&
1066 763 : 'Action: check psp design or input file.',ch10,&
1067 1526 : 'Assume experienced user. Execution will continue.'
1068 763 : ABI_WARNING(msg)
1069 : end if
1070 :
1071 2717 : if (lloc>lmax .and. pspcod/=4 .and. pspcod/=8 .and. pspcod/=10) then
1072 : write(msg, '(a,2i12,a,a,a,a)' )&
1073 0 : 'lloc,lmax=',lloc,lmax,ch10,&
1074 0 : 'chosen l of local psp exceeds range from input data.',ch10,&
1075 0 : 'Action: check pseudopotential input file.'
1076 0 : ABI_ERROR(msg)
1077 : end if
1078 :
1079 : ! Does the pspcod agree with type of calculation (paw or not)?
1080 2717 : if (((pspcod/=7.and.pspcod/=17).and.psps%usepaw==1).or.((pspcod==7.or.pspcod==17).and.psps%usepaw==0)) then
1081 : write(msg, '(a,i0,a,a,i0,a)' )&
1082 0 : 'In reading atomic psp file, finds pspcod= ',pspcod,ch10,&
1083 0 : 'This is not an allowed value with usepaw= ',psps%usepaw,'.'
1084 0 : ABI_BUG(msg)
1085 : end if
1086 :
1087 2717 : if (.not.psps%vlspl_recipSpace .and. (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10 .and. pspcod /= 7)) then
1088 : ! The following "if" statement can substitute the one just before once libBigDFT
1089 : ! has been upgraded to include pspcod 10
1090 : ! if (.not.psps%vlspl_recipSpace .and. (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10)) then
1091 : write(msg, '(a,i0,2a)' )&
1092 0 : 'In reading atomic psp file, finds pspcod=',pspcod,ch10,&
1093 0 : 'This is not an allowed value with real space computation.'
1094 0 : ABI_BUG(msg)
1095 : end if
1096 :
1097 : ! MJV 16/6/2009 added pspcod 11 for upf format
1098 : !if( pspcod<1 .or. (pspcod>11.and.pspcod/=17) ) then
1099 2717 : if( pspcod<1) then
1100 : write(msg, '(a,i0,4a)' )&
1101 0 : 'In reading atomic psp file, finds pspcod= ',pspcod,ch10,&
1102 0 : 'This is not an allowed value. Allowed values are 1-12 or 17 .',ch10,&
1103 0 : 'Action: check pseudopotential input file.'
1104 0 : ABI_ERROR(msg)
1105 : end if
1106 :
1107 : ! -----------------------------------------------------------------------
1108 : ! Set various terms to 0 in case not defined below
1109 8151 : ABI_MALLOC(e990,(psps%mpssoang))
1110 5434 : ABI_MALLOC(e999,(psps%mpssoang))
1111 5434 : ABI_MALLOC(rcpsp,(psps%mpssoang))
1112 5434 : ABI_MALLOC(rms,(psps%mpssoang))
1113 5434 : ABI_MALLOC(epspsp,(psps%mpssoang))
1114 5434 : ABI_MALLOC(ekb1,(psps%mpssoang))
1115 5434 : ABI_MALLOC(ekb2,(psps%mpssoang))
1116 18585 : e990(:)=zero ;e999(:)=zero
1117 18585 : rcpsp(:)=zero;rms(:)=zero
1118 18585 : ekb1(:)=zero ;ekb2(:)=zero
1119 10651 : epspsp(:)=zero
1120 2717 : qchrg=zero
1121 :
1122 : ! ----------------------------------------------------------------------
1123 2717 : if (pspcod==1 .or. pspcod==4)then
1124 :
1125 : ! Teter pseudopotential (pspcod=1 or 4)
1126 : call psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,&
1127 : e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,&
1128 : mmax,psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,pspcod,qchrg,psps%qgrid_ff,&
1129 694 : rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp))
1130 :
1131 328 : else if (pspcod==2)then
1132 :
1133 : ! GTH pseudopotential
1134 328 : call psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion)
1135 62320 : xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0
1136 :
1137 318 : else if (pspcod==3)then
1138 :
1139 : ! HGH pseudopotential
1140 318 : call psp3in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), vlspl,zion)
1141 17232 : xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0
1142 :
1143 145 : else if (pspcod==5)then
1144 :
1145 : ! Old phoney pseudopotentials
1146 : call psp5in(ekb,ekb1,ekb2,epsatm,epspsp,&
1147 : e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,&
1148 : mmax,psps%mpsang,psps%mpssoang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%pspso(ipsp),qchrg,psps%qgrid_ff,&
1149 145 : rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp))
1150 :
1151 179 : else if (pspcod==6)then
1152 : ! FHI pseudopotentials
1153 : call psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,&
1154 : psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%optnlxccc,psps%positron,qchrg,psps%qgrid_ff,psps%useylm,vlspl,&
1155 179 : xcccrc,xccc1d,zion,psps%znuclpsp(ipsp))
1156 :
1157 : else if (pspcod==7)then
1158 : ! PAW "pseudopotentials"
1159 : call pawpsp_7in(epsatm,ffspl,dtset%icoulomb,ABS(dtset%hyb_mixing),dtset%ixc,&
1160 : & lmax,psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,&
1161 : & pawrad,pawtab,dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,&
1162 : & dtset%usewvl,dtset%usexcnhat_orig,vlspl,xcccrc,dtset%xclevel,&
1163 : & dtset%xc_denpos,zion,psps%znuclpsp(ipsp),&
1164 231 : & xc_taupos=dtset%xc_taupos,el_temp=el_temp)
1165 :
1166 479 : else if (pspcod==8)then
1167 :
1168 : ! DRH pseudopotentials
1169 : call psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,&
1170 : psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc,psps%pspso(ipsp),&
1171 479 : qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad)
1172 :
1173 : if (nc_debug) then
1174 : call psp_dump_outputs("PSP8",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, &
1175 : psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, &
1176 : indlmn,nproj,ekb,ffspl,vlspl,xccc1d)
1177 : end if
1178 :
1179 15 : else if (pspcod==9)then
1180 :
1181 : #if defined HAVE_LIBPSML
1182 : call psp9in(psps%filpsp(ipsp),ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,&
1183 : psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc, &
1184 : psps%pspso(ipsp),qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,&
1185 15 : xcccrc,xccc1d,xcctau1d,zion,psps%znuclpsp(ipsp),nctab,maxrad)
1186 :
1187 : if (nc_debug) then
1188 : call psp_dump_outputs("PSML",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, &
1189 : psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, &
1190 : indlmn,nproj,ekb,ffspl,vlspl,xccc1d,xcctau1d)
1191 : end if
1192 : #else
1193 : write(msg,'(2a)') &
1194 : 'ABINIT is not compiled with XML support for reading this type of pseudopotential ', trim(psps%filpsp(ipsp))
1195 : ABI_ERROR(msg)
1196 : #endif
1197 :
1198 12 : else if (pspcod==10)then
1199 :
1200 : ! HGH pseudopotential, full h/k matrix read
1201 12 : call psp10in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), vlspl,zion)
1202 84 : xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0
1203 :
1204 : ! NB for pspcod 11 the reading has already been done above.
1205 : else if (pspcod==17)then
1206 : ! PAW XML pseudopotentials
1207 : call pawpsp_17in(epsatm,ffspl,dtset%icoulomb,ipsp,ABS(dtset%hyb_mixing),dtset%ixc,lmax,&
1208 : & psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,pawpsp_header,pawrad,pawtab,&
1209 : & dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,dtset%usewvl,&
1210 : & dtset%usexcnhat_orig,vlspl,xcccrc,&
1211 : & dtset%xclevel,dtset%xc_denpos,pspheads_tmp%zionpsp,psps%znuclpsp(ipsp),&
1212 305 : & xc_taupos=dtset%xc_taupos,el_temp=el_temp)
1213 305 : call paw_setup_free(paw_setuploc)
1214 : end if
1215 :
1216 2717 : close (unit=tmp_unit)
1217 :
1218 : ! ----------------------------------------------------------------------
1219 2717 : if (pspcod==2 .or. pspcod==3 .or. pspcod==10)then
1220 658 : write(msg, '(a,a,a,a,a,a,a,a,a,a)' )ch10,&
1221 658 : ' pspatm : COMMENT -',ch10,&
1222 658 : ' the projectors are not normalized,',ch10,&
1223 658 : ' so that the KB energies are not consistent with ',ch10,&
1224 658 : ' definition in PRB44, 8503 (1991). ',ch10,& ! [[cite:Gonze1991]]
1225 1316 : ' However, this does not influence the results obtained hereafter.'
1226 1974 : call wrtout([std_out, ab_out], msg)
1227 : ! The following lines are added to keep backward compatibilty
1228 658 : maxrad=zero
1229 : #if defined HAVE_BIGDFT
1230 : do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0
1231 : if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp))
1232 : end do
1233 : if (abs(maxrad)<=tol12) then
1234 : psps%gth_params%radii_cf(ipsp,3)=zero
1235 : else
1236 : psps%gth_params%radii_cf(ipsp,3)=max( &
1237 : & min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, &
1238 : & psps%gth_params%radii_cf(ipsp,2))
1239 : end if
1240 : #endif
1241 : end if
1242 :
1243 2717 : if (pspcod/=7.and.pspcod/=17) then
1244 2181 : write(msg, '(a,f14.8,a,a)' ) ' pspatm : epsatm=',epsatm,ch10,' --- l ekb(1:nproj) -->'
1245 6543 : call wrtout([std_out, ab_out], msg)
1246 : !print *, "nproj", nproj
1247 2181 : iln0=0
1248 12782 : do ilmn=1,psps%lmnmax
1249 10601 : iln=indlmn(5,ilmn)
1250 12782 : if (iln>iln0) then
1251 4706 : il=indlmn(1,ilmn)
1252 4706 : if (indlmn(6,ilmn)==1) then
1253 4012 : iln0=iln0+nproj(il+1)
1254 10427 : write(msg, '(13x,i1,4f12.6)' ) il,(ekb(iln+ii),ii=0,nproj(il+1)-1)
1255 : else
1256 : ! Note il = ll here i.e. the s channel in the SOC part is not included in nproj.
1257 694 : iln0=iln0+nproj(il+psps%mpsang)
1258 694 : if (dtset%spnorbscl /= one) then
1259 : call wrtout([std_out, ab_out], &
1260 6 : sjoin(" Rescaling spin-orbit KB energies using spnorbscl.", ftoa(dtset%spnorbscl)))
1261 10 : ekb(iln:iln + nproj(il + psps%mpsang) - 1) = ekb(iln:iln + nproj(il + psps%mpsang) - 1) * dtset%spnorbscl
1262 : end if
1263 694 : if (dtset%so_psp(ipsp) == 0) cycle
1264 1074 : write(msg, '(2x,a,i1,4f12.6)' ) 'spin-orbit ',il,(ekb(iln+ii),ii=0,nproj(il+psps%mpsang)-1)
1265 : end if
1266 13950 : call wrtout([std_out, ab_out], msg)
1267 : end if
1268 : end do
1269 : end if
1270 :
1271 : ! NC: Evalute spline-fit of the model core charge in reciprocal space.
1272 : ! idem for kinetic energy density, which might be 0
1273 : ! TODO: Be careful, because we will be using the PAW routines in which tcore is always available!
1274 : ! Should add a test with 2 NC pseudos: one with NLCC and the other without!
1275 2717 : if (psps%usepaw == 0) then
1276 2181 : call nctab_eval_tcorespl(nctab, psps%n1xccc, xcccrc, xccc1d, xcctau1d, psps%mqgrid_vl, psps%qgrid_vl)
1277 : end if
1278 :
1279 2717 : write(msg,'(3a)') ' pspatm: atomic psp has been read ',' and splines computed',ch10
1280 8151 : call wrtout([std_out, ab_out], msg)
1281 :
1282 2717 : ABI_FREE(e990)
1283 2717 : ABI_FREE(e999)
1284 2717 : ABI_FREE(rcpsp)
1285 2717 : ABI_FREE(rms)
1286 2717 : ABI_FREE(ekb1)
1287 2717 : ABI_FREE(ekb2)
1288 2717 : ABI_FREE(epspsp)
1289 2717 : ABI_FREE(nproj)
1290 :
1291 5434 : if (dtset%prtvol > 9 .and. psps%usepaw==0 .and. psps%lmnmax>3) then
1292 23 : write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_pspdata), ipsp, ".dat"
1293 23 : if (open_file(filnam, msg, newunit=unt) /= 0) then
1294 0 : ABI_ERROR(msg)
1295 : end if
1296 23 : write (unt,*) '# Pseudopotential data in reciprocal space as used by ABINIT'
1297 23 : write (unt,'(a)', ADVANCE='NO') '# index vlocal '
1298 23 : if (psps%lnmax > 0) write (unt,'(a,I3)', ADVANCE='NO') ' 1st proj(l=', indlmn(1,1)
1299 23 : if (psps%lnmax > 1) write (unt,'(a,I3)', ADVANCE='NO') ') 2nd(l=', indlmn(1,2)
1300 23 : if (psps%lnmax > 2) write (unt,'(a,I3,a)', ADVANCE='NO') ') 3rd(l=', indlmn(1,3), ')'
1301 23 : write (unt,*)
1302 :
1303 66046 : do ii = 1, psps%mqgrid_vl
1304 66023 : write(unt, '(I5,E24.16)', ADVANCE='NO') ii, vlspl(ii,1)
1305 66023 : if (psps%lnmax > 0) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,1)
1306 66023 : if (psps%lnmax > 1) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,2)
1307 66023 : if (psps%lnmax > 2) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,3)
1308 66046 : write(unt, *)
1309 : end do
1310 23 : close(unt)
1311 :
1312 23 : write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_nlcc_derivs), ipsp, ".dat"
1313 23 : if (open_file(filnam, msg, newunit=unt) /= 0) then
1314 0 : ABI_ERROR(msg)
1315 : end if
1316 23 : write (unt,*) '# Non-linear core corrections'
1317 23 : write (unt,*) '# r, pseudocharge, 1st, 2nd, 3rd, 4th, 5th derivatives'
1318 25033 : do ii = 1, psps%n1xccc
1319 25010 : write (unt,*) xcccrc*(ii-1)/(psps%n1xccc-1), xccc1d(ii,1), xccc1d(ii,2), &
1320 50043 : xccc1d(ii,3), xccc1d(ii,4), xccc1d(ii,5), xccc1d(ii,6)
1321 : end do
1322 23 : write (unt,*) '# r, pseudochg kinE, 1st, 2nd, 3rd, 4th, 5th derivatives'
1323 25033 : do ii = 1, psps%n1xccc
1324 25010 : write (unt,*) xcccrc*(ii-1)/(psps%n1xccc-1), xcctau1d(ii,1), xcctau1d(ii,2), &
1325 50043 : xcctau1d(ii,3), xcctau1d(ii,4), xcctau1d(ii,5), xcctau1d(ii,6)
1326 : end do
1327 23 : close(unt)
1328 : end if
1329 :
1330 : end if ! me=0
1331 :
1332 2923 : if (paral_mode==1) then
1333 276 : call timab(48,1,tsec)
1334 276 : call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc)
1335 276 : call timab(48,2,tsec)
1336 : end if
1337 :
1338 2923 : if (psps%usepaw==1) then
1339 60823 : indlmn(:,:)=0
1340 57722 : indlmn(1:6,1:pawtab%lmn_size)=pawtab%indlmn(1:6,1:pawtab%lmn_size)
1341 : end if
1342 :
1343 : !--------------------------------------------------------------------
1344 : !WVL+PAW:
1345 2923 : if (dtset%usepaw==1 .and. (dtset%icoulomb /= 0 .or. dtset%usewvl==1)) then
1346 : #if defined HAVE_BIGDFT
1347 : psps%gth_params%psppar(:,:,ipsp) = UNINITIALIZED(1._dp)
1348 : psps%gth_params%radii_cf(ipsp,:) = UNINITIALIZED(1._dp)
1349 : call wvl_descr_psp_fill(psps%gth_params, ipsp, psps%pspxc(1), int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), 0)
1350 : #endif
1351 :
1352 : ! The following lines are added to keep backward compatibilty
1353 0 : maxrad=zero
1354 : #if defined HAVE_BIGDFT
1355 : do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0
1356 : if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp))
1357 : end do
1358 : if (abs(maxrad)<=tol12) then
1359 : !== MT COMMENT
1360 : ! Damien wants to activate this (in order to directly compare to bigDFT):
1361 : psps%gth_params%radii_cf(ipsp,3)= psps%gth_params%radii_cf(ipsp,2)
1362 : ! But, this changes strongly file references.
1363 : ! So, I keep this, waiting for Tonatiuh s validation
1364 : psps%gth_params%radii_cf(ipsp,3) = (psps%gth_params%radii_cf(ipsp,1)+psps%gth_params%radii_cf(ipsp,2))*half
1365 : !== MT COMMENT
1366 : else
1367 : psps%gth_params%radii_cf(ipsp,3)=max( &
1368 : min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, &
1369 : psps%gth_params%radii_cf(ipsp,2))
1370 : end if
1371 : if(present(comm_mpi)) then
1372 : call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss,comm_mpi)
1373 : else
1374 : call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss)
1375 : end if
1376 : #endif
1377 : end if
1378 :
1379 : !end of WVL+PAW section
1380 : !----------------------------------------------------
1381 :
1382 : return
1383 :
1384 : ! Handle IO error
1385 : 10 continue
1386 0 : ABI_ERROR(errmsg)
1387 :
1388 2923 : end subroutine pspatm
1389 : !!***
1390 :
1391 : !!****f* ABINIT/psp_dump_outputs
1392 : !! NAME
1393 : !! psp_dump_outputs
1394 : !!
1395 : !! FUNCTION
1396 : !! Debugging routines used to dumo PSP data in Yaml format.
1397 : !!
1398 : !! SOURCE
1399 :
1400 : subroutine psp_dump_outputs(pfx,pspcod,lmnmax,lnmax,mpssoang, &
1401 : mqgrid,n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, &
1402 : indlmn,nproj,ekb,ffspl,vlspl,xccc1d,xcctau1d)
1403 :
1404 : !Arguments ------------------------------------
1405 : !scalars
1406 : character(len=*), intent(in) :: pfx
1407 : integer,intent(in) :: pspcod,lmnmax,lnmax,mpssoang,mqgrid,n1xccc
1408 : integer,intent(in) :: mmax
1409 : real(dp),intent(in) :: maxrad,epsatm,qchrg,xcccrc
1410 : type(nctab_t),intent(in) :: nctab
1411 : !arrays
1412 : integer,intent(in) :: indlmn(6,lmnmax),nproj(mpssoang)
1413 : real(dp),intent(in) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
1414 : real(dp),intent(in) :: xccc1d(n1xccc,6)
1415 : real(dp),intent(in),optional :: xcctau1d(n1xccc,6)
1416 :
1417 : !Local variables ------------------------------
1418 : !scalars
1419 : integer, parameter :: dump = 64
1420 : integer :: ierr, i, j ,k
1421 : character(len=500) :: msg
1422 : ! *********************************************************************
1423 :
1424 : open(unit=dump, file=trim(pfx)//"_psp_info.yaml", status='REPLACE', err=10, iostat=ierr)
1425 :
1426 : write(dump,'(3a)') "%YAML 1.2", ch10, "---"
1427 :
1428 : write(dump, '(2a)') ch10, "# Pseudopotential info"
1429 : write(dump, '(a,1x,i8)') "pspcod:", pspcod
1430 :
1431 : write(dump, '(2a)') ch10, "# Array dimensions"
1432 : write(dump, '(a)') "dims:"
1433 : write(dump, '(4x,a,1x,i8)') "lmnmax:", lmnmax
1434 : write(dump, '(4x,a,1x,i8)') "lnmax:", lnmax
1435 : write(dump, '(4x,a,1x,i8)') "mpssoang:", mpssoang
1436 : write(dump, '(4x,a,1x,i8)') "mqgrid:", mqgrid
1437 : write(dump, '(4x,a,1x,i8)') "n1xccc:", n1xccc
1438 : write(dump, '(4x,a,1x,i8)') "mmax:", mmax
1439 :
1440 : write(dump, '(2a)') ch10, "# Quantities"
1441 : write(dump, '(a,1x,e12.5)') "maxrad:", maxrad
1442 : write(dump, '(a,1x,e12.5)') "epsatm:", epsatm
1443 : write(dump, '(a,1x,e12.5)') "qchrg:", qchrg
1444 : write(dump, '(a,1x,e12.5)') "xcccrc:", xcccrc
1445 :
1446 : write(dump, '(2a)') ch10, "# Structure: nctab"
1447 : write(dump, '(a)') "nctab:"
1448 : write(dump,'(4x,a,":",1x,i4)') "mqgrid_vl", nctab%mqgrid_vl
1449 : write(dump,'(4x,a,":",1x,l4)') "has_tvale", nctab%has_tvale
1450 : write(dump,'(4x,a,":",1x,l4)') "has_tcore", nctab%has_tcore
1451 : write(dump,'(4x,a,":",1x,e12.5)') "dncdq0", nctab%dncdq0
1452 : write(dump,'(4x,a,":",1x,e12.5)') "d2ncdq0", nctab%d2ncdq0
1453 : write(dump,'(4x,a,":",1x,e12.5)') "dnvdq0", nctab%dnvdq0
1454 :
1455 : if ( nctab%has_tvale ) then
1456 : write(dump, '(2a)') ch10, "# Array: nctab_tvalespl(mqgrid_vl,2)"
1457 : write(dump, '(a)') "nctab_tvalespl:"
1458 : do j=1,2
1459 : do i=1,nctab%mqgrid_vl
1460 : if ( i == 1 ) then
1461 : write(dump,'(4x,a,1x,e12.5)') "- -", nctab%tvalespl(i,j)
1462 : else
1463 : write(dump,'(4x,a,1x,e12.5)') " -", nctab%tvalespl(i,j)
1464 : end if
1465 : end do
1466 : end do
1467 : end if
1468 :
1469 : if ( nctab%has_tcore ) then
1470 : write(dump, '(2a)') ch10, "# Array: nctab_tcorespl(mqgrid_vl,2)"
1471 : write(dump, '(a)') "nctab_tcorespl:"
1472 : do j=1,2
1473 : do i=1,nctab%mqgrid_vl
1474 : if ( i == 1 ) then
1475 : write(dump,'(4x,a,1x,e12.5)') "- -", nctab%tcorespl(i,j)
1476 : else
1477 : write(dump,'(4x,a,1x,e12.5)') " -", nctab%tcorespl(i,j)
1478 : end if
1479 : end do
1480 : end do
1481 : end if
1482 :
1483 : write(dump, '(2a)') ch10, "# Array: integer indlmn(6,lmnmax)"
1484 : write(dump, '(a)') "indlmn:"
1485 : do i=1,lmnmax
1486 : write(dump,'(4x,a,i4,5(",",i4),a)') "- [", indlmn(:,i), "]"
1487 : end do
1488 :
1489 : write(dump, '(2a)') ch10, "# Array: integer nproj(mpssoang)"
1490 : write(dump, '(a)') "nproj:"
1491 : do i=1,mpssoang
1492 : write(dump,'(4x,"-",1x,i4)') nproj(i)
1493 : end do
1494 :
1495 : write(dump, '(2a)') ch10, "# Array: double ekb(lnmax)"
1496 : write(dump, '(a)') "ekb:"
1497 : do i=1,lnmax
1498 : write(dump,'(4x,"-",1x,e12.5)') ekb(i)
1499 : end do
1500 :
1501 : write(dump, '(2a)') ch10, "# Array: ffspl(mqgrid,2,lnmax)"
1502 : write(dump, '(a)') "ffspl:"
1503 : do k=1,lnmax
1504 : do j=1,2
1505 : do i=1,mqgrid
1506 : if ( (i == 1) .and. (j == 1) ) then
1507 : write(dump,'(4x,a,1x,e12.5)') "- - -", ffspl(i,j,k)
1508 : else if ( i == 1 ) then
1509 : write(dump,'(4x,a,1x,e12.5)') " - -", ffspl(i,j,k)
1510 : else
1511 : write(dump,'(4x,a,1x,e12.5)') " -", ffspl(i,j,k)
1512 : end if
1513 : end do
1514 : end do
1515 : end do
1516 :
1517 : write(dump, '(2a)') ch10, "# Array: vlspl(mqgrid,2)"
1518 : write(dump, '(a)') "vlspl:"
1519 : do j=1,2
1520 : do i=1,mqgrid
1521 : if ( i == 1 ) then
1522 : write(dump,'(4x,a,1x,e12.5)') "- -", vlspl(i,j)
1523 : else
1524 : write(dump,'(4x,a,1x,e12.5)') " -", vlspl(i,j)
1525 : end if
1526 : end do
1527 : end do
1528 :
1529 : write(dump, '(2a)') ch10, "# Array: xccc1d(n1xccc,6)"
1530 : write(dump, '(a)') "xccc1d:"
1531 : do j=1,6
1532 : do i=1,n1xccc
1533 : if ( i == 1 ) then
1534 : write(dump,'(4x,a,1x,e12.5)') "- -", xccc1d(i,j)
1535 : else
1536 : write(dump,'(4x,a,1x,e12.5)') " -", xccc1d(i,j)
1537 : end if
1538 : end do
1539 : end do
1540 :
1541 : if (present(xcctau1d)) then
1542 : write(dump, '(2a)') ch10, "# Array: xcctau1d(n1xccc,6)"
1543 : write(dump, '(a)') "xcctau1d:"
1544 : do j=1,6
1545 : do i=1,n1xccc
1546 : if ( i == 1 ) then
1547 : write(dump,'(4x,a,1x,e12.5)') "- -", xcctau1d(i,j)
1548 : else
1549 : write(dump,'(4x,a,1x,e12.5)') " -", xcctau1d(i,j)
1550 : end if
1551 : end do
1552 : end do
1553 : end if
1554 :
1555 : write (dump,'(2a)') ch10, "..."
1556 :
1557 : close(dump)
1558 :
1559 : return
1560 : 10 continue
1561 :
1562 : if (ierr /= 0) then
1563 : write(msg,'(3a,i0)') "Error writing pseudopotential information", ch10, "IOSTAT=", ierr
1564 : ABI_WARNING(msg)
1565 : end if
1566 :
1567 : end subroutine psp_dump_outputs
1568 : !!***
1569 :
1570 : end module m_pspini
1571 : !!***
|
