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