Branch data 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 : 6988 : 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 : 6988 : type(pawrad_type) :: pawrad_dum
150 : 6988 : type(pawtab_type) :: pawtab_dum
151 : 6988 : 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 : 6988 : integer,allocatable :: indlmn_alch(:,:,:),new_pspso(:)
158 : 6988 : integer,pointer :: indlmn(:,:)
159 : : real(dp),save :: epsatm(npspmax)
160 : : real(dp) :: tsec(2)
161 : 6988 : real(dp),allocatable :: dvlspl(:,:),dvlspl_alch(:,:,:),ekb(:),ekb_alch(:,:)
162 : 6988 : real(dp),allocatable :: epsatm_alch(:),ffspl(:,:,:),ffspl_alch(:,:,:,:)
163 : 6988 : real(dp),allocatable :: vlspl(:,:),vlspl_alch(:,:,:),xccc1d(:,:)
164 : 6988 : real(dp),allocatable :: xcctau1d(:,:)
165 : 6988 : real(dp),allocatable :: xccc1d_alch(:,:,:),xcccrc_alch(:)
166 : 6988 : real(dp),allocatable :: xcctau1d_alch(:,:,:)
167 : 6988 : type(nctab_t),target,allocatable :: nctab_alch(:)
168 : : ! *************************************************************************
169 : :
170 : : DBG_ENTER("COLL")
171 : :
172 : : ! Keep track of time spent in this subroutine
173 : 6988 : call timab(15,1,tsec)
174 : :
175 : : !-------------------------------------------------------------
176 : : ! Some initializations
177 : : !-------------------------------------------------------------
178 : :
179 : : !Useful sizes
180 : 6988 : ntypat=psps%ntypat
181 : 6988 : mtypalch=psps%mtypalch
182 : 6988 : npsp=psps%npsp
183 [ - + ]: 6988 : 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 : 6988 : qmax=1.2d0 * sqrt(gsqcut)
190 : : !ffnl is always computed in reciprocal space
191 : 6988 : dq=qmax/(one*(psps%mqgrid_ff-1))
192 [ + + ]: 21126347 : do ii=1,psps%mqgrid_ff
193 : 21126347 : psps%qgrid_ff(ii)=(ii-1)*dq
194 : : end do
195 [ + + ]: 6988 : 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 [ + - ]: 6988 : if (psps%vlspl_recipSpace) then
199 : 6988 : 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 [ + + ]: 21175078 : do ii=1,psps%mqgrid_vl
208 : 21175078 : psps%qgrid_vl(ii)=(ii-1)*dq
209 : : end do
210 : :
211 : : !Determine whether new optional data requests have changed
212 : 6988 : paw_options=0;paw_size=0
213 [ + + ]: 6988 : if (psps%usepaw==1) then
214 : 1475 : paw_size=size(pawtab)
215 [ + + - + : 1475 : has_kij=(dtset%positron/=0.or.abs(dtset%effmass_free-one)>tol8.or.dtset%orbmag>0.or.dtset%use_rcpaw==1)
+ + + + ]
216 : 1475 : has_tvale=.true. ! Will be modified later (depending on PAW dataset format)
217 : 1475 : has_nabla=.false.
218 [ + + + - ]: 1475 : has_shapefncg=(dtset%optdriver==RUNL_GSTATE.and.((dtset%iprcel>=20.and.dtset%iprcel<70).or.dtset%iprcel>=80))
219 [ + - - + ]: 1475 : has_wvl=(dtset%usewvl==1.or.dtset%icoulomb/=0)
220 [ + - + + ]: 1475 : has_tproj=(dtset%usewvl==1.or.dtset%use_rcpaw==1) ! projectors will be free at the end of the psp reading
221 [ + + ]: 3595 : has_vminushalf=(maxval(dtset%ldaminushalf)==1)
222 : 1475 : has_coretau=(dtset%usekden>=1)
223 : 1475 : paw_add_core=(dtset%paw_add_core==1)
224 [ + + ]: 1475 : if (has_kij) paw_options(1)=1
225 : 1475 : if (has_tvale) paw_options(2)=1
226 : : if (has_nabla) paw_options(5)=1
227 [ - + ]: 1475 : if (has_shapefncg) paw_options(6)=1
228 [ - + ]: 1475 : if (has_wvl) paw_options(7)=1
229 [ + + ]: 1475 : if (has_tproj) paw_options(8)=1
230 [ + + ]: 1475 : if (has_vminushalf)paw_options(9)=1
231 [ + + ]: 1475 : if (has_coretau) paw_options(10)=1
232 [ + + ]: 1475 : if (paw_add_core) paw_options(11)=1
233 : : !if (dtset%prtvclmb /= 0) then
234 : 1475 : paw_options(3) = 1
235 : 1475 : 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 [ - + + - : 20964 : ABI_MALLOC(new_pspso,(npsp))
+ - + - -
+ - + ]
242 [ + + + + ]: 77238 : if (ifirst==1) pspso_old(:)=-1
243 [ + + + + ]: 77238 : if (ifirst==1) pspso_zero(:)=-1
244 [ + + ]: 16899 : do ipsp=1,npsp
245 : 9911 : 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 [ + + + + : 9911 : if (psps%pspso(ipsp)==pspso_old(ipsp).or. &
+ + ]
249 : : & (psps%pspso(ipsp)==1.and.pspso_old(ipsp)==pspso_zero(ipsp))) then
250 : 7958 : new_pspso(ipsp)=0
251 : : end if
252 : : ! No new characteristics if PAW
253 [ + + ]: 9911 : if (psps%usepaw==1) new_pspso(ipsp)=0
254 : : ! Prepare the saving of the intrinsic pseudopotential characteristics
255 [ + + ]: 16899 : if(psps%pspso(ipsp)==1) pspso_zero(ipsp)=0
256 : : end do
257 : :
258 : : !Compute the general condition for new computation of pseudopotentials
259 : 6988 : 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 [ + + + + : 85000 : & .or. (use_rcpaw_old==1) &
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + - +
- + - + +
+ + + + +
+ + + + -
+ - + - +
+ + + + +
- + - - +
+ ]
285 : 2282 : & ) gencond=1
286 : :
287 [ + - + + ]: 6988 : if (present(comm_mpi).and.psps%usepaw==1) then
288 [ + + ]: 1475 : if(xmpi_comm_size(comm_mpi)>1) call xmpi_sum(gencond,comm_mpi,ierr)
289 [ + + ]: 1475 : if (gencond/=0) gencond=1
290 : : end if
291 [ - + ]: 6988 : ABI_FREE(new_pspso)
292 : :
293 : : !-------------------------------------------------------------
294 : : ! Following section is only reached when new computation
295 : : ! of pseudopotentials is needed
296 : : !-------------------------------------------------------------
297 : :
298 [ + + ]: 6988 : if (gencond==1) then
299 : :
300 : 2282 : write(msg, '(a,a)' ) ch10,'--- Pseudopotential description ------------------------------------------------'
301 : 2282 : call wrtout(ab_out,msg)
302 : :
303 [ + + + - : 6846 : ABI_MALLOC(ekb,(psps%dimekb*(1-psps%usepaw)))
+ - + + -
+ - + ]
304 [ + + + - : 6846 : ABI_MALLOC(xccc1d,(psps%n1xccc*(1-psps%usepaw),6))
+ - + + -
+ - + ]
305 [ + + + - : 4564 : ABI_MALLOC(xcctau1d,(psps%n1xccc*(1-psps%usepaw),6))
- + - + ]
306 [ - + + - : 9128 : ABI_MALLOC(ffspl,(psps%mqgrid_ff,2,psps%lnmax))
- + + - +
- + - - +
- + ]
307 [ - + + - : 6846 : ABI_MALLOC(vlspl,(psps%mqgrid_vl,2))
+ - + - -
+ - + ]
308 [ - + ]: 2282 : if (.not.psps%vlspl_recipSpace) then
309 [ # # # # : 0 : ABI_MALLOC(dvlspl,(psps%mqgrid_vl,2))
# # # # ]
310 : : else
311 [ - + ]: 2282 : ABI_MALLOC(dvlspl,(0,0))
312 : : end if
313 : :
314 : : ! PAW: reset flags for optional data
315 [ + + ]: 2282 : 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 : 536 : & 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 [ + + ]: 2282 : if(mtypalch==0)then
328 : : ! No alchemical pseudoatom, in all datasets, npsp=ntypat
329 : :
330 [ + + ]: 5138 : do ipsp=1,npsp
331 : :
332 : 2882 : xcccrc=zero
333 [ + + + + : 93615791 : ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero
+ + + + +
+ + + ]
334 [ - + - - : 2882 : if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero
- - ]
335 [ + + ]: 2882 : if (psps%usepaw==0) then
336 [ + + + + ]: 17106681 : xccc1d(:,:)=zero
337 [ + + + + ]: 17106681 : xcctau1d(:,:)=zero
338 : : end if
339 : 2882 : indlmn=>psps%indlmn(:,:,ipsp)
340 [ + + + + ]: 136246 : indlmn(:,:)=0
341 : :
342 [ + - ]: 2882 : write(msg, '(a,i4,a,t38,a)' )'- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp))
343 [ + + ]: 8646 : 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 [ + + ]: 2882 : if(psps%usepaw==0)then
348 : : call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,&
349 [ - + ]: 2121 : pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,psps%nctab(ipsp))
350 [ + + ]: 11094 : psps%ekb(:,ipsp)=ekb(:)
351 [ + + + + ]: 17106681 : psps%xccc1d(:,:,ipsp)=xccc1d(:,:)
352 [ + + + + ]: 17106681 : psps%xcctau1d(:,:,ipsp)=xcctau1d(:,:)
353 : : else
354 [ + - ]: 761 : 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 [ - + ]: 761 : pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,xcctau1d,nctab_dum,comm_mpi=comm_mpi_)
357 [ + + - + ]: 761 : 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 : 2882 : psps%xcccrc(ipsp)=xcccrc
364 : 2882 : psps%znucltypat(ipsp)=psps%znuclpsp(ipsp)
365 [ + + + + : 75969514 : psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:)
+ + ]
366 [ + + + + ]: 17640186 : psps%vlspl(:,:,ipsp)=vlspl(:,:)
367 [ - + - - : 5138 : 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 [ - + ]: 2282 : ABI_FREE(ekb)
583 [ - + ]: 2282 : ABI_FREE(ffspl)
584 [ - + ]: 2282 : ABI_FREE(vlspl)
585 [ - + ]: 2282 : ABI_FREE(xccc1d)
586 [ - + ]: 2282 : ABI_FREE(xcctau1d)
587 [ - + ]: 2282 : 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 [ + + ]: 16899 : do ipsp=1,npsp
598 : 16899 : psps%epsatm(ipsp)=epsatm(ipsp)
599 : : enddo
600 : 6988 : call pspcor(ecore,epsatm,dtset%natom,ntypat,dtset%typat,psps%ziontypat)
601 [ + + ]: 6988 : if(abs(ecore_old-ecore)>tol8*abs(ecore_old+ecore))then
602 : 1480 : 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 [ + + ]: 1480 : if(dtset%iscf>=0) call wrtout(ab_out,msg)
605 : 1480 : call wrtout(std_out,msg)
606 : : end if
607 : :
608 : : !End of pseudopotential output section
609 : 6988 : write(msg, '(2a)' )'--------------------------------------------------------------------------------',ch10
610 : 6988 : call wrtout(ab_out,msg)
611 : :
612 : : !-------------------------------------------------------------
613 : : ! Keep track of this call to the routine
614 : : !-------------------------------------------------------------
615 : :
616 [ + + ]: 6988 : if (ifirst==1) ifirst=0
617 : :
618 : 6988 : mqgridff_old=psps%mqgrid_ff
619 : 6988 : mqgridvl_old=psps%mqgrid_vl
620 : 6988 : mpssoang_old=psps%mpssoang
621 : 6988 : spnorbscl_old = dtset%spnorbscl
622 : 6988 : ixc_old=dtset%ixc
623 : 6988 : hyb_mixing_old=dtset%hyb_mixing
624 [ + + ]: 6988 : gsqcut_old=gsqcut;if (psps%usepaw==1) gsqcutdg_old=gsqcutdg
625 : 6988 : lmnmax_old=psps%lmnmax
626 : 6988 : lnmax_old=psps%lnmax
627 : 6988 : optnlxccc_old=psps%optnlxccc
628 : 6988 : usepaw_old=psps%usepaw
629 : 6988 : dimekb_old=psps%dimekb
630 : 6988 : useylm_old=psps%useylm
631 : 6988 : pawxcdev_old=dtset%pawxcdev
632 : 6988 : positron_old=dtset%positron
633 : 6988 : usewvl_old = dtset%usewvl
634 : 6988 : usekden_old = dtset%usekden
635 : 6988 : usexcnhat_old=dtset%usexcnhat_orig
636 : 6988 : paw_size_old=paw_size
637 : 6988 : use_rcpaw_old=dtset%use_rcpaw
638 : 6988 : ecore_old=ecore
639 : 6988 : paw_options_old(:)=paw_options(:)
640 : :
641 [ + + ]: 16899 : do ipsp=1,npsp
642 : 9911 : pspso_old(ipsp)=psps%pspso(ipsp)
643 [ + + ]: 16899 : if(pspso_zero(ipsp)==0)pspso_zero(ipsp)=psps%pspso(ipsp)
644 : : end do
645 [ + + + + ]: 63378 : psps%mproj = maxval(psps%indlmn(3,:,:))
646 : :
647 [ + + + + ]: 9270 : 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 [ + + ]: 6988 : 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 : 6988 : call timab(15,2,tsec)
658 : :
659 : : DBG_EXIT("COLL")
660 : :
661 [ - + - + : 20964 : 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 : 7006 : 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 : 7006 : charge = 0.d0
700 : 7006 : esum = 0.d0
701 [ + + ]: 23250 : do ia=1,natom
702 : : ! compute pseudocharge:
703 : 16244 : charge=charge+zion(typat(ia))
704 : : ! add pseudocore energies together:
705 : 23250 : esum = esum + epsatm(typat(ia))
706 : : end do
707 : :
708 : 7006 : ecore=charge*esum
709 : :
710 : 7006 : 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 : 2950 : subroutine pspatm(dq,dtset,dtfil,ekb,epsatm,ffspl,indlmn,ipsp,pawrad,pawtab,&
814 : 2950 : 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 : 2950 : type(pspheader_type) :: pspheads_tmp
848 : : !arrays
849 : 2950 : integer,allocatable :: nproj(:)
850 : : real(dp) :: tsec(2),ecut_tmp(3,2)
851 : 2950 : real(dp),allocatable :: e990(:),e999(:),ekb1(:),ekb2(:),epspsp(:),rcpsp(:)
852 : 2950 : 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 : 2950 : 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 : 2950 : paral_mode=0
865 [ + + ]: 2950 : if (present(comm_mpi)) then
866 [ - + + + ]: 761 : if (psps%usepaw==1.and.xmpi_comm_size(comm_mpi)>1) paral_mode=1
867 : : end if
868 : 284 : me=0;if (paral_mode==1) me=xmpi_comm_rank(comm_mpi)
869 : :
870 [ + + ]: 2473 : if (paral_mode == 1) then
871 [ - + ]: 284 : ABI_CHECK(psps%usepaw==1, "paral_mode==1 is only compatible with PAW, see call to pawpsp_bcast below")
872 : : end if
873 : :
874 : 2950 : nctab%has_tvale = .False.; nctab%has_tcore = .False.
875 : 2950 : pspcod = -1
876 : : !Get electronic temperature from dtset
877 [ + + - + : 2950 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
- + ]
878 : :
879 [ + + ]: 2950 : if (me==0) then
880 : : ! Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case
881 [ + + - + ]: 2738 : 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 [ + - ]: 2738 : write(msg, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(psps%filpsp(ipsp))
890 [ + + ]: 8214 : call wrtout([std_out, ab_out], msg)
891 : :
892 : : ! Check if the file pseudopotential file is written in (XML| XML-PAW | UPF)
893 : 2738 : 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 [ - + + - : 8214 : ABI_MALLOC(nproj,(psps%mpssoang))
+ - + - -
+ - + ]
898 [ + + ]: 10743 : nproj(:)=0
899 : :
900 [ + + + + ]: 2738 : 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 [ - + ]: 2396 : 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 [ + - ]: 2396 : rewind (unit=tmp_unit,err=10,iomsg=errmsg)
909 : :
910 : : ! Read and write some description of file from first line (character data)
911 [ + - ]: 2396 : read (tmp_unit,'(a)',err=10,iomsg=errmsg) title
912 : 2396 : write(msg, '(a,a)' ) '- ',trim(title)
913 [ + + ]: 7188 : call wrtout([std_out, ab_out], msg)
914 : :
915 : : ! Read and write more data describing psp parameters
916 [ + - ]: 2396 : read (tmp_unit,*,err=10,iomsg=errmsg) znucl,zion,pspdat
917 : 2396 : write(msg,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat'
918 [ + + ]: 7188 : call wrtout([std_out, ab_out], msg)
919 : :
920 [ + - ]: 2396 : read (tmp_unit,*,err=10,iomsg=errmsg) pspcod,pspxc,lmax,lloc,mmax,r2well
921 [ + + ]: 2396 : if(pspxc<0) then
922 : : write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) &
923 : 194 : pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
924 : : else
925 : : write(msg, '(4i5,i10,f10.5,t47,a)' ) &
926 : 2202 : pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
927 : : end if
928 [ + + ]: 7188 : call wrtout([std_out, ab_out], msg)
929 : :
930 [ + + + + ]: 342 : 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 [ + + + - ]: 327 : else if (usexml == 1 .and. xmlpaw == 1) then
967 : : write(msg,'(a,a)') &
968 [ + - ]: 316 : '- pspatm : Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp))
969 [ + + ]: 948 : call wrtout([std_out, ab_out], msg)
970 : :
971 : : ! Return header information
972 : 316 : call pawpsxml2ab(psps%filpsp(ipsp),ecut_tmp, pspheads_tmp,0)
973 [ - + ]: 316 : ABI_FREE(pspheads_tmp%nproj)
974 [ - + ]: 316 : ABI_FREE(pspheads_tmp%nprojso)
975 : 316 : lmax=pspheads_tmp%lmax
976 : 316 : pspxc=pspheads_tmp%pspxc
977 : 316 : znucl=pspheads_tmp%znuclpsp
978 : 316 : pawpsp_header%basis_size=pspheads_tmp%pawheader%basis_size
979 : 316 : pawpsp_header%l_size=pspheads_tmp%pawheader%l_size
980 : 316 : pawpsp_header%lmn_size=pspheads_tmp%pawheader%lmn_size
981 : 316 : pawpsp_header%mesh_size=pspheads_tmp%pawheader%mesh_size
982 : 316 : pawpsp_header%pawver=pspheads_tmp%pawheader%pawver
983 : 316 : pawpsp_header%shape_type=pspheads_tmp%pawheader%shape_type
984 : 316 : pawpsp_header%rpaw=pspheads_tmp%pawheader%rpaw
985 : 316 : pawpsp_header%rshp=pspheads_tmp%pawheader%rshp
986 : 316 : 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 [ + + ]: 2738 : (pspcod/=7).and.(pspcod/=17))then
1027 : : ! If pspso requires internal characteristics, set it to 1 for non-HGH psps
1028 [ - + ]: 1219 : if(psps%pspso(ipsp)==1) psps%pspso(ipsp)=0
1029 [ - + ]: 1219 : 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 [ - + ]: 2738 : 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 [ - + ]: 2738 : 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 [ + + ]: 2738 : if (dtset%ixc==0) then
1061 : 15 : ABI_WARNING('Note that input ixc=0 => no xc is being used.')
1062 [ + + ]: 2723 : else if(dtset%ixc/=pspxc) then
1063 : : write(msg, '(a,i0,a,i0,8a)' )&
1064 : 765 : 'Pseudopotential file pspxc: ',pspxc,', not equal to input ixc: ',dtset%ixc,'.',ch10,&
1065 : 765 : 'These parameters must agree to get the same xc in ABINIT code as in psp construction.',ch10,&
1066 : 765 : 'Action: check psp design or input file.',ch10,&
1067 : 1530 : 'Assume experienced user. Execution will continue.'
1068 : 765 : ABI_WARNING(msg)
1069 : : end if
1070 : :
1071 [ + + ]: 2738 : 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 [ + + + - : 2738 : 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 [ - + - - : 2738 : 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 [ - + ]: 2738 : 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 [ - + + - : 8214 : ABI_MALLOC(e990,(psps%mpssoang))
+ - + - -
+ - + ]
1110 [ - + + - : 5476 : ABI_MALLOC(e999,(psps%mpssoang))
- + - + ]
1111 [ - + + - : 5476 : ABI_MALLOC(rcpsp,(psps%mpssoang))
- + - + ]
1112 [ - + + - : 5476 : ABI_MALLOC(rms,(psps%mpssoang))
- + - + ]
1113 [ - + + - : 5476 : ABI_MALLOC(epspsp,(psps%mpssoang))
- + - + ]
1114 [ - + + - : 5476 : ABI_MALLOC(ekb1,(psps%mpssoang))
- + - + ]
1115 [ - + + - : 5476 : ABI_MALLOC(ekb2,(psps%mpssoang))
- + - + ]
1116 [ + + + + ]: 18748 : e990(:)=zero ;e999(:)=zero
1117 [ + + + + ]: 18748 : rcpsp(:)=zero;rms(:)=zero
1118 [ + + + + ]: 18748 : ekb1(:)=zero ;ekb2(:)=zero
1119 [ + + ]: 10743 : epspsp(:)=zero
1120 : 2738 : qchrg=zero
1121 : :
1122 : : ! ----------------------------------------------------------------------
1123 [ + + + + : 2738 : 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 : 696 : 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 : 233 : & xc_taupos=dtset%xc_taupos,el_temp=el_temp)
1165 : :
1166 : 485 : 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 : 485 : 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 : 316 : & xc_taupos=dtset%xc_taupos,el_temp=el_temp)
1213 : 316 : call paw_setup_free(paw_setuploc)
1214 : : end if
1215 : :
1216 : 2738 : close (unit=tmp_unit)
1217 : :
1218 : : ! ----------------------------------------------------------------------
1219 [ + + ]: 2738 : 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 [ + + ]: 2738 : if (pspcod/=7.and.pspcod/=17) then
1244 : 2189 : write(msg, '(a,f14.8,a,a)' ) ' pspatm : epsatm=',epsatm,ch10,' --- l ekb(1:nproj) -->'
1245 [ + + ]: 6567 : call wrtout([std_out, ab_out], msg)
1246 : : !print *, "nproj", nproj
1247 : 2189 : iln0=0
1248 [ + + ]: 12869 : do ilmn=1,psps%lmnmax
1249 : 10680 : iln=indlmn(5,ilmn)
1250 [ + + ]: 12869 : if (iln>iln0) then
1251 : 4738 : il=indlmn(1,ilmn)
1252 [ + + ]: 4738 : if (indlmn(6,ilmn)==1) then
1253 : 4033 : iln0=iln0+nproj(il+1)
1254 [ + - + + ]: 10493 : 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 : 705 : iln0=iln0+nproj(il+psps%mpsang)
1258 [ + + ]: 705 : 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 [ + + ]: 705 : if (dtset%so_psp(ipsp) == 0) cycle
1264 [ + - + + ]: 1092 : write(msg, '(2x,a,i1,4f12.6)' ) 'spin-orbit ',il,(ekb(iln+ii),ii=0,nproj(il+psps%mpsang)-1)
1265 : : end if
1266 [ + + ]: 14025 : 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 [ + + ]: 2738 : if (psps%usepaw == 0) then
1276 : 2189 : call nctab_eval_tcorespl(nctab, psps%n1xccc, xcccrc, xccc1d, xcctau1d, psps%mqgrid_vl, psps%qgrid_vl)
1277 : : end if
1278 : :
1279 : 2738 : write(msg,'(3a)') ' pspatm: atomic psp has been read ',' and splines computed',ch10
1280 [ + + ]: 8214 : call wrtout([std_out, ab_out], msg)
1281 : :
1282 [ - + ]: 2738 : ABI_FREE(e990)
1283 [ - + ]: 2738 : ABI_FREE(e999)
1284 [ - + ]: 2738 : ABI_FREE(rcpsp)
1285 [ - + ]: 2738 : ABI_FREE(rms)
1286 [ - + ]: 2738 : ABI_FREE(ekb1)
1287 [ - + ]: 2738 : ABI_FREE(ekb2)
1288 [ - + ]: 2738 : ABI_FREE(epspsp)
1289 [ - + ]: 2738 : ABI_FREE(nproj)
1290 : :
1291 [ + + + + : 5476 : 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 [ + + ]: 2950 : if (paral_mode==1) then
1333 : 284 : call timab(48,1,tsec)
1334 : 284 : call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc)
1335 : 284 : call timab(48,2,tsec)
1336 : : end if
1337 : :
1338 [ + + ]: 2950 : if (psps%usepaw==1) then
1339 [ + + + + ]: 62116 : indlmn(:,:)=0
1340 [ + + + + ]: 58910 : indlmn(1:6,1:pawtab%lmn_size)=pawtab%indlmn(1:6,1:pawtab%lmn_size)
1341 : : end if
1342 : :
1343 : : !--------------------------------------------------------------------
1344 : : !WVL+PAW:
1345 [ + + + - : 2950 : 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 [ - + - + : 2950 : 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 : : !!***
|
