1!
2! Copyright (C) 2004-2008 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!---------------------------------------------------------------------
9subroutine import_upf ( )
10  !---------------------------------------------------------------------
11  !
12  !   read "is"-th pseudopotential in the Unified Pseudopotential Format
13  !   from unit "iunps" - convert and copy to internal PWscf variables
14  !   return error code in "ierr" (success: ierr=0)
15  !
16  ! PWSCF modules
17  !
18  use constants, only : fpi
19  use kinds, only : dp
20  use radial_grids, only : ndmx, radial_grid_type, allocate_radial_grid, &
21                           nullify_radial_grid, deallocate_radial_grid, &
22                           radial_grid_copy
23  use ld1inc, only : file_pseudo, zval, nlcc, pseudotype, etots, lmax, lsave_wfc,&
24                     zed, nbeta, betas, lls, jjs, ikk, els, rcut, rcutus, &
25                     lloc, vpsloc, grid, nwfs, bmat, qq, qvan, qvanl, rhoc, &
26                     rhos, phis, which_augfun, lpaw, rmatch_augfun, pawsetup, psipaw
27  use funct, only: set_dft_from_name
28  !
29  use pseudo_types, only : pseudo_upf, pseudo_config, deallocate_pseudo_upf
30  use paw_type
31#if defined (__use_fox)
32  use upf_module, only : read_upf_new
33#else
34  use read_upf_new_module, only : read_upf_new
35#endif
36  !
37  implicit none
38  !
39  integer :: iunps, ierr, ibeta, jbeta, kbeta, l, l1, l2
40  !
41  !     Local variables
42  !
43  integer :: nb, ios
44  TYPE (pseudo_upf) :: upf
45  !
46  ierr = 1
47  call read_upf_new ( file_pseudo, upf, ierr)
48  !
49  if (ierr>0) &
50     call errore('import_upf','reading pseudo upf',abs(ierr))
51  !
52  zval  = upf%zp
53  nlcc = upf%nlcc
54  call set_dft_from_name (upf%dft)
55
56  if (upf%typ.eq.'NC'.OR.upf%typ.eq.'SL') then
57     pseudotype=2
58  else
59     pseudotype=3
60  endif
61  lpaw = upf%tpawp
62
63  etots=upf%etotps
64  lmax = upf%lmax
65  grid%mesh = upf%mesh
66  call allocate_radial_grid(grid, grid%mesh)
67  grid%r  (1:grid%mesh) = upf%r  (1:upf%mesh)
68  grid%rab(1:grid%mesh) = upf%rab(1:upf%mesh)
69  grid%r2 (1:grid%mesh) = grid%r(1:grid%mesh)**2
70  grid%sqr(1:grid%mesh) = sqrt(grid%r(1:grid%mesh))
71  if (.not.upf%has_so) then
72     if (grid%r(1) > 0.0_dp) then
73        !
74        ! r(i+1) = exp(xmin)/zmesh * exp(i*dx)
75        !
76        grid%dx=log(grid%r(grid%mesh)/grid%r(1))/(grid%mesh-1)
77        grid%rmax=grid%r(grid%mesh)
78        grid%xmin=log(zed*grid%r(1))
79        grid%zmesh=zed
80     else
81        !
82        ! r(i+1) = exp(xmin)/zmesh * ( exp(i*dx) - 1 )
83        !
84        grid%dx=log( (grid%r(3)-grid%r(2)) / grid%r(2) )
85        grid%rmax=grid%r(grid%mesh)
86        grid%zmesh=zed
87        grid%xmin=log(zed*grid%r(2)**2/(grid%r(3)-2.0_dp*grid%r(2)))
88     end if
89  else
90     grid%dx=upf%dx
91     grid%xmin=upf%xmin
92     grid%zmesh=upf%zmesh
93     grid%rmax=exp(grid%xmin+(grid%mesh-1)*grid%dx)/grid%zmesh
94  endif
95  if (abs(exp(grid%xmin+(grid%mesh-1)*grid%dx)/zed-grid%rmax).gt.1.e-6_dp) &
96       call errore('read_pseudoup','mesh not supported',1)
97
98  nwfs = upf%nwfc
99
100  nbeta= upf%nbeta
101  lls(1:nbeta)=upf%lll(1:nbeta)
102
103  if (upf%has_so) then
104     jjs(1:nbeta)=upf%jjj(1:nbeta)
105  else
106     jjs=0.0_dp
107  endif
108  !
109  !
110  do nb=1,nbeta
111     ikk(nb)=upf%kbeta(nb)
112     els(nb)=upf%els_beta(nb)
113     rcut(nb)=upf%rcut(nb)
114     rcutus(nb)=upf%rcutus(nb)
115  end do
116  betas(1:grid%mesh, 1:nbeta) = upf%beta(1:upf%mesh, 1:upf%nbeta)
117  bmat(1:nbeta, 1:nbeta) = upf%dion(1:upf%nbeta, 1:upf%nbeta)
118  !
119  if (pseudotype.eq.3) then
120     qq(1:nbeta,1:nbeta) = upf%qqq(1:upf%nbeta,1:upf%nbeta)
121     do ibeta=1,nbeta
122        do jbeta=ibeta,nbeta
123           kbeta = jbeta * (jbeta-1) / 2 + ibeta
124           if (upf%q_with_l .or. lpaw) then
125              l1=upf%lll(ibeta)
126              l2=upf%lll(jbeta)
127              do l=abs(l1-l2), l1+l2
128                 qvanl(1:upf%mesh,ibeta,jbeta,l)=upf%qfuncl(1:upf%mesh,kbeta,l)
129                 if (ibeta /= jbeta) qvanl (1:grid%mesh, jbeta, ibeta, l)= &
130                                    upf%qfuncl(1:upf%mesh,kbeta,l)
131              enddo
132              qvan(1:upf%mesh,ibeta,jbeta)=upf%qfuncl(1:upf%mesh,kbeta,0)
133              if (ibeta /= jbeta) qvan(1:grid%mesh, jbeta, ibeta)= &
134                                    upf%qfuncl(1:upf%mesh,kbeta,0)
135              which_augfun='PSQ'
136           else
137              qvan (1:grid%mesh, ibeta, jbeta) = upf%qfunc(1:upf%mesh,kbeta)
138              if (ibeta /= jbeta) qvan (1:grid%mesh, jbeta, ibeta)= &
139                                       upf%qfunc(1:upf%mesh,kbeta)
140              which_augfun='AE'
141           endif
142        enddo
143     enddo
144  else
145     qq=0.0_dp
146     qvan=0.0_dp
147  endif
148  !
149  if (upf%nlcc) then
150     rhoc(1:grid%mesh) = upf%rho_atc(1:upf%mesh)*fpi*grid%r2(1:upf%mesh)
151  else
152     rhoc(:) = 0.0_dp
153  end if
154  rhos=0.0_dp
155  rhos (1:grid%mesh,1) = upf%rho_at (1:upf%mesh)
156
157  !phis(1:grid%mesh,1:nwfs)=upf%chi(1:grid%mesh,1:nwfs)
158  if(upf%has_wfc) then
159      lsave_wfc = .true.
160      phis(1:grid%mesh,1:nbeta)=upf%pswfc(1:grid%mesh,1:nbeta)
161      psipaw(1:grid%mesh,1:nbeta)=upf%aewfc(1:grid%mesh,1:nbeta)
162  endif
163  !!! TEMP
164  lloc = -1
165  vpsloc(1:grid%mesh) = upf%vloc(1:upf%mesh)
166  !!!
167  ! paw:
168  if (lpaw) then
169    which_augfun = upf%paw%augshape
170    rmatch_augfun = upf%paw%raug
171    call allocate_pseudo_paw( pawsetup, grid%mesh, nbeta, lmax )
172    CALL nullify_radial_grid( pawsetup%grid )
173    call allocate_radial_grid(pawsetup%grid,grid%mesh)
174    call set_pawsetup( pawsetup, upf )
175    CALL radial_grid_copy(grid, pawsetup%grid)
176  endif
177
178  CALL deallocate_pseudo_upf( upf )
179
180
181end subroutine import_upf
182
183SUBROUTINE set_pawsetup(pawset_, upf_)
184USE kinds, ONLY : DP
185USE constants, ONLY : fpi
186USE paw_type, ONLY : paw_t
187USE pseudo_types, ONLY: pseudo_upf
188USE ld1_parameters,   ONLY: nwfsx
189USE atomic_paw,    ONLY : compute_nonlocal_coeff_ion
190IMPLICIT NONE
191TYPE(paw_t), INTENT(INOUT) :: pawset_
192TYPE(pseudo_upf), INTENT(IN) :: upf_
193REAL(DP), ALLOCATABLE :: ddd_(:,:)
194INTEGER :: mesh, nbeta,ih,jh,ijh
195
196   nbeta=upf_%nbeta
197   mesh=upf_%mesh
198   pawset_%augfun=0.0_DP
199   pawset_%augmom=0.0_DP
200   pawset_%enl(:) = 0.0_DP
201   if (upf_%has_so) then
202      pawset_%jj(1:nbeta) = upf_%jjj(1:nbeta)
203      pawset_%rel=2
204   else
205      pawset_%jj(:) = 0.0_DP
206      pawset_%rel=1
207   endif
208   pawset_%l(1:nbeta) = upf_%lll(1:nbeta)
209   pawset_%ikk(1:nbeta) = upf_%kbeta(1:nbeta)
210   pawset_%oc(1:nbeta) = upf_%paw%oc(1:nbeta)
211   pawset_%aewfc(1:mesh,1:nbeta) = upf_%aewfc(1:mesh,1:nbeta)
212   pawset_%pswfc(1:mesh,1:nbeta) = upf_%pswfc(1:mesh,1:nbeta)
213   IF (upf_%has_so) &
214   pawset_%aewfc_rel(1:mesh,1:nbeta) = upf_%paw%aewfc_rel(1:mesh,1:nbeta)
215   pawset_%proj(1:mesh,1:nbeta) = upf_%beta(1:mesh,1:nbeta)
216
217   DO ih = 1,nbeta
218   DO jh = ih,nbeta
219      ijh = jh * (jh-1) / 2 + ih
220      pawset_%augfun(1:mesh,ih,jh,0:upf_%paw%lmax_aug) = &
221                           upf_%qfuncl(1:mesh,ijh,0:upf_%paw%lmax_aug)
222      IF ( ih /= jh ) &
223      pawset_%augfun(1:mesh,jh,ih,0:upf_%paw%lmax_aug) = &
224                           upf_%qfuncl(1:mesh,ijh,0:upf_%paw%lmax_aug)
225   ENDDO
226   ENDDO
227
228   pawset_%augmom(1:nbeta,1:nbeta,0:upf_%paw%lmax_aug) = &
229                  upf_%paw%augmom(1:nbeta,1:nbeta,0:upf_%paw%lmax_aug)
230   pawset_%aeccharge(1:mesh) = upf_%paw%ae_rho_atc(1:mesh)*fpi*upf_%r(1:mesh)**2
231   pawset_%psccharge(1:mesh) = upf_%rho_atc(1:mesh)*fpi*upf_%r(1:mesh)**2
232   pawset_%pscharge(1:mesh) = upf_%rho_at(1:mesh)
233   pawset_%aeloc(1:mesh) = upf_%paw%ae_vloc(1:mesh)
234   pawset_%psloc(1:mesh) = upf_%vloc(1:mesh)
235   pawset_%kdiff(1:nbeta,1:nbeta) = 0.0_DP
236   pawset_%dion (1:nbeta,1:nbeta) = upf_%dion(1:nbeta,1:nbeta)
237   pawset_%symbol=upf_%psd
238   pawset_%zval=upf_%zp
239   pawset_%z=upf_%zmesh
240   pawset_%nlcc=upf_%nlcc
241   pawset_%nwfc=upf_%nbeta
242   pawset_%irc=upf_%kkbeta
243   pawset_%lmax=upf_%lmax
244   pawset_%rmatch_augfun=upf_%paw%raug
245!
246!  The kinetic energy must be recalculated
247!
248   ALLOCATE(ddd_(nwfsx,nwfsx))
249
250   CALL compute_nonlocal_coeff_ion(ddd_, pawset_)
251   pawset_%kdiff(1:nbeta,1:nbeta) = upf_%dion(1:nbeta,1:nbeta)- &
252                                    ddd_(1:nbeta,1:nbeta)
253
254   DEALLOCATE(ddd_)
255
256END SUBROUTINE set_pawsetup
257