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 1164 : 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 1164 : integer,allocatable :: nproj(:),nprojso(:)
107 1164 : 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 1164 : test_paw=0
124 :
125 2736 : do ipsp=1,npsp
126 :
127 1572 : pspheads(ipsp)%filpsp=trim(filnam(ipsp))
128 :
129 : ! Check if the file is written in XML
130 1572 : usexml = 0
131 1572 : if (open_file(filnam(ipsp), msg, newunit=unt, form="formatted", status="old") /= 0) then
132 0 : ABI_ERROR(msg)
133 : end if
134 :
135 1572 : rewind(unit=unt, err=10, iomsg=errmsg)
136 1572 : read(unt, "(a)", err=10, iomsg=errmsg) testxml
137 :
138 1572 : if(testxml(1:5)=='<?xml')then
139 199 : usexml = 1
140 199 : read(unt,*, err=10, iomsg=errmsg) testxml
141 199 : 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 1572 : useupf = 0
153 1572 : 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 1572 : close(unit=unt, err=10, iomsg=errmsg)
163 :
164 : ! Check if pseudopotential file is a QE UPF1 file
165 1572 : if (useupf == 0) then
166 1566 : if (open_file(filnam(ipsp), msg, newunit=unt, form="formatted", status="old") /= 0) then
167 0 : ABI_ERROR(msg)
168 : end if
169 :
170 1566 : rewind(unit=unt, err=10, iomsg=errmsg)
171 1566 : read(unt,*,err=10,iomsg=errmsg) testxml ! just a string, no relation to xml.
172 1566 : if (testxml(1:9)=='<PP_INFO>') then
173 : useupf = 1
174 : else
175 1565 : useupf = 0
176 : end if
177 1566 : close(unit=unt,err=10,iomsg=errmsg)
178 : end if
179 1572 : if(test_paw==0) then
180 6615 : ABI_ICALLOC(nproj,(0:3))
181 5292 : ABI_ICALLOC(nprojso,(3))
182 : endif
183 :
184 : ! Read the header of the pseudopotential file
185 1572 : if (usexml /= 1 .and. useupf == 0) then
186 : ! Open the psp file and read a normal abinit style header
187 1366 : if (open_file(filnam(ipsp), msg, newunit=unt, form='formatted', status='old') /= 0) then
188 0 : ABI_ERROR(msg)
189 : end if
190 1366 : rewind (unit=unt, err=10, iomsg=errmsg)
191 :
192 : ! Read the three first lines
193 1366 : read(unt, '(a)', err=10, iomsg=errmsg) pspheads(ipsp)%title
194 1366 : read(unt,*, err=10, iomsg=errmsg)pspheads(ipsp)%znuclpsp,pspheads(ipsp)%zionpsp,pspheads(ipsp)%pspdat
195 1366 : read(unt,*, err=10, iomsg=errmsg)pspheads(ipsp)%pspcod,pspheads(ipsp)%pspxc,pspheads(ipsp)%lmax,idum,mmax
196 :
197 1366 : pspcod=pspheads(ipsp)%pspcod
198 1366 : lmax=pspheads(ipsp)%lmax
199 1366 : write(msg,'(a,f5.1,a,i4,a,i4)')' read the values zionpsp=',pspheads(ipsp)%zionpsp,' , pspcod=',pspcod,' , lmax=',lmax
200 1366 : call wrtout(std_out,msg,'PERS')
201 :
202 1366 : 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 10440 : nproj(0:3)=0 ; nprojso(1:3)=0
209 : endif
210 :
211 1366 : pspheads(ipsp)%xccc=0
212 1366 : pspheads(ipsp)%pspso=0
213 :
214 206 : 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 195 : else if(usexml==1.and.test_paw==1)then
237 :
238 : write(msg,'(4a)') &
239 188 : '- inpspheads : Reading pseudopotential header in XML form from ',ch10,&
240 376 : '- ',trim(filnam(ipsp))
241 564 : call wrtout([std_out, ab_out], msg)
242 :
243 188 : call pawpsxml2ab(filnam(ipsp),ecut_tmp(:,:,ipsp), pspheads(ipsp),1)
244 188 : 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 9432 : pspheads(ipsp)%GTHradii = zero
278 1572 : if(test_paw/=1.and.lmax>3) then
279 0 : ABI_ERROR('NC with l>3 not implmented')
280 : endif
281 1572 : if(.not.allocated(pspheads(ipsp)%nproj)) then
282 4119 : ABI_MALLOC(pspheads(ipsp)%nproj,(0:lmax))
283 : endif
284 1572 : if(.not.allocated(pspheads(ipsp)%nprojso)) then
285 4119 : ABI_MALLOC(pspheads(ipsp)%nprojso,(lmax))
286 : endif
287 1572 : if(.not.allocated(nproj)) then
288 564 : ABI_MALLOC(nproj,(0:lmax))
289 : endif
290 1572 : if(.not.allocated(nprojso)) then
291 564 : ABI_MALLOC(nprojso,(lmax))
292 : endif
293 :
294 1572 : 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 285 : read(unt,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
461 285 : if (fchrg>1.d-15) pspheads(ipsp)%xccc=n1xccc_default
462 285 : read(unt,*, err=10, iomsg=errmsg) nproj(0:lmax)
463 285 : read(unt,*, err=10, iomsg=errmsg) extension_switch
464 285 : 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 221 : pspso=0
469 : end if
470 285 : 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 195 : 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 1572 : if (pspcod /= 17) then
512 4770 : pspheads(ipsp)%nproj(0:lmax)=nproj(0:lmax)
513 3388 : 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 1572 : close(unt)
519 :
520 : ! Compute md5 checksum
521 1572 : pspheads(ipsp)%md5_checksum = md5_sum_from_file(filnam(ipsp))
522 1572 : ABI_FREE(nproj)
523 2736 : 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 1164 : mpsang=1
528 1164 : n1xccc=pspheads(1)%xccc
529 2736 : do ii=1,npsp
530 1572 : mpsang=max(pspheads(ii)%lmax+1,mpsang)
531 2736 : n1xccc=max(pspheads(ii)%xccc,n1xccc)
532 : end do
533 :
534 1164 : write(msg,'(2a,i0,a,i0,a)')ch10,' inpspheads: deduce mpsang = ',mpsang,', n1xccc = ',n1xccc,'.'
535 1164 : call wrtout(std_out,msg,'PERS')
536 :
537 : ! Test: if one psp is PAW, all must be
538 1164 : if (test_paw==1) then
539 641 : do ipsp=1,npsp
540 641 : 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 1442 : 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
589 : !arrays
590 1442 : integer,allocatable :: list_int(:),lmax_(:)
591 : real(dp) :: tsec(2)
592 1442 : real(dp),allocatable :: list_dpr(:)
593 1442 : character(len=fnlen),allocatable :: list_char(:)
594 : #endif
595 : !*************************************************************************
596 :
597 : #if defined HAVE_MPI
598 1442 : call timab(48,1,tsec)
599 :
600 1442 : comm = xmpi_world
601 :
602 : ! Broadcast the characters (file names and titles)
603 4326 : ABI_MALLOC(list_char,(3*npsp))
604 3394 : list_char(1:npsp)=pspheads(1:npsp)%filpsp
605 3394 : list_char(npsp+1:2*npsp)=pspheads(1:npsp)%title
606 3394 : list_char(2*npsp+1:3*npsp)=pspheads(1:npsp)%md5_checksum
607 :
608 1442 : call xmpi_bcast(list_char,master,comm,ierr)
609 :
610 3394 : pspheads(1:npsp)%filpsp=list_char(1:npsp)
611 3394 : pspheads(1:npsp)%title=list_char(npsp+1:2*npsp)
612 3394 : pspheads(1:npsp)%md5_checksum=list_char(2*npsp+1:3*npsp)(1:md5_slen)
613 1442 : ABI_FREE(list_char)
614 :
615 : ! Brodcast the integers
616 1442 : list_size=6*npsp+1
617 4326 : ABI_MALLOC(list_int,(list_size))
618 3394 : list_int(1 : npsp) = pspheads(1:npsp)%lmax
619 3394 : list_int(1+ npsp: 2*npsp) = pspheads(1:npsp)%xccc
620 3394 : list_int(1+ 2*npsp: 3*npsp) = pspheads(1:npsp)%pspxc
621 3394 : list_int(1+ 3*npsp: 4*npsp) = pspheads(1:npsp)%pspdat
622 3394 : list_int(1+ 4*npsp: 5*npsp) = pspheads(1:npsp)%pspcod
623 3394 : list_int(1+ 5*npsp: 6*npsp) = pspheads(1:npsp)%pspso
624 1442 : list_int(6*npsp+1) = test_paw
625 :
626 1442 : call xmpi_bcast(list_int,master,comm,ierr)
627 :
628 3394 : pspheads(1:npsp)%lmax = list_int(1 : npsp)
629 3394 : pspheads(1:npsp)%xccc = list_int(1+ npsp: 2*npsp)
630 3394 : pspheads(1:npsp)%pspxc = list_int(1+ 2*npsp: 3*npsp)
631 3394 : pspheads(1:npsp)%pspdat = list_int(1+ 3*npsp: 4*npsp)
632 3394 : pspheads(1:npsp)%pspcod = list_int(1+ 4*npsp: 5*npsp)
633 3394 : pspheads(1:npsp)%pspso = list_int(1+ 5*npsp: 6*npsp)
634 1442 : test_paw = list_int(6*npsp+1)
635 :
636 1442 : ABI_FREE(list_int)
637 :
638 1442 : list_size=0
639 4326 : ABI_MALLOC(lmax_,(npsp))
640 3394 : do ipsp=1,npsp
641 1952 : lmax_(ipsp)=max(pspheads(ipsp)%lmax,0)
642 1952 : if(.not.allocated(pspheads(ipsp)%nproj)) then
643 1140 : ABI_MALLOC(pspheads(ipsp)%nproj,(0:lmax_(ipsp)))
644 : endif
645 1952 : if(.not.allocated(pspheads(ipsp)%nprojso)) then
646 1140 : ABI_MALLOC(pspheads(ipsp)%nprojso,(lmax_(ipsp)))
647 : endif
648 3394 : list_size=list_size+2*lmax_(ipsp)+1
649 : enddo
650 4326 : ABI_MALLOC(list_int,(list_size))
651 9002 : list_int=0
652 : ii=0
653 3394 : do ipsp=1,npsp
654 8150 : do il=0,lmax_(ipsp)
655 4756 : ii=ii+1
656 6708 : list_int(ii) = pspheads(ipsp)%nproj(il)
657 : enddo
658 : enddo
659 3394 : do ipsp=1,npsp
660 6198 : do il=1,lmax_(ipsp)
661 2804 : ii=ii+1
662 4756 : list_int(ii) = pspheads(ipsp)%nprojso(il)
663 : enddo
664 : enddo
665 :
666 1442 : call xmpi_bcast(list_int,master,comm,ierr)
667 :
668 1442 : ii=0
669 3394 : do ipsp=1,npsp
670 8150 : do il=0,lmax_(ipsp)
671 4756 : ii=ii+1
672 6708 : pspheads(ipsp)%nproj(il) = list_int(ii)
673 : enddo
674 : enddo
675 3394 : do ipsp=1,npsp
676 6198 : do il=1,lmax_(ipsp)
677 2804 : ii=ii+1
678 4756 : pspheads(ipsp)%nprojso(il) = list_int(ii)
679 : enddo
680 : enddo
681 1442 : ABI_FREE(list_int)
682 1442 : ABI_FREE(lmax_)
683 :
684 : ! Unbeliveable, this cannot be sent with the others, for woopy
685 2884 : ABI_MALLOC(list_int,(npsp))
686 3394 : list_int(1:npsp) = pspheads(1:npsp)%usewvl
687 1442 : call xmpi_bcast(list_int,master,comm,ierr)
688 3394 : pspheads(1:npsp)%usewvl = list_int(1:npsp)
689 1442 : ABI_FREE(list_int)
690 :
691 : ! Broadcast zionpsp and znuclpsp
692 4326 : ABI_MALLOC(list_dpr,(7*npsp))
693 3394 : list_dpr(1 : npsp) = pspheads(1:npsp)%zionpsp
694 3394 : list_dpr(1+ npsp:2*npsp) = pspheads(1:npsp)%znuclpsp
695 3394 : list_dpr(1+2*npsp:3*npsp) = pspheads(1:npsp)%GTHradii(0)
696 3394 : list_dpr(1+3*npsp:4*npsp) = pspheads(1:npsp)%GTHradii(1)
697 3394 : list_dpr(1+4*npsp:5*npsp) = pspheads(1:npsp)%GTHradii(2)
698 3394 : list_dpr(1+5*npsp:6*npsp) = pspheads(1:npsp)%GTHradii(3)
699 3394 : list_dpr(1+6*npsp:7*npsp) = pspheads(1:npsp)%GTHradii(4)
700 :
701 1442 : call xmpi_bcast(list_dpr,master,comm,ierr)
702 :
703 3394 : pspheads(1:npsp)%zionpsp = list_dpr(1 : npsp)
704 3394 : pspheads(1:npsp)%znuclpsp = list_dpr(1+ npsp:2*npsp)
705 3394 : pspheads(1:npsp)%GTHradii(0) = list_dpr(1+2*npsp:3*npsp)
706 3394 : pspheads(1:npsp)%GTHradii(1) = list_dpr(1+3*npsp:4*npsp)
707 3394 : pspheads(1:npsp)%GTHradii(2) = list_dpr(1+4*npsp:5*npsp)
708 3394 : pspheads(1:npsp)%GTHradii(3) = list_dpr(1+5*npsp:6*npsp)
709 3394 : pspheads(1:npsp)%GTHradii(4) = list_dpr(1+6*npsp:7*npsp)
710 1442 : ABI_FREE(list_dpr)
711 :
712 : ! Broadcast additional integers for PAW psps (testpaw was sent, previously)
713 1442 : if (test_paw==1) then
714 1548 : ABI_MALLOC(list_int,(6*npsp))
715 955 : list_int(1 : npsp)=pspheads(1:npsp)%pawheader%basis_size
716 955 : list_int(1+ npsp:2*npsp)=pspheads(1:npsp)%pawheader%l_size
717 955 : list_int(1+2*npsp:3*npsp)=pspheads(1:npsp)%pawheader%lmn_size
718 955 : list_int(1+3*npsp:4*npsp)=pspheads(1:npsp)%pawheader%mesh_size
719 955 : list_int(1+4*npsp:5*npsp)=pspheads(1:npsp)%pawheader%pawver
720 955 : list_int(1+5*npsp:6*npsp)=pspheads(1:npsp)%pawheader%shape_type
721 :
722 387 : call xmpi_bcast(list_int,master,comm,ierr)
723 :
724 955 : pspheads(1:npsp)%pawheader%basis_size=list_int(1 : npsp)
725 955 : pspheads(1:npsp)%pawheader%l_size =list_int(1+ npsp:2*npsp)
726 955 : pspheads(1:npsp)%pawheader%lmn_size =list_int(1+2*npsp:3*npsp)
727 955 : pspheads(1:npsp)%pawheader%mesh_size =list_int(1+3*npsp:4*npsp)
728 955 : pspheads(1:npsp)%pawheader%pawver =list_int(1+4*npsp:5*npsp)
729 955 : pspheads(1:npsp)%pawheader%shape_type=list_int(1+5*npsp:6*npsp)
730 387 : ABI_FREE(list_int)
731 :
732 : ! broadcast rpaw values
733 1548 : ABI_MALLOC(list_dpr,(2*npsp))
734 :
735 955 : list_dpr(1 : npsp) = pspheads(1:npsp)%pawheader%rpaw
736 955 : list_dpr(1+1*npsp:2*npsp) = pspheads(1:npsp)%pawheader%rshp
737 :
738 387 : call xmpi_bcast(list_dpr,master,comm,ierr)
739 :
740 955 : pspheads(1:npsp)%pawheader%rpaw = list_dpr(1 : npsp)
741 955 : pspheads(1:npsp)%pawheader%rshp = list_dpr(1+ npsp:2*npsp)
742 :
743 387 : ABI_FREE(list_dpr)
744 : end if
745 :
746 1442 : call timab(48,2,tsec)
747 :
748 : #else
749 : ! Code to use unused dummy arguments
750 : if(pspheads(1)%lmax == -10) pspheads(1)%lmax=-10
751 : if(test_paw == -1) test_paw = -1
752 : #endif
753 :
754 1442 : end subroutine pspheads_comm
755 : !!***
756 :
757 : !!****f* m_pspheads/pawpsxml2ab
758 : !! NAME
759 : !! pawpsxml2ab
760 : !!
761 : !! FUNCTION
762 : !! From a XML format pseudopotential file which has already been read in,
763 : !! convert to abinit internal datastructures.
764 : !!
765 : !! INPUTS
766 : !! ecut_tmp(3,2)= possible ecut values as read in psp files
767 : !! filenam= input file name (atomicdata XML)
768 : !! option= 1 if header only is read; 0 if the whole data are read
769 : !!
770 : !! OUTPUT
771 : !! pspheads data structure is filled
772 : !!
773 : !! SOURCE
774 :
775 494 : subroutine pawpsxml2ab(filnam, ecut_tmp, pspheads, option)
776 :
777 : !Arguments ------------------------------------
778 : !scalars
779 : integer, intent(in) :: option
780 : character(len=fnlen), intent(in) :: filnam
781 : type(pspheader_type),intent(inout) :: pspheads !vz_i
782 : !arrays
783 : real(dp),intent(inout) :: ecut_tmp(3,2)
784 :
785 : !Local variables-------------------------------
786 : integer :: ii,il,lloc,lmax,pspcod,pspxc
787 : real(dp) :: r2well,zionpsp,znuclpsp
788 : ! character(len=100) :: xclibxc, msg
789 : ! *********************************************************************
790 :
791 494 : if (option==1) then
792 188 : call rdpawpsxml_header(ecut_tmp,filnam,paw_setuploc)
793 188 : paw_setuploc%idgrid= paw_setuploc%radial_grid(1)%id
794 : else
795 306 : call rdpawpsxml(filnam,paw_setuploc)
796 : end if
797 :
798 494 : call pawpsp_read_header_xml(lloc,lmax,pspcod, pspxc,paw_setuploc,r2well,zionpsp,znuclpsp)
799 :
800 494 : pspheads%lmax=lmax
801 494 : pspheads%pspxc=pspxc
802 494 : pspheads%zionpsp=zionpsp
803 494 : pspheads%znuclpsp=znuclpsp
804 :
805 : call pawpsp_read_pawheader(pspheads%pawheader%basis_size,&
806 : pspheads%lmax,pspheads%pawheader%lmn_size,&
807 : pspheads%pawheader%l_size,pspheads%pawheader%mesh_size,&
808 : pspheads%pawheader%pawver,paw_setuploc,pspheads%pawheader%rpaw,&
809 494 : pspheads%pawheader%rshp,pspheads%pawheader%shape_type)
810 :
811 1482 : ABI_MALLOC(pspheads%nproj,(0:lmax))
812 1482 : ABI_MALLOC(pspheads%nprojso,(lmax))
813 1587 : pspheads%nproj=0
814 1587 : do il=0,pspheads%lmax
815 6716 : do ii=1,pspheads%pawheader%basis_size
816 6222 : if(paw_setuploc%valence_states%state(ii)%ll==il) pspheads%nproj(il)=pspheads%nproj(il)+1
817 : end do
818 : end do
819 :
820 1093 : pspheads%nprojso=0
821 494 : pspheads%pspdat=27061961
822 494 : pspheads%pspso=1
823 494 : pspheads%xccc=1
824 494 : pspheads%title=paw_setuploc%atom%symbol
825 :
826 494 : if (option==1) call paw_setup_free(paw_setuploc)
827 :
828 494 : end subroutine pawpsxml2ab
829 : !!***
830 :
831 : !!****f* m_pspheads/upf1_to_psphead
832 : !! NAME
833 : !! upf1_to_psphead
834 : !!
835 : !! FUNCTION
836 : !! This routine wraps a call to a PWSCF module, which reads in
837 : !! a UPF1 (PWSCF / Espresso) format pseudopotential, then transfers
838 : !! data for the HEADER of abinit psps only!
839 : !!
840 : !! INPUTS
841 : !! filpsp = name of file with UPF1 data
842 : !!
843 : !! OUTPUT
844 : !! pspxc = index of xc functional for this pseudo
845 : !! lmax_ = maximal angular momentum
846 : !! znucl = charge of species nucleus
847 : !! zion = valence charge
848 : !! n1xccc = default number of points. Set to 0 if no nlcc is present
849 : !! nproj_l= number of projectors for each channel
850 : !! nprojso_l= number of projectors for each channel for SO correction projectors
851 : !!
852 : !! SOURCE
853 :
854 3 : subroutine upf1_to_psphead(filpsp, znucl, zion, pspxc, lmax_, n1xccc, nproj_l, nprojso_l)
855 :
856 : use m_read_upf_pwscf, only : read_pseudo
857 : use pseudo_pwscf ! pwscf module with all data explicit!
858 :
859 : !Arguments -------------------------------
860 : character(len=fnlen), intent(in) :: filpsp
861 : integer,intent(inout) :: n1xccc
862 : integer,intent(out) :: pspxc, lmax_
863 : real(dp),intent(out) :: znucl, zion
864 : !arrays
865 : integer,intent(out) :: nproj_l(0:3)
866 : integer,intent(out) :: nprojso_l(1:3)
867 :
868 : !Local variables -------------------------
869 : integer :: iproj, ll, iunit
870 : character(len=500) :: msg
871 : type(atomdata_t) :: atom
872 :
873 : ! *********************************************************************
874 :
875 1 : if (open_file(filpsp, msg, newunit=iunit, status='old',form='formatted') /= 0) then
876 0 : ABI_ERROR(msg)
877 : end if
878 :
879 : ! read in psp data to static data in pseudo module, for ipsx == 1
880 1 : call read_pseudo(1,iunit)
881 1 : close (iunit)
882 :
883 : ! copy over to abinit internal arrays and vars
884 : ! FIXME: The API is broken. It does not recognize PBEsol
885 : ! should use upfdft_to_ixc
886 1 : call upfxc2abi(dft(1), pspxc)
887 1 : lmax_ = lmax(1)
888 1 : call atomdata_from_symbol(atom,psd(1))
889 1 : znucl = atom%znucl
890 1 : zion = zp(1)
891 :
892 1 : nproj_l = 0
893 4 : do iproj = 1, nbeta(1)
894 3 : ll = lll(iproj,1)
895 4 : nproj_l(ll) = nproj_l(ll) + 1
896 : end do
897 :
898 1 : nprojso_l = 0 !FIXME deal with so
899 : !do iproj = 1, nbeta(1)
900 : !nprojso_l(ll+1) = nprojso_l(ll+1) + 1
901 : !end do
902 :
903 1 : if (.not. nlcc(1)) n1xccc = 0
904 :
905 1 : end subroutine upf1_to_psphead
906 : !!***
907 :
908 : !!****f* m_pspheads/upf2_to_psphead
909 : !! NAME
910 : !! upf2_to_psphead
911 : !!
912 : !! FUNCTION
913 : !! This routine wraps a call to a PWSCF module, which reads in
914 : !! a UPF2 (PWSCF / Espresso) format pseudopotential, then transfers
915 : !! data for the HEADER of abinit psps only!
916 : !!
917 : !! INPUTS
918 : !! filpsp = name of file with UPF1 data
919 : !!
920 : !! OUTPUT
921 : !! pspxc = index of xc functional for this pseudo
922 : !! lmax_ = maximal angular momentum
923 : !! znucl = charge of species nucleus
924 : !! zion = valence charge
925 : !! n1xccc = default number of points. Set to 0 if no nlcc is present
926 : !! nproj_l= number of projectors for each channel
927 : !! nprojso_l= number of projectors for each channel for SO correction projectors
928 : !!
929 : !! SOURCE
930 :
931 6 : subroutine upf2_to_psphead(filpsp, znucl, zion, pspxc, lmax, n1xccc, nproj_l, nprojso_l)
932 :
933 : !Arguments -------------------------------
934 : character(len=fnlen), intent(in) :: filpsp
935 : integer,intent(inout) :: n1xccc
936 : integer,intent(out) :: pspxc, lmax
937 : real(dp),intent(out) :: znucl, zion
938 : !arrays
939 : integer,intent(out) :: nproj_l(0:3), nprojso_l(1:3)
940 :
941 : !Local variables -------------------------
942 : integer :: ierr , iprj, ll, mmax, irad
943 : real(dp) :: amesh, damesh
944 : character(len=500) :: msg
945 : logical :: linear_mesh
946 6 : type(pseudo_upf) :: upf
947 : type(atomdata_t) :: atom
948 : ! arrays
949 6 : real(dp),allocatable :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:)
950 : ! *********************************************************************
951 :
952 : ! See also https://github.com/QEF/qeschemas/blob/master/UPF/qe_pp-0.99.xsd
953 6 : call read_upf_new(filpsp, upf, ierr)
954 6 : ABI_CHECK(ierr == 0, sjoin("read_upf_new returned ierr:", itoa(ierr)))
955 :
956 6 : call atomdata_from_symbol(atom, upf%psd)
957 6 : znucl = atom%znucl
958 6 : zion = upf%zp
959 6 : lmax = upf%lmax
960 6 : mmax = upf%mesh
961 :
962 : ! Consistency check
963 6 : ABI_CHECK(upf%typ == "NC", sjoin("Only NC pseudos in UPF2 format are supported while type is:", upf%typ))
964 6 : ABI_CHECK(upfdft_to_ixc(upf%dft, pspxc, msg) == 0, msg)
965 6 : if (.not. upf%nlcc) n1xccc = 0
966 :
967 : ! Check that rad grid is linear starting at zero
968 6 : linear_mesh = .True.
969 6 : amesh = upf%r(2) - upf%r(1); damesh = zero
970 8730 : do irad=2,mmax-1
971 8730 : damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1)))
972 : end do
973 6 : linear_mesh = damesh < tol8
974 :
975 6 : if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then
976 : write(msg,'(3a)')&
977 0 : 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
978 0 : 'Action: check your pseudopotential file.'
979 0 : ABI_ERROR(msg)
980 : end if
981 :
982 6 : nproj_l = 0; nprojso_l = 0
983 :
984 6 : if (.not. upf%has_so) then
985 : ! Scalar case
986 27 : do iprj=1,upf%nbeta
987 23 : ll = upf%lll(iprj)
988 27 : nproj_l(ll) = nproj_l(ll) + 1
989 : end do
990 :
991 : else
992 : ! Pseudo in j = l + s representation.
993 2 : call upf2_jl2srso(upf, nproj_l, nprojso_l, vsr, esr, vso, eso)
994 :
995 2 : ABI_FREE(vsr)
996 2 : ABI_FREE(esr)
997 2 : ABI_FREE(vso)
998 2 : ABI_FREE(eso)
999 : end if
1000 :
1001 6 : call deallocate_pseudo_upf(upf)
1002 :
1003 6 : end subroutine upf2_to_psphead
1004 : !!***
1005 :
1006 : !!****f* m_pspheads/upf2_jl2srso
1007 : !! NAME
1008 : !! upf2_jl2srso
1009 : !!
1010 : !! FUNCTION
1011 : !!
1012 : !! INPUTS
1013 : !!
1014 : !! OUTPUT
1015 : !!
1016 : !! SOURCE
1017 :
1018 6 : subroutine upf2_jl2srso(upf, nproj_l, nprojso_l, vsr, esr, vso, eso)
1019 :
1020 : !Arguments -------------------------------
1021 : type(pseudo_upf),intent(in) :: upf
1022 : !arrays
1023 : integer,intent(out) :: nproj_l(0:3), nprojso_l(1:3)
1024 : real(dp),allocatable,intent(out) :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:)
1025 :
1026 : !Local variables -------------------------
1027 : integer :: iprj, ii, ll, l1, il, ik, lmax, mmax, mxprj
1028 : real(dp) :: jtot, eprmin !eps_srso,
1029 : !character(len=500) :: msg
1030 : ! arrays
1031 : integer :: irc6(6),nproj6(6), done_ilk(6,2)
1032 6 : real(dp),allocatable :: vkb(:,:,:,:), evkb(:,:,:)
1033 : ! *********************************************************************
1034 :
1035 6 : lmax = upf%lmax; mmax = upf%mesh
1036 6 : nproj_l = 0; nprojso_l = 0
1037 :
1038 : ! Pseudo in j = l + s representation.
1039 6 : irc6 = zero; nproj6 = zero
1040 66 : do iprj=1,upf%nbeta
1041 60 : ll = upf%lll(iprj)
1042 60 : nproj6(ll+1) = nproj6(ll+1) + 1
1043 : !irc6(ll+1) = max(upf%kbeta(iprj), irc6(ll+1))
1044 66 : irc6(ll+1) = mmax
1045 : end do
1046 :
1047 : ! Divide by two for l > 0 as this is sr_so_r expects.
1048 36 : nproj6(2:) = nproj6(2:) / 2
1049 42 : mxprj = maxval(nproj6)
1050 :
1051 36 : ABI_MALLOC(vkb, (mmax,mxprj,4,2))
1052 24 : ABI_MALLOC(evkb, (mxprj,4,2))
1053 30 : ABI_MALLOC(vsr, (mmax,2*mxprj,4))
1054 24 : ABI_MALLOC(esr, (2*mxprj,4))
1055 18 : ABI_MALLOC(vso, (mmax,2*mxprj,4))
1056 12 : ABI_MALLOC(eso, (2*mxprj,4))
1057 :
1058 6 : done_ilk = 0
1059 66 : do iprj=1,upf%nbeta
1060 60 : jtot = upf%jjj(iprj)
1061 60 : ll = upf%lll(iprj)
1062 60 : il = ll + 1
1063 60 : if (ll == 0) then
1064 : ik = 1
1065 : else
1066 : ! l+1/2 --> ik 1, l-1/2 --> ik 2
1067 48 : if (abs(jtot - (ll + half)) < tol6) then
1068 : ik = 1
1069 24 : else if (abs(jtot - (ll - half)) < tol6) then
1070 : ik = 2
1071 : else
1072 0 : ABI_ERROR(sjoin("Cannot detect ik index from jtot:", ftoa(jtot)))
1073 : end if
1074 : end if
1075 :
1076 60 : done_ilk(il, ik) = done_ilk(il, ik) + 1
1077 60 : ii = done_ilk(il, ik)
1078 60 : evkb(ii,il,ik) = upf%dion(iprj,iprj) * half ! convert from Rydberg to Ha
1079 98586 : vkb(:,ii,il,ik) = upf%beta(:,iprj)
1080 : end do
1081 :
1082 6 : call sr_so_r(lmax, irc6, nproj6, upf%r, mmax, mxprj, evkb, vkb, vsr, esr, vso, eso)
1083 :
1084 : ! MG: This is done in oncvpsp 3.3 but not in oncvpsp4
1085 : ! drop sr, so orthonormal projectors with neglibible coefficients
1086 : ! modify cutoff if desired
1087 :
1088 6 : eprmin=2.0d-5
1089 6 : write(std_out,'(/a,1p,e10.2,a)') 'Orthonormal projectors with coefficients <', &
1090 12 : eprmin,' Ha will be dropped'
1091 :
1092 24 : do l1=1,lmax+1
1093 18 : if(abs(esr(3,l1))<eprmin) esr(3,l1)=0.0d0
1094 18 : if(abs(esr(4,l1))<eprmin) esr(4,l1)=0.0d0
1095 18 : if(abs(eso(3,l1))<eprmin) eso(3,l1)=0.0d0
1096 24 : if(abs(eso(4,l1))<eprmin) eso(4,l1)=0.0d0
1097 : end do
1098 :
1099 : #if 0
1100 : ! MG: This is done in oncvpsp 4 but not in oncvpsp 3.3
1101 : ! set smallest components to zero (following the approach used in oncvpsp)
1102 : eps_srso=1.0d-3
1103 : do l1=1,lmax+1
1104 : if (nproj6(l1) > 0) then
1105 : do iprj=2,2*nproj6(l1)
1106 : if (abs(esr(iprj,l1)) < eps_srso*abs(esr(1,l1))) esr(iprj,l1) = 0.0d0
1107 : if (l1 == 1) cycle
1108 : if (abs(eso(iprj,l1)) < eps_srso*abs(eso(1,l1))) eso(iprj,l1) = 0.0d0
1109 : end do
1110 : end if
1111 : end do
1112 : #endif
1113 :
1114 : ! set up projector number for sr_so calculations based on non-zero coefficients
1115 : ! note that energies and projectors have been sorted sr_so_r
1116 : ! so the relevant projectors are packed in the first positions.
1117 24 : do l1=1,lmax+1
1118 18 : ll = l1 - 1
1119 96 : do ii=1,2*nproj6(l1)
1120 72 : if (abs(esr(ii,l1)) > 0.0d0) nproj_l(ll) = nproj_l(ll) + 1
1121 90 : if (abs(eso(ii,l1)) > 0.0d0) nprojso_l(ll) = nprojso_l(ll) + 1
1122 : end do
1123 : !write(std_out, '(a,3(i0,1x))')' ll, nproj_l, nprojso_l',ll, nproj_l(ll), nprojso_l(ll)
1124 : end do
1125 :
1126 6 : ABI_FREE(vkb)
1127 6 : ABI_FREE(evkb)
1128 :
1129 6 : end subroutine upf2_jl2srso
1130 : !!***
1131 :
1132 : !!****f* m_pspheads/upfxc2abi
1133 : !! NAME
1134 : !! upfxc2abi
1135 : !!
1136 : !! FUNCTION
1137 : !! This routine wraps a call to an OCTOPUS module, which reformats
1138 : !! a UPF (PWSCF / Espresso) string describing XC functionals,
1139 : !! and returns the abinit internal code pspxc
1140 : !!
1141 : !! INPUTS
1142 : !! dft = string with x/c functionals from PWSCF format
1143 : !!
1144 : !! OUTPUT
1145 : !! pspxc = index of xc functional for this pseudo
1146 : !!
1147 : !! NOTES
1148 : !! FIXME: extend to more functionals with libxc
1149 : !! Could be included in separate module, eg read_upf_pwscf or funct_pwscf
1150 : !! Left without defs_basis or calls to abinit routines ON PURPOSE
1151 : !!
1152 : !! SOURCE
1153 :
1154 2 : subroutine upfxc2abi(dft, pspxc)
1155 :
1156 : use funct_pwscf ! pwscf module for naming xc functionals
1157 :
1158 : !Arguments -------------------------------
1159 : character(len=*), intent(in) :: dft
1160 : integer, intent(out) :: pspxc
1161 :
1162 : !Local variables -------------------------
1163 : integer :: iexch,icorr,igcx,igcc
1164 : integer :: totalindex, offset
1165 : ! *********************************************************************
1166 :
1167 : !extract from char*20 :: dft(:)
1168 : !### The following has been copied from pwscf src/Modules/upf_to_internal.f90:
1169 : !workaround for rrkj format - it contains the indices, not the name
1170 2 : if ( dft(1:6)=='INDEX:') then
1171 0 : read( dft(7:10), '(4i1)') iexch,icorr,igcx,igcc
1172 0 : call set_dft_from_indices(iexch,icorr,igcx,igcc)
1173 : else
1174 2 : call set_dft_from_name( dft )
1175 2 : iexch = get_iexch()
1176 2 : icorr = get_icorr()
1177 2 : igcx = get_igcx()
1178 2 : igcc = get_igcc()
1179 : end if
1180 :
1181 : !reset dft string to avoid stray spaces
1182 2 : call set_dft_from_indices(iexch,icorr,igcx,igcc)
1183 2 : write(std_out,'(a)') ' upf2abinit: XC string from pseudopotential is :'
1184 2 : write(std_out,'(3a)') '>', dft, '<'
1185 2 : ABI_WARNING("upfxc2abi is not guaranteed to return the right ixc from QE XC string e.g. PBEsol. Please crosscheck!")
1186 :
1187 2 : offset = 100
1188 2 : totalindex = offset*offset*offset*iexch + offset*offset*icorr + offset*igcx + igcc
1189 0 : select case (totalindex)
1190 : case (00000000) !(" NOX NOC NOGX NOGC") ! no xc
1191 0 : pspxc = 0
1192 : case (01010000) !(" SLA PZ NOGX NOGC") ! slater exchange + Perdew Zunger
1193 0 : pspxc = 2
1194 : case (01050000) !(" SLA WIG NOGX NOGC") ! slater exchange + Wigner corr
1195 0 : pspxc = 4
1196 : case (01060000) !(" SLA HL NOGX NOGC") ! Hedin + Lundqvist
1197 0 : pspxc = 5
1198 : case (02000000) !(" SL1 NOC NOGX NOGC") ! full slater exchange
1199 0 : pspxc = 6
1200 : case (01040000) !(" SLA PW NOGX NOGC") ! slater exchange + Perdew Wang
1201 2 : pspxc = 7
1202 : case (01000000) !(" SLA NOC NOGX NOGC") ! Perdew Wang + no corr
1203 0 : pspxc = 8
1204 : case (01040304) !(" SLA PW PBX PBC") ! LDA + PBE GGA
1205 0 : pspxc = 11 ! PBE
1206 : case (01000300) !(" SLA NOC PBX NOGC") ! exchange part of PBE GGA
1207 0 : pspxc = 12
1208 : case (01040404) !(" SLA PW RPB PBC") ! rev PBE
1209 0 : pspxc = 14
1210 : case (00000505) !(" NOX NOC HTCH HTCH") ! HTCH 120
1211 0 : pspxc = 17
1212 : case (01030103) !(" SLA LYP B88 BLYP") ! BLYP
1213 0 : pspxc = -106131
1214 : case (01040101) !(" SLA PW B88 P86") ! BP86
1215 0 : pspxc = -106132
1216 : case (00030603) !(" NOX LYP OPTX BLYP") ! OLYP
1217 0 : pspxc = -110131
1218 : ! case () !("R2SCAN01") ! R2SCAN01
1219 : ! pspxc = -645642
1220 : ! FIXME: important cases left to be patched with libxc:
1221 : ! vosko wilkins nusair
1222 : ! ortiz ballone
1223 : ! pbe0
1224 : ! Gunnarson-Lunqvist
1225 : ! make general approach: check gradient parts first, then lda.
1226 : ! event. check if they are consistent.
1227 : case default
1228 2 : ABI_ERROR('upf2abinit: XC functional not recognized')
1229 : end select
1230 :
1231 2 : end subroutine upfxc2abi
1232 : !!***
1233 :
1234 : !!****f* m_pspheads/updft_to_ixc
1235 : !! NAME
1236 : !! updft_to_ixc
1237 : !!
1238 : !! FUNCTION
1239 : !! Returns the abinit internal `ixc` from `dft` string with XC functional in QE format.
1240 : !!
1241 : !! SOURCE
1242 :
1243 16 : integer function upfdft_to_ixc(dft, ixc, msg) result(ierr)
1244 :
1245 : !Arguments ------------------------------------
1246 : character(len=*),intent(in) :: dft
1247 : character(len=*),intent(out) :: msg
1248 : integer,intent(out) :: ixc
1249 :
1250 : !Local variables-------------------------------
1251 : integer :: start !, ii
1252 : character(len=500) :: x_name, c_name, gcx_name, gcc_name
1253 : !*************************************************************************
1254 :
1255 : ! This list taken from oncvpsp/src/upfout.f90
1256 : ! It should be OK as long as the UPF2 NC pseudos are generated with oncvpsp
1257 : ! but it does not cover all QE possibilities.
1258 16 : ierr = 0; msg = ""
1259 16 : ixc = 0
1260 0 : select case (dft)
1261 : case ("PZ")
1262 0 : ixc = -001009
1263 : case ("PBE")
1264 12 : ixc = -101130 !; ixc = 11
1265 : case ("PW91")
1266 0 : ixc = -109134
1267 : case ("PBESOL")
1268 0 : ixc = -116133
1269 : case ("REVPBE")
1270 0 : ixc = -102130
1271 : case ("BP")
1272 0 : ixc = -106132
1273 : case ("BLYP")
1274 0 : ixc = -106131
1275 : case ("WC")
1276 0 : ixc = -118130
1277 : case ('SLA PW NOGX NOGC', "SLA+PW") ! string produced by oncvpsp3
1278 0 : ixc = -1012
1279 : case ('R2SCAN01')
1280 4 : ixc = -645642
1281 : case default
1282 16 : ierr = 1
1283 : end select
1284 :
1285 : ! Extract substrings with
1286 : ! 1) exchange
1287 : ! 2) correlation
1288 : ! 3) gradient correction, exchange
1289 : ! 4) gradient correction, correlation
1290 : if (ierr == 1) then
1291 0 : ierr = 0; start = 1
1292 0 : ABI_CHECK(next_token(dft, start, x_name) == 0 , "Error reading x_name")
1293 0 : ABI_CHECK(next_token(dft, start, c_name) == 0 , "Error reading c_name")
1294 0 : ABI_CHECK(next_token(dft, start, gcx_name) == 0 , "Error reading gcx_name")
1295 0 : ABI_CHECK(next_token(dft, start, gcc_name) == 0 , "Error reading gcc_name")
1296 : !call remove_non_ascii(gcc_name)
1297 : !print *, "dft: `", trim(dft), "`"
1298 : !print *, "x_name: `", trim(x_name), "`, c_name: `", trim(c_name), &
1299 : ! "`, gcx_name: `", trim(gcx_name), "`, gcc_name: `", trim(gcc_name), "`"
1300 :
1301 0 : if (x_name == "SLA" .and. c_name == "PW") then
1302 : !print *, "in first if", gcx_name == "NOGX", trim(gcc_name) == "NOGC"
1303 : !print *, "len_trim(gcc_name)", len_trim(gcc_name)
1304 0 : if (gcx_name == "NOGX" .and. gcc_name == "NOGC") then
1305 0 : ixc = -1012
1306 : else
1307 : ierr = 1
1308 : !print *, "in second ierr"
1309 : end if
1310 : else
1311 : ierr = 1
1312 : end if
1313 : end if
1314 :
1315 : if (ierr == 1) then
1316 : write(msg, "(5a)") &
1317 0 : "Cannot find ABINIT ixc value corresponding to QE dft string: `", trim(dft), "`", ch10, &
1318 0 : "Please update mapping in m_pspheads/upfdft_to_ixc."
1319 : end if
1320 :
1321 16 : end function upfdft_to_ixc
1322 : !!***
1323 :
1324 : ! Copyright (c) 1989-2019 by D. R. Hamann, Mat-Sim Research LLC and Rutgers
1325 : ! University
1326 : !
1327 : !
1328 : ! This program is free software: you can redistribute it and/or modify
1329 : ! it under the terms of the GNU General Public License as published by
1330 : ! the Free Software Foundation, either version 3 of the License, or
1331 : ! (at your option) any later version.
1332 : !
1333 : ! This program is distributed in the hope that it will be useful,
1334 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
1335 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1336 : ! GNU General Public License for more details.
1337 : !
1338 : ! You should have received a copy of the GNU General Public License
1339 : ! along with this program. If not, see <http://www.gnu.org/licenses/>.
1340 : !
1341 6 : subroutine sr_so_r(lmax,irc,nproj,rr,mmax,mxprj,evkb,vkb, &
1342 6 : & vsr,esr,vso,eso)
1343 :
1344 : ! reformulates non-local potentials based on j = l +/- 1/2 to scalar-
1345 : ! relativistic and L dot S projectors
1346 : ! uses relationship <L dot S> = (J^2 - L^2 - S^2)/2
1347 : ! so L dot S = +/- l/2 for j = l +/- 1/2
1348 :
1349 : !lmax maximum angular momentum
1350 : !irc core radii indices
1351 : !nproj number of projectors for each l
1352 : !rr log radial grid
1353 : !mmax size of radial grid
1354 : !mmax dimension of log grid
1355 : !mxprj dimension of number of projectors
1356 : !vkb vkb projectors
1357 : !evkb coefficients of BKB projectors
1358 : !vsr normalized scalar projectors
1359 : !esr energy coefficients of vscal
1360 : !vso normalized spin-orbig projectors
1361 : !esol energy coefficients of vso
1362 :
1363 : !implicit none
1364 : !integer, parameter :: dp=kind(1.0d0)
1365 :
1366 : !Input variables
1367 : integer,intent(in) :: lmax,mmax,mxprj
1368 : integer,intent(in) :: irc(6),nproj(6)
1369 : real(dp),intent(in) :: rr(mmax),vkb(mmax,mxprj,4,2),evkb(mxprj,4,2)
1370 :
1371 : !Output variables
1372 : real(dp),intent(out) :: vsr(mmax,2*mxprj,4),esr(2*mxprj,4)
1373 : real(dp),intent(out) :: vso(mmax,2*mxprj,4),eso(2*mxprj,4)
1374 :
1375 : !Local variables
1376 : integer :: ii,jj,kk,ik1,ik2,ip1,ip2,ipk,ll,l1,info,nn
1377 : real(dp) :: amesh
1378 : real(dp) :: apk,tt
1379 12 : real(dp) :: sovl(2*mxprj,2*mxprj),sovlev(2*mxprj),ascl(2*mxprj,2*mxprj),aso(2*mxprj,2*mxprj)
1380 12 : real(dp) :: asclst(2*mxprj,2*mxprj),wsclst(2*mxprj),asost(2*mxprj,2*mxprj),wsost(2*mxprj)
1381 12 : real(dp) :: asclt(2*mxprj,2*mxprj),asot(2*mxprj,2*mxprj)
1382 12 : real(dp) :: sphalf(2*mxprj,2*mxprj),smhalf(2*mxprj,2*mxprj)
1383 12 : real(dp) :: fscl(mxprj),fso(mxprj),work(10*mxprj)
1384 6 : real(dp), allocatable :: vkbt(:,:),vkbst(:,:)
1385 : logical :: sorted
1386 : character(len=500) :: msg
1387 :
1388 : ABI_UNUSED(irc(1))
1389 :
1390 : ! Check that rad grid is linear starting at zero
1391 : !linear_mesh = .True.
1392 6 : amesh = rr(2) - rr(1) !; damesh = zero
1393 : !do irad=2,mmax-1
1394 : ! damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1)))
1395 : !end do
1396 : !linear_mesh = damesh < tol8
1397 :
1398 : !if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then
1399 : ! write(msg,'(3a)')&
1400 : ! 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,&
1401 : ! 'Action: check your pseudopotential file.'
1402 : ! ABI_ERROR(msg)
1403 : !end if
1404 :
1405 : !allocate(vkbt(mmax,2*mxprj),vkbst(mmax,2*mxprj))
1406 30 : ABI_MALLOC(vkbt, (mmax,2*mxprj))
1407 18 : ABI_MALLOC(vkbst, (mmax,2*mxprj))
1408 :
1409 24 : do l1=1,lmax+1
1410 18 : ll=l1-1
1411 :
1412 18 : if(ll==0) then
1413 39438 : vsr(:,:,l1)=0.0d0
1414 39438 : vso(:,:,l1)=0.0d0
1415 30 : esr(:,l1)=0.0d0
1416 30 : eso(:,l1)=0.0d0
1417 6 : if(nproj(l1)>=1) then
1418 18 : do ii=1,nproj(l1)
1419 19716 : vsr(:,ii,l1)=vkb(:,ii,l1,1)
1420 18 : esr(ii,l1)=evkb(ii,l1,1)
1421 : end do
1422 : end if
1423 : cycle
1424 : end if
1425 :
1426 12 : nn=2*nproj(l1)
1427 :
1428 12 : fscl(1)=(ll+1)/dble(2*ll+1)
1429 12 : fscl(2)=ll/dble(2*ll+1)
1430 12 : fso(1)=2/dble(2*ll+1)
1431 12 : fso(2)=-2/dble(2*ll+1)
1432 :
1433 : ! construct overlap matrix and diagonal energy matrices
1434 :
1435 252 : sovl(:,:)=0.0d0
1436 252 : ascl(:,:)=0.0d0
1437 252 : aso(:,:)=0.0d0
1438 78876 : vkbt(:,:)=0.0d0
1439 :
1440 36 : do ik1=1,2
1441 84 : do ip1=1,nproj(l1)
1442 48 : ii=ip1+(ik1-1)*nproj(l1)
1443 :
1444 48 : ascl(ii,ii)=fscl(ik1)*evkb(ip1,l1,ik1)
1445 48 : aso(ii,ii)=fso(ik1)*evkb(ip1,l1,ik1)
1446 :
1447 78864 : vkbt(:,ii)=vkb(:,ip1,l1,ik1)
1448 :
1449 168 : do ik2=1,2
1450 336 : do ip2=1,nproj(l1)
1451 192 : jj=ip2+(ik2-1)*nproj(l1)
1452 :
1453 : ! MG: This routine cannot be used as it assumes log mesh.
1454 : !call vpinteg(vkb(1,ip1,l1,ik1),vkb(1,ip2,l1,ik2),irc(l1),2*l1, &
1455 : ! sovl(ii,jj),rr)
1456 :
1457 : ! So we replace it with simpson integration
1458 315552 : sovl(ii,jj) = simpson(amesh, vkb(:,ip1,l1,ik1) * vkb(:,ip2,l1,ik2))
1459 : end do
1460 : end do
1461 : end do
1462 : end do
1463 :
1464 12 : call dsyev( 'V', 'U', nn, sovl, 2*mxprj, sovlev, work, 10*mxprj, info )
1465 :
1466 12 : if(info .ne. 0) then
1467 0 : write(msg,'(a,i4)') 'sr_so_r: S matrix eigenvalue ERROR, info=',info
1468 0 : ABI_ERROR(msg)
1469 : end if
1470 :
1471 : ! construct S^(-1/2) AND s^(1/2)
1472 :
1473 60 : do jj=1,nn
1474 48 : tt=sqrt(sovlev(jj))
1475 252 : do ii=1,nn
1476 192 : sphalf(ii,jj)=tt*sovl(ii,jj)
1477 240 : smhalf(ii,jj)=sovl(ii,jj)/tt
1478 : end do
1479 : end do
1480 :
1481 : ! take linear combinations to form orthonormal basis functions
1482 :
1483 78876 : vkbst(:,:)=0.0d0
1484 :
1485 60 : do jj=1,nn
1486 252 : do ii=1,nn
1487 315504 : vkbst(:,jj)=vkbst(:,jj) + smhalf(ii,jj)*vkbt(:,ii)
1488 : end do
1489 : end do
1490 :
1491 : ! construct A^(-1)* = S^(1/2)^T A^(-1) S^(1/2)
1492 :
1493 252 : asclt(:,:)=0.0d0
1494 252 : asclst(:,:)=0.0d0
1495 252 : asot(:,:)=0.0d0
1496 252 : asost(:,:)=0.0d0
1497 :
1498 60 : do ii=1,nn
1499 252 : do jj=1,nn
1500 1008 : do kk=1,nn
1501 768 : asclt(ii,jj)=asclt(ii,jj)+ascl(ii,kk)*sphalf(kk,jj)
1502 960 : asot(ii,jj) =asot(ii,jj) + aso(ii,kk)*sphalf(kk,jj)
1503 : end do
1504 : end do
1505 : end do
1506 :
1507 60 : do ii=1,nn
1508 252 : do jj=1,nn
1509 1008 : do kk=1,nn
1510 768 : asclst(ii,jj)=asclst(ii,jj)+asclt(kk,jj)*sphalf(kk,ii)
1511 960 : asost(ii,jj) =asost(ii,jj) + asot(kk,jj)*sphalf(kk,ii)
1512 : end do
1513 : end do
1514 : end do
1515 :
1516 : ! find eigenvalues and eigenvectors of the A* matrices
1517 :
1518 : ! SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
1519 :
1520 12 : call dsyev( 'V', 'U', nn, asclst, 2*mxprj, wsclst, work, 10*mxprj, info )
1521 :
1522 12 : if(info .ne. 0) then
1523 0 : write(msg,'(a,i4)') 'sr_so_r: A* matrix eigenvalue ERROR, info=',info
1524 0 : ABI_ERROR(msg)
1525 : end if
1526 :
1527 12 : call dsyev( 'V', 'U', nn, asost, 2*mxprj, wsost, work, 10*mxprj, info )
1528 :
1529 :
1530 12 : if(info .ne. 0) then
1531 0 : write(msg,'(a,i4)') 'sr_so_r: A* matrix eigenvalue ERROR, info=',info
1532 0 : ABI_ERROR(msg)
1533 : end if
1534 :
1535 : ! take linear combinations to form orthonormal projectors
1536 :
1537 78876 : vsr(:,:,l1)=0.0d0
1538 78876 : vso(:,:,l1)=0.0d0
1539 60 : esr(:,l1)=0.0d0
1540 60 : eso(:,l1)=0.0d0
1541 :
1542 60 : do ii=1,nn
1543 48 : esr(ii,l1)=wsclst(ii)
1544 48 : eso(ii,l1)= wsost(ii)
1545 252 : do jj=1,nn
1546 315456 : vsr(:,ii,l1)=vsr(:,ii,l1) + asclst(jj,ii)*vkbst(:,jj)
1547 315504 : vso(:,ii,l1)= vso(:,ii,l1) + asost(jj,ii)*vkbst(:,jj)
1548 : end do
1549 : end do
1550 :
1551 : ! bubble-sort on coefficient magnitudes for scalar and then s-o
1552 : ! (Yes, I know bubble-sort is the least-efficient sorting algorithm.)
1553 :
1554 30 : do ii=1,100
1555 : sorted=.true.
1556 120 : do jj=2,nn
1557 120 : if(abs(esr(jj-1,l1))<abs(esr(jj,l1))) then
1558 36 : tt=esr(jj,l1)
1559 59148 : vkbt(:,1)=vsr(:,jj,l1)
1560 36 : esr(jj,l1)=esr(jj-1,l1)
1561 59148 : vsr(:,jj,l1)=vsr(:,jj-1,l1)
1562 36 : esr(jj-1,l1)=tt
1563 59148 : vsr(:,jj-1,l1)=vkbt(:,1)
1564 : sorted=.false.
1565 : end if
1566 : end do
1567 30 : if(sorted) exit
1568 : end do
1569 :
1570 48 : do ii=1,100
1571 : sorted=.true.
1572 192 : do jj=2,nn
1573 192 : if(abs(eso(jj-1,l1))<abs(eso(jj,l1))) then
1574 54 : tt=eso(jj,l1)
1575 88722 : vkbt(:,1)=vso(:,jj,l1)
1576 54 : eso(jj,l1)=eso(jj-1,l1)
1577 88722 : vso(:,jj,l1)=vso(:,jj-1,l1)
1578 54 : eso(jj-1,l1)=tt
1579 88722 : vso(:,jj-1,l1)=vkbt(:,1)
1580 : sorted=.false.
1581 : end if
1582 : end do
1583 48 : if(sorted) exit
1584 : end do
1585 :
1586 : write(std_out,'(/a,i2)') &
1587 12 : & ' Orthonormal scalar projector coefficients, l = ',ll
1588 12 : write(std_out,'(1p,6e12.4)') (esr(jj,l1),jj=1,nn)
1589 : write(std_out,'(/a,i2)') &
1590 12 : & ' Orthonormal spin-orbit projector coefficients, l = ',ll
1591 12 : write(std_out,'(1p,6e12.4)') (eso(jj,l1),jj=1,nn)
1592 :
1593 : ! Set sign of projectors (physically irrelevant) so that they are positive
1594 : ! at their peak (needed for compaisons apparently)
1595 :
1596 66 : do jj=1,nn
1597 : apk=0.0d0
1598 78864 : do ii=1,mmax
1599 78864 : if(abs(vso(ii,jj,l1))>apk) then
1600 3102 : apk=abs(vso(ii,jj,l1))
1601 3102 : ipk=ii
1602 : end if
1603 : end do
1604 48 : if(vso(ipk,jj,l1)<0.0d0) then
1605 31278 : vso(:,jj,l1)=-vso(:,jj,l1)
1606 : end if
1607 48 : apk=0.0d0
1608 78864 : do ii=1,mmax
1609 78864 : if(abs(vsr(ii,jj,l1))>apk) then
1610 2910 : apk=abs(vsr(ii,jj,l1))
1611 2910 : ipk=ii
1612 : end if
1613 : end do
1614 60 : if(vsr(ipk,jj,l1)<0.0d0) then
1615 41136 : vsr(:,jj,l1)=-vsr(:,jj,l1)
1616 : end if
1617 : end do
1618 :
1619 : end do ! l1
1620 :
1621 6 : ABI_FREE(vkbt)
1622 6 : ABI_FREE(vkbst)
1623 6 : return
1624 : end subroutine sr_so_r
1625 : !!***
1626 :
1627 : ! Copyright (c) 1989-2019 by D. R. Hamann, Mat-Sim Research LLC and Rutgers
1628 : ! University
1629 : !
1630 : !
1631 : ! This program is free software: you can redistribute it and/or modify
1632 : ! it under the terms of the GNU General Public License as published by
1633 : ! the Free Software Foundation, either version 3 of the License, or
1634 : ! (at your option) any later version.
1635 : !
1636 : ! This program is distributed in the hope that it will be useful,
1637 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
1638 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1639 : ! GNU General Public License for more details.
1640 : !
1641 : ! You should have received a copy of the GNU General Public License
1642 : ! along with this program. If not, see <http://www.gnu.org/licenses/>.
1643 : !
1644 : ! subroutine vpinteg(gg,hh,nn,mm,ss,rr)
1645 : !
1646 : !! integrals that go into construction of Vanderbilt separable pseudopotential
1647 : !
1648 : !! product of functions gg*hh goes like rr**mm at rr -> 0
1649 : !! integral on usual log mesh from 1 to nn
1650 : !
1651 : !!Input variables
1652 : ! integer,intent(in) :: nn,mm
1653 : ! real(dp),intent(in) :: gg(nn),hh(nn),rr(nn)
1654 : !
1655 : !!Output variable
1656 : ! real(dp),intent(out) :: ss
1657 : !
1658 : !!Local variables
1659 : ! real(dp) :: r0,amesh,al
1660 : ! integer :: ii
1661 : !
1662 : ! al = 0.01d0 * dlog(rr(101)/rr(1))
1663 : ! amesh = exp(al)
1664 : !
1665 : ! r0=rr(1)/dsqrt(amesh)
1666 : ! ss=r0**(mm+1)*(gg(1)*hh(1)/rr(1)**mm)/dfloat(mm+1)
1667 : !
1668 : ! do ii = 4, nn - 3
1669 : ! ss = ss + al*gg(ii)*hh(ii)*rr(ii)
1670 : ! end do
1671 : !
1672 : ! ss=ss + al*(23.d0*rr(nn-2)*gg(nn-2)*hh(nn-2) &
1673 : !& + 28.d0*rr(nn-1)*gg(nn-1)*hh(nn-1) &
1674 : !& + 9.d0*rr(nn )*gg(nn )*hh(nn ))/24.d0
1675 : !
1676 : !
1677 : ! return
1678 : ! end subroutine vpinteg
1679 :
1680 : end module m_pspheads
1681 : !!***
|
