1!
2! Copyright (C) 2001-2020 Quantum ESPRESSO 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!
9!----------------------------------------------------------------------------
10MODULE io_rho_xml
11  !----------------------------------------------------------------------------
12  !
13  USE kinds,       ONLY : DP
14  USE io_files,    ONLY : create_directory
15  USE io_base,     ONLY : write_rhog, read_rhog
16  !
17  PRIVATE
18  PUBLIC :: write_scf, read_scf
19  !
20  ! {read|write}_rho: read or write the charge density
21  ! {read|write}_scf: as above, plus ldaU ns, PAW becsum, meta-GGA
22  !
23  CONTAINS
24
25    SUBROUTINE write_scf ( rho, nspin )
26      !
27      USE paw_variables,    ONLY : okpaw
28      USE ldaU,             ONLY : lda_plus_u, hub_back, lda_plus_u_kind, nsg
29      USE funct,            ONLY : dft_is_meta
30      USE noncollin_module, ONLY : noncolin
31      USE spin_orb,         ONLY : domag
32      USE scf,              ONLY : scf_type
33      !
34      USE cell_base,        ONLY : bg, tpiba
35      USE gvect,            ONLY : ig_l2g, mill
36      USE control_flags,    ONLY : gamma_only
37      USE io_files,         ONLY : restart_dir
38      USE io_global,        ONLY : ionode, ionode_id, stdout
39      USE mp_pools,         ONLY : my_pool_id
40      USE mp_bands,         ONLY : my_bgrp_id, root_bgrp_id, &
41                                   root_bgrp, intra_bgrp_comm
42      USE mp_images,        ONLY : intra_image_comm
43      USE mp,               ONLY : mp_bcast
44
45      !
46      IMPLICIT NONE
47      TYPE(scf_type),   INTENT(IN)           :: rho
48      INTEGER,          INTENT(IN)           :: nspin
49      !
50      CHARACTER (LEN=256) :: dirname
51      INTEGER :: nspin_, iunocc, iunpaw, ierr
52      INTEGER, EXTERNAL :: find_free_unit
53
54      dirname = restart_dir ( )
55      CALL create_directory( dirname )
56      ! in the following case do not read or write polarization
57      IF ( noncolin .AND. .NOT.domag ) THEN
58         nspin_ = 1
59      ELSE
60         nspin_ = nspin
61      ENDIF
62      ! Write G-space density
63      IF ( my_pool_id == 0 .AND. my_bgrp_id == root_bgrp_id ) &
64           CALL write_rhog( TRIM(dirname) // "charge-density", &
65           root_bgrp, intra_bgrp_comm, &
66           bg(:,1)*tpiba, bg(:,2)*tpiba, bg(:,3)*tpiba, &
67           gamma_only, mill, ig_l2g, rho%of_g(:,1:nspin_) )
68      !
69      ! Write kinetic energy density density
70      IF ( dft_is_meta() ) THEN
71         IF ( my_pool_id == 0 .AND. my_bgrp_id == root_bgrp_id ) &
72              CALL write_rhog( TRIM(dirname) // "ekin-density", &
73              root_bgrp, intra_bgrp_comm, &
74              bg(:,1)*tpiba, bg(:,2)*tpiba, bg(:,3)*tpiba, &
75              gamma_only, mill, ig_l2g, rho%kin_g(:,1:nspin_) )
76         WRITE(stdout,'(5x,"Writing meta-gga kinetic term")')
77      ENDIF
78
79      ! Then write the other terms to separate files
80
81      IF ( lda_plus_u ) THEN
82         !
83         iunocc = find_free_unit ()
84         IF ( ionode ) THEN
85            OPEN ( UNIT=iunocc, FILE = TRIM(dirname) // 'occup.txt', &
86                 FORM='formatted', STATUS='unknown' )
87            IF (lda_plus_u_kind.EQ.0) THEN
88               WRITE( iunocc, * , iostat = ierr) rho%ns
89               IF (hub_back) WRITE( iunocc, * , iostat = ierr) rho%nsb
90            ELSEIF (lda_plus_u_kind.EQ.1) THEN
91               IF (noncolin) THEN
92                  WRITE( iunocc, * , iostat = ierr) rho%ns_nc
93               ELSE
94                  WRITE( iunocc, * , iostat = ierr) rho%ns
95               ENDIF
96            ELSEIF (lda_plus_u_kind.EQ.2) THEN
97               WRITE( iunocc, * , iostat = ierr) nsg
98               ! Write Hubbard_V to file
99               CALL write_V
100            ENDIF
101         ENDIF
102         CALL mp_bcast( ierr, ionode_id, intra_image_comm )
103         IF ( ierr/=0 ) CALL errore('write_scf', 'Writing ldaU ns', 1)
104         IF ( ionode ) THEN
105            CLOSE( UNIT = iunocc, STATUS = 'KEEP' )
106         ENDIF
107         !
108      END IF
109      !
110      IF ( okpaw ) THEN
111         !
112         iunpaw = find_free_unit ()
113         IF ( ionode ) THEN
114            OPEN ( UNIT=iunpaw, FILE = TRIM(dirname) // 'paw.txt', &
115                 FORM='formatted', STATUS='unknown' )
116            WRITE( iunpaw, * , iostat = ierr) rho%bec
117         END IF
118         CALL mp_bcast( ierr, ionode_id, intra_image_comm )
119         IF ( ierr/=0 ) CALL errore('write_scf', 'Writing PAW becsum',1)
120         IF ( ionode ) THEN
121            CLOSE( UNIT = iunpaw, STATUS = 'KEEP' )
122         ENDIF
123         !
124      END IF
125
126      RETURN
127    END SUBROUTINE write_scf
128
129    SUBROUTINE read_scf ( rho, nspin, gamma_only )
130      !
131      USE scf,              ONLY : scf_type
132      USE paw_variables,    ONLY : okpaw
133      USE ldaU,             ONLY : lda_plus_u, starting_ns, hub_back, &
134                                   lda_plus_u_kind, nsg
135      USE noncollin_module, ONLY : noncolin
136      USE spin_orb,         ONLY : domag
137      USE gvect,            ONLY : ig_l2g
138      USE funct,            ONLY : dft_is_meta
139      USE io_files,         ONLY : restart_dir
140      USE io_global,        ONLY : ionode, ionode_id, stdout
141      USE mp_bands,         ONLY : root_bgrp, intra_bgrp_comm
142      USE mp_images,        ONLY : intra_image_comm
143      USE mp,               ONLY : mp_bcast, mp_sum
144      !
145      IMPLICIT NONE
146      TYPE(scf_type),   INTENT(INOUT)        :: rho
147      INTEGER,          INTENT(IN)           :: nspin
148      LOGICAL, OPTIONAL,INTENT(IN)           :: gamma_only
149      !
150      CHARACTER(LEN=256) :: dirname
151      LOGICAL :: lexist
152      INTEGER :: nspin_, iunocc, iunpaw, ierr
153      INTEGER, EXTERNAL :: find_free_unit
154
155      dirname = restart_dir ( )
156      ! in the following case do not read or write polarization
157      IF ( noncolin .AND. .NOT.domag ) THEN
158         nspin_=1
159      ELSE
160         nspin_=nspin
161      ENDIF
162      ! read charge density
163      CALL read_rhog( TRIM(dirname) // "charge-density", &
164           root_bgrp, intra_bgrp_comm, &
165           ig_l2g, nspin_, rho%of_g, gamma_only )
166      IF ( nspin > nspin_) rho%of_g(:,nspin_+1:nspin) = (0.0_dp, 0.0_dp)
167      !
168      ! read kinetic energy density
169      IF ( dft_is_meta() ) THEN
170         CALL read_rhog( TRIM(dirname) // "ekin-density", &
171              root_bgrp, intra_bgrp_comm, &
172              ig_l2g, nspin_, rho%kin_g, gamma_only, ierr )
173         IF ( ierr == 0 ) THEN
174            WRITE(stdout,'(5x,"Reading meta-gga kinetic term")')
175         ELSE
176            rho%kin_g(:,:) = (0.0_dp, 0.0_dp)
177            WRITE(stdout,'(5x,"BEWARE: kinetic-energy density file not found,",&
178                    & " Kinetic-energy density set to 0")')
179         ENDIF
180      END IF
181
182      IF ( lda_plus_u ) THEN
183         !
184         ! The occupations ns also need to be read in order to build up
185         ! the potential
186         !
187         iunocc = find_free_unit ()
188         IF ( ionode ) THEN
189            OPEN ( UNIT=iunocc, FILE = TRIM(dirname) // 'occup.txt', &
190                 FORM='formatted', STATUS='old', IOSTAT=ierr )
191            IF (lda_plus_u_kind.EQ.0) THEN
192               READ( UNIT = iunocc, FMT = *, iostat = ierr ) rho%ns
193               IF (hub_back) READ( UNIT = iunocc, FMT = * , iostat = ierr) rho%nsb
194            ELSEIF (lda_plus_u_kind.EQ.1) THEN
195               IF (noncolin) THEN
196                  READ( UNIT = iunocc, FMT = *, iostat = ierr ) rho%ns_nc
197               ELSE
198                  READ( UNIT = iunocc, FMT = *, iostat = ierr ) rho%ns
199               ENDIF
200            ELSEIF (lda_plus_u_kind.EQ.2) THEN
201               READ( UNIT = iunocc, FMT = * , iostat = ierr) nsg
202            ENDIF
203         ENDIF
204         !
205         CALL mp_bcast( ierr, ionode_id, intra_image_comm )
206         IF ( ierr/=0 ) CALL errore('read_scf', 'Reading ldaU ns', 1)
207         !
208         IF ( ionode ) THEN
209            CLOSE( UNIT = iunocc, STATUS = 'KEEP')
210         ELSE
211            IF (lda_plus_u_kind.EQ.0) THEN
212               rho%ns(:,:,:,:) = 0.D0
213               IF (hub_back) rho%nsb(:,:,:,:) = 0.D0
214            ELSEIF (lda_plus_u_kind.EQ.1) THEN
215               IF (noncolin) THEN
216                  rho%ns_nc(:,:,:,:) = 0.D0
217               ELSE
218                  rho%ns(:,:,:,:) = 0.D0
219               ENDIF
220            ELSEIF (lda_plus_u_kind.EQ.2) THEN
221               nsg(:,:,:,:,:) = (0.d0, 0.d0)
222            ENDIF
223         ENDIF
224         !
225         IF (lda_plus_u_kind.EQ.0) THEN
226            CALL mp_sum(rho%ns, intra_image_comm)
227            IF (hub_back) CALL mp_sum(rho%nsb, intra_image_comm)
228         ELSEIF (lda_plus_u_kind.EQ.1) THEN
229            IF (noncolin) THEN
230               CALL mp_sum(rho%ns_nc, intra_image_comm)
231            ELSE
232               CALL mp_sum(rho%ns, intra_image_comm)
233            ENDIF
234         ELSEIF (lda_plus_u_kind.EQ.2) THEN
235            CALL mp_sum(nsg, intra_image_comm)
236         ENDIF
237         !
238         ! If projections on Hubbard manifold are read from file, there is no
239         ! need to set starting values: reset them to prevent any problem
240         starting_ns = -1.0_dp
241         !
242      ENDIF
243      !
244      IF ( okpaw ) THEN
245         !
246         ! Also the PAW coefficients are needed:
247         !
248         iunpaw = find_free_unit ()
249         IF ( ionode ) THEN
250            OPEN ( UNIT=iunpaw, FILE = TRIM(dirname) // 'paw.txt', &
251                 FORM='formatted', STATUS='old', IOSTAT=ierr )
252            READ( UNIT = iunpaw, FMT = *, iostat=ierr ) rho%bec
253         END IF
254         CALL mp_bcast( ierr, ionode_id, intra_image_comm )
255         IF ( ierr/=0 ) CALL errore('read_scf', 'Reading PAW becsum',1)
256         IF ( ionode ) THEN
257            CLOSE( UNIT = iunpaw, STATUS = 'KEEP')
258         ELSE
259            rho%bec(:,:,:) = 0.D0
260         END IF
261         CALL mp_sum(rho%bec, intra_image_comm)
262         !
263      END IF
264      !
265      RETURN
266    END SUBROUTINE read_scf
267    !
268END MODULE io_rho_xml
269