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 1162 : 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 1162 : integer,allocatable :: nproj(:),nprojso(:)
107 1162 : 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 1162 : test_paw=0
124 :
125 2732 : do ipsp=1,npsp
126 :
127 1570 : pspheads(ipsp)%filpsp=trim(filnam(ipsp))
128 :
129 : ! Check if the file is written in XML
130 1570 : usexml = 0
131 1570 : if (open_file(filnam(ipsp), msg, newunit=unt, form="formatted", status="old") /= 0) then
132 0 : ABI_ERROR(msg)
133 : end if
134 :
135 1570 : rewind(unit=unt, err=10, iomsg=errmsg)
136 1570 : read(unt, "(a)", err=10, iomsg=errmsg) testxml
137 :
138 1570 : if(testxml(1:5)=='<?xml')then
139 198 : usexml = 1
140 198 : read(unt,*, err=10, iomsg=errmsg) testxml
141 198 : 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 1570 : useupf = 0
153 1570 : 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 1570 : close(unit=unt, err=10, iomsg=errmsg)
163 :
164 : ! Check if pseudopotential file is a QE UPF1 file
165 1570 : if (useupf == 0) then
166 1564 : if (open_file(filnam(ipsp), msg, newunit=unt, form="formatted", status="old") /= 0) then
167 0 : ABI_ERROR(msg)
168 : end if
169 :
170 1564 : rewind(unit=unt, err=10, iomsg=errmsg)
171 1564 : read(unt,*,err=10,iomsg=errmsg) testxml ! just a string, no relation to xml.
172 1564 : if (testxml(1:9)=='<PP_INFO>') then
173 : useupf = 1
174 : else
175 1563 : useupf = 0
176 : end if
177 1564 : close(unit=unt,err=10,iomsg=errmsg)
178 : end if
179 1570 : if(test_paw==0) then
180 6610 : ABI_ICALLOC(nproj,(0:3))
181 5288 : ABI_ICALLOC(nprojso,(3))
182 : endif
183 :
184 : ! Read the header of the pseudopotential file
185 1570 : if (usexml /= 1 .and. useupf == 0) then
186 : ! Open the psp file and read a normal abinit style header
187 1365 : if (open_file(filnam(ipsp), msg, newunit=unt, form='formatted', status='old') /= 0) then
188 0 : ABI_ERROR(msg)
189 : end if
190 1365 : rewind (unit=unt, err=10, iomsg=errmsg)
191 :
192 : ! Read the three first lines
193 1365 : read(unt, '(a)', err=10, iomsg=errmsg) pspheads(ipsp)%title
194 1365 : read(unt,*, err=10, iomsg=errmsg)pspheads(ipsp)%znuclpsp,pspheads(ipsp)%zionpsp,pspheads(ipsp)%pspdat
195 1365 : read(unt,*, err=10, iomsg=errmsg)pspheads(ipsp)%pspcod,pspheads(ipsp)%pspxc,pspheads(ipsp)%lmax,idum,mmax
196 :
197 1365 : pspcod=pspheads(ipsp)%pspcod
198 1365 : lmax=pspheads(ipsp)%lmax
199 1365 : write(msg,'(a,f5.1,a,i4,a,i4)')' read the values zionpsp=',pspheads(ipsp)%zionpsp,' , pspcod=',pspcod,' , lmax=',lmax
200 1365 : call wrtout(std_out,msg,'PERS')
201 :
202 1365 : if(test_paw==1) then
203 183 : ABI_MALLOC(nproj,(0:lmax))
204 183 : ABI_MALLOC(nprojso,(lmax))
205 205 : nproj(0:lmax)=0
206 144 : nprojso(1:lmax)=0
207 : else
208 10432 : nproj(0:3)=0 ; nprojso(1:3)=0
209 : endif
210 :
211 1365 : pspheads(ipsp)%xccc=0
212 1365 : pspheads(ipsp)%pspso=0
213 :
214 205 : 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 194 : else if(usexml==1.and.test_paw==1)then
237 :
238 : write(msg,'(4a)') &
239 187 : '- inpspheads : Reading pseudopotential header in XML form from ',ch10,&
240 374 : '- ',trim(filnam(ipsp))
241 561 : call wrtout([std_out, ab_out], msg)
242 :
243 187 : call pawpsxml2ab(filnam(ipsp),ecut_tmp(:,:,ipsp), pspheads(ipsp),1)
244 187 : 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 9420 : pspheads(ipsp)%GTHradii = zero
278 1570 : if(test_paw/=1.and.lmax>3) then
279 0 : ABI_ERROR('NC with l>3 not implmented')
280 : endif
281 1570 : if(.not.allocated(pspheads(ipsp)%nproj)) then
282 4116 : ABI_MALLOC(pspheads(ipsp)%nproj,(0:lmax))
283 : endif
284 1570 : if(.not.allocated(pspheads(ipsp)%nprojso)) then
285 4116 : ABI_MALLOC(pspheads(ipsp)%nprojso,(lmax))
286 : endif
287 1570 : if(.not.allocated(nproj)) then
288 561 : ABI_MALLOC(nproj,(0:lmax))
289 : endif
290 1570 : if(.not.allocated(nprojso)) then
291 561 : ABI_MALLOC(nprojso,(lmax))
292 : endif
293 :
294 1570 : if (pspcod==1 .or. pspcod==4)then
295 :
296 : ! Teter format
297 1583 : do ilmax=0,lmax
298 1119 : read(unt,*, err=10, iomsg=errmsg) lang,e990,e999,nproj(ilmax)
299 1583 : read(unt,*, err=10, iomsg=errmsg)
300 : end do
301 464 : read(unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
302 464 : 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 195 : test_paw=1;pspheads(ipsp)%pawheader%pawver=1
381 195 : read (unt,'(a80)', err=10, iomsg=errmsg) pspline
382 195 : pspline=adjustl(pspline)
383 195 : if (pspline(1:3)=="paw".or.pspline(1:3)=="PAW") &
384 195 : & read(unit=pspline(4:80),fmt=*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%pawver
385 195 : 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 195 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%basis_size,pspheads(ipsp)%pawheader%lmn_size
399 585 : ABI_MALLOC(orb,(pspheads(ipsp)%pawheader%basis_size))
400 1149 : orb(:)=0
401 195 : read (unt,*, err=10, iomsg=errmsg) (orb(ii), ii=1,pspheads(ipsp)%pawheader%basis_size)
402 195 : pspheads(ipsp)%pawheader%mesh_size=mmax
403 195 : read (unt,*, err=10, iomsg=errmsg) nmesh
404 928 : do ii=1,nmesh
405 928 : read(unt,*, err=10, iomsg=errmsg)
406 : end do
407 195 : read (unt,*, err=10, iomsg=errmsg) pspheads(ipsp)%pawheader%rpaw
408 195 : pspheads(ipsp)%pawheader%rshp=pspheads(ipsp)%pawheader%rpaw
409 195 : read (unt,'(a80)', err=10, iomsg=errmsg) pspline
410 195 : pspline=adjustl(pspline); write(std_out,*) pspline
411 195 : read(unit=pspline,fmt=*) pspheads(ipsp)%pawheader%shape_type
412 195 : if (pspheads(ipsp)%pawheader%pawver==2.and.&
413 0 : pspheads(ipsp)%pawheader%shape_type==3) pspheads(ipsp)%pawheader%shape_type=-1
414 195 : 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 682 : do ilmax=0,lmax
421 3233 : do ii=1,pspheads(ipsp)%pawheader%basis_size
422 3038 : if(orb(ii)==ilmax) nproj(ilmax)=nproj(ilmax)+1
423 : end do
424 : end do
425 1149 : pspheads(ipsp)%pawheader%l_size=2*maxval(orb)+1
426 195 : pspheads(ipsp)%xccc=1 ! We suppose apriori that cc is used (but n1xccc is not used in PAW)
427 195 : 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 284 : read(unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
461 284 : if (fchrg>1.d-15) pspheads(ipsp)%xccc=n1xccc_default
462 284 : read(unt,*, err=10, iomsg=errmsg) nproj(0:lmax)
463 284 : read(unt,*, err=10, iomsg=errmsg) extension_switch
464 284 : if(any(extension_switch == [2, 3])) then
465 64 : pspso=2
466 64 : read(unt,*,err=10,iomsg=errmsg) nprojso(1:lmax)
467 : else
468 220 : pspso=0
469 : end if
470 284 : 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 194 : 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 1570 : if (pspcod /= 17) then
512 4766 : pspheads(ipsp)%nproj(0:lmax)=nproj(0:lmax)
513 3385 : 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 1570 : close(unt)
519 :
520 : ! Compute md5 checksum
521 1570 : pspheads(ipsp)%md5_checksum = md5_sum_from_file(filnam(ipsp))
522 1570 : ABI_FREE(nproj)
523 2732 : 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 1162 : mpsang=1
528 1162 : n1xccc=pspheads(1)%xccc
529 2732 : do ii=1,npsp
530 1570 : mpsang=max(pspheads(ii)%lmax+1,mpsang)
531 2732 : n1xccc=max(pspheads(ii)%xccc,n1xccc)
532 : end do
533 :
534 1162 : write(msg,'(2a,i0,a,i0,a)')ch10,' inpspheads: deduce mpsang = ',mpsang,', n1xccc = ',n1xccc,'.'
535 1162 : call wrtout(std_out,msg,'PERS')
536 :
537 : ! Test: if one psp is PAW, all must be
538 1162 : if (test_paw==1) then
539 639 : do ipsp=1,npsp
540 639 : 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 1440 : 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 1440 : integer,allocatable :: list_int(:)
591 : real(dp) :: tsec(2)
592 1440 : real(dp),allocatable :: list_dpr(:)
593 1440 : character(len=fnlen),allocatable :: list_char(:)
594 : #endif
595 : !*************************************************************************
596 :
597 : #if defined HAVE_MPI
598 1440 : call timab(48,1,tsec)
599 :
600 1440 : comm = xmpi_world
601 :
602 : ! Broadcast the characters (file names and titles)
603 4320 : ABI_MALLOC(list_char,(3*npsp))
604 3390 : list_char(1:npsp)=pspheads(1:npsp)%filpsp
605 3390 : list_char(npsp+1:2*npsp)=pspheads(1:npsp)%title
606 3390 : list_char(2*npsp+1:3*npsp)=pspheads(1:npsp)%md5_checksum
607 :
608 1440 : call xmpi_bcast(list_char,master,comm,ierr)
609 :
610 3390 : pspheads(1:npsp)%filpsp=list_char(1:npsp)
611 3390 : pspheads(1:npsp)%title=list_char(npsp+1:2*npsp)
612 3390 : pspheads(1:npsp)%md5_checksum=list_char(2*npsp+1:3*npsp)(1:md5_slen)
613 1440 : ABI_FREE(list_char)
614 :
615 : ! Brodcast the integers
616 1440 : list_size=6*npsp+1
617 4320 : ABI_MALLOC(list_int,(list_size))
618 3390 : list_int(1 : npsp) = pspheads(1:npsp)%lmax
619 3390 : list_int(1+ npsp: 2*npsp) = pspheads(1:npsp)%xccc
620 3390 : list_int(1+ 2*npsp: 3*npsp) = pspheads(1:npsp)%pspxc
621 3390 : list_int(1+ 3*npsp: 4*npsp) = pspheads(1:npsp)%pspdat
622 3390 : list_int(1+ 4*npsp: 5*npsp) = pspheads(1:npsp)%pspcod
623 3390 : list_int(1+ 5*npsp: 6*npsp) = pspheads(1:npsp)%pspso
624 1440 : list_int(6*npsp+1) = test_paw
625 :
626 1440 : call xmpi_bcast(list_int,master,comm,ierr)
627 :
628 3390 : pspheads(1:npsp)%lmax = list_int(1 : npsp)
629 3390 : pspheads(1:npsp)%xccc = list_int(1+ npsp: 2*npsp)
630 3390 : pspheads(1:npsp)%pspxc = list_int(1+ 2*npsp: 3*npsp)
631 3390 : pspheads(1:npsp)%pspdat = list_int(1+ 3*npsp: 4*npsp)
632 3390 : pspheads(1:npsp)%pspcod = list_int(1+ 4*npsp: 5*npsp)
633 3390 : pspheads(1:npsp)%pspso = list_int(1+ 5*npsp: 6*npsp)
634 1440 : test_paw = list_int(6*npsp+1)
635 :
636 1440 : ABI_FREE(list_int)
637 :
638 1440 : list_size=0
639 3390 : do ipsp=1,npsp
640 1950 : if(.not.allocated(pspheads(ipsp)%nproj)) then
641 1140 : ABI_MALLOC(pspheads(ipsp)%nproj,(0:pspheads(ipsp)%lmax))
642 : endif
643 1950 : if(.not.allocated(pspheads(ipsp)%nprojso)) then
644 1140 : ABI_MALLOC(pspheads(ipsp)%nprojso,(pspheads(ipsp)%lmax))
645 : endif
646 1950 : lmax_=max(pspheads(ipsp)%lmax,0)
647 3390 : list_size=list_size+2*lmax_+1
648 : enddo
649 4320 : ABI_MALLOC(list_int,(list_size))
650 8992 : list_int=0
651 : ii=0
652 3390 : do ipsp=1,npsp
653 8139 : do il=0,pspheads(ipsp)%lmax
654 4749 : ii=ii+1
655 6699 : list_int(ii) = pspheads(ipsp)%nproj(il)
656 : enddo
657 : enddo
658 3390 : do ipsp=1,npsp
659 6191 : do il=1,pspheads(ipsp)%lmax
660 2801 : ii=ii+1
661 4751 : list_int(ii) = pspheads(ipsp)%nprojso(il)
662 : enddo
663 : enddo
664 :
665 1440 : call xmpi_bcast(list_int,master,comm,ierr)
666 :
667 1440 : ii=0
668 3390 : do ipsp=1,npsp
669 8139 : do il=0,pspheads(ipsp)%lmax
670 4749 : ii=ii+1
671 6699 : pspheads(ipsp)%nproj(il) = list_int(ii)
672 : enddo
673 : enddo
674 3390 : do ipsp=1,npsp
675 6191 : do il=1,pspheads(ipsp)%lmax
676 2801 : ii=ii+1
677 4751 : pspheads(ipsp)%nprojso(il) = list_int(ii)
678 : enddo
679 : enddo
680 1440 : ABI_FREE(list_int)
681 :
682 : ! Unbeliveable, this cannot be sent with the others, for woopy
683 4320 : ABI_MALLOC(list_int,(npsp))
684 3390 : list_int(1:npsp) = pspheads(1:npsp)%usewvl
685 1440 : call xmpi_bcast(list_int,master,comm,ierr)
686 3390 : pspheads(1:npsp)%usewvl = list_int(1:npsp)
687 1440 : ABI_FREE(list_int)
688 :
689 : ! Broadcast zionpsp and znuclpsp
690 4320 : ABI_MALLOC(list_dpr,(7*npsp))
691 3390 : list_dpr(1 : npsp) = pspheads(1:npsp)%zionpsp
692 3390 : list_dpr(1+ npsp:2*npsp) = pspheads(1:npsp)%znuclpsp
693 3390 : list_dpr(1+2*npsp:3*npsp) = pspheads(1:npsp)%GTHradii(0)
694 3390 : list_dpr(1+3*npsp:4*npsp) = pspheads(1:npsp)%GTHradii(1)
695 3390 : list_dpr(1+4*npsp:5*npsp) = pspheads(1:npsp)%GTHradii(2)
696 3390 : list_dpr(1+5*npsp:6*npsp) = pspheads(1:npsp)%GTHradii(3)
697 3390 : list_dpr(1+6*npsp:7*npsp) = pspheads(1:npsp)%GTHradii(4)
698 :
699 1440 : call xmpi_bcast(list_dpr,master,comm,ierr)
700 :
701 3390 : pspheads(1:npsp)%zionpsp = list_dpr(1 : npsp)
702 3390 : pspheads(1:npsp)%znuclpsp = list_dpr(1+ npsp:2*npsp)
703 3390 : pspheads(1:npsp)%GTHradii(0) = list_dpr(1+2*npsp:3*npsp)
704 3390 : pspheads(1:npsp)%GTHradii(1) = list_dpr(1+3*npsp:4*npsp)
705 3390 : pspheads(1:npsp)%GTHradii(2) = list_dpr(1+4*npsp:5*npsp)
706 3390 : pspheads(1:npsp)%GTHradii(3) = list_dpr(1+5*npsp:6*npsp)
707 3390 : pspheads(1:npsp)%GTHradii(4) = list_dpr(1+6*npsp:7*npsp)
708 1440 : ABI_FREE(list_dpr)
709 :
710 : ! Broadcast additional integers for PAW psps (testpaw was sent, previously)
711 1440 : if (test_paw==1) then
712 1544 : ABI_MALLOC(list_int,(6*npsp))
713 953 : list_int(1 : npsp)=pspheads(1:npsp)%pawheader%basis_size
714 953 : list_int(1+ npsp:2*npsp)=pspheads(1:npsp)%pawheader%l_size
715 953 : list_int(1+2*npsp:3*npsp)=pspheads(1:npsp)%pawheader%lmn_size
716 953 : list_int(1+3*npsp:4*npsp)=pspheads(1:npsp)%pawheader%mesh_size
717 953 : list_int(1+4*npsp:5*npsp)=pspheads(1:npsp)%pawheader%pawver
718 953 : list_int(1+5*npsp:6*npsp)=pspheads(1:npsp)%pawheader%shape_type
719 :
720 386 : call xmpi_bcast(list_int,master,comm,ierr)
721 :
722 953 : pspheads(1:npsp)%pawheader%basis_size=list_int(1 : npsp)
723 953 : pspheads(1:npsp)%pawheader%l_size =list_int(1+ npsp:2*npsp)
724 953 : pspheads(1:npsp)%pawheader%lmn_size =list_int(1+2*npsp:3*npsp)
725 953 : pspheads(1:npsp)%pawheader%mesh_size =list_int(1+3*npsp:4*npsp)
726 953 : pspheads(1:npsp)%pawheader%pawver =list_int(1+4*npsp:5*npsp)
727 953 : pspheads(1:npsp)%pawheader%shape_type=list_int(1+5*npsp:6*npsp)
728 386 : ABI_FREE(list_int)
729 :
730 : ! broadcast rpaw values
731 1544 : ABI_MALLOC(list_dpr,(2*npsp))
732 :
733 953 : list_dpr(1 : npsp) = pspheads(1:npsp)%pawheader%rpaw
734 953 : list_dpr(1+1*npsp:2*npsp) = pspheads(1:npsp)%pawheader%rshp
735 :
736 386 : call xmpi_bcast(list_dpr,master,comm,ierr)
737 :
738 953 : pspheads(1:npsp)%pawheader%rpaw = list_dpr(1 : npsp)
739 953 : pspheads(1:npsp)%pawheader%rshp = list_dpr(1+ npsp:2*npsp)
740 :
741 386 : ABI_FREE(list_dpr)
742 : end if
743 :
744 1440 : 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 1440 : 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 492 : 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 492 : if (option==1) then
790 187 : call rdpawpsxml_header(ecut_tmp,filnam,paw_setuploc)
791 187 : paw_setuploc%idgrid= paw_setuploc%radial_grid(1)%id
792 : else
793 305 : call rdpawpsxml(filnam,paw_setuploc)
794 : end if
795 :
796 492 : call pawpsp_read_header_xml(lloc,lmax,pspcod, pspxc,paw_setuploc,r2well,zionpsp,znuclpsp)
797 :
798 492 : pspheads%lmax=lmax
799 492 : pspheads%pspxc=pspxc
800 492 : pspheads%zionpsp=zionpsp
801 492 : 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 492 : pspheads%pawheader%rshp,pspheads%pawheader%shape_type)
808 :
809 1476 : ABI_MALLOC(pspheads%nproj,(0:lmax))
810 1476 : ABI_MALLOC(pspheads%nprojso,(lmax))
811 1581 : pspheads%nproj=0
812 1581 : do il=0,pspheads%lmax
813 6694 : do ii=1,pspheads%pawheader%basis_size
814 6202 : if(paw_setuploc%valence_states%state(ii)%ll==il) pspheads%nproj(il)=pspheads%nproj(il)+1
815 : end do
816 : end do
817 :
818 1089 : pspheads%nprojso=0
819 492 : pspheads%pspdat=27061961
820 492 : pspheads%pspso=1
821 492 : pspheads%xccc=1
822 492 : pspheads%title=paw_setuploc%atom%symbol
823 :
824 492 : if (option==1) call paw_setup_free(paw_setuploc)
825 :
826 492 : 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 : !!***
|
