1! 2! Copyright (C) 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 export_upf(filename, unit_loc) 10 !--------------------------------------------------------------------- 11 ! 12 use constants, only : fpi 13 use kinds, only : dp 14 use ld1inc, only : author, nlcc, zval, lpaw, write_coulomb, & 15 etots, rel, ecutwfc, ecutrho, iswitch, & 16 nwfts, nbeta, lmax, which_augfun, elts, octs, llts, & 17 nnts, rcutusts, rcutts, rcut, rcutus, els, ikk, nwfs, & 18 lls, nns, ocs, beta, bmat, qq, qvan, qvanl, rcloc, lloc, & 19 betas, grid, rhos, phits, psipaw, vpsloc, phis, & 20 rmatch_augfun, etot, etots, jjs, pawsetup, nn, & 21 core_state, ll, el, nwf, psi, vpot, nconf, zed, & 22 jjts, vpstot, lltsc, rcuttsc, rcutustsc, eltsc, & 23 lsave_wfc, wfc_ae_recon, wfc_ps_recon, tm, enlts, & 24 nstoaets, pseudotype, enls, rhoc, vnl, vpsloc, & 25 lgipaw_reconstruction, use_paw_as_gipaw, use_xsd 26 use funct, only: get_dft_name 27 use global_version, only: version_number 28 ! 29 use pseudo_types, only : pseudo_upf, pseudo_config, & 30 deallocate_pseudo_upf, deallocate_pseudo_config 31#if defined (__use_fox) 32 use write_upf_module, only: write_upf 33#else 34 use write_upf_new, only: write_upf 35#endif 36 ! 37 implicit none 38 ! 39 CHARACTER(len=*),INTENT(IN) :: filename 40 INTEGER,INTENT(IN):: unit_loc 41 ! 42 integer :: ibeta, jbeta, kbeta, l, ind, l1, l2 43 ! 44 ! Local variables 45 ! 46 integer :: nb, mesh 47 TYPE (pseudo_upf) :: upf 48 TYPE (pseudo_config) :: at_conf 49 CHARACTER(len=2), external :: atom_name 50 CHARACTER(len=9) :: day, hour 51 52 call date_and_tim(day,hour) 53 ! 54 IF (iswitch < 4 ) THEN 55 upf%generated='Generated using "atomic" code by A. Dal Corso & 56 & v.' // TRIM (version_number) 57 58 ELSE IF (iswitch==4) THEN 59 upf%generated='Generated using LDA-1/2 implemented by Leonardo& 60 & Matheus Marion Jorge' 61 ENDIF 62 upf%author=trim(author) 63 upf%date=trim(day) 64 upf%nv = "2.0.1" ! format version 65 ! 66 upf%zp = zval 67 upf%nlcc = nlcc 68 upf%dft = get_dft_name() 69 upf%psd = atom_name(nint(zed)) 70 71 if( pseudotype == 3) then 72 upf%tvanp = .true. 73 upf%typ='USPP' 74 else 75 upf%tvanp = .false. 76 upf%typ='NC' 77 endif 78 if(lpaw) upf%typ='PAW' 79 if(write_coulomb) upf%typ='1/r' 80 81 upf%tpawp = lpaw 82 upf%tcoulombp = write_coulomb 83 upf%has_gipaw = lgipaw_reconstruction 84 upf%paw_as_gipaw = use_paw_as_gipaw 85 upf%etotps = etots 86 upf%has_so = (rel == 2) 87 IF (rel == 2) THEN 88 upf%rel='full' 89 ELSE IF (rel == 1) THEN 90 upf%rel='scalar' 91 ELSE IF (rel < 1) THEN 92 upf%rel='no' 93 ELSE 94 call errore('export_upf', 'Unknown relativistic',1) 95 ENDIF 96 ! 97 upf%ecutwfc = ecutwfc 98 upf%ecutrho = max(ecutrho, ecutwfc*4._dp) 99 ! 100 upf%nwfc = nwfts 101 upf%nbeta = nbeta 102 ! 103 if (.not. lpaw) then 104 upf%lmax = lmax 105 upf%q_with_l = (which_augfun == 'PSQ') 106 else 107 upf%lmax = pawsetup%lmax 108 upf%q_with_l = .true. 109 endif 110 upf%lmax_rho = 2*upf%lmax 111 upf%nqlc = 2* upf%lmax+1 112 113 upf%mesh = grid%mesh 114 upf%dx = grid%dx 115 upf%xmin = grid%xmin 116 upf%zmesh = grid%zmesh 117 upf%rmax = grid%rmax 118 ! 119 upf%r = grid%r 120 upf%rab = grid%rab 121 ! 122 ! when possible, write semilocal PP's in the UPF file - may be 123 ! useful if one wants to use PPs in the UPF format in other codes 124 ! 125 if( pseudotype == 1 ) then 126 if ( rel == 2 ) then 127 allocate(upf%vnl(1:grid%mesh, 0:upf%lmax,2)) 128 else 129 allocate(upf%vnl(1:grid%mesh, 0:upf%lmax,1)) 130 end if 131 do nb=1, nbeta 132 l=lls(nb) 133 if ( rel < 2 .or. l == 0 .or. & 134 abs(jjs(nb)-l+0.5_dp) < 0.001_dp) then 135 ind = 1 136 else if ( rel == 2 .and. l > 0 .and. & 137 abs(jjs(nb)-l-0.5_dp) < 0.001_dp) then 138 ind = 2 139 endif 140 upf%vnl(1:grid%mesh,l,ind) = vnl(1:grid%mesh,l,ind) + & 141 vpsloc(1:grid%mesh) 142 end do 143 end if 144 ! 145 allocate(upf%lll(nbeta)) 146 upf%lll(1:nbeta) = lls(1:nbeta) 147 ! 148 ! *initial* wavefunctions indexes and parameters 149 allocate(upf%els(upf%nwfc), upf%oc(upf%nwfc), & 150 upf%nchi(upf%nwfc), upf%lchi(upf%nwfc), & 151 upf%epseu(upf%nwfc), upf%rcut_chi(upf%nwfc), & 152 upf%rcutus_chi(upf%nwfc) ) 153 upf%els(1:upf%nwfc) = elts(1:upf%nwfc) 154 upf%oc(1:upf%nwfc) = octs(1:upf%nwfc) 155 upf%lchi(1:upf%nwfc) = llts(1:upf%nwfc) 156 upf%nchi(1:upf%nwfc) = nnts(1:upf%nwfc) 157 upf%epseu(1:upf%nwfc) = enlts(1:upf%nwfc) 158 upf%rcut_chi(1:upf%nwfc) = rcutts(1:upf%nwfc) 159 upf%rcutus_chi(1:upf%nwfc) = rcutusts(1:upf%nwfc) 160 ! 161 ! projectors indexes and parameters 162 ! 163 allocate(upf%kbeta(nbeta), upf%els_beta(nbeta),& 164 upf%rcut(nbeta), upf%rcutus(nbeta)) 165 do nb=1,nbeta 166 upf%kbeta(nb) = ikk(nb) 167 upf%els_beta(nb)= els(nb) 168 upf%rcut(nb) = rcut(nb) 169 upf%rcutus(nb) = rcutus(nb) 170 end do 171 upf%kkbeta = maxval(upf%kbeta(1:nbeta)) 172 ! 173 ! Save GENERATION configuration: not needed to use the pseudopotential, 174 ! but must be saved for reference and for re-generating the pseudo 175 ! 176 at_conf%nwfs = nwfs 177 if (tm) then 178 at_conf%pseud = 'troullier-martins' 179 else 180 at_conf%pseud = 'rrkj' 181 endif 182 183 allocate(at_conf%els (nwfs),& 184 at_conf%nns (nwfs),& 185 at_conf%lls (nwfs),& 186 at_conf%ocs (nwfs),& 187 at_conf%rcut (nwfs),& 188 at_conf%rcutus(nwfs),& 189 at_conf%enls (nwfs)) 190 at_conf%els (1:nwfs) = els (1:nwfs) ! label (char*2) 191 at_conf%nns (1:nwfs) = nns (1:nwfs) ! n 192 at_conf%lls (1:nwfs) = lls (1:nwfs) ! l 193 at_conf%ocs (1:nwfs) = ocs (1:nwfs) ! occupation 194 at_conf%rcut (1:nwfs) = rcut (1:nwfs) ! inner cutoff radius 195 at_conf%rcutus(1:nwfs) = rcutus(1:nwfs) ! outer cutoff radius 196 at_conf%enls (1:nwfs) = enls (1:nwfs) ! one-particle energy 197 198 199 ! projectors 200 allocate(upf%beta(grid%mesh, upf%nbeta)) 201 upf%beta(1:grid%mesh, 1:upf%nbeta) = betas(1:grid%mesh, 1:nbeta) 202 ! 203 ! hamiltonian terms 204 allocate(upf%dion(upf%nbeta, upf%nbeta)) 205 upf%dion(1:upf%nbeta, 1:upf%nbeta) = bmat(1:nbeta, 1:nbeta) 206 ! 207 if (pseudotype.eq.3) then 208 allocate(upf%qqq(upf%nbeta, upf%nbeta)) 209 upf%qqq(1:upf%nbeta,1:upf%nbeta) = qq(1:nbeta,1:nbeta) 210 ! 211 upf%qqq_eps = 1.e-12_dp ! (hardcoded) 212 upf%nqf = 0 ! polinomial expansion of aug.charge is not supported by atomic 213 ! 214 if (upf%q_with_l .or. lpaw) then 215 allocate(upf%qfuncl(upf%mesh, upf%nbeta*(upf%nbeta+1)/2, 0:2*upf%lmax)) 216 else 217 allocate(upf%qfunc(upf%mesh, upf%nbeta*(upf%nbeta+1)/2)) 218 endif 219 ! 220 if(lpaw) qvanl(1:grid%mesh,:,:,:) = pawsetup%augfun(1:grid%mesh,:,:,:) 221 do ibeta=1,nbeta 222 do jbeta=ibeta,nbeta 223 kbeta = jbeta * (jbeta-1) / 2 + ibeta 224 if (upf%q_with_l .or. lpaw) then 225 l1=upf%lll(ibeta) 226 l2=upf%lll(jbeta) 227 do l=abs(l1-l2), l1+l2 228 upf%qfuncl(1:grid%mesh,kbeta,l) = qvanl(1:grid%mesh,ibeta,jbeta,l) 229 enddo 230 else 231 upf%qfunc(1:grid%mesh,kbeta) = qvan (1:grid%mesh, ibeta, jbeta) 232 endif 233 enddo 234 enddo 235 ! 236 endif 237 ! 238 allocate(upf%rho_atc(upf%mesh)) 239 if (upf%nlcc) then 240 upf%rho_atc(1:grid%mesh) = rhoc(1:grid%mesh)/fpi/grid%r2(1:grid%mesh) 241 else 242 upf%rho_atc(:) = 0.0_dp 243 end if 244 245 allocate(upf%rho_at(upf%mesh)) 246 upf%rho_at (1:grid%mesh) = rhos (1:grid%mesh,1) 247 ! 248 allocate(upf%chi(upf%mesh,upf%nwfc)) 249 upf%chi(1:grid%mesh,1:upf%nwfc) = phits(1:grid%mesh,1:upf%nwfc) 250 ! 251 allocate(upf%vloc(upf%mesh)) 252 upf%vloc(1:grid%mesh) = vpsloc(1:grid%mesh) 253 upf%lloc = lloc 254 upf%rcloc = rcloc 255 ! 256 ! 257 if (upf%has_so) CALL export_upf_so() 258 if (upf%tpawp) CALL export_upf_paw() 259 if (upf%has_gipaw) CALL export_upf_gipaw() 260 upf%has_wfc = lsave_wfc 261 if (upf%has_wfc) CALL export_upf_wfc() 262 ! 263 if (use_xsd) then 264 CALL write_upf( FILENAME = TRIM(filename) , UPF=upf, SCHEMA = 'qe_pp', CONF = at_conf, U_INPUT = unit_loc) 265 else 266 CALL write_upf( FILENAME = TRIM(filename), UPF= upf, SCHEMA = 'v2', CONF = at_conf, U_INPUT = unit_loc) 267 endif 268 ! 269 CALL deallocate_pseudo_upf( upf ) 270 CALL deallocate_pseudo_config( at_conf ) 271 272 RETURN 273 274 CONTAINS 275 276 SUBROUTINE export_upf_wfc 277 ALLOCATE( upf%aewfc(upf%mesh, upf%nbeta), upf%pswfc(upf%mesh, upf%nbeta) ) 278 upf%aewfc(1:upf%mesh,1:upf%nbeta) = psipaw(1:upf%mesh,1:upf%nbeta) 279 upf%pswfc(1:upf%mesh,1:upf%nbeta) = phis(1:upf%mesh,1:upf%nbeta) 280 END SUBROUTINE export_upf_wfc 281 282 SUBROUTINE export_upf_so 283 ALLOCATE( upf%nn(upf%nwfc), upf%jchi(upf%nwfc), upf%jjj(upf%nbeta) ) 284 285 upf%els(1:upf%nwfc) = elts(1:upf%nwfc) 286 upf%nn(1:upf%nwfc) = nnts(1:upf%nwfc) 287 upf%lchi(1:upf%nwfc) = llts(1:upf%nwfc) 288 upf%jchi(1:upf%nwfc) = jjts(1:upf%nwfc) 289 ! 290 upf%lll(1:upf%nbeta) = lls(1:upf%nbeta) 291 upf%jjj(1:upf%nbeta) = jjs(1:upf%nbeta) 292 293 END SUBROUTINE export_upf_so 294 ! 295 SUBROUTINE export_upf_paw 296 INTEGER :: co,n !EMINE 297 upf%paw_data_format = 2 298 ! 299 upf%paw%core_energy = etot -etots 300 upf%paw%lmax_aug = 2*upf%lmax 301 upf%paw%augshape = which_augfun 302 upf%paw%raug = rmatch_augfun 303 upf%paw%iraug = pawsetup%irc 304 305 allocate(upf%paw%ae_rho_atc(upf%mesh)) 306 upf%paw%ae_rho_atc(1:upf%mesh) = pawsetup%aeccharge(1:upf%mesh)/fpi/grid%r2(1:grid%mesh) 307 ! 308 allocate(upf%paw%ae_vloc(upf%mesh)) 309 upf%paw%ae_vloc(1:upf%mesh) = pawsetup%aeloc(1:upf%mesh) 310 ! 311 allocate(upf%paw%oc(upf%nbeta)) 312 do nb = 1,upf%nbeta 313 upf%paw%oc(nb) = max(pawsetup%oc(nb),0._dp) 314 enddo 315 ! 316 allocate(upf%paw%augmom(upf%nbeta, upf%nbeta, 0:2*upf%lmax)) 317 upf%paw%augmom(1:upf%nbeta,1:upf%nbeta,0:2*upf%lmax) & 318 = pawsetup%augmom(1:upf%nbeta,1:upf%nbeta,0:2*upf%lmax) 319 ! 320 upf%kkbeta = max(upf%kkbeta, upf%paw%iraug) 321 322 IF (upf%has_so) THEN 323 ALLOCATE( upf%paw%aewfc_rel(upf%mesh, upf%nbeta) ) 324 upf%paw%aewfc_rel(1:upf%mesh,1:upf%nbeta) = & 325 pawsetup%aewfc_rel(1:upf%mesh,1:upf%nbeta) 326 ENDIF 327 ! 328 !upf%paw%pfunc(:) = not used when writing, reconstructed from upf%aewfc 329 !upf%paw%ptfunc(:) = not used when writing, reconstructed from upf%pswfc 330 !=============================================================== 331 !For PAW pseudopotentials, now we also include core information: 332 !even when lgipaw_reconstruction = .false. 333 !EMINE 334 upf%gipaw_ncore_orbitals = COUNT(core_state(1:nwf)) 335 co = upf%gipaw_ncore_orbitals 336 ALLOCATE ( & 337 upf%gipaw_core_orbital_n(co), & 338 upf%gipaw_core_orbital_l(co), & 339 upf%gipaw_core_orbital_el(co), & 340 upf%gipaw_core_orbital(upf%mesh,co)) 341 upf%gipaw_core_orbital_n(1:co) = nn(1:co) 342 upf%gipaw_core_orbital_l(1:co) = ll(1:co) 343 upf%gipaw_core_orbital_el(1:co) = el(1:co) 344 DO n = 1,co 345 upf%gipaw_core_orbital(1:upf%mesh,n) & 346 = psi(1:upf%mesh,1,n) 347 ENDDO 348 !================================================================ 349 RETURN 350 END SUBROUTINE export_upf_paw 351 SUBROUTINE export_upf_gipaw 352 INTEGER :: co,nw,n,nb 353 354 IF ( nconf /= 1 ) CALL errore ( "write_gipaw_orbitals", & 355 "The (GI)PAW reconstruction requires one test configuration", abs(nconf) ) 356 357 upf%gipaw_data_format = 2 ! The version of the format 358 359 upf%gipaw_ncore_orbitals = COUNT(core_state(1:nwf)) 360 co = upf%gipaw_ncore_orbitals 361 upf%gipaw_wfs_nchannels = nwfts 362 nw = upf%gipaw_wfs_nchannels 363 364 ALLOCATE ( & 365 upf%gipaw_core_orbital_n(co), & 366 upf%gipaw_core_orbital_l(co), & 367 upf%gipaw_core_orbital_el(co), & 368 upf%gipaw_core_orbital(upf%mesh,co), & 369 upf%gipaw_wfs_el(nw), & 370 upf%gipaw_wfs_ll(nw), & 371 upf%gipaw_wfs_rcut(nw), & 372 upf%gipaw_wfs_rcutus(nw), & 373 upf%gipaw_wfs_ae(upf%mesh,nw), & 374 upf%gipaw_wfs_ps(upf%mesh,nw), & 375 upf%gipaw_vlocal_ae(upf%mesh), & 376 upf%gipaw_vlocal_ps(upf%mesh) & 377 ) 378 379 upf%gipaw_core_orbital_n(1:co) = nn(1:co) 380 upf%gipaw_core_orbital_l(1:co) = ll(1:co) 381 upf%gipaw_core_orbital_el(1:co) = el(1:co) 382 383 DO n = 1,co 384 upf%gipaw_core_orbital(1:upf%mesh,n) & 385 = psi(1:upf%mesh,1,n) 386 ENDDO 387 388 upf%gipaw_vlocal_ae(1:upf%mesh) & 389 = grid%r(1:upf%mesh) * vpot(1:upf%mesh,1) 390 upf%gipaw_vlocal_ps(1:upf%mesh) & 391 = grid%r(1:upf%mesh) * vpstot(1:upf%mesh,1) 392 393 upf%gipaw_wfs_el(1:nw) = elts(1:nw) 394 upf%gipaw_wfs_ll(1:nw) = lltsc(1:nw,1) 395 396 ! Find the suitable core radii to be written out 397 !*apsi WARNING: DOES NOT WORK WITH VANDERBILT PP YET 398 DO nb = 1,nw 399 upf%gipaw_wfs_rcut(nb) = -1.0_dp 400 DO n = 1, nwfs 401 IF ( els(n) == eltsc(nb,1) ) THEN 402 upf%gipaw_wfs_rcut(nb) = rcuttsc(nb,1) 403 upf%gipaw_wfs_rcutus(nb) = rcutustsc(nb,1) 404 END IF 405 END DO 406 ! 407 IF ( upf%gipaw_wfs_rcut(nb) < 0.0_dp ) THEN 408 DO n = 1, nwfs 409 ! If there is one with the same l... 410 IF ( lltsc(nb,1) == lls(n) ) THEN 411 upf%gipaw_wfs_rcut(nb) = rcuttsc(nb,1) 412 upf%gipaw_wfs_rcutus(nb) = rcutustsc(nb,1) 413 END IF 414 END DO 415 END IF 416 ENDDO 417 418 DO n = 1,nw 419 ! 420 upf%gipaw_wfs_ae(1:upf%mesh,n) = wfc_ae_recon(1:upf%mesh,nstoaets(n)) 421 ! 422 upf%gipaw_wfs_ps(1:upf%mesh,n) = wfc_ps_recon(1:upf%mesh,n) 423 ENDDO 424 425 426 RETURN 427 END SUBROUTINE export_upf_gipaw 428 ! 429END SUBROUTINE export_upf 430