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