Branch data Line data Source code
1 : : !!****m* ABINIT/m_pspheads
2 : : !! NAME
3 : : !! m_pspheads
4 : : !!
5 : : !! FUNCTION
6 : : !! Functions used to read the pseudopotential header of each psp file,
7 : : !! in order to initialize pspheads(1:npsp).
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, FrD, AF, MT, FJ, MJV, MG, DRH)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : MODULE m_pspheads
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_errors
28 : : use m_xmpi
29 : : use m_atomdata
30 : : use m_hash_md5
31 : : use m_psxml2ab
32 : : #if defined HAVE_LIBPSML
33 : : use m_psml
34 : : #endif
35 : : #if defined HAVE_BIGDFT
36 : : use BigDFT_API, only: atomic_info, psp_from_data
37 : : #endif
38 : :
39 : : use defs_datatypes, only : pspheader_type
40 : : use m_time, only : timab
41 : : use m_io_tools, only : open_file
42 : : use m_numeric_tools,only : simpson
43 : : use m_fstrings, only : basename, lstrip, sjoin, startswith, atoi, itoa, ftoa, toupper, next_token
44 : : use m_pawpsp, only : pawpsp_read_header_xml,pawpsp_read_pawheader
45 : : use m_pawxmlps, only : rdpawpsxml,rdpawpsxml_header, paw_setup_free,paw_setuploc
46 : : use pseudo_types, only : pseudo_upf, deallocate_pseudo_upf !, pseudo_config
47 : : use read_upf_new_module, only : read_upf_new
48 : :
49 : : implicit none
50 : :
51 : : private
52 : : !!***
53 : :
54 : : public :: inpspheads ! Initialize pspheads(1:npsp).
55 : : public :: pspheads_comm ! Communicate pspheads to all processors
56 : : public :: pawpsxml2ab
57 : :
58 : : public :: upf2_jl2srso
59 : : public :: upfxc2abi ! UPF XcC to Abinit pspxc (DEPRECATED. Only used for UPF1)
60 : : public :: upfdft_to_ixc ! UPF2 dft to Abinit pspxc.
61 : :
62 : : contains
63 : : !!***
64 : :
65 : : !!****f* m_pspheads/inpspheads
66 : : !! NAME
67 : : !! inpspheads
68 : : !!
69 : : !! FUNCTION
70 : : !! Read the pseudopotential header of each psp file, in order to initialize pspheads(1:npsp).
71 : : !!
72 : : !! INPUTS
73 : : !! npsp=number of pseudopotentials
74 : : !!
75 : : !! OUTPUT
76 : : !! pspheads(npsp)=<type pspheader_type>=all the important information from the
77 : : !! pseudopotential file headers, as well as the psp file names
78 : : !! ecut_tmp(3,2,npsp)= possible ecut values as read in psp files
79 : : !!
80 : : !! SOURCE
81 : :
82 : 1176 : subroutine inpspheads(filnam, npsp, pspheads, ecut_tmp)
83 : :
84 : : !Arguments ------------------------------------
85 : : !scalars
86 : : integer,intent(in) :: npsp
87 : : !arrays
88 : : real(dp),intent(inout) :: ecut_tmp(3,2,10)
89 : : character(len=fnlen), intent(in) :: filnam(npsp)
90 : : type(pspheader_type),intent(inout) :: pspheads(npsp) !vz_i
91 : :
92 : : !Local variables-------------------------------
93 : : !In case a xc core correction is to be taken into account,
94 : : !the n1xccc value will be given by n1xccc_default. Otherwise it is set to 0.
95 : : !scalars
96 : : integer,parameter :: n1xccc_default=2501
97 : : integer :: extension_switch
98 : : integer :: idum,ii,ilmax,ipsp,lang,lmax,mmax,mpsang,n1xccc,nmesh
99 : : integer :: pspcod,pspso,test_paw,usexml,unt,useupf
100 : : real(dp) :: al,e990,e999,fchrg,qchrg,r1,rchrg,rr
101 : : character(len=3) :: testxc
102 : : character(len=500) :: msg,errmsg
103 : : character(len=70) :: testxml
104 : : character(len=80) :: pspline
105 : : !arrays
106 : 1176 : integer,allocatable :: nproj(:),nprojso(:)
107 : 1176 : integer,allocatable :: orb(:)
108 : : real(dp) :: hdum(3)
109 : : #if defined HAVE_BIGDFT
110 : : !new variables for wvl+paw
111 : : character(len=2) :: symbol
112 : : integer :: iasctype,nzatom, nelpsp, npspcode_,ixc_
113 : : real(dp) :: rcov,ehomo
114 : : real(dp) :: psppar(0:4,0:6)
115 : : logical :: exists
116 : : #endif
117 : : #if defined HAVE_LIBPSML
118 : : character(len=3) :: atmsymb
119 : : character(len=30) :: creator
120 : : #endif
121 : : !*************************************************************************
122 : :
123 : 1176 : test_paw=0
124 : :
125 [ + + ]: 2767 : do ipsp=1,npsp
126 : :
127 : 1591 : pspheads(ipsp)%filpsp=trim(filnam(ipsp))
128 : :
129 : : ! Check if the file is written in XML
130 : 1591 : usexml = 0
131 [ - + ]: 1591 : if (open_file(filnam(ipsp), msg, newunit=unt, form="formatted", status="old") /= 0) then
132 : 0 : ABI_ERROR(msg)
133 : : end if
134 : :
135 [ + - ]: 1591 : rewind(unit=unt, err=10, iomsg=errmsg)
136 [ + - ]: 1591 : read(unt, "(a)", err=10, iomsg=errmsg) testxml
137 : :
138 [ + + ]: 1591 : if(testxml(1:5)=='<?xml')then
139 : 208 : usexml = 1
140 [ + - ]: 208 : read(unt,*, err=10, iomsg=errmsg) testxml
141 [ + + ]: 208 : if(testxml(1:4)=='<paw')then
142 : : test_paw = 1
143 : : else
144 : 11 : test_paw = 0
145 : : end if
146 : : else
147 : : usexml = 0
148 : : end if
149 : :
150 : : ! Check if pseudopotential file is a QE UPF2 file
151 : : ! "<UPF version="2.0.1">
152 : 1591 : useupf = 0
153 [ + + ]: 1591 : if (testxml(1:4) == '<UPF') then
154 : 6 : ii = index(testxml, '"')
155 [ + - ]: 6 : if (ii /= 0) then
156 : 6 : useupf = atoi(testxml(ii+1:ii+1))
157 : : else
158 : 0 : ABI_ERROR(sjoin("Cannot find version attribute in UPF2 file:", filnam(ipsp)))
159 : : end if
160 : : end if
161 : :
162 [ + - ]: 1591 : close(unit=unt, err=10, iomsg=errmsg)
163 : :
164 : : ! Check if pseudopotential file is a QE UPF1 file
165 [ + + ]: 1591 : if (useupf == 0) then
166 [ - + ]: 1585 : if (open_file(filnam(ipsp), msg, newunit=unt, form="formatted", status="old") /= 0) then
167 : 0 : ABI_ERROR(msg)
168 : : end if
169 : :
170 [ + - ]: 1585 : rewind(unit=unt, err=10, iomsg=errmsg)
171 [ + - ]: 1585 : read(unt,*,err=10,iomsg=errmsg) testxml ! just a string, no relation to xml.
172 [ + + ]: 1585 : if (testxml(1:9)=='<PP_INFO>') then
173 : : useupf = 1
174 : : else
175 : 1584 : useupf = 0
176 : : end if
177 [ + - ]: 1585 : close(unit=unt,err=10,iomsg=errmsg)
178 : : end if
179 [ + + ]: 1591 : if(test_paw==0) then
180 [ - + + + ]: 6655 : ABI_ICALLOC(nproj,(0:3))
181 [ - + + + ]: 5324 : ABI_ICALLOC(nprojso,(3))
182 : : endif
183 : :
184 : : ! Read the header of the pseudopotential file
185 [ + + ]: 1591 : if (usexml /= 1 .and. useupf == 0) then
186 : : ! Open the psp file and read a normal abinit style header
187 [ - + ]: 1376 : if (open_file(filnam(ipsp), msg, newunit=unt, form='formatted', status='old') /= 0) then
188 : 0 : ABI_ERROR(msg)
189 : : end if
190 [ + - ]: 1376 : rewind (unit=unt, err=10, iomsg=errmsg)
191 : :
192 : : ! Read the three first lines
193 [ + - ]: 1376 : read(unt, '(a)', err=10, iomsg=errmsg) pspheads(ipsp)%title
194 [ + - ]: 1376 : read(unt,*, err=10, iomsg=errmsg)pspheads(ipsp)%znuclpsp,pspheads(ipsp)%zionpsp,pspheads(ipsp)%pspdat
195 [ + - ]: 1376 : read(unt,*, err=10, iomsg=errmsg)pspheads(ipsp)%pspcod,pspheads(ipsp)%pspxc,pspheads(ipsp)%lmax,idum,mmax
196 : :
197 : 1376 : pspcod=pspheads(ipsp)%pspcod
198 : 1376 : lmax=pspheads(ipsp)%lmax
199 : 1376 : write(msg,'(a,f5.1,a,i4,a,i4)')' read the values zionpsp=',pspheads(ipsp)%zionpsp,' , pspcod=',pspcod,' , lmax=',lmax
200 : 1376 : call wrtout(std_out,msg,'PERS')
201 : :
202 [ + + ]: 1376 : if(test_paw==1) then
203 [ - + + - : 189 : ABI_MALLOC(nproj,(0:lmax))
+ - + - -
+ - + -
+ ]
204 [ + + + - : 189 : ABI_MALLOC(nprojso,(lmax))
+ - + + -
+ - + -
+ ]
205 [ + + ]: 213 : nproj(0:lmax)=0
206 [ + + ]: 150 : nprojso(1:lmax)=0
207 : : else
208 [ + + + + ]: 10504 : nproj(0:3)=0 ; nprojso(1:3)=0
209 : : endif
210 : :
211 : 1376 : pspheads(ipsp)%xccc=0
212 : 1376 : pspheads(ipsp)%pspso=0
213 : :
214 [ + + ]: 215 : else if (usexml==1 .and. test_paw==0) then
215 : : #if defined HAVE_LIBPSML
216 : : write(msg,'(4a)') &
217 : 11 : '- inpspheads : Reading pseudopotential header in XML form from ',ch10,&
218 [ + - ]: 22 : & '- ',trim(filnam(ipsp))
219 [ + + ]: 33 : call wrtout([std_out, ab_out], msg)
220 : :
221 : : ! could pass pspheads(ipsp) directly and fill all of it in psxml2ab
222 : 11 : call psxml2abheader( filnam(ipsp), pspheads(ipsp), atmsymb, creator, 1)
223 : :
224 : : ! save some stuff locally for this ipsp
225 : 11 : pspcod = pspheads(ipsp)%pspcod
226 : 11 : lmax = pspheads(ipsp)%lmax
227 [ + + ]: 50 : nproj(0:min(lmax,3)) = pspheads(ipsp)%nproj(0:min(lmax,3))
228 [ + + ]: 39 : nprojso(1:min(lmax,3)) = pspheads(ipsp)%nprojso(1:min(lmax,3))
229 : :
230 : : #else
231 : : write(msg, '(2a)') "XML norm-conserving pseudopotential has been input,", &
232 : : " but abinit is not compiled with libPSML support. Reconfigure and recompile."
233 : : ABI_ERROR(msg)
234 : : #endif
235 : :
236 [ + + ]: 204 : else if(usexml==1.and.test_paw==1)then
237 : :
238 : : write(msg,'(4a)') &
239 : 197 : '- inpspheads : Reading pseudopotential header in XML form from ',ch10,&
240 [ + - ]: 394 : '- ',trim(filnam(ipsp))
241 [ + + ]: 591 : call wrtout([std_out, ab_out], msg)
242 : :
243 : 197 : call pawpsxml2ab(filnam(ipsp),ecut_tmp(:,:,ipsp), pspheads(ipsp),1)
244 : 197 : pspcod=17; pspheads(ipsp)%pspcod=pspcod
245 : :
246 [ - + ]: 7 : else if (useupf > 0) then
247 : 7 : pspheads(ipsp)%xccc = n1xccc_default ! will be set to 0 if no nlcc
248 : :
249 [ + + ]: 7 : if (useupf == 1) then
250 : 1 : pspheads(ipsp)%pspcod = 11
251 : : call upf1_to_psphead(filnam(ipsp), pspheads(ipsp)%znuclpsp, pspheads(ipsp)%zionpsp, pspheads(ipsp)%pspxc, &
252 : 1 : pspheads(ipsp)%lmax, pspheads(ipsp)%xccc, nproj, nprojso)
253 : :
254 : : ! FIXME: generalize for SO pseudos
255 : 1 : pspheads(ipsp)%pspso = 0
256 : :
257 : : else
258 : : ! UPF2 format
259 : 6 : pspheads(ipsp)%pspcod = 12
260 : : call upf2_to_psphead(filnam(ipsp), pspheads(ipsp)%znuclpsp, pspheads(ipsp)%zionpsp, pspheads(ipsp)%pspxc, &
261 : 6 : pspheads(ipsp)%lmax, pspheads(ipsp)%xccc, nproj, nprojso)
262 : :
263 [ + + + + : 22 : pspheads(ipsp)%pspso = merge(2, 0, any(nprojso > 0))
+ + ]
264 : : end if
265 : :
266 : 7 : pspcod = pspheads(ipsp)%pspcod
267 : 7 : lmax = pspheads(ipsp)%lmax
268 : : end if
269 : :
270 : : !write(std_out,*) "pspheads(ipsp)%znuclpsp", pspheads(ipsp)%znuclpsp
271 : : !write(std_out,*) "pspheads(ipsp)%zionpsp", pspheads(ipsp)%zionpsp
272 : : !write(std_out,*) "pspheads(ipsp)%pspcod", pspheads(ipsp)%pspcod
273 : : !write(std_out,*) "pspheads(ipsp)%pspxc", pspheads(ipsp)%pspxc
274 : : !write(std_out,*) "pspheads(ipsp)%lmax", pspheads(ipsp)%lmax
275 : :
276 : : ! Initialize nproj, nprojso, pspso, as well as xccc, for each type of psp
277 [ + + ]: 9546 : pspheads(ipsp)%GTHradii = zero
278 [ + + - + ]: 1591 : if(test_paw/=1.and.lmax>3) then
279 : 0 : ABI_ERROR('NC with l>3 not implmented')
280 : : endif
281 [ + + ]: 1591 : if(.not.allocated(pspheads(ipsp)%nproj)) then
282 [ + + + - : 4149 : ABI_MALLOC(pspheads(ipsp)%nproj,(0:lmax))
+ - + + -
+ - + -
+ ]
283 : : endif
284 [ + + ]: 1591 : if(.not.allocated(pspheads(ipsp)%nprojso)) then
285 [ + + + - : 4149 : ABI_MALLOC(pspheads(ipsp)%nprojso,(lmax))
+ - + + -
+ - + -
+ ]
286 : : endif
287 [ + + ]: 1591 : if(.not.allocated(nproj)) then
288 [ - + + - : 591 : ABI_MALLOC(nproj,(0:lmax))
+ - + - -
+ - + -
+ ]
289 : : endif
290 [ + + ]: 1591 : if(.not.allocated(nprojso)) then
291 [ - + + - : 591 : ABI_MALLOC(nprojso,(lmax))
+ - + - -
+ - + -
+ ]
292 : : endif
293 : :
294 [ + + + + : 1591 : if (pspcod==1 .or. pspcod==4)then
+ + + + +
+ ]
295 : :
296 : : ! Teter format
297 [ + + ]: 1590 : do ilmax=0,lmax
298 [ + - ]: 1124 : read(unt,*, err=10, iomsg=errmsg) lang,e990,e999,nproj(ilmax)
299 [ + - ]: 1590 : read(unt,*, err=10, iomsg=errmsg)
300 : : end do
301 [ + - ]: 466 : read(unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
302 [ + + ]: 466 : if (fchrg>1.d-15) pspheads(ipsp)%xccc=n1xccc_default
303 : :
304 : : else if (pspcod==2)then
305 : :
306 : : ! GTH pseudopotentials
307 [ + - ]: 132 : read(unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%GTHradii(0) !rloc
308 [ + - ]: 132 : read(unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%GTHradii(1),hdum(1),hdum(2)
309 [ + + ]: 132 : if(abs(hdum(1))>1.d-9) nproj(0)=1
310 [ + + ]: 132 : if(abs(hdum(2))>1.d-9) nproj(0)=2
311 [ + - ]: 132 : read(unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%GTHradii(2),hdum(3)
312 [ + + ]: 132 : if(abs(hdum(3))>1.d-9) nproj(1)=1
313 : :
314 : : else if (pspcod==3)then
315 : :
316 : : ! HGH pseudopotentials
317 [ + - ]: 125 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%GTHradii(0) !rloc
318 [ + + ]: 434 : do ilmax=0,lmax
319 [ + - ]: 309 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%GTHradii(ilmax + 1),hdum(1),hdum(2),hdum(3)
320 [ + + ]: 309 : if (abs(hdum(1))>1.d-9)nproj(ilmax)=1
321 [ + + ]: 309 : if (abs(hdum(2))>1.d-9)nproj(ilmax)=2
322 [ + + ]: 309 : if (abs(hdum(3))>1.d-9)nproj(ilmax)=3
323 [ + + ]: 309 : if (ilmax>0.and.ilmax<3) then
324 [ + - ]: 182 : read (unt,*, err=10, iomsg=errmsg) hdum(1),hdum(2),hdum(3)
325 [ + - ]: 182 : if (abs(hdum(1))>1.d-9)nprojso(ilmax)=1
326 [ + + ]: 182 : if (abs(hdum(2))>1.d-9)nprojso(ilmax)=2
327 [ + + ]: 182 : if (abs(hdum(3))>1.d-9)nprojso(ilmax)=3
328 [ + - ]: 182 : if(nprojso(ilmax)>0)pspheads(ipsp)%pspso=2
329 : : end if
330 [ + + ]: 434 : if (ilmax==3) then
331 [ + - ]: 2 : read (unt,*, err=10, iomsg=errmsg) hdum(1)
332 [ + - ]: 2 : if (abs(hdum(1))>1.d-9)nprojso(3)=1
333 [ - + ]: 2 : if(nprojso(3)>0)pspheads(ipsp)%pspso=2
334 : : end if
335 : : end do
336 : :
337 : : else if (pspcod==5)then
338 : :
339 : : ! PHONEY pseudopotentials
340 : : ! read parameter for Hamman grid
341 : 65 : pspso=1
342 [ + - ]: 65 : read (unt,fmt=*,err=50,end=50) r1,al,pspso
343 : : 50 continue
344 [ + + ]: 246 : do ilmax=0,lmax
345 [ + - ]: 181 : read (unt,*, err=10, iomsg=errmsg) lang,e990,e999,nproj(ilmax)
346 [ + - ]: 181 : read (unt,*, err=10, iomsg=errmsg)
347 [ + + + + ]: 246 : if (ilmax>0.and.pspso/=1) then
348 [ + - ]: 2 : read (unt,*, err=10, iomsg=errmsg) lang,e990,e999,nprojso(ilmax)
349 [ + - ]: 2 : read (unt,*, err=10, iomsg=errmsg)
350 : 2 : pspheads(ipsp)%pspso=pspso
351 : : ! Meaning of pspso internally to ABINIT has been changed in v5.4
352 : : ! So file must contain pspso 1, but ABINIT will have pspso 0.
353 : : if(pspso==1)pspheads(ipsp)%pspso=0
354 : : end if
355 : : end do
356 [ + - ]: 65 : read (unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
357 [ + + ]: 65 : if (fchrg>1.d-15) pspheads(ipsp)%xccc=n1xccc_default
358 : :
359 : : else if (pspcod==6)then
360 : :
361 : : ! FHI pseudopotentials
362 : 94 : read (unt, '(a3)') testxc
363 : : ! Note: prior to version 2.2, this 4th line started with 4-- ,
364 : : ! and no core-correction was available.
365 [ + + ]: 94 : if(testxc/='4--')then
366 [ + - ]: 55 : backspace(unt, err=10, iomsg=errmsg)
367 [ + - ]: 55 : read (unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
368 : : else
369 : 39 : fchrg=0.0_dp
370 : : end if
371 [ + + ]: 94 : if (fchrg>1.d-15) pspheads(ipsp)%xccc=n1xccc_default
372 : : ! XG020728 : Should take lloc into account ??
373 [ + + ]: 345 : do ilmax=0,lmax
374 : 345 : nproj(ilmax)=1
375 : : end do
376 : :
377 : : else if (pspcod==7)then
378 : :
379 : : ! PAW pseudopotentials
380 : 199 : test_paw=1;pspheads(ipsp)%pawheader%pawver=1
381 [ + - ]: 199 : read (unt,'(a80)', err=10, iomsg=errmsg) pspline
382 : 199 : pspline=adjustl(pspline)
383 [ - + - - ]: 199 : if (pspline(1:3)=="paw".or.pspline(1:3)=="PAW") &
384 [ + - ]: 199 : & read(unit=pspline(4:80),fmt=*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%pawver
385 [ - + ]: 199 : if (pspheads(ipsp)%pawheader%pawver==1) then ! Compatibility with Abinit v4.2.x
386 : 0 : read (unit=pspline,fmt=*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%basis_size,&
387 [ # # ]: 0 : pspheads(ipsp)%pawheader%lmn_size
388 [ # # # # : 0 : ABI_MALLOC(orb,(pspheads(ipsp)%pawheader%basis_size))
# # # # #
# # # ]
389 [ # # ]: 0 : orb(:)=0
390 [ # # ]: 0 : read (unt,*, err=10, iomsg=errmsg) (orb(ii), ii=1,pspheads(ipsp)%pawheader%basis_size)
391 [ # # ]: 0 : read (unt,*, err=10, iomsg=errmsg)
392 [ # # ]: 0 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%rpaw
393 : 0 : pspheads(ipsp)%pawheader%rshp=pspheads(ipsp)%pawheader%rpaw
394 [ # # ]: 0 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%mesh_size
395 [ # # ]: 0 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%shape_type
396 [ # # ]: 0 : if (pspheads(ipsp)%pawheader%shape_type==3) pspheads(ipsp)%pawheader%shape_type=-1
397 : : else
398 [ + - ]: 199 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%basis_size,pspheads(ipsp)%pawheader%lmn_size
399 [ - + + - : 597 : ABI_MALLOC(orb,(pspheads(ipsp)%pawheader%basis_size))
+ - + - -
+ - + ]
400 [ + + ]: 1176 : orb(:)=0
401 [ + - ]: 199 : read (unt,*, err=10, iomsg=errmsg) (orb(ii), ii=1,pspheads(ipsp)%pawheader%basis_size)
402 : 199 : pspheads(ipsp)%pawheader%mesh_size=mmax
403 [ + - ]: 199 : read (unt,*, err=10, iomsg=errmsg) nmesh
404 [ + + ]: 945 : do ii=1,nmesh
405 [ + - ]: 945 : read(unt,*, err=10, iomsg=errmsg)
406 : : end do
407 [ + - ]: 199 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%rpaw
408 : 199 : pspheads(ipsp)%pawheader%rshp=pspheads(ipsp)%pawheader%rpaw
409 [ + - ]: 199 : read (unt,'(a80)', err=10, iomsg=errmsg) pspline
410 : 199 : pspline=adjustl(pspline); write(std_out,*) pspline
411 : 199 : read(unit=pspline,fmt=*) pspheads(ipsp)%pawheader%shape_type
412 [ + + - + ]: 199 : if (pspheads(ipsp)%pawheader%pawver==2.and.&
413 : 0 : pspheads(ipsp)%pawheader%shape_type==3) pspheads(ipsp)%pawheader%shape_type=-1
414 [ + + + + ]: 199 : if (pspheads(ipsp)%pawheader%pawver>=3.and.pspheads(ipsp)%pawheader%shape_type==-1) then
415 [ + - ]: 2 : rr=zero;read(unit=pspline,fmt=*,err=20,end=20) ii,rr
416 : : 20 continue
417 [ + - ]: 2 : if (rr>=tol8) pspheads(ipsp)%pawheader%rshp=rr
418 : : end if
419 : : end if
420 [ + + ]: 698 : do ilmax=0,lmax
421 [ + + ]: 3322 : do ii=1,pspheads(ipsp)%pawheader%basis_size
422 [ + + ]: 3123 : if(orb(ii)==ilmax) nproj(ilmax)=nproj(ilmax)+1
423 : : end do
424 : : end do
425 [ + + ]: 1176 : pspheads(ipsp)%pawheader%l_size=2*maxval(orb)+1
426 : 199 : pspheads(ipsp)%xccc=1 ! We suppose apriori that cc is used (but n1xccc is not used in PAW)
427 [ - + ]: 199 : ABI_FREE(orb)
428 : :
429 : : #if defined HAVE_BIGDFT
430 : : ! WVL+PAW case, need to define GTHradii
431 : : if(pspheads(ipsp)%usewvl==1) then
432 : : ! Obtain the HGH parameters by default from BigDFT
433 : :
434 : : call atomic_info(int(pspheads(ipsp)%znuclpsp), int(pspheads(ipsp)%zionpsp), &
435 : : symbol = symbol, ehomo = ehomo, rcov = rcov, nsccode = iasctype)
436 : :
437 : : ! I use the XC: Perdew, Burke & Ernzerhof as default, since
438 : : ! other XC potentials may not be in the BigDFT table.
439 : : ixc_=1
440 : : call psp_from_data(symbol, nzatom, nelpsp, npspcode_, ixc_, psppar, exists)
441 : : if(.not. exists) then
442 : : write(msg,'(4a)')ch10,&
443 : : "Chemical element not found in BigDFT table",ch10,&
444 : : "Action: upgrade BigDFT table"
445 : : ABI_BUG(msg)
446 : : end if
447 : : !
448 : : ! pspheads(ipsp)%pawheader%rpaw/4.0d0
449 : : pspheads(ipsp)%GTHradii(0)=psppar(0,0) !rloc
450 : : pspheads(ipsp)%GTHradii(1)=psppar(1,0) !rrs
451 : : pspheads(ipsp)%GTHradii(2)=psppar(2,0) !rrp
452 : : ! pspheads(ipsp)%GTHradii(1) = one / sqrt(abs(two * ehomo))
453 : : ! write(*,*)pspheads(ipsp)%GTHradii(:)
454 : : end if
455 : : #endif
456 : :
457 : : else if (pspcod==8)then
458 : :
459 : : ! DRH pseudopotentials
460 [ + - ]: 289 : read(unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
461 [ + + ]: 289 : if (fchrg>1.d-15) pspheads(ipsp)%xccc=n1xccc_default
462 [ + - ]: 289 : read(unt,*, err=10, iomsg=errmsg) nproj(0:lmax)
463 [ + - ]: 289 : read(unt,*, err=10, iomsg=errmsg) extension_switch
464 [ + + ]: 289 : if(any(extension_switch == [2, 3])) then
465 : 66 : pspso=2
466 [ + - ]: 66 : read(unt,*,err=10,iomsg=errmsg) nprojso(1:lmax)
467 : : else
468 : 223 : pspso=0
469 : : end if
470 : 289 : pspheads(ipsp)%pspso=pspso
471 : :
472 : : else if (pspcod==9)then
473 : : ! placeholder: nothing to do everything is read above
474 : :
475 : : else if (pspcod==10)then
476 : :
477 : : ! HGH pseudopotentials, full h/k matrices
478 [ + - ]: 6 : read (unt,*,err=10,iomsg=errmsg) pspheads(ipsp)%GTHradii(0) !rloc
479 [ + - ]: 6 : read (unt,*,err=10,iomsg=errmsg) idum
480 [ - + ]: 6 : if(idum-1/=lmax) then
481 : 0 : ABI_ERROR("in inpspheads: nnonloc-1 /= lmax")
482 : : end if
483 [ + + ]: 19 : do ilmax=0,lmax
484 : : read (unt,*,err=10,iomsg=errmsg) &
485 [ + - ]: 13 : pspheads(ipsp)%GTHradii(ilmax + 1),nproj(ilmax),(hdum(idum),idum=1,nproj(ilmax))
486 [ + + ]: 17 : do idum=2,nproj(ilmax) !skip the rest of h_ij
487 [ + - ]: 17 : read (unt,*, err=10, iomsg=errmsg)
488 : : end do
489 [ + + ]: 13 : if (ilmax==0) cycle
490 : 7 : nprojso(ilmax)=nproj(ilmax)
491 [ - + ]: 13 : if(nprojso(ilmax)>0)then
492 : 7 : pspheads(ipsp)%pspso=2
493 [ + + ]: 15 : do idum=1,nprojso(ilmax) !skip the rest of k_ij
494 [ + - ]: 15 : read (unt,*, err=10, iomsg=errmsg)
495 : : end do
496 : : end if
497 : : end do
498 : :
499 [ + - ]: 204 : else if (any(pspcod == [11, 12, 17])) then
500 : : ! already done above
501 : :
502 : : else
503 : : write(msg, '(a,i0,4a)' )&
504 : 0 : 'The pseudopotential code (pspcod) read from file is ',pspcod,ch10,&
505 : 0 : 'This value is not allowed.',ch10,&
506 : 0 : 'Action: use a correct pseudopotential file.'
507 : 0 : ABI_ERROR(msg)
508 : : end if ! pspcod
509 : :
510 : : ! Store in pspheads
511 [ + + ]: 1591 : if (pspcod /= 17) then
512 [ + + ]: 4809 : pspheads(ipsp)%nproj(0:lmax)=nproj(0:lmax)
513 [ + + ]: 3417 : pspheads(ipsp)%nprojso(1:lmax)=nprojso(1:lmax)
514 : : end if
515 : : !write(std_out,'(a,*(i0,1x))') 'nproj = ', pspheads(ipsp)%nproj(:)
516 : : !write(std_out,'(a,*(i0,1x))') 'nprojso = ', pspheads(ipsp)%nprojso(:)
517 : :
518 : 1591 : close(unt)
519 : :
520 : : ! Compute md5 checksum
521 : 1591 : pspheads(ipsp)%md5_checksum = md5_sum_from_file(filnam(ipsp))
522 [ - + ]: 1591 : ABI_FREE(nproj)
523 [ - + ]: 2767 : ABI_FREE(nprojso)
524 : : end do ! ipsp=1,npsp
525 : :
526 : : ! Note that mpsang is the max of 1+lmax, with minimal value 1 (even for local psps, at present)
527 : 1176 : mpsang=1
528 : 1176 : n1xccc=pspheads(1)%xccc
529 [ + + ]: 2767 : do ii=1,npsp
530 : 1591 : mpsang=max(pspheads(ii)%lmax+1,mpsang)
531 : 2767 : n1xccc=max(pspheads(ii)%xccc,n1xccc)
532 : : end do
533 : :
534 : 1176 : write(msg,'(2a,i0,a,i0,a)')ch10,' inpspheads: deduce mpsang = ',mpsang,', n1xccc = ',n1xccc,'.'
535 : 1176 : call wrtout(std_out,msg,'PERS')
536 : :
537 : : ! Test: if one psp is PAW, all must be
538 [ + + ]: 1176 : if (test_paw==1) then
539 [ + + ]: 661 : do ipsp=1,npsp
540 [ - + ]: 661 : if (all(pspheads(ipsp)%pspcod /= [7, 17])) then
541 : : write(msg, '(5a)' )&
542 : 0 : 'One pseudopotential is PAW (pspcod=7 or 17) !',ch10,&
543 : 0 : 'All pseudopotentials must be PAW (this is not the case here) !',ch10,&
544 : 0 : 'Action: use only PAW pseudopotential files.'
545 : 0 : ABI_ERROR(msg)
546 : : end if
547 : : end do
548 : : end if
549 : :
550 : : return
551 : :
552 : : ! Handle IO error
553 : : 10 continue
554 : 0 : ABI_ERROR(errmsg)
555 : :
556 [ # # # # : 0 : end subroutine inpspheads
# # ]
557 : : !!***
558 : :
559 : : !!****f* ABINIT/pspheads_comm
560 : : !! NAME
561 : : !! pspheads_comm
562 : : !!
563 : : !! FUNCTION
564 : : !! Communicate pspheads to all processors
565 : : !!
566 : : !! INPUTS
567 : : !! npsp=number of pseudopotentials
568 : : !! test_paw=0 if no PAW, 1 if PAW
569 : : !!
570 : : !! SIDE EFFECTS
571 : : !! pspheads(npsp)=<type pspheader_type>=all the important information from the
572 : : !! pseudopotential file headers, as well as the psp file names. On one processor at input,
573 : : !! on all processors at output
574 : : !!
575 : : !! SOURCE
576 : :
577 : 1460 : subroutine pspheads_comm(npsp,pspheads,test_paw)
578 : :
579 : : !Arguments ------------------------------------
580 : : integer,intent(in) :: npsp
581 : : integer,intent(inout) :: test_paw
582 : : type(pspheader_type),intent(inout) :: pspheads(npsp)
583 : :
584 : : !Local variables-------------------------------
585 : : #if defined HAVE_MPI
586 : : !scalars
587 : : integer,parameter :: master=0
588 : : integer :: ierr,comm,ii,ipsp,il,list_size,lmax_
589 : : !arrays
590 : 1460 : integer,allocatable :: list_int(:)
591 : : real(dp) :: tsec(2)
592 : 1460 : real(dp),allocatable :: list_dpr(:)
593 : 1460 : character(len=fnlen),allocatable :: list_char(:)
594 : : #endif
595 : : !*************************************************************************
596 : :
597 : : #if defined HAVE_MPI
598 : 1460 : call timab(48,1,tsec)
599 : :
600 : 1460 : comm = xmpi_world
601 : :
602 : : ! Broadcast the characters (file names and titles)
603 [ + - - + : 4380 : ABI_MALLOC(list_char,(3*npsp))
+ - + - +
- - + - +
- + ]
604 [ + + ]: 3440 : list_char(1:npsp)=pspheads(1:npsp)%filpsp
605 [ + + ]: 3440 : list_char(npsp+1:2*npsp)=pspheads(1:npsp)%title
606 [ + + ]: 3440 : list_char(2*npsp+1:3*npsp)=pspheads(1:npsp)%md5_checksum
607 : :
608 : 1460 : call xmpi_bcast(list_char,master,comm,ierr)
609 : :
610 [ + + ]: 3440 : pspheads(1:npsp)%filpsp=list_char(1:npsp)
611 [ + + ]: 3440 : pspheads(1:npsp)%title=list_char(npsp+1:2*npsp)
612 [ + + ]: 3440 : pspheads(1:npsp)%md5_checksum=list_char(2*npsp+1:3*npsp)(1:md5_slen)
613 [ - + ]: 1460 : ABI_FREE(list_char)
614 : :
615 : : ! Brodcast the integers
616 : 1460 : list_size=6*npsp+1
617 [ - + + - : 4380 : ABI_MALLOC(list_int,(list_size))
+ - + - -
+ - + -
+ ]
618 [ + + ]: 3440 : list_int(1 : npsp) = pspheads(1:npsp)%lmax
619 [ + + ]: 3440 : list_int(1+ npsp: 2*npsp) = pspheads(1:npsp)%xccc
620 [ + + ]: 3440 : list_int(1+ 2*npsp: 3*npsp) = pspheads(1:npsp)%pspxc
621 [ + + ]: 3440 : list_int(1+ 3*npsp: 4*npsp) = pspheads(1:npsp)%pspdat
622 [ + + ]: 3440 : list_int(1+ 4*npsp: 5*npsp) = pspheads(1:npsp)%pspcod
623 [ + + ]: 3440 : list_int(1+ 5*npsp: 6*npsp) = pspheads(1:npsp)%pspso
624 : 1460 : list_int(6*npsp+1) = test_paw
625 : :
626 : 1460 : call xmpi_bcast(list_int,master,comm,ierr)
627 : :
628 [ + + ]: 3440 : pspheads(1:npsp)%lmax = list_int(1 : npsp)
629 [ + + ]: 3440 : pspheads(1:npsp)%xccc = list_int(1+ npsp: 2*npsp)
630 [ + + ]: 3440 : pspheads(1:npsp)%pspxc = list_int(1+ 2*npsp: 3*npsp)
631 [ + + ]: 3440 : pspheads(1:npsp)%pspdat = list_int(1+ 3*npsp: 4*npsp)
632 [ + + ]: 3440 : pspheads(1:npsp)%pspcod = list_int(1+ 4*npsp: 5*npsp)
633 [ + + ]: 3440 : pspheads(1:npsp)%pspso = list_int(1+ 5*npsp: 6*npsp)
634 : 1460 : test_paw = list_int(6*npsp+1)
635 : :
636 : 1460 : ABI_FREE(list_int)
637 : :
638 : 1460 : list_size=0
639 [ + + ]: 3440 : do ipsp=1,npsp
640 [ + + ]: 1980 : if(.not.allocated(pspheads(ipsp)%nproj)) then
641 [ - + + - : 1167 : ABI_MALLOC(pspheads(ipsp)%nproj,(0:pspheads(ipsp)%lmax))
+ - + - -
+ - + -
+ ]
642 : : endif
643 [ + + ]: 1980 : if(.not.allocated(pspheads(ipsp)%nprojso)) then
644 [ + + + - : 1167 : ABI_MALLOC(pspheads(ipsp)%nprojso,(pspheads(ipsp)%lmax))
+ - + + -
+ - + -
+ ]
645 : : endif
646 : 1980 : lmax_=max(pspheads(ipsp)%lmax,0)
647 : 3440 : list_size=list_size+2*lmax_+1
648 : : enddo
649 [ - + + - : 4380 : ABI_MALLOC(list_int,(list_size))
+ - + - -
+ - + ]
650 [ + + ]: 9142 : list_int=0
651 : : ii=0
652 [ + + ]: 3440 : do ipsp=1,npsp
653 [ + + ]: 8269 : do il=0,pspheads(ipsp)%lmax
654 : 4829 : ii=ii+1
655 : 6809 : list_int(ii) = pspheads(ipsp)%nproj(il)
656 : : enddo
657 : : enddo
658 [ + + ]: 3440 : do ipsp=1,npsp
659 [ + + ]: 6291 : do il=1,pspheads(ipsp)%lmax
660 : 2851 : ii=ii+1
661 : 4831 : list_int(ii) = pspheads(ipsp)%nprojso(il)
662 : : enddo
663 : : enddo
664 : :
665 : 1460 : call xmpi_bcast(list_int,master,comm,ierr)
666 : :
667 : 1460 : ii=0
668 [ + + ]: 3440 : do ipsp=1,npsp
669 [ + + ]: 8269 : do il=0,pspheads(ipsp)%lmax
670 : 4829 : ii=ii+1
671 : 6809 : pspheads(ipsp)%nproj(il) = list_int(ii)
672 : : enddo
673 : : enddo
674 [ + + ]: 3440 : do ipsp=1,npsp
675 [ + + ]: 6291 : do il=1,pspheads(ipsp)%lmax
676 : 2851 : ii=ii+1
677 : 4831 : pspheads(ipsp)%nprojso(il) = list_int(ii)
678 : : enddo
679 : : enddo
680 [ - + ]: 1460 : ABI_FREE(list_int)
681 : :
682 : : ! Unbeliveable, this cannot be sent with the others, for woopy
683 [ - + + - : 4380 : ABI_MALLOC(list_int,(npsp))
+ - + - -
+ - + ]
684 [ + + ]: 3440 : list_int(1:npsp) = pspheads(1:npsp)%usewvl
685 : 1460 : call xmpi_bcast(list_int,master,comm,ierr)
686 [ + + ]: 3440 : pspheads(1:npsp)%usewvl = list_int(1:npsp)
687 [ - + ]: 1460 : ABI_FREE(list_int)
688 : :
689 : : ! Broadcast zionpsp and znuclpsp
690 [ + - - + : 4380 : ABI_MALLOC(list_dpr,(7*npsp))
+ - + - +
- - + - +
- + ]
691 [ + + ]: 3440 : list_dpr(1 : npsp) = pspheads(1:npsp)%zionpsp
692 [ + + ]: 3440 : list_dpr(1+ npsp:2*npsp) = pspheads(1:npsp)%znuclpsp
693 [ + + ]: 3440 : list_dpr(1+2*npsp:3*npsp) = pspheads(1:npsp)%GTHradii(0)
694 [ + + ]: 3440 : list_dpr(1+3*npsp:4*npsp) = pspheads(1:npsp)%GTHradii(1)
695 [ + + ]: 3440 : list_dpr(1+4*npsp:5*npsp) = pspheads(1:npsp)%GTHradii(2)
696 [ + + ]: 3440 : list_dpr(1+5*npsp:6*npsp) = pspheads(1:npsp)%GTHradii(3)
697 [ + + ]: 3440 : list_dpr(1+6*npsp:7*npsp) = pspheads(1:npsp)%GTHradii(4)
698 : :
699 : 1460 : call xmpi_bcast(list_dpr,master,comm,ierr)
700 : :
701 [ + + ]: 3440 : pspheads(1:npsp)%zionpsp = list_dpr(1 : npsp)
702 [ + + ]: 3440 : pspheads(1:npsp)%znuclpsp = list_dpr(1+ npsp:2*npsp)
703 [ + + ]: 3440 : pspheads(1:npsp)%GTHradii(0) = list_dpr(1+2*npsp:3*npsp)
704 [ + + ]: 3440 : pspheads(1:npsp)%GTHradii(1) = list_dpr(1+3*npsp:4*npsp)
705 [ + + ]: 3440 : pspheads(1:npsp)%GTHradii(2) = list_dpr(1+4*npsp:5*npsp)
706 [ + + ]: 3440 : pspheads(1:npsp)%GTHradii(3) = list_dpr(1+5*npsp:6*npsp)
707 [ + + ]: 3440 : pspheads(1:npsp)%GTHradii(4) = list_dpr(1+6*npsp:7*npsp)
708 [ - + ]: 1460 : ABI_FREE(list_dpr)
709 : :
710 : : ! Broadcast additional integers for PAW psps (testpaw was sent, previously)
711 [ + + ]: 1460 : if (test_paw==1) then
712 [ + - - + : 1600 : ABI_MALLOC(list_int,(6*npsp))
+ - + - +
- - + - +
- + ]
713 [ + + ]: 990 : list_int(1 : npsp)=pspheads(1:npsp)%pawheader%basis_size
714 [ + + ]: 990 : list_int(1+ npsp:2*npsp)=pspheads(1:npsp)%pawheader%l_size
715 [ + + ]: 990 : list_int(1+2*npsp:3*npsp)=pspheads(1:npsp)%pawheader%lmn_size
716 [ + + ]: 990 : list_int(1+3*npsp:4*npsp)=pspheads(1:npsp)%pawheader%mesh_size
717 [ + + ]: 990 : list_int(1+4*npsp:5*npsp)=pspheads(1:npsp)%pawheader%pawver
718 [ + + ]: 990 : list_int(1+5*npsp:6*npsp)=pspheads(1:npsp)%pawheader%shape_type
719 : :
720 : 400 : call xmpi_bcast(list_int,master,comm,ierr)
721 : :
722 [ + + ]: 990 : pspheads(1:npsp)%pawheader%basis_size=list_int(1 : npsp)
723 [ + + ]: 990 : pspheads(1:npsp)%pawheader%l_size =list_int(1+ npsp:2*npsp)
724 [ + + ]: 990 : pspheads(1:npsp)%pawheader%lmn_size =list_int(1+2*npsp:3*npsp)
725 [ + + ]: 990 : pspheads(1:npsp)%pawheader%mesh_size =list_int(1+3*npsp:4*npsp)
726 [ + + ]: 990 : pspheads(1:npsp)%pawheader%pawver =list_int(1+4*npsp:5*npsp)
727 [ + + ]: 990 : pspheads(1:npsp)%pawheader%shape_type=list_int(1+5*npsp:6*npsp)
728 [ - + ]: 400 : ABI_FREE(list_int)
729 : :
730 : : ! broadcast rpaw values
731 [ + - - + : 1600 : ABI_MALLOC(list_dpr,(2*npsp))
+ - + - +
- - + - +
- + ]
732 : :
733 [ + + ]: 990 : list_dpr(1 : npsp) = pspheads(1:npsp)%pawheader%rpaw
734 [ + + ]: 990 : list_dpr(1+1*npsp:2*npsp) = pspheads(1:npsp)%pawheader%rshp
735 : :
736 : 400 : call xmpi_bcast(list_dpr,master,comm,ierr)
737 : :
738 [ + + ]: 990 : pspheads(1:npsp)%pawheader%rpaw = list_dpr(1 : npsp)
739 [ + + ]: 990 : pspheads(1:npsp)%pawheader%rshp = list_dpr(1+ npsp:2*npsp)
740 : :
741 [ - + ]: 400 : ABI_FREE(list_dpr)
742 : : end if
743 : :
744 : 1460 : call timab(48,2,tsec)
745 : :
746 : : #else
747 : : ! Code to use unused dummy arguments
748 : : if(pspheads(1)%lmax == -10) pspheads(1)%lmax=-10
749 : : if(test_paw == -1) test_paw = -1
750 : : #endif
751 : :
752 [ - + - + : 1460 : end subroutine pspheads_comm
- + ]
753 : : !!***
754 : :
755 : : !!****f* m_pspheads/pawpsxml2ab
756 : : !! NAME
757 : : !! pawpsxml2ab
758 : : !!
759 : : !! FUNCTION
760 : : !! From a XML format pseudopotential file which has already been read in,
761 : : !! convert to abinit internal datastructures.
762 : : !!
763 : : !! INPUTS
764 : : !! ecut_tmp(3,2)= possible ecut values as read in psp files
765 : : !! filenam= input file name (atomicdata XML)
766 : : !! option= 1 if header only is read; 0 if the whole data are read
767 : : !!
768 : : !! OUTPUT
769 : : !! pspheads data structure is filled
770 : : !!
771 : : !! SOURCE
772 : :
773 : 515 : subroutine pawpsxml2ab(filnam, ecut_tmp, pspheads, option)
774 : :
775 : : !Arguments ------------------------------------
776 : : !scalars
777 : : integer, intent(in) :: option
778 : : character(len=fnlen), intent(in) :: filnam
779 : : type(pspheader_type),intent(inout) :: pspheads !vz_i
780 : : !arrays
781 : : real(dp),intent(inout) :: ecut_tmp(3,2)
782 : :
783 : : !Local variables-------------------------------
784 : : integer :: ii,il,lloc,lmax,pspcod,pspxc
785 : : real(dp) :: r2well,zionpsp,znuclpsp
786 : : ! character(len=100) :: xclibxc, msg
787 : : ! *********************************************************************
788 : :
789 [ + + ]: 515 : if (option==1) then
790 : 197 : call rdpawpsxml_header(ecut_tmp,filnam,paw_setuploc)
791 : 197 : paw_setuploc%idgrid= paw_setuploc%radial_grid(1)%id
792 : : else
793 : 318 : call rdpawpsxml(filnam,paw_setuploc)
794 : : end if
795 : :
796 : 515 : call pawpsp_read_header_xml(lloc,lmax,pspcod, pspxc,paw_setuploc,r2well,zionpsp,znuclpsp)
797 : :
798 : 515 : pspheads%lmax=lmax
799 : 515 : pspheads%pspxc=pspxc
800 : 515 : pspheads%zionpsp=zionpsp
801 : 515 : pspheads%znuclpsp=znuclpsp
802 : :
803 : : call pawpsp_read_pawheader(pspheads%pawheader%basis_size,&
804 : : pspheads%lmax,pspheads%pawheader%lmn_size,&
805 : : pspheads%pawheader%l_size,pspheads%pawheader%mesh_size,&
806 : : pspheads%pawheader%pawver,paw_setuploc,pspheads%pawheader%rpaw,&
807 : 515 : pspheads%pawheader%rshp,pspheads%pawheader%shape_type)
808 : :
809 [ - + + - : 1545 : ABI_MALLOC(pspheads%nproj,(0:lmax))
+ - + - -
+ - + -
+ ]
810 [ + + + - : 1545 : ABI_MALLOC(pspheads%nprojso,(lmax))
+ - + + -
+ - + -
+ ]
811 [ + + ]: 1652 : pspheads%nproj=0
812 [ + + ]: 1652 : do il=0,pspheads%lmax
813 [ + + ]: 6969 : do ii=1,pspheads%pawheader%basis_size
814 [ + + ]: 6454 : if(paw_setuploc%valence_states%state(ii)%ll==il) pspheads%nproj(il)=pspheads%nproj(il)+1
815 : : end do
816 : : end do
817 : :
818 [ + + ]: 1137 : pspheads%nprojso=0
819 : 515 : pspheads%pspdat=27061961
820 : 515 : pspheads%pspso=1
821 : 515 : pspheads%xccc=1
822 : 515 : pspheads%title=paw_setuploc%atom%symbol
823 : :
824 [ + + ]: 515 : if (option==1) call paw_setup_free(paw_setuploc)
825 : :
826 : 515 : end subroutine pawpsxml2ab
827 : : !!***
828 : :
829 : : !!****f* m_pspheads/upf1_to_psphead
830 : : !! NAME
831 : : !! upf1_to_psphead
832 : : !!
833 : : !! FUNCTION
834 : : !! This routine wraps a call to a PWSCF module, which reads in
835 : : !! a UPF1 (PWSCF / Espresso) format pseudopotential, then transfers
836 : : !! data for the HEADER of abinit psps only!
837 : : !!
838 : : !! INPUTS
839 : : !! filpsp = name of file with UPF1 data
840 : : !!
841 : : !! OUTPUT
842 : : !! pspxc = index of xc functional for this pseudo
843 : : !! lmax_ = maximal angular momentum
844 : : !! znucl = charge of species nucleus
845 : : !! zion = valence charge
846 : : !! n1xccc = default number of points. Set to 0 if no nlcc is present
847 : : !! nproj_l= number of projectors for each channel
848 : : !! nprojso_l= number of projectors for each channel for SO correction projectors
849 : : !!
850 : : !! SOURCE
851 : :
852 : 3 : subroutine upf1_to_psphead(filpsp, znucl, zion, pspxc, lmax_, n1xccc, nproj_l, nprojso_l)
853 : :
854 : : use m_read_upf_pwscf, only : read_pseudo
855 : : use pseudo_pwscf ! pwscf module with all data explicit!
856 : :
857 : : !Arguments -------------------------------
858 : : character(len=fnlen), intent(in) :: filpsp
859 : : integer,intent(inout) :: n1xccc
860 : : integer,intent(out) :: pspxc, lmax_
861 : : real(dp),intent(out) :: znucl, zion
862 : : !arrays
863 : : integer,intent(out) :: nproj_l(0:3)
864 : : integer,intent(out) :: nprojso_l(1:3)
865 : :
866 : : !Local variables -------------------------
867 : : integer :: iproj, ll, iunit
868 : : character(len=500) :: msg
869 : : type(atomdata_t) :: atom
870 : :
871 : : ! *********************************************************************
872 : :
873 [ - + ]: 1 : if (open_file(filpsp, msg, newunit=iunit, status='old',form='formatted') /= 0) then
874 : 0 : ABI_ERROR(msg)
875 : : end if
876 : :
877 : : ! read in psp data to static data in pseudo module, for ipsx == 1
878 : 1 : call read_pseudo(1,iunit)
879 : 1 : close (iunit)
880 : :
881 : : ! copy over to abinit internal arrays and vars
882 : : ! FIXME: The API is broken. It does not recognize PBEsol
883 : : ! should use upfdft_to_ixc
884 : 1 : call upfxc2abi(dft(1), pspxc)
885 : 1 : lmax_ = lmax(1)
886 : 1 : call atomdata_from_symbol(atom,psd(1))
887 : 1 : znucl = atom%znucl
888 : 1 : zion = zp(1)
889 : :
890 : 1 : nproj_l = 0
891 [ + + ]: 4 : do iproj = 1, nbeta(1)
892 : 3 : ll = lll(iproj,1)
893 : 4 : nproj_l(ll) = nproj_l(ll) + 1
894 : : end do
895 : :
896 : 1 : nprojso_l = 0 !FIXME deal with so
897 : : !do iproj = 1, nbeta(1)
898 : : !nprojso_l(ll+1) = nprojso_l(ll+1) + 1
899 : : !end do
900 : :
901 [ - + ]: 1 : if (.not. nlcc(1)) n1xccc = 0
902 : :
903 : 1 : end subroutine upf1_to_psphead
904 : : !!***
905 : :
906 : : !!****f* m_pspheads/upf2_to_psphead
907 : : !! NAME
908 : : !! upf2_to_psphead
909 : : !!
910 : : !! FUNCTION
911 : : !! This routine wraps a call to a PWSCF module, which reads in
912 : : !! a UPF2 (PWSCF / Espresso) format pseudopotential, then transfers
913 : : !! data for the HEADER of abinit psps only!
914 : : !!
915 : : !! INPUTS
916 : : !! filpsp = name of file with UPF1 data
917 : : !!
918 : : !! OUTPUT
919 : : !! pspxc = index of xc functional for this pseudo
920 : : !! lmax_ = maximal angular momentum
921 : : !! znucl = charge of species nucleus
922 : : !! zion = valence charge
923 : : !! n1xccc = default number of points. Set to 0 if no nlcc is present
924 : : !! nproj_l= number of projectors for each channel
925 : : !! nprojso_l= number of projectors for each channel for SO correction projectors
926 : : !!
927 : : !! SOURCE
928 : :
929 : 6 : subroutine upf2_to_psphead(filpsp, znucl, zion, pspxc, lmax, n1xccc, nproj_l, nprojso_l)
930 : :
931 : : !Arguments -------------------------------
932 : : character(len=fnlen), intent(in) :: filpsp
933 : : integer,intent(inout) :: n1xccc
934 : : integer,intent(out) :: pspxc, lmax
935 : : real(dp),intent(out) :: znucl, zion
936 : : !arrays
937 : : integer,intent(out) :: nproj_l(0:3), nprojso_l(1:3)
938 : :
939 : : !Local variables -------------------------
940 : : integer :: ierr , iprj, ll, mmax, irad
941 : : real(dp) :: amesh, damesh
942 : : character(len=500) :: msg
943 : : logical :: linear_mesh
944 : 6 : type(pseudo_upf) :: upf
945 : : type(atomdata_t) :: atom
946 : : ! arrays
947 : 6 : real(dp),allocatable :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:)
948 : : ! *********************************************************************
949 : :
950 : : ! See also https://github.com/QEF/qeschemas/blob/master/UPF/qe_pp-0.99.xsd
951 : 6 : call read_upf_new(filpsp, upf, ierr)
952 [ - + ]: 6 : ABI_CHECK(ierr == 0, sjoin("read_upf_new returned ierr:", itoa(ierr)))
953 : :
954 : 6 : call atomdata_from_symbol(atom, upf%psd)
955 : 6 : znucl = atom%znucl
956 : 6 : zion = upf%zp
957 : 6 : lmax = upf%lmax
958 : 6 : mmax = upf%mesh
959 : :
960 : : ! Consistency check
961 [ - + ]: 6 : ABI_CHECK(upf%typ == "NC", sjoin("Only NC pseudos in UPF2 format are supported while type is:", upf%typ))
962 [ - + ]: 6 : ABI_CHECK(upfdft_to_ixc(upf%dft, pspxc, msg) == 0, msg)
963 [ - + ]: 6 : if (.not. upf%nlcc) n1xccc = 0
964 : :
965 : : ! Check that rad grid is linear starting at zero
966 : 6 : linear_mesh = .True.
967 : 6 : amesh = upf%r(2) - upf%r(1); damesh = zero
968 [ + + ]: 8730 : do irad=2,mmax-1
969 : 8730 : damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1)))
970 : : end do
971 : 6 : linear_mesh = damesh < tol8
972 : :
973 [ + - - + ]: 6 : if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then
974 : : write(msg,'(3a)')&
975 : 0 : 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
976 : 0 : 'Action: check your pseudopotential file.'
977 : 0 : ABI_ERROR(msg)
978 : : end if
979 : :
980 : 6 : nproj_l = 0; nprojso_l = 0
981 : :
982 [ + + ]: 6 : if (.not. upf%has_so) then
983 : : ! Scalar case
984 [ + + ]: 27 : do iprj=1,upf%nbeta
985 : 23 : ll = upf%lll(iprj)
986 : 27 : nproj_l(ll) = nproj_l(ll) + 1
987 : : end do
988 : :
989 : : else
990 : : ! Pseudo in j = l + s representation.
991 : 2 : call upf2_jl2srso(upf, nproj_l, nprojso_l, vsr, esr, vso, eso)
992 : :
993 [ - + ]: 2 : ABI_FREE(vsr)
994 [ - + ]: 2 : ABI_FREE(esr)
995 [ - + ]: 2 : ABI_FREE(vso)
996 [ - + ]: 2 : ABI_FREE(eso)
997 : : end if
998 : :
999 : 6 : call deallocate_pseudo_upf(upf)
1000 : :
1001 [ - + - + : 6 : end subroutine upf2_to_psphead
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
1002 : : !!***
1003 : :
1004 : : !!****f* m_pspheads/upf2_jl2srso
1005 : : !! NAME
1006 : : !! upf2_jl2srso
1007 : : !!
1008 : : !! FUNCTION
1009 : : !!
1010 : : !! INPUTS
1011 : : !!
1012 : : !! OUTPUT
1013 : : !!
1014 : : !! SOURCE
1015 : :
1016 : 6 : subroutine upf2_jl2srso(upf, nproj_l, nprojso_l, vsr, esr, vso, eso)
1017 : :
1018 : : !Arguments -------------------------------
1019 : : type(pseudo_upf),intent(in) :: upf
1020 : : !arrays
1021 : : integer,intent(out) :: nproj_l(0:3), nprojso_l(1:3)
1022 : : real(dp),allocatable,intent(out) :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:)
1023 : :
1024 : : !Local variables -------------------------
1025 : : integer :: iprj, ii, ll, l1, il, ik, lmax, mmax, mxprj
1026 : : real(dp) :: jtot, eprmin !eps_srso,
1027 : : !character(len=500) :: msg
1028 : : ! arrays
1029 : : integer :: irc6(6),nproj6(6), done_ilk(6,2)
1030 : 6 : real(dp),allocatable :: vkb(:,:,:,:), evkb(:,:,:)
1031 : : ! *********************************************************************
1032 : :
1033 : 6 : lmax = upf%lmax; mmax = upf%mesh
1034 : 6 : nproj_l = 0; nprojso_l = 0
1035 : :
1036 : : ! Pseudo in j = l + s representation.
1037 : 6 : irc6 = zero; nproj6 = zero
1038 [ + + ]: 66 : do iprj=1,upf%nbeta
1039 : 60 : ll = upf%lll(iprj)
1040 : 60 : nproj6(ll+1) = nproj6(ll+1) + 1
1041 : : !irc6(ll+1) = max(upf%kbeta(iprj), irc6(ll+1))
1042 : 66 : irc6(ll+1) = mmax
1043 : : end do
1044 : :
1045 : : ! Divide by two for l > 0 as this is sr_so_r expects.
1046 [ + + ]: 36 : nproj6(2:) = nproj6(2:) / 2
1047 [ + + ]: 42 : mxprj = maxval(nproj6)
1048 : :
1049 [ - + + - : 36 : ABI_MALLOC(vkb, (mmax,mxprj,4,2))
- + + - +
- + - + -
+ - - + -
+ ]
1050 [ - + + - : 24 : ABI_MALLOC(evkb, (mxprj,4,2))
+ - + - +
- - + -
+ ]
1051 [ - + + - : 30 : ABI_MALLOC(vsr, (mmax,2*mxprj,4))
+ - - + +
- + - + -
+ - - + -
+ - + ]
1052 [ - + + - : 24 : ABI_MALLOC(esr, (2*mxprj,4))
+ - + - +
- - + - +
- + ]
1053 [ - + + - : 18 : ABI_MALLOC(vso, (mmax,2*mxprj,4))
- + + - -
+ - + -
+ ]
1054 [ - + + - : 12 : ABI_MALLOC(eso, (2*mxprj,4))
- + - + -
+ ]
1055 : :
1056 : 6 : done_ilk = 0
1057 [ + + ]: 66 : do iprj=1,upf%nbeta
1058 : 60 : jtot = upf%jjj(iprj)
1059 : 60 : ll = upf%lll(iprj)
1060 : 60 : il = ll + 1
1061 [ + + ]: 60 : if (ll == 0) then
1062 : : ik = 1
1063 : : else
1064 : : ! l+1/2 --> ik 1, l-1/2 --> ik 2
1065 [ + + ]: 48 : if (abs(jtot - (ll + half)) < tol6) then
1066 : : ik = 1
1067 [ - + ]: 24 : else if (abs(jtot - (ll - half)) < tol6) then
1068 : : ik = 2
1069 : : else
1070 : 0 : ABI_ERROR(sjoin("Cannot detect ik index from jtot:", ftoa(jtot)))
1071 : : end if
1072 : : end if
1073 : :
1074 : 60 : done_ilk(il, ik) = done_ilk(il, ik) + 1
1075 : 60 : ii = done_ilk(il, ik)
1076 : 60 : evkb(ii,il,ik) = upf%dion(iprj,iprj) * half ! convert from Rydberg to Ha
1077 [ + + ]: 98586 : vkb(:,ii,il,ik) = upf%beta(:,iprj)
1078 : : end do
1079 : :
1080 : 6 : call sr_so_r(lmax, irc6, nproj6, upf%r, mmax, mxprj, evkb, vkb, vsr, esr, vso, eso)
1081 : :
1082 : : ! MG: This is done in oncvpsp 3.3 but not in oncvpsp4
1083 : : ! drop sr, so orthonormal projectors with neglibible coefficients
1084 : : ! modify cutoff if desired
1085 : :
1086 : 6 : eprmin=2.0d-5
1087 : 6 : write(std_out,'(/a,1p,e10.2,a)') 'Orthonormal projectors with coefficients <', &
1088 : 12 : eprmin,' Ha will be dropped'
1089 : :
1090 [ + + ]: 24 : do l1=1,lmax+1
1091 [ + + ]: 18 : if(abs(esr(3,l1))<eprmin) esr(3,l1)=0.0d0
1092 [ + + ]: 18 : if(abs(esr(4,l1))<eprmin) esr(4,l1)=0.0d0
1093 [ + + ]: 18 : if(abs(eso(3,l1))<eprmin) eso(3,l1)=0.0d0
1094 [ + + ]: 24 : if(abs(eso(4,l1))<eprmin) eso(4,l1)=0.0d0
1095 : : end do
1096 : :
1097 : : #if 0
1098 : : ! MG: This is done in oncvpsp 4 but not in oncvpsp 3.3
1099 : : ! set smallest components to zero (following the approach used in oncvpsp)
1100 : : eps_srso=1.0d-3
1101 : : do l1=1,lmax+1
1102 : : if (nproj6(l1) > 0) then
1103 : : do iprj=2,2*nproj6(l1)
1104 : : if (abs(esr(iprj,l1)) < eps_srso*abs(esr(1,l1))) esr(iprj,l1) = 0.0d0
1105 : : if (l1 == 1) cycle
1106 : : if (abs(eso(iprj,l1)) < eps_srso*abs(eso(1,l1))) eso(iprj,l1) = 0.0d0
1107 : : end do
1108 : : end if
1109 : : end do
1110 : : #endif
1111 : :
1112 : : ! set up projector number for sr_so calculations based on non-zero coefficients
1113 : : ! note that energies and projectors have been sorted sr_so_r
1114 : : ! so the relevant projectors are packed in the first positions.
1115 [ + + ]: 24 : do l1=1,lmax+1
1116 : 18 : ll = l1 - 1
1117 [ + + ]: 96 : do ii=1,2*nproj6(l1)
1118 [ + + ]: 72 : if (abs(esr(ii,l1)) > 0.0d0) nproj_l(ll) = nproj_l(ll) + 1
1119 [ + + ]: 90 : if (abs(eso(ii,l1)) > 0.0d0) nprojso_l(ll) = nprojso_l(ll) + 1
1120 : : end do
1121 : : !write(std_out, '(a,3(i0,1x))')' ll, nproj_l, nprojso_l',ll, nproj_l(ll), nprojso_l(ll)
1122 : : end do
1123 : :
1124 [ - + ]: 6 : ABI_FREE(vkb)
1125 [ - + ]: 6 : ABI_FREE(evkb)
1126 : :
1127 : 6 : end subroutine upf2_jl2srso
1128 : : !!***
1129 : :
1130 : : !!****f* m_pspheads/upfxc2abi
1131 : : !! NAME
1132 : : !! upfxc2abi
1133 : : !!
1134 : : !! FUNCTION
1135 : : !! This routine wraps a call to an OCTOPUS module, which reformats
1136 : : !! a UPF (PWSCF / Espresso) string describing XC functionals,
1137 : : !! and returns the abinit internal code pspxc
1138 : : !!
1139 : : !! INPUTS
1140 : : !! dft = string with x/c functionals from PWSCF format
1141 : : !!
1142 : : !! OUTPUT
1143 : : !! pspxc = index of xc functional for this pseudo
1144 : : !!
1145 : : !! NOTES
1146 : : !! FIXME: extend to more functionals with libxc
1147 : : !! Could be included in separate module, eg read_upf_pwscf or funct_pwscf
1148 : : !! Left without defs_basis or calls to abinit routines ON PURPOSE
1149 : : !!
1150 : : !! SOURCE
1151 : :
1152 : 2 : subroutine upfxc2abi(dft, pspxc)
1153 : :
1154 : : use funct_pwscf ! pwscf module for naming xc functionals
1155 : :
1156 : : !Arguments -------------------------------
1157 : : character(len=*), intent(in) :: dft
1158 : : integer, intent(out) :: pspxc
1159 : :
1160 : : !Local variables -------------------------
1161 : : integer :: iexch,icorr,igcx,igcc
1162 : : integer :: totalindex, offset
1163 : : ! *********************************************************************
1164 : :
1165 : : !extract from char*20 :: dft(:)
1166 : : !### The following has been copied from pwscf src/Modules/upf_to_internal.f90:
1167 : : !workaround for rrkj format - it contains the indices, not the name
1168 [ - + ]: 2 : if ( dft(1:6)=='INDEX:') then
1169 : 0 : read( dft(7:10), '(4i1)') iexch,icorr,igcx,igcc
1170 : 0 : call set_dft_from_indices(iexch,icorr,igcx,igcc)
1171 : : else
1172 : 2 : call set_dft_from_name( dft )
1173 : 2 : iexch = get_iexch()
1174 : 2 : icorr = get_icorr()
1175 : 2 : igcx = get_igcx()
1176 : 2 : igcc = get_igcc()
1177 : : end if
1178 : :
1179 : : !reset dft string to avoid stray spaces
1180 : 2 : call set_dft_from_indices(iexch,icorr,igcx,igcc)
1181 : 2 : write(std_out,'(a)') ' upf2abinit: XC string from pseudopotential is :'
1182 : 2 : write(std_out,'(3a)') '>', dft, '<'
1183 : 2 : ABI_WARNING("upfxc2abi is not guaranteed to return the right ixc from QE XC string e.g. PBEsol. Please crosscheck!")
1184 : :
1185 : 2 : offset = 100
1186 : 2 : totalindex = offset*offset*offset*iexch + offset*offset*icorr + offset*igcx + igcc
1187 : 0 : select case (totalindex)
1188 : : case (00000000) !(" NOX NOC NOGX NOGC") ! no xc
1189 : 0 : pspxc = 0
1190 : : case (01010000) !(" SLA PZ NOGX NOGC") ! slater exchange + Perdew Zunger
1191 : 0 : pspxc = 2
1192 : : case (01050000) !(" SLA WIG NOGX NOGC") ! slater exchange + Wigner corr
1193 : 0 : pspxc = 4
1194 : : case (01060000) !(" SLA HL NOGX NOGC") ! Hedin + Lundqvist
1195 : 0 : pspxc = 5
1196 : : case (02000000) !(" SL1 NOC NOGX NOGC") ! full slater exchange
1197 : 0 : pspxc = 6
1198 : : case (01040000) !(" SLA PW NOGX NOGC") ! slater exchange + Perdew Wang
1199 : 2 : pspxc = 7
1200 : : case (01000000) !(" SLA NOC NOGX NOGC") ! Perdew Wang + no corr
1201 : 0 : pspxc = 8
1202 : : case (01040304) !(" SLA PW PBX PBC") ! LDA + PBE GGA
1203 : 0 : pspxc = 11 ! PBE
1204 : : case (01000300) !(" SLA NOC PBX NOGC") ! exchange part of PBE GGA
1205 : 0 : pspxc = 12
1206 : : case (01040404) !(" SLA PW RPB PBC") ! rev PBE
1207 : 0 : pspxc = 14
1208 : : case (00000505) !(" NOX NOC HTCH HTCH") ! HTCH 120
1209 : 0 : pspxc = 17
1210 : : case (01030103) !(" SLA LYP B88 BLYP") ! BLYP
1211 : 0 : pspxc = -106131
1212 : : case (01040101) !(" SLA PW B88 P86") ! BP86
1213 : 0 : pspxc = -106132
1214 : : case (00030603) !(" NOX LYP OPTX BLYP") ! OLYP
1215 : 0 : pspxc = -110131
1216 : : ! case () !("R2SCAN01") ! R2SCAN01
1217 : : ! pspxc = -645642
1218 : : ! FIXME: important cases left to be patched with libxc:
1219 : : ! vosko wilkins nusair
1220 : : ! ortiz ballone
1221 : : ! pbe0
1222 : : ! Gunnarson-Lunqvist
1223 : : ! make general approach: check gradient parts first, then lda.
1224 : : ! event. check if they are consistent.
1225 : : case default
1226 [ - - - - : 2 : ABI_ERROR('upf2abinit: XC functional not recognized')
- + - - -
- - - - -
- ]
1227 : : end select
1228 : :
1229 : 2 : end subroutine upfxc2abi
1230 : : !!***
1231 : :
1232 : : !!****f* m_pspheads/updft_to_ixc
1233 : : !! NAME
1234 : : !! updft_to_ixc
1235 : : !!
1236 : : !! FUNCTION
1237 : : !! Returns the abinit internal `ixc` from `dft` string with XC functional in QE format.
1238 : : !!
1239 : : !! SOURCE
1240 : :
1241 : 16 : integer function upfdft_to_ixc(dft, ixc, msg) result(ierr)
1242 : :
1243 : : !Arguments ------------------------------------
1244 : : character(len=*),intent(in) :: dft
1245 : : character(len=*),intent(out) :: msg
1246 : : integer,intent(out) :: ixc
1247 : :
1248 : : !Local variables-------------------------------
1249 : : integer :: start !, ii
1250 : : character(len=500) :: x_name, c_name, gcx_name, gcc_name
1251 : : !*************************************************************************
1252 : :
1253 : : ! This list taken from oncvpsp/src/upfout.f90
1254 : : ! It should be OK as long as the UPF2 NC pseudos are generated with oncvpsp
1255 : : ! but it does not cover all QE possibilities.
1256 : 16 : ierr = 0; msg = ""
1257 : 16 : ixc = 0
1258 : 0 : select case (dft)
1259 : : case ("PZ")
1260 : 0 : ixc = -001009
1261 : : case ("PBE")
1262 : 12 : ixc = -101130 !; ixc = 11
1263 : : case ("PW91")
1264 : 0 : ixc = -109134
1265 : : case ("PBESOL")
1266 : 0 : ixc = -116133
1267 : : case ("REVPBE")
1268 : 0 : ixc = -102130
1269 : : case ("BP")
1270 : 0 : ixc = -106132
1271 : : case ("BLYP")
1272 : 0 : ixc = -106131
1273 : : case ("WC")
1274 : 0 : ixc = -118130
1275 : : case ('SLA PW NOGX NOGC', "SLA+PW") ! string produced by oncvpsp3
1276 : 0 : ixc = -1012
1277 : : case ('R2SCAN01')
1278 : 4 : ixc = -645642
1279 : : case default
1280 [ - - + - : 16 : ierr = 1
- - - - -
- + ]
1281 : : end select
1282 : :
1283 : : ! Extract substrings with
1284 : : ! 1) exchange
1285 : : ! 2) correlation
1286 : : ! 3) gradient correction, exchange
1287 : : ! 4) gradient correction, correlation
1288 : : if (ierr == 1) then
1289 : 0 : ierr = 0; start = 1
1290 [ # # ]: 0 : ABI_CHECK(next_token(dft, start, x_name) == 0 , "Error reading x_name")
1291 [ # # ]: 0 : ABI_CHECK(next_token(dft, start, c_name) == 0 , "Error reading c_name")
1292 [ # # ]: 0 : ABI_CHECK(next_token(dft, start, gcx_name) == 0 , "Error reading gcx_name")
1293 [ # # ]: 0 : ABI_CHECK(next_token(dft, start, gcc_name) == 0 , "Error reading gcc_name")
1294 : : !call remove_non_ascii(gcc_name)
1295 : : !print *, "dft: `", trim(dft), "`"
1296 : : !print *, "x_name: `", trim(x_name), "`, c_name: `", trim(c_name), &
1297 : : ! "`, gcx_name: `", trim(gcx_name), "`, gcc_name: `", trim(gcc_name), "`"
1298 : :
1299 [ # # # # ]: 0 : if (x_name == "SLA" .and. c_name == "PW") then
1300 : : !print *, "in first if", gcx_name == "NOGX", trim(gcc_name) == "NOGC"
1301 : : !print *, "len_trim(gcc_name)", len_trim(gcc_name)
1302 [ # # # # ]: 0 : if (gcx_name == "NOGX" .and. gcc_name == "NOGC") then
1303 : 0 : ixc = -1012
1304 : : else
1305 : : ierr = 1
1306 : : !print *, "in second ierr"
1307 : : end if
1308 : : else
1309 : : ierr = 1
1310 : : end if
1311 : : end if
1312 : :
1313 : : if (ierr == 1) then
1314 : : write(msg, "(5a)") &
1315 : 0 : "Cannot find ABINIT ixc value corresponding to QE dft string: `", trim(dft), "`", ch10, &
1316 : 0 : "Please update mapping in m_pspheads/upfdft_to_ixc."
1317 : : end if
1318 : :
1319 [ + - ]: 16 : end function upfdft_to_ixc
1320 : : !!***
1321 : :
1322 : : ! Copyright (c) 1989-2019 by D. R. Hamann, Mat-Sim Research LLC and Rutgers
1323 : : ! University
1324 : : !
1325 : : !
1326 : : ! This program is free software: you can redistribute it and/or modify
1327 : : ! it under the terms of the GNU General Public License as published by
1328 : : ! the Free Software Foundation, either version 3 of the License, or
1329 : : ! (at your option) any later version.
1330 : : !
1331 : : ! This program is distributed in the hope that it will be useful,
1332 : : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
1333 : : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1334 : : ! GNU General Public License for more details.
1335 : : !
1336 : : ! You should have received a copy of the GNU General Public License
1337 : : ! along with this program. If not, see <http://www.gnu.org/licenses/>.
1338 : : !
1339 : 6 : subroutine sr_so_r(lmax,irc,nproj,rr,mmax,mxprj,evkb,vkb, &
1340 : 6 : & vsr,esr,vso,eso)
1341 : :
1342 : : ! reformulates non-local potentials based on j = l +/- 1/2 to scalar-
1343 : : ! relativistic and L dot S projectors
1344 : : ! uses relationship <L dot S> = (J^2 - L^2 - S^2)/2
1345 : : ! so L dot S = +/- l/2 for j = l +/- 1/2
1346 : :
1347 : : !lmax maximum angular momentum
1348 : : !irc core radii indices
1349 : : !nproj number of projectors for each l
1350 : : !rr log radial grid
1351 : : !mmax size of radial grid
1352 : : !mmax dimension of log grid
1353 : : !mxprj dimension of number of projectors
1354 : : !vkb vkb projectors
1355 : : !evkb coefficients of BKB projectors
1356 : : !vsr normalized scalar projectors
1357 : : !esr energy coefficients of vscal
1358 : : !vso normalized spin-orbig projectors
1359 : : !esol energy coefficients of vso
1360 : :
1361 : : !implicit none
1362 : : !integer, parameter :: dp=kind(1.0d0)
1363 : :
1364 : : !Input variables
1365 : : integer,intent(in) :: lmax,mmax,mxprj
1366 : : integer,intent(in) :: irc(6),nproj(6)
1367 : : real(dp),intent(in) :: rr(mmax),vkb(mmax,mxprj,4,2),evkb(mxprj,4,2)
1368 : :
1369 : : !Output variables
1370 : : real(dp),intent(out) :: vsr(mmax,2*mxprj,4),esr(2*mxprj,4)
1371 : : real(dp),intent(out) :: vso(mmax,2*mxprj,4),eso(2*mxprj,4)
1372 : :
1373 : : !Local variables
1374 : : integer :: ii,jj,kk,ik1,ik2,ip1,ip2,ipk,ll,l1,info,nn
1375 : : real(dp) :: amesh
1376 : : real(dp) :: apk,tt
1377 : 12 : real(dp) :: sovl(2*mxprj,2*mxprj),sovlev(2*mxprj),ascl(2*mxprj,2*mxprj),aso(2*mxprj,2*mxprj)
1378 : 12 : real(dp) :: asclst(2*mxprj,2*mxprj),wsclst(2*mxprj),asost(2*mxprj,2*mxprj),wsost(2*mxprj)
1379 : 12 : real(dp) :: asclt(2*mxprj,2*mxprj),asot(2*mxprj,2*mxprj)
1380 : 12 : real(dp) :: sphalf(2*mxprj,2*mxprj),smhalf(2*mxprj,2*mxprj)
1381 : 12 : real(dp) :: fscl(mxprj),fso(mxprj),work(10*mxprj)
1382 : 6 : real(dp), allocatable :: vkbt(:,:),vkbst(:,:)
1383 : : logical :: sorted
1384 : : character(len=500) :: msg
1385 : :
1386 : : ABI_UNUSED(irc(1))
1387 : :
1388 : : ! Check that rad grid is linear starting at zero
1389 : : !linear_mesh = .True.
1390 : 6 : amesh = rr(2) - rr(1) !; damesh = zero
1391 : : !do irad=2,mmax-1
1392 : : ! damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1)))
1393 : : !end do
1394 : : !linear_mesh = damesh < tol8
1395 : :
1396 : : !if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then
1397 : : ! write(msg,'(3a)')&
1398 : : ! 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
1399 : : ! 'Action: check your pseudopotential file.'
1400 : : ! ABI_ERROR(msg)
1401 : : !end if
1402 : :
1403 : : !allocate(vkbt(mmax,2*mxprj),vkbst(mmax,2*mxprj))
1404 [ - + + - : 30 : ABI_MALLOC(vkbt, (mmax,2*mxprj))
+ - - + +
- + - + -
- + - + ]
1405 [ - + + - : 18 : ABI_MALLOC(vkbst, (mmax,2*mxprj))
- + + - -
+ - + ]
1406 : :
1407 [ + + ]: 24 : do l1=1,lmax+1
1408 : 18 : ll=l1-1
1409 : :
1410 [ + + ]: 18 : if(ll==0) then
1411 [ + + + + ]: 39438 : vsr(:,:,l1)=0.0d0
1412 [ + + + + ]: 39438 : vso(:,:,l1)=0.0d0
1413 [ + + ]: 30 : esr(:,l1)=0.0d0
1414 [ + + ]: 30 : eso(:,l1)=0.0d0
1415 [ + - ]: 6 : if(nproj(l1)>=1) then
1416 [ + + ]: 18 : do ii=1,nproj(l1)
1417 [ + + ]: 19716 : vsr(:,ii,l1)=vkb(:,ii,l1,1)
1418 : 18 : esr(ii,l1)=evkb(ii,l1,1)
1419 : : end do
1420 : : end if
1421 : : cycle
1422 : : end if
1423 : :
1424 : 12 : nn=2*nproj(l1)
1425 : :
1426 : 12 : fscl(1)=(ll+1)/dble(2*ll+1)
1427 : 12 : fscl(2)=ll/dble(2*ll+1)
1428 : 12 : fso(1)=2/dble(2*ll+1)
1429 : 12 : fso(2)=-2/dble(2*ll+1)
1430 : :
1431 : : ! construct overlap matrix and diagonal energy matrices
1432 : :
1433 [ + + + + ]: 252 : sovl(:,:)=0.0d0
1434 [ + + + + ]: 252 : ascl(:,:)=0.0d0
1435 [ + + + + ]: 252 : aso(:,:)=0.0d0
1436 [ + + + + ]: 78876 : vkbt(:,:)=0.0d0
1437 : :
1438 [ + + ]: 36 : do ik1=1,2
1439 [ + + ]: 84 : do ip1=1,nproj(l1)
1440 : 48 : ii=ip1+(ik1-1)*nproj(l1)
1441 : :
1442 : 48 : ascl(ii,ii)=fscl(ik1)*evkb(ip1,l1,ik1)
1443 : 48 : aso(ii,ii)=fso(ik1)*evkb(ip1,l1,ik1)
1444 : :
1445 [ + + ]: 78864 : vkbt(:,ii)=vkb(:,ip1,l1,ik1)
1446 : :
1447 [ + + ]: 168 : do ik2=1,2
1448 [ + + ]: 336 : do ip2=1,nproj(l1)
1449 : 192 : jj=ip2+(ik2-1)*nproj(l1)
1450 : :
1451 : : ! MG: This routine cannot be used as it assumes log mesh.
1452 : : !call vpinteg(vkb(1,ip1,l1,ik1),vkb(1,ip2,l1,ik2),irc(l1),2*l1, &
1453 : : ! sovl(ii,jj),rr)
1454 : :
1455 : : ! So we replace it with simpson integration
1456 [ + - + + ]: 315552 : sovl(ii,jj) = simpson(amesh, vkb(:,ip1,l1,ik1) * vkb(:,ip2,l1,ik2))
1457 : : end do
1458 : : end do
1459 : : end do
1460 : : end do
1461 : :
1462 : 12 : call dsyev( 'V', 'U', nn, sovl, 2*mxprj, sovlev, work, 10*mxprj, info )
1463 : :
1464 [ - + ]: 12 : if(info .ne. 0) then
1465 : 0 : write(msg,'(a,i4)') 'sr_so_r: S matrix eigenvalue ERROR, info=',info
1466 : 0 : ABI_ERROR(msg)
1467 : : end if
1468 : :
1469 : : ! construct S^(-1/2) AND s^(1/2)
1470 : :
1471 [ + + ]: 60 : do jj=1,nn
1472 : 48 : tt=sqrt(sovlev(jj))
1473 [ + + ]: 252 : do ii=1,nn
1474 : 192 : sphalf(ii,jj)=tt*sovl(ii,jj)
1475 : 240 : smhalf(ii,jj)=sovl(ii,jj)/tt
1476 : : end do
1477 : : end do
1478 : :
1479 : : ! take linear combinations to form orthonormal basis functions
1480 : :
1481 [ + + + + ]: 78876 : vkbst(:,:)=0.0d0
1482 : :
1483 [ + + ]: 60 : do jj=1,nn
1484 [ + + ]: 252 : do ii=1,nn
1485 [ + + ]: 315504 : vkbst(:,jj)=vkbst(:,jj) + smhalf(ii,jj)*vkbt(:,ii)
1486 : : end do
1487 : : end do
1488 : :
1489 : : ! construct A^(-1)* = S^(1/2)^T A^(-1) S^(1/2)
1490 : :
1491 [ + + + + ]: 252 : asclt(:,:)=0.0d0
1492 [ + + + + ]: 252 : asclst(:,:)=0.0d0
1493 [ + + + + ]: 252 : asot(:,:)=0.0d0
1494 [ + + + + ]: 252 : asost(:,:)=0.0d0
1495 : :
1496 [ + + ]: 60 : do ii=1,nn
1497 [ + + ]: 252 : do jj=1,nn
1498 [ + + ]: 1008 : do kk=1,nn
1499 : 768 : asclt(ii,jj)=asclt(ii,jj)+ascl(ii,kk)*sphalf(kk,jj)
1500 : 960 : asot(ii,jj) =asot(ii,jj) + aso(ii,kk)*sphalf(kk,jj)
1501 : : end do
1502 : : end do
1503 : : end do
1504 : :
1505 [ + + ]: 60 : do ii=1,nn
1506 [ + + ]: 252 : do jj=1,nn
1507 [ + + ]: 1008 : do kk=1,nn
1508 : 768 : asclst(ii,jj)=asclst(ii,jj)+asclt(kk,jj)*sphalf(kk,ii)
1509 : 960 : asost(ii,jj) =asost(ii,jj) + asot(kk,jj)*sphalf(kk,ii)
1510 : : end do
1511 : : end do
1512 : : end do
1513 : :
1514 : : ! find eigenvalues and eigenvectors of the A* matrices
1515 : :
1516 : : ! SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
1517 : :
1518 : 12 : call dsyev( 'V', 'U', nn, asclst, 2*mxprj, wsclst, work, 10*mxprj, info )
1519 : :
1520 [ - + ]: 12 : if(info .ne. 0) then
1521 : 0 : write(msg,'(a,i4)') 'sr_so_r: A* matrix eigenvalue ERROR, info=',info
1522 : 0 : ABI_ERROR(msg)
1523 : : end if
1524 : :
1525 : 12 : call dsyev( 'V', 'U', nn, asost, 2*mxprj, wsost, work, 10*mxprj, info )
1526 : :
1527 : :
1528 [ - + ]: 12 : if(info .ne. 0) then
1529 : 0 : write(msg,'(a,i4)') 'sr_so_r: A* matrix eigenvalue ERROR, info=',info
1530 : 0 : ABI_ERROR(msg)
1531 : : end if
1532 : :
1533 : : ! take linear combinations to form orthonormal projectors
1534 : :
1535 [ + + + + ]: 78876 : vsr(:,:,l1)=0.0d0
1536 [ + + + + ]: 78876 : vso(:,:,l1)=0.0d0
1537 [ + + ]: 60 : esr(:,l1)=0.0d0
1538 [ + + ]: 60 : eso(:,l1)=0.0d0
1539 : :
1540 [ + + ]: 60 : do ii=1,nn
1541 : 48 : esr(ii,l1)=wsclst(ii)
1542 : 48 : eso(ii,l1)= wsost(ii)
1543 [ + + ]: 252 : do jj=1,nn
1544 [ + + ]: 315456 : vsr(:,ii,l1)=vsr(:,ii,l1) + asclst(jj,ii)*vkbst(:,jj)
1545 [ + + ]: 315504 : vso(:,ii,l1)= vso(:,ii,l1) + asost(jj,ii)*vkbst(:,jj)
1546 : : end do
1547 : : end do
1548 : :
1549 : : ! bubble-sort on coefficient magnitudes for scalar and then s-o
1550 : : ! (Yes, I know bubble-sort is the least-efficient sorting algorithm.)
1551 : :
1552 [ - + ]: 30 : do ii=1,100
1553 : : sorted=.true.
1554 [ + + ]: 120 : do jj=2,nn
1555 [ + + ]: 120 : if(abs(esr(jj-1,l1))<abs(esr(jj,l1))) then
1556 : 36 : tt=esr(jj,l1)
1557 [ + + ]: 59148 : vkbt(:,1)=vsr(:,jj,l1)
1558 : 36 : esr(jj,l1)=esr(jj-1,l1)
1559 [ + + ]: 59148 : vsr(:,jj,l1)=vsr(:,jj-1,l1)
1560 : 36 : esr(jj-1,l1)=tt
1561 [ + + ]: 59148 : vsr(:,jj-1,l1)=vkbt(:,1)
1562 : : sorted=.false.
1563 : : end if
1564 : : end do
1565 [ + + ]: 30 : if(sorted) exit
1566 : : end do
1567 : :
1568 [ + - ]: 48 : do ii=1,100
1569 : : sorted=.true.
1570 [ + + ]: 192 : do jj=2,nn
1571 [ + + ]: 192 : if(abs(eso(jj-1,l1))<abs(eso(jj,l1))) then
1572 : 54 : tt=eso(jj,l1)
1573 [ + + ]: 88722 : vkbt(:,1)=vso(:,jj,l1)
1574 : 54 : eso(jj,l1)=eso(jj-1,l1)
1575 [ + + ]: 88722 : vso(:,jj,l1)=vso(:,jj-1,l1)
1576 : 54 : eso(jj-1,l1)=tt
1577 [ + + ]: 88722 : vso(:,jj-1,l1)=vkbt(:,1)
1578 : : sorted=.false.
1579 : : end if
1580 : : end do
1581 [ + + ]: 48 : if(sorted) exit
1582 : : end do
1583 : :
1584 : : write(std_out,'(/a,i2)') &
1585 : 12 : & ' Orthonormal scalar projector coefficients, l = ',ll
1586 : 12 : write(std_out,'(1p,6e12.4)') (esr(jj,l1),jj=1,nn)
1587 : : write(std_out,'(/a,i2)') &
1588 : 12 : & ' Orthonormal spin-orbit projector coefficients, l = ',ll
1589 : 12 : write(std_out,'(1p,6e12.4)') (eso(jj,l1),jj=1,nn)
1590 : :
1591 : : ! Set sign of projectors (physically irrelevant) so that they are positive
1592 : : ! at their peak (needed for compaisons apparently)
1593 : :
1594 [ + + ]: 66 : do jj=1,nn
1595 : : apk=0.0d0
1596 [ + + ]: 78864 : do ii=1,mmax
1597 [ + + ]: 78864 : if(abs(vso(ii,jj,l1))>apk) then
1598 : 3102 : apk=abs(vso(ii,jj,l1))
1599 : 3102 : ipk=ii
1600 : : end if
1601 : : end do
1602 [ + + ]: 48 : if(vso(ipk,jj,l1)<0.0d0) then
1603 [ + + ]: 31278 : vso(:,jj,l1)=-vso(:,jj,l1)
1604 : : end if
1605 : 48 : apk=0.0d0
1606 [ + + ]: 78864 : do ii=1,mmax
1607 [ + + ]: 78864 : if(abs(vsr(ii,jj,l1))>apk) then
1608 : 2910 : apk=abs(vsr(ii,jj,l1))
1609 : 2910 : ipk=ii
1610 : : end if
1611 : : end do
1612 [ + + ]: 60 : if(vsr(ipk,jj,l1)<0.0d0) then
1613 [ + + ]: 41136 : vsr(:,jj,l1)=-vsr(:,jj,l1)
1614 : : end if
1615 : : end do
1616 : :
1617 : : end do ! l1
1618 : :
1619 [ - + ]: 6 : ABI_FREE(vkbt)
1620 [ - + ]: 6 : ABI_FREE(vkbst)
1621 : 6 : return
1622 : : end subroutine sr_so_r
1623 : : !!***
1624 : :
1625 : : ! Copyright (c) 1989-2019 by D. R. Hamann, Mat-Sim Research LLC and Rutgers
1626 : : ! University
1627 : : !
1628 : : !
1629 : : ! This program is free software: you can redistribute it and/or modify
1630 : : ! it under the terms of the GNU General Public License as published by
1631 : : ! the Free Software Foundation, either version 3 of the License, or
1632 : : ! (at your option) any later version.
1633 : : !
1634 : : ! This program is distributed in the hope that it will be useful,
1635 : : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
1636 : : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1637 : : ! GNU General Public License for more details.
1638 : : !
1639 : : ! You should have received a copy of the GNU General Public License
1640 : : ! along with this program. If not, see <http://www.gnu.org/licenses/>.
1641 : : !
1642 : : ! subroutine vpinteg(gg,hh,nn,mm,ss,rr)
1643 : : !
1644 : : !! integrals that go into construction of Vanderbilt separable pseudopotential
1645 : : !
1646 : : !! product of functions gg*hh goes like rr**mm at rr -> 0
1647 : : !! integral on usual log mesh from 1 to nn
1648 : : !
1649 : : !!Input variables
1650 : : ! integer,intent(in) :: nn,mm
1651 : : ! real(dp),intent(in) :: gg(nn),hh(nn),rr(nn)
1652 : : !
1653 : : !!Output variable
1654 : : ! real(dp),intent(out) :: ss
1655 : : !
1656 : : !!Local variables
1657 : : ! real(dp) :: r0,amesh,al
1658 : : ! integer :: ii
1659 : : !
1660 : : ! al = 0.01d0 * dlog(rr(101)/rr(1))
1661 : : ! amesh = exp(al)
1662 : : !
1663 : : ! r0=rr(1)/dsqrt(amesh)
1664 : : ! ss=r0**(mm+1)*(gg(1)*hh(1)/rr(1)**mm)/dfloat(mm+1)
1665 : : !
1666 : : ! do ii = 4, nn - 3
1667 : : ! ss = ss + al*gg(ii)*hh(ii)*rr(ii)
1668 : : ! end do
1669 : : !
1670 : : ! ss=ss + al*(23.d0*rr(nn-2)*gg(nn-2)*hh(nn-2) &
1671 : : !& + 28.d0*rr(nn-1)*gg(nn-1)*hh(nn-1) &
1672 : : !& + 9.d0*rr(nn )*gg(nn )*hh(nn ))/24.d0
1673 : : !
1674 : : !
1675 : : ! return
1676 : : ! end subroutine vpinteg
1677 : :
1678 : : end module m_pspheads
1679 : : !!***
|
