Branch data Line data Source code
1 : : !!****m* ABINIT/m_ioarr
2 : : !! NAME
3 : : !! m_ioarr
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides routines to read/write arrays given on the FFT mesh (densities, potentials ...).
7 : : !! The code supports both Fortran files as well as netcdf files in a transparent way.
8 : : !! The appropriate IO layer is selected from files extensions: netcdf primitives are used if the
9 : : !! file ends with `.nc`. If all the other cases we read/write files in Fortran format.
10 : : !! MPI-IO primitives are used when the FFT arrays are MPI distributed.
11 : : !!
12 : : !! COPYRIGHT
13 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MVer, MT, MG)
14 : : !! This file is distributed under the terms of the
15 : : !! GNU General Public License, see ~abinit/COPYING
16 : : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : MODULE m_ioarr
27 : :
28 : : use defs_basis
29 : : use m_abicore
30 : : USE_MPI
31 : : use m_xmpi
32 : : use m_wffile
33 : : use m_errors
34 : : use m_nctk
35 : : use m_dtset
36 : : use m_crystal
37 : : use m_ebands
38 : : use m_hdr
39 : : use m_pawrhoij
40 : : use netcdf
41 : :
42 : : use defs_abitypes, only : mpi_type
43 : : use defs_wvltypes, only : wvl_denspot_type
44 : : use m_time, only : cwtime, cwtime_report, timab
45 : : use m_io_tools, only : iomode_from_fname, iomode2str, open_file, get_unit
46 : : use m_fstrings, only : sjoin, itoa, endswith, ltoa
47 : : use m_numeric_tools, only : interpolate_denpot
48 : : use m_geometry, only : metric
49 : : use m_mpinfo, only : destroy_mpi_enreg, ptabs_fourdp, initmpi_seq
50 : : use m_fourier_interpol,only : fourier_interpol
51 : :
52 : : implicit none
53 : :
54 : : #ifdef HAVE_MPI1
55 : : include 'mpif.h'
56 : : #endif
57 : :
58 : : private
59 : :
60 : : public :: ioarr ! Read or write rho(r) or v(r), either ground-state or response-functions.
61 : : public :: fftdatar_write ! Write an array in real space. IO library is automatically selected
62 : : ! from the file extension and the number of FFT processors:
63 : : public :: fftdatar_write_from_hdr ! Write an array in real-space to file plus crystal_t and ebands_t
64 : : public :: read_rhor ! Read rhor from DEN file.
65 : : public :: fort_denpot_skip ! Skip the header and the DEN/POT records (Fortran format)
66 : :
67 : : private :: denpot_spin_convert ! Convert a density/potential from a spin representation to another
68 : :
69 : : CONTAINS !====================================================================================================
70 : : !!***
71 : :
72 : : !----------------------------------------------------------------------
73 : :
74 : : !!****f* m_ioarr/ioarr
75 : : !!
76 : : !! NAME
77 : : !! ioarr
78 : : !!
79 : : !! FUNCTION
80 : : !! Read or write rho(r) or v(r), either ground-state or response-functions.
81 : : !! If ground-state, these arrays are real, if response-functions, these arrays are complex.
82 : : !! (in general, an array stored in unformatted form on a real space fft grid).
83 : : !! rdwr=1 to read, 2 to write
84 : : !!
85 : : !! This subroutine should be called only by one processor in the writing mode
86 : : !!
87 : : !! INPUTS
88 : : !! (some may be output)
89 : : !! accessfil=
90 : : !! 0 for FORTRAN_IO
91 : : !! 3 for ETSF_IO
92 : : !! 4 for MPI_IO
93 : : !! cplex=1 for real array, 2 for complex
94 : : !! nfft=Number of FFT points treated by this node.
95 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
96 : : !! dtset <type(dataset_type)>=all input variables for this dataset
97 : : !! fform=integer specification for data type:
98 : : !! 2 for wf; 52 for density; 102 for potential
99 : : !! old format (prior to ABINITv2.0): 1, 51 and 101.
100 : : !! fildata=file name
101 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
102 : : !! if rdwr=1 , used to compare with the hdr of the read disk file
103 : : !! if rdwr=2 , used as the header of the written disk file
104 : : !! mpi_enreg=information about MPI parallelization
105 : : !! rdwr=choice parameter, see above
106 : : !! rdwrpaw=1 only if rhoij PAW quantities have to be read (if rdwr=1)
107 : : !! [single_proc]=True if only ONE MPI process is calling this routine. This usually happens when
108 : : !! master calls ioarr to read data that is then broadcasted in the caller. Default: False.
109 : : !! Note that singleproc is not compatible with FFT parallelism because nfft is assumed to be
110 : : !! the total number of points in the FFT mesh.
111 : : !!
112 : : !! OUTPUT
113 : : !! (see side effects)
114 : : !!
115 : : !! SIDE EFFECTS
116 : : !! Input/Output
117 : : !! arr(cplex*nfft,nspden)=array on real space grid, returned for rdwr=1, input for rdwr=2
118 : : !! etotal=total energy (Ha), returned for rdwr=1
119 : : !! === if rdwrpaw/=0 ===
120 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
121 : : !!
122 : : !! SOURCE
123 : :
124 : 0 : subroutine ioarr(accessfil,arr,dtset,etotal,fform,fildata,hdr,mpi_enreg, &
125 : 0 : & ngfft,cplex,nfft,pawrhoij,rdwr,rdwrpaw,wvl_den,single_proc)
126 : :
127 : : !Arguments ------------------------------------
128 : : !scalars
129 : : integer,intent(in) :: accessfil,cplex,nfft,rdwr,rdwrpaw
130 : : integer,intent(inout) :: fform
131 : : real(dp),intent(inout) :: etotal
132 : : character(len=*),intent(in) :: fildata
133 : : logical,optional,intent(in) :: single_proc
134 : : type(MPI_type),intent(inout) :: mpi_enreg
135 : : type(dataset_type),intent(in) :: dtset
136 : : type(hdr_type),intent(inout) :: hdr
137 : : type(wvl_denspot_type),optional, intent(in) :: wvl_den
138 : : !arrays
139 : : integer,intent(in) :: ngfft(18)
140 : : real(dp),intent(inout),target :: arr(cplex*nfft,dtset%nspden)
141 : : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
142 : :
143 : : !Local variables-------------------------------
144 : : integer :: ncid,ncerr
145 : : character(len=fnlen) :: file_etsf
146 : : #ifdef HAVE_BIGDFT
147 : : integer :: i,i1,i2,i3,ia,ind,n1,n2,n3
148 : : integer :: zindex,zstart,zstop
149 : : #endif
150 : : !scalars
151 : : integer,parameter :: master=0
152 : : logical,parameter :: ALLOW_FFTINTERP=.True.
153 : : logical :: need_fftinterp,icheck_fft,qeq0
154 : : integer :: in_unt,out_unt,nfftot_in,nfftot_out,nspden,ncplxfft
155 : : integer :: iomode,fform_dum,iarr,ierr,ispden,me,me_fft,comm_fft
156 : : integer :: comm_cell,usewvl,unt
157 : : integer :: restart,restartpaw,spaceComm,spaceComm_io
158 : : real(dp) :: cputime,walltime,gflops
159 : : character(len=500) :: msg,errmsg
160 : : character(len=fnlen) :: my_fildata
161 : : character(len=nctk_slen) :: varname
162 : 0 : type(hdr_type),target :: hdr0
163 : : type(wffile_type) :: wff
164 : 0 : type(MPI_type) :: MPI_enreg_seq
165 : : !arrays
166 : : integer :: ngfft_in(18),ngfft_out(18)
167 : 0 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
168 : 0 : real(dp), contiguous, pointer :: arr_file(:,:),my_density(:,:)
169 : 0 : real(dp),allocatable :: rhor_file(:,:),rhog_in(:,:),rhor_out(:,:),rhog_out(:,:)
170 : 0 : type(pawrhoij_type),pointer:: pawrhoij__(:)
171 : : ! *************************************************************************
172 : :
173 : : DBG_ENTER("COLL")
174 : :
175 : 0 : ncplxfft = cplex*nfft
176 : :
177 : 0 : restartpaw=0
178 : 0 : my_fildata = fildata
179 : 0 : nspden = dtset%nspden; usewvl = dtset%usewvl
180 : :
181 : : ! Check validity of arguments--only rho(r) (51,52) and V(r) (101,102) are presently supported
182 [ # # ]: 0 : if ( (fform-1)/2 /=25 .and. (fform-1)/2 /=50 ) then
183 : 0 : write(msg,'(a,i0,a)')' Input fform= ',fform,' not allowed.'
184 : 0 : ABI_BUG(msg)
185 : : end if
186 : :
187 : : ! Print input fform
188 [ # # # # ]: 0 : if ( (fform-1)/2==25 .and. rdwr==1) then
189 : 0 : msg = ' ioarr: reading density data '
190 [ # # # # ]: 0 : else if ( (fform-1)/2==25 .and. rdwr==2) then
191 : 0 : msg = ' ioarr: writing density data'
192 [ # # # # ]: 0 : else if ( (fform-1)/2==50 .and. rdwr==1) then
193 : 0 : msg = ' ioarr: reading potential data'
194 [ # # # # ]: 0 : else if ( (fform-1)/2==50 .and. rdwr==2) then
195 : 0 : msg = ' ioarr: writing potential data'
196 : : end if
197 : 0 : call wrtout(std_out,msg)
198 : :
199 [ # # ]: 0 : call wrtout(std_out, 'ioarr: file name is: '//TRIM(fildata))
200 : :
201 [ # # ]: 0 : if (accessfil == IO_MODE_ETSF) then ! Initialize filename in case of ETSF file.
202 : 0 : file_etsf = nctk_ncify(fildata)
203 : 0 : call wrtout(std_out,sjoin('file name for ETSF access: ', file_etsf))
204 : : end if
205 : :
206 : : !Some definitions for MPI-IO access
207 : 0 : spaceComm = mpi_enreg%comm_cell
208 : 0 : comm_cell = mpi_enreg%comm_cell
209 : 0 : comm_fft = mpi_enreg%comm_fft
210 : :
211 [ # # ]: 0 : if (accessfil == 4) then
212 : 0 : iomode=IO_MODE_MPI
213 [ # # ]: 0 : if (rdwr==1) then
214 : : spaceComm=mpi_enreg%comm_cell
215 : : else
216 : 0 : spaceComm=mpi_enreg%comm_fft
217 : : end if
218 : 0 : me=xmpi_comm_rank(spaceComm)
219 [ # # ]: 0 : if (mpi_enreg%nproc_fft>1) then
220 : 0 : me_fft=mpi_enreg%me_fft
221 : 0 : spaceComm_io=mpi_enreg%comm_fft
222 : : else
223 : 0 : me_fft=0
224 : 0 : spaceComm_io=xmpi_comm_self
225 : : end if
226 : : end if
227 [ # # ]: 0 : if (usewvl==1) then
228 : 0 : spaceComm=mpi_enreg%comm_cell
229 : 0 : me=xmpi_comm_rank(spaceComm)
230 : : end if
231 : :
232 : : ! Change communicators and ranks if we are calling ioarr with one single processor.
233 [ # # ]: 0 : if (present(single_proc)) then
234 [ # # ]: 0 : if (single_proc) then
235 : 0 : spaceComm = xmpi_comm_self
236 : 0 : spaceComm_io = xmpi_comm_self
237 [ # # ]: 0 : ABI_CHECK(mpi_enreg%nproc_fft == 1, "single_proc cannot be used when nproc_fft > 1")
238 : 0 : comm_cell = xmpi_comm_self
239 : 0 : comm_fft = xmpi_comm_self
240 : 0 : me = 0
241 : : end if
242 : : end if
243 : :
244 : : !=======================================
245 : : !Handle input from disk file
246 : : !=======================================
247 : :
248 : 0 : call cwtime(cputime, walltime, gflops, "start")
249 : :
250 [ # # ]: 0 : if (rdwr==1) then
251 [ # # ]: 0 : if (accessfil == 0 .or. accessfil == 4) then
252 : :
253 : : ! Here master checks if the input rho(r) is given on a FFT mesh that quals
254 : : ! the one used in the run. If not, we perform a Fourier interpolation, we write the
255 : : ! interpolated rho(r) to a temporary file and we use this file to restart.
256 [ # # ]: 0 : if (ALLOW_FFTINTERP .and. usewvl==0) then
257 : 0 : need_fftinterp = .False.; icheck_fft = .True.
258 : : ! only master checks the FFT mesh if MPI-IO. All processors read ngfft if Fortran-IO
259 : : ! Note that, when Fortran-IO is used, we don't know if the routine is called
260 : : ! by a single processor or by all procs in comm_cell hence we cannot broadcast my_fildata
261 : : ! inside spaceComm as done if accessfil == 4
262 [ # # ]: 0 : if (accessfil == 4) icheck_fft = (xmpi_comm_rank(spaceComm)==master)
263 : :
264 [ # # ]: 0 : if (icheck_fft) then
265 [ # # ]: 0 : if (open_file(fildata,msg,newunit=in_unt,form='unformatted',status='old') /= 0) then
266 : 0 : ABI_ERROR(msg)
267 : : end if
268 : :
269 : 0 : call hdr_io(fform_dum,hdr0,rdwr,in_unt)
270 [ # # # # ]: 0 : need_fftinterp = (ANY(hdr%ngfft/=hdr0%ngfft) )
271 : 0 : qeq0=(hdr%qptn(1)**2+hdr%qptn(2)**2+hdr%qptn(3)**2<1.d-14)
272 : : ! FIXME: SHould handle double-grid if PAW
273 [ # # ]: 0 : nfftot_in = product(hdr0%ngfft(1:3))
274 [ # # ]: 0 : nfftot_out = product(hdr%ngfft(1:3))
275 : :
276 [ # # ]: 0 : if (need_fftinterp) then
277 : : write(msg, "(2a,2(a,3(i0,1x)))")&
278 : 0 : "Will perform Fourier interpolation since in and out ngfft differ",ch10,&
279 : 0 : "ngfft in file: ",hdr0%ngfft,", expected ngfft: ",hdr%ngfft
280 : 0 : ABI_WARNING(msg)
281 : :
282 : : ! Read rho(r) from file, interpolate it, write data and change fildata
283 [ # # # # : 0 : ABI_MALLOC(rhor_file, (cplex*nfftot_in, hdr0%nspden))
# # # # #
# # # # #
# # ]
284 [ # # # # : 0 : ABI_MALLOC(rhog_in, (2, nfftot_in))
# # # # #
# # # ]
285 [ # # # # : 0 : ABI_MALLOC(rhor_out, (cplex*nfftot_out, hdr0%nspden))
# # # # #
# # # # #
# # ]
286 [ # # # # : 0 : ABI_MALLOC(rhog_out, (2, nfftot_out))
# # # # #
# # # ]
287 : :
288 [ # # ]: 0 : do ispden=1,hdr0%nspden
289 [ # # ]: 0 : read(in_unt, err=10, iomsg=errmsg) (rhor_file(iarr,ispden), iarr=1,cplex*nfftot_in)
290 : : end do
291 : :
292 [ # # # # ]: 0 : ngfft_in = dtset%ngfft; ngfft_out = dtset%ngfft
293 [ # # # # ]: 0 : ngfft_in(1:3) = hdr0%ngfft(1:3); ngfft_out(1:3) = hdr%ngfft(1:3)
294 [ # # # # ]: 0 : ngfft_in(4:6) = hdr0%ngfft(1:3); ngfft_out(4:6) = hdr%ngfft(1:3)
295 [ # # # # ]: 0 : ngfft_in(9:18) = 0; ngfft_out(9:18) = 0
296 : 0 : ngfft_in(10) = 1; ngfft_out(10) = 1
297 : :
298 : 0 : call initmpi_seq(MPI_enreg_seq)
299 : : ! Which one is coarse? Note that this part is not very robust and can fail!
300 [ # # ]: 0 : if (ngfft_in(2) * ngfft_in(3) < ngfft_out(2) * ngfft_out(3)) then
301 : 0 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft_in(2),ngfft_in(3),'all')
302 : 0 : call MPI_enreg_seq%distribfft%init_seq('f',ngfft_out(2),ngfft_out(3),'all')
303 : : else
304 : 0 : call MPI_enreg_seq%distribfft%init_seq('f',ngfft_in(2),ngfft_in(3),'all')
305 : 0 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft_out(2),ngfft_out(3),'all')
306 : : end if
307 : :
308 : : call fourier_interpol(cplex,hdr0%nspden,0,0,nfftot_in,ngfft_in,nfftot_out,ngfft_out,&
309 : 0 : MPI_enreg_seq,rhor_file,rhor_out,rhog_in,rhog_out)
310 : :
311 : 0 : call destroy_mpi_enreg(MPI_enreg_seq)
312 : :
313 : : ! MG Hack: Change fildata so that we will use this file to read the correct rho(r)
314 : : ! FIXME: This should be done in a cleaner way!
315 [ # # ]: 0 : my_fildata = trim(fildata)//"__fftinterp_rhor__"
316 [ # # ]: 0 : if (my_fildata == fildata) my_fildata = "__fftinterp_rhor__"
317 [ # # ]: 0 : if (open_file(my_fildata,msg,newunit=out_unt,form='unformatted',status='unknown') /= 0) then
318 : 0 : ABI_ERROR(msg)
319 : : end if
320 : 0 : call hdr_io(fform_dum,hdr,2,out_unt)
321 [ # # ]: 0 : do ispden=1,hdr0%nspden
322 [ # # ]: 0 : write(out_unt, err=10, iomsg=errmsg) (rhor_out(iarr,ispden),iarr=1,cplex*nfftot_out)
323 : : end do
324 : 0 : close(out_unt)
325 : :
326 [ # # ]: 0 : ABI_FREE(rhor_file)
327 [ # # ]: 0 : ABI_FREE(rhog_in)
328 [ # # ]: 0 : ABI_FREE(rhor_out)
329 [ # # ]: 0 : ABI_FREE(rhog_out)
330 : : end if ! need_fftinterp
331 : :
332 : 0 : call hdr0%free()
333 [ # # ]: 0 : close(in_unt, err=10, iomsg=errmsg)
334 : : end if ! master
335 [ # # ]: 0 : if (accessfil == 4) call xmpi_bcast(my_fildata,master,spaceComm,ierr)
336 : : end if
337 : :
338 [ # # ]: 0 : if (accessfil == 4) then
339 : 0 : unt = get_unit()
340 : 0 : call WffOpen(iomode,spaceComm,my_fildata,ierr,wff,0,me,unt,spaceComm_io)
341 : 0 : call hdr_io(fform_dum,hdr0,rdwr,wff)
342 : : ! Compare the internal header and the header from the file
343 : 0 : call hdr_check(fform,fform_dum,hdr,hdr0,'COLL',restart,restartpaw)
344 : :
345 : : else
346 [ # # ]: 0 : if (open_file(my_fildata, msg, newunit=unt, form="unformatted", status="old", action="read") /= 0) then
347 : 0 : ABI_ERROR(msg)
348 : : end if
349 : : ! Initialize hdr0, thanks to reading of unwff1
350 : 0 : call hdr_io(fform_dum,hdr0,rdwr,unt)
351 : : ! Compare the internal header and the header from the file
352 : 0 : call hdr_check(fform,fform_dum,hdr,hdr0,'COLL',restart,restartpaw)
353 : : end if
354 : 0 : etotal=hdr0%etot
355 : :
356 : : ! NOTE: should check that restart is possible !!
357 : : !call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
358 : :
359 : : ! If nspden[file] /= nspden, need a temporary array
360 [ # # ]: 0 : if (hdr0%nspden/=nspden) then
361 [ # # # # : 0 : ABI_MALLOC(arr_file,(cplex*nfft,hdr0%nspden))
# # # # #
# # # # #
# # ]
362 : : else
363 : 0 : arr_file => arr
364 : : end if
365 : :
366 : : ! Read data
367 [ # # ]: 0 : do ispden=1,hdr0%nspden
368 [ # # ]: 0 : if(accessfil == 4) then
369 : 0 : call xderiveRRecInit(wff,ierr)
370 : 0 : call xderiveRead(wff,arr(1:ncplxfft,ispden),ncplxfft,spaceComm_io,ierr)
371 : 0 : call xderiveRRecEnd(wff,ierr)
372 : :
373 : : !call xmpio_read_dp(mpi_fh,offset,sc_mode,ncount,buf,fmarker,mpierr,advance)
374 : : !do idat=1,ndat
375 : : ! do i3=1,n3
376 : : ! if( fftn3_distrib(i3) == me_fft) then
377 : : ! i3_local = ffti3_local(i3)
378 : : ! i3_ldat = i3_local + (idat - 1) * nd3proc
379 : : ! do i2=1,n2
380 : : ! frbase=n1*(i2-1+n2*(i3_local-1)) + (idat - 1) * nfft
381 : : ! do i1=1,n1
382 : : ! fofr(i1+frbase)=workr(1,i1,i2,i3_ldat)
383 : : ! end do
384 : : ! end do
385 : : ! end if
386 : : ! end do
387 : : !end do
388 : :
389 : : else
390 [ # # ]: 0 : read(unt, err=10, iomsg=errmsg) (arr_file(iarr,ispden),iarr=1,ncplxfft)
391 : : end if
392 : : end do
393 : :
394 [ # # ]: 0 : if (accessfil == 4) then
395 : 0 : call wffclose(wff,ierr)
396 : : else
397 [ # # ]: 0 : close (unit=unt, err=10, iomsg=errmsg)
398 : : end if
399 : :
400 [ # # ]: 0 : else if (accessfil == 3) then
401 : :
402 : : ! Read the header and broadcast it in comm_cell
403 : : ! FIXME: Use xmpi_comm_self for the time-being because, in loper, ioarr
404 : : ! is called by me==0
405 : 0 : call hdr0%from_fname(file_etsf, fform_dum, comm_cell)
406 [ # # ]: 0 : ABI_CHECK(fform_dum/=0, "hdr_read_from_fname returned fform 0")
407 : :
408 : : ! Compare the internal header and the header from the file
409 : 0 : call hdr_check(fform, fform_dum, hdr, hdr0, 'COLL', restart, restartpaw)
410 : :
411 : : ! If nspden[file] /= nspden, need a temporary array
412 [ # # ]: 0 : if (hdr0%nspden /= nspden) then
413 [ # # # # : 0 : ABI_MALLOC(arr_file,(cplex*nfft,hdr0%nspden))
# # # # #
# # # # #
# # ]
414 : : else
415 : 0 : arr_file => arr
416 : : end if
417 : :
418 [ # # ]: 0 : if (usewvl == 1) then
419 : : ! Read the array
420 [ # # ]: 0 : if (fform==52) then ! density
421 : 0 : varname = "density"
422 [ # # ]: 0 : else if (fform==102) then ! all potential forms!!!!
423 : 0 : varname = "exchange_correlation_potential"
424 : : end if
425 : :
426 : : ! Open the file
427 [ # # ]: 0 : NCF_CHECK(nctk_open_read(ncid, file_etsf, xmpi_comm_self))
428 [ # # ]: 0 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, varname), arr_file))
429 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
430 : : else
431 : : ! Get MPI-FFT tables from input ngfft
432 : 0 : call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
433 : :
434 : : ! Get the name of the netcdf variable from the ABINIT extension and read data.
435 : 0 : varname = varname_from_fname(file_etsf)
436 : 0 : ncerr = nctk_read_datar(file_etsf,varname,ngfft,cplex,nfft,hdr0%nspden,comm_fft,fftn3_distrib,ffti3_local,arr)
437 [ # # ]: 0 : NCF_CHECK(ncerr)
438 : : end if
439 : :
440 : : else
441 : 0 : write(msg,'(a,i0,a)')'Bad value for accessfil', accessfil, ' on read '
442 : 0 : ABI_BUG(msg)
443 : : end if
444 : :
445 : 0 : call wrtout(std_out,sjoin("data read from disk file: ", fildata))
446 : :
447 : 0 : etotal=hdr0%etot
448 : :
449 : : ! Possibly need to convert the potential/density spin components
450 [ # # ]: 0 : if (hdr0%nspden/=nspden) then
451 : 0 : call denpot_spin_convert(arr_file,hdr0%nspden,arr,nspden,fform)
452 [ # # ]: 0 : ABI_FREE(arr_file)
453 : : end if
454 : :
455 : : ! Eventually copy (or distribute) PAW data
456 [ # # # # ]: 0 : if (rdwrpaw==1.and.restartpaw/=0) then
457 : 0 : pawrhoij__ => hdr0%pawrhoij ! Trick needed by nvhpc 23.9
458 [ # # ]: 0 : if (size(pawrhoij__) /= size(pawrhoij)) then
459 : : call pawrhoij_copy(hdr0%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
460 : 0 : keep_nspden=.true.)
461 : : else
462 : 0 : call pawrhoij_copy(hdr0%pawrhoij,pawrhoij,keep_nspden=.true.)
463 : : end if
464 : : end if
465 : :
466 [ # # # # ]: 0 : if (accessfil == 0 .or. accessfil == 3 .or. accessfil == 4) call hdr0%free()
467 : :
468 : : ! =======================================
469 : : ! Set up for writing data
470 : : ! =======================================
471 [ # # ]: 0 : else if (rdwr==2) then
472 : :
473 : : ! In the wavelet case (isolated boundary counditions), the
474 : : ! arr array has a buffer that we need to remove.
475 [ # # ]: 0 : if (usewvl == 1) then
476 : : #ifdef HAVE_BIGDFT
477 : : zindex = wvl_den%denspot%dpbox%nscatterarr(me, 3)
478 : : if (wvl_den%denspot%rhod%geocode == 'F') then
479 : : n1 = (wvl_den%denspot%dpbox%ndims(1) - 31) / 2
480 : : n2 = (wvl_den%denspot%dpbox%ndims(2) - 31) / 2
481 : : n3 = (wvl_den%denspot%dpbox%ndims(3) - 31) / 2
482 : : zstart = max(15 - zindex, 0)
483 : : zstop = wvl_den%denspot%dpbox%nscatterarr(me, 2) + &
484 : : & wvl_den%denspot%dpbox%nscatterarr(me, 4) - &
485 : : & max(zindex + wvl_den%denspot%dpbox%nscatterarr(me, 2) &
486 : : & - 2 * n3 - 15, 0)
487 : : else
488 : : ABI_ERROR('ioarr: WVL not implemented yet.')
489 : : end if
490 : : if (zstop - zstart + 1 > 0) then
491 : : ! Our slab contains (zstop - zstart + 1) elements
492 : : ABI_MALLOC(my_density,((n1*2)*(n2*2)*(zstop-zstart),nspden))
493 : : ! We copy the data except the buffer to my_density
494 : : ind = 0
495 : :
496 : : do i3 = zstart, zstop - 1, 1
497 : : ia = (i3 - 1) * dtset%ngfft(1) * dtset%ngfft(2)
498 : : do i2 = 0, 2 * n2 - 1, 1
499 : : i = ia + (i2 + 14) * dtset%ngfft(1) + 14
500 : : do i1 = 0, 2 * n1 - 1, 1
501 : : i = i + 1
502 : : ind = ind + 1
503 : : my_density(ind, :) = arr(i, :)
504 : : end do
505 : : end do
506 : : end do
507 : : else
508 : : nullify(my_density)
509 : : end if
510 : : #else
511 : 0 : BIGDFT_NOTENABLED_ERROR()
512 : : if(.false. .and. present(wvl_den))then
513 : : write(std_out,*)' One should not be here'
514 : : endif
515 : : #endif
516 : : end if
517 : :
518 : : ! Make sure ngfft agrees with hdr%ngfft.
519 [ # # ]: 0 : if (usewvl == 0) then
520 [ # # # # : 0 : if (any(ngfft(:3) /= hdr%ngfft(:3))) then
# # ]
521 : 0 : write(msg,"(2(a,3(1x,i0)))")"input ngfft: ",ngfft(:3),"differs from hdr%ngfft: ",hdr%ngfft(:3)
522 : 0 : ABI_ERROR(msg)
523 : : end if
524 : : end if
525 : :
526 [ # # ]: 0 : if (accessfil == 0 .or. accessfil == 4) then
527 [ # # ]: 0 : if(accessfil == 4) then
528 : 0 : unt = get_unit()
529 : 0 : call WffOpen(iomode,spaceComm,fildata,ierr,wff,0,me,unt)
530 : 0 : call hdr_io(fform,hdr,rdwr,wff)
531 : : else
532 [ # # ]: 0 : if (open_file(fildata, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then
533 : 0 : ABI_ERROR(msg)
534 : : end if
535 : :
536 : : ! Write header
537 : 0 : call hdr_io(fform,hdr,rdwr,unt)
538 : : end if
539 : :
540 : : ! Write actual data
541 [ # # ]: 0 : do ispden=1,nspden
542 [ # # ]: 0 : if(accessfil == 4) then
543 : 0 : call xderiveWRecInit(wff,ierr,me_fft)
544 : 0 : call xderiveWrite(wff,arr(1:ncplxfft,ispden),ncplxfft,spaceComm_io,ierr)
545 : 0 : call xderiveWRecEnd(wff,ierr,me_fft)
546 : : else
547 [ # # ]: 0 : if (usewvl == 0) then
548 [ # # ]: 0 : write(unt, err=10, iomsg=errmsg) (arr(iarr,ispden),iarr=1,ncplxfft)
549 : : else
550 [ # # ]: 0 : write(unt, err=10, iomsg=errmsg) (my_density(iarr,ispden),iarr=1,size(my_density, 1))
551 : : end if
552 : : end if
553 : : end do
554 : :
555 [ # # ]: 0 : if(accessfil == 4) then
556 : 0 : call WffClose(wff,ierr)
557 : : else
558 [ # # ]: 0 : close(unt, err=10, iomsg=errmsg)
559 : : end if
560 : :
561 [ # # ]: 0 : else if ( accessfil == 3 ) then
562 : :
563 : : ! Master in comm_fft creates the file and writes the header.
564 [ # # ]: 0 : if (xmpi_comm_rank(comm_fft) == 0) then
565 : 0 : call hdr%write_to_fname(file_etsf, fform)
566 : : end if
567 : 0 : call xmpi_barrier(comm_fft)
568 : :
569 : : ! Write the array
570 [ # # ]: 0 : if (usewvl == 0) then
571 : : ! Get MPI-FFT tables from input ngfft
572 : 0 : call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
573 : :
574 : 0 : varname = varname_from_fname(file_etsf)
575 : 0 : ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden,comm_fft,fftn3_distrib,ffti3_local,arr)
576 [ # # ]: 0 : NCF_CHECK(ncerr)
577 : : else
578 [ # # ]: 0 : NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self))
579 : :
580 [ # # ]: 0 : if (fform==52) then ! density
581 : 0 : varname = "density"
582 : : if (usewvl == 0) then
583 : : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), arr))
584 : : else
585 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), my_density))
586 : : end if
587 [ # # ]: 0 : else if (fform==102) then ! all potential forms!!!!
588 : 0 : varname = "exchange_correlation_potential"
589 [ # # ]: 0 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, varname), arr))
590 : : end if
591 : :
592 [ # # ]: 0 : NCF_CHECK(nf90_close(ncid))
593 : : end if
594 : :
595 : : else
596 : 0 : write(msg,'(a,i0,a)')'Bad value for accessfil', accessfil, ' on write '
597 : 0 : ABI_ERROR(msg)
598 : : end if
599 : :
600 [ # # # # ]: 0 : if (usewvl == 1 .and. associated(my_density)) then
601 : 0 : ABI_FREE(my_density)
602 : : end if
603 : :
604 : 0 : call wrtout(std_out,sjoin(' Data written to disk file:', fildata))
605 : :
606 : : else
607 : 0 : write(msg,'(a,i0,a)')'Called with rdwr = ',rdwr,' not allowed.'
608 : 0 : ABI_BUG(msg)
609 : : end if
610 : :
611 : 0 : call cwtime_report(" IO operation", cputime, walltime, gflops)
612 : :
613 : : DBG_EXIT("COLL")
614 : :
615 : 0 : return
616 : :
617 : : ! Handle Fortran IO error
618 : : 10 continue
619 : 0 : ABI_ERROR(errmsg)
620 : :
621 [ # # # # : 0 : end subroutine ioarr
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
622 : : !!***
623 : :
624 : : !----------------------------------------------------------------------
625 : :
626 : : !!****f* m_ioarr/fftdatar_write
627 : : !! NAME
628 : : !! fftdatar_write
629 : : !!
630 : : !! FUNCTION
631 : : !! Write an array in real space on the FFT box to file.
632 : : !! The array can be real or complex depending on cplex
633 : : !! IO library is automatically selected from the file extension and the number of FFT processors:
634 : : !!
635 : : !! 1) If path ends with ".nc", the netcdf library is used else Fortran format.
636 : : !!
637 : : !! 2) If nproc_fft > 1, parallel IO is used (if available)
638 : : !!
639 : : !! INPUTS
640 : : !! varname=Name of the variable to write (used if ETSF-IO).
641 : : !! path=File name
642 : : !! iomode=
643 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
644 : : !! crystal<crystal_t>= data type gathering info on symmetries and unit cell (used if etsf_io)
645 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
646 : : !! cplex=1 for real array, 2 for complex
647 : : !! nfft=Number of FFT points treated by this node.
648 : : !! nspden=Number of spin-density components.
649 : : !! datar(cplex*nfft,nspden)=array on the real space FFT grid.
650 : : !! mpi_enreg=information about MPI parallelization
651 : : !! [ebands]<ebands_t>=data type with energies and occupations (used if etsf_io)
652 : : !!
653 : : !! OUTPUT
654 : : !! Only writing
655 : : !!
656 : : !! NOTES
657 : : !! The string passed to fftdatar_write (first argument) gives the name used to store the data in the netcdf file
658 : : !! The function varname_from_fname defined in the module m_hdr.F90 gives the mapping between the Abinit
659 : : !! file extension and the netcdf name e.g. foo_VHXC.nc --> vxc
660 : : !! This function is used in cut3d so that we can immediately select the data to analyze without having
661 : : !! to prompt the user
662 : : !! Remember to update varname_from_fname if you add a new file or if you change the name of the variable.
663 : : !!
664 : : !! fform i.e. the integer specification for data type is automatically initialized from varname.
665 : : !!
666 : : !! SOURCE
667 : :
668 : 7842 : subroutine fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands)
669 : :
670 : : !Arguments ------------------------------------
671 : : !scalars
672 : : integer,intent(in) :: iomode,cplex,nfft,nspden
673 : : character(len=*),intent(in) :: varname,path
674 : : type(hdr_type),intent(inout) :: hdr
675 : : type(crystal_t),intent(in) :: crystal
676 : : type(ebands_t),optional,intent(in) :: ebands
677 : : type(MPI_type),intent(in) :: mpi_enreg
678 : : !arrays
679 : : integer,intent(in) :: ngfft(18)
680 : : real(dp),intent(inout) :: datar(cplex*nfft,nspden)
681 : : !type(pawrhoij_type),optional,intent(inout) :: pawrhoij_all(hdr%usepaw*crystal%natom)
682 : :
683 : : !Local variables-------------------------------
684 : : !!scalars
685 : : integer,parameter :: master=0
686 : : integer :: n1,n2,n3,comm_fft,nproc_fft,me_fft,iarr,ierr,ii,ispden,unt,mpierr,fform
687 : : integer :: i3_glob,my_iomode
688 : : integer(kind=XMPI_OFFSET_KIND) :: hdr_offset,my_offset,nfft_tot
689 : : integer :: ncid,ncerr
690 : : character(len=fnlen) :: file_etsf
691 : : real(dp) :: cputime,walltime,gflops
692 : : character(len=500) :: msg,errmsg
693 : : type(abifile_t) :: abifile
694 : : !arrays
695 : 7842 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
696 : 7842 : integer(XMPI_OFFSET_KIND) :: bsize_frecord(nspden)
697 : : ! *************************************************************************
698 : :
699 : 7842 : abifile = abifile_from_varname(varname)
700 : 7842 : if (abifile%fform == 0) then
701 : 0 : ABI_ERROR(sjoin("Cannot find any abifile object associated to varname:", varname))
702 : : end if
703 : : ! Get fform from abifile. TODO: check file extension
704 : 7842 : fform = abifile%fform
705 : :
706 : 7842 : comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
707 : 7842 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = n1*n2*n3
708 : :
709 : : ! Select iomode
710 : : ! Use Fortran IO if nproc_fft 1, in principle this is not needed because the
711 : : ! MPI-IO code should produce binary files that are readable with Fortran-IO
712 : : ! but it seems that NAG uses its own binary format
713 : 7842 : my_iomode = iomode
714 [ + + ]: 7842 : if (my_iomode /= IO_MODE_ETSF .and. nproc_fft == 1) my_iomode = IO_MODE_FORTRAN
715 [ + + ]: 7842 : if (nproc_fft > 1 .and. my_iomode == IO_MODE_FORTRAN) my_iomode = IO_MODE_MPI
716 : :
717 : 7842 : call wrtout(std_out, sjoin(ch10, "fftdatar_write: About to write data to:", path, "with iomode:",iomode2str(my_iomode)))
718 : 7842 : call cwtime(cputime, walltime, gflops, "start")
719 : :
720 : : ! Get MPI-FFT tables from input ngfft
721 : 7842 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
722 : :
723 : 28 : select case (my_iomode)
724 : : case (IO_MODE_FORTRAN)
725 [ - + ]: 28 : ABI_CHECK(nproc_fft == 1, "MPI-IO must be enabled when FFT parallelism is used")
726 [ - + ]: 28 : if (open_file(path, msg, newunit=unt, form='unformatted', status='unknown', action="write") /= 0) then
727 : 0 : ABI_ERROR(msg)
728 : : end if
729 : 28 : call hdr%fort_write(unt, fform, ierr)
730 [ - + ]: 28 : ABI_CHECK(ierr==0, "ierr !=0")
731 [ + + ]: 59 : do ii=1,nspden
732 [ + - ]: 59 : write(unt, err=10, iomsg=errmsg) (datar(iarr,ii), iarr=1,cplex * nfft)
733 : : end do
734 [ + - ]: 28 : close(unt, err=10, iomsg=errmsg)
735 : :
736 : : ! Write PAW rhoij
737 : : !call pawrhoij_io(hdr%pawrhoij,unit,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,headform,"Write")
738 : : !call pawrhoij_io(rhoij_ptr,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,&
739 : : ! HDR_LATEST_HEADFORM,"Write",form="netcdf")
740 : :
741 : : #ifdef HAVE_MPI_IO
742 : : case (IO_MODE_MPI)
743 : : ! Find the first z-plane treated by this node.
744 : : ! WARNING: Here I assume that the z-planes in real space
745 : : ! are distributed in contiguous blocks (as usually done in MPI-FFT)
746 [ + - ]: 1188 : do i3_glob=1,n3
747 [ + + ]: 1188 : if (me_fft == fftn3_distrib(i3_glob)) exit
748 : : end do
749 [ - + ]: 134 : ABI_CHECK(i3_glob /= n3 +1, "This processor does not have z-planes!")
750 : :
751 : : ! Master writes the header.
752 [ + + ]: 134 : if (me_fft == master) call hdr%write_to_fname(path, fform)
753 : 134 : call xmpi_barrier(comm_fft) ! TODO: Non-blocking barrier.
754 : :
755 : 134 : call MPI_FILE_OPEN(comm_fft, path, MPI_MODE_RDWR, xmpio_info, unt, mpierr)
756 : 134 : ABI_CHECK_MPI(mpierr,"MPI_FILE_OPEN")
757 : :
758 : : ! Skip the header and get the offset of the header
759 : 134 : call hdr_mpio_skip(unt,fform,hdr_offset)
760 : : !write(std_out,*)"i3_glob, nfft, hdr_offset,",i3_glob,nfft,hdr_offset,fftn3_distrib == me_fft
761 : :
762 : : ! Each proc writes a contiguous slice of the nspden records.
763 : : ! my_offset is the position inside the Fortran record.
764 [ + + ]: 356 : do ispden=1,nspden
765 : : my_offset = hdr_offset + xmpio_bsize_frm + ((ispden - 1) * 2 * xmpio_bsize_frm) + &
766 : 222 : ((i3_glob-1) * cplex * n1 * n2 * xmpi_bsize_dp) + ((ispden-1) * cplex * nfft_tot * xmpi_bsize_dp)
767 : 222 : call MPI_FILE_WRITE_AT_ALL(unt,my_offset,datar(:,ispden),cplex*nfft,MPI_DOUBLE_PRECISION,MPI_STATUS_IGNORE,mpierr)
768 : 356 : ABI_CHECK_MPI(mpierr,"MPI_FILE_WRITE_AT_ALL")
769 : : end do
770 : :
771 : : ! master writes the fortran record markers.
772 [ + + ]: 134 : if (me_fft == master) then
773 [ + + ]: 178 : bsize_frecord = cplex * nfft_tot * xmpi_bsize_dp
774 : : #if 1
775 : 67 : my_offset = hdr_offset
776 [ + + ]: 178 : do ispden=1,nspden
777 : 111 : call xmpio_write_frm(unt,my_offset,xmpio_single,bsize_frecord(ispden),mpierr)
778 : 178 : ABI_CHECK_MPI(mpierr,"xmpio_write_frm")
779 : : end do
780 : : #else
781 : : ! TODO: Understand why this code does not work!
782 : : call xmpio_write_frmarkers(unt,hdr_offset,xmpio_single,nspden,bsize_frecord,ierr)
783 : : ABI_CHECK(ierr==0, "xmpio_write_frmarkers")
784 : : #endif
785 : : end if
786 : :
787 : 134 : call MPI_FILE_CLOSE(unt,mpierr)
788 : 134 : ABI_CHECK_MPI(mpierr,"FILE_CLOSE!")
789 : :
790 : : ! Add full pawrhoij datastructure at the end of the file.
791 : : !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then
792 : : ! if (open_file(path, msg, newunit=unt, form='unformatted', status='old', action="write", access="append") /= 0) then
793 : : ! ABI_ERROR(msg)
794 : : ! end if
795 : : ! call pawrhoij_io(pawrhoij_all,un,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write")
796 : : ! close(unt)
797 : : !end if
798 : : #endif
799 : :
800 : : case (IO_MODE_ETSF)
801 : 7680 : file_etsf = nctk_ncify(path)
802 : :
803 : : ! Write datar.
804 : : ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden, &
805 : 7680 : comm_fft,fftn3_distrib,ffti3_local,datar,action="create")
806 [ - + ]: 7680 : NCF_CHECK(ncerr)
807 : 7680 : call xmpi_barrier(comm_fft)
808 : :
809 : : ! Master writes the header.
810 [ + + ]: 7680 : if (xmpi_comm_rank(comm_fft) == master) then
811 [ - + ]: 7522 : NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self))
812 [ - + ]: 7522 : NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.))
813 : : ! Add information on the crystalline structure.
814 [ - + ]: 7522 : NCF_CHECK(crystal%ncwrite(ncid))
815 [ + + ]: 7522 : if (present(ebands)) then
816 [ - + ]: 3674 : NCF_CHECK(ebands%ncwrite(ncid))
817 : : end if
818 : :
819 : : ! Add full pawrhoij datastructure.
820 : : !if (present(pawrhoij_all) .and. me_fft == master .and. hdr%usepaw == 1) then
821 : : ! call pawrhoij_io(pawrhoij_all,ncid,hdr%nsppol,hdr%nspinor,hdr%nspden,hdr%lmn_size,hdr%typat,hdr%headform,"Write", form="netcdf")
822 : : !end if
823 : :
824 [ - + ]: 7522 : NCF_CHECK(nf90_close(ncid))
825 : : end if
826 : :
827 : : case default
828 [ + + + - ]: 7976 : ABI_ERROR(sjoin("Wrong iomode:",itoa(my_iomode)))
829 : : end select
830 : :
831 : 7842 : call cwtime_report(" IO operation", cputime, walltime, gflops)
832 : :
833 : : return
834 : :
835 : : ! Handle Fortran IO error
836 : : 10 continue
837 : 0 : ABI_ERROR(errmsg)
838 : :
839 [ - + ]: 7842 : end subroutine fftdatar_write
840 : : !!***
841 : :
842 : : !----------------------------------------------------------------------
843 : :
844 : : !!****f* m_ioarr/fftdatar_write_from_hdr
845 : : !! NAME
846 : : !! fftdatar_write_from_hdr
847 : : !!
848 : : !! FUNCTION
849 : : !! Write an array in real space on the FFT box to file.
850 : : !! crystal and ebands are constructed from the Abinit header.
851 : : !!
852 : : !! TODO
853 : : !! This routine will be removed when crystal_t and ebands_t will become standard objects
854 : : !! available in the GS/DFPT part.
855 : : !!
856 : : !! INPUTS
857 : : !! [eigen](mband*hdr%nkpt*hdr%nsppol)=GS eigenvalues
858 : : !! See fftdatar_write for the meaning of the other variables.
859 : : !!
860 : : !! OUTPUT
861 : : !!
862 : : !! SOURCE
863 : :
864 [ + + + - : 3883 : subroutine fftdatar_write_from_hdr(varname,path,iomode,hdr,ngfft,cplex,nfft,nspden,datar,mpi_enreg,eigen)
- + ]
865 : :
866 : : !Arguments ------------------------------------
867 : : !scalars
868 : : integer,intent(in) :: iomode,cplex,nfft,nspden
869 : : character(len=*),intent(in) :: varname,path
870 : : type(hdr_type),intent(inout) :: hdr
871 : : type(MPI_type),intent(in) :: mpi_enreg
872 : : !arrays
873 : : integer,intent(in) :: ngfft(18)
874 : : real(dp),intent(inout) :: datar(cplex*nfft,nspden)
875 : : real(dp),optional,intent(in) :: eigen(:)
876 : :
877 : : !Local variables-------------------------------
878 : : !!scalars
879 : : integer :: mband
880 : : type(crystal_t) :: crystal
881 : 3883 : type(ebands_t) :: ebands
882 : : !arrays
883 : 3883 : real(dp),allocatable :: ene3d(:,:,:)
884 : : ! *************************************************************************
885 : :
886 : 3883 : crystal = hdr%get_crystal()
887 : :
888 [ + + ]: 3883 : if (present(eigen)) then
889 [ + + ]: 1476 : mband = maxval(hdr%nband)
890 [ - + ]: 36 : ABI_CHECK(size(eigen) == mband * hdr%nkpt * hdr%nsppol, "Wrong size(eigen)")
891 [ - + + - : 180 : ABI_MALLOC(ene3d, (mband, hdr%nkpt, hdr%nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
892 : 36 : call unpack_eneocc(hdr%nkpt, hdr%nsppol, mband, hdr%nband, eigen, ene3d)
893 : 36 : call ebands%from_hdr(hdr, mband, ene3d)
894 [ - + ]: 36 : ABI_FREE(ene3d)
895 : :
896 : 36 : call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg,ebands=ebands)
897 : 36 : call ebands%free()
898 : : else
899 : 3847 : call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg)
900 : : end if
901 : :
902 : 3883 : call crystal%free()
903 : :
904 [ - + - + : 3883 : end subroutine fftdatar_write_from_hdr
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
905 : : !!***
906 : :
907 : : !----------------------------------------------------------------------
908 : :
909 : : !!****f* m_ioarr/read_rhor
910 : : !! NAME
911 : : !! read_rhor
912 : : !!
913 : : !! FUNCTION
914 : : !! Read the DEN file with name fname reporting the density on the real FFT mesh
915 : : !! specified through the input variable ngfft. If the FFT mesh asked in input and that found
916 : : !! on file differ, the routine performs a FFT interpolation and renormalize the density so that it
917 : : !! integrates to the correct number of electrons. The interpolation is done only for NC.
918 : : !! For PAW, this is not possible because one should include the onsite contribution so this task
919 : : !! is delegated to the caller.
920 : : !!
921 : : !! INPUTS
922 : : !! fname=Name of the file
923 : : !! cplex=1 if array is real, 2 if complex e.g. DFPT density.
924 : : !! nspden=Number of spin density components.
925 : : !! nfft=Number of FFT points (treated by this processor)
926 : : !! ngfft(18)=Info on the FFT mesh.
927 : : !! pawread= 1 if pawrhoij should be read from file, 0 otherwise. Meaningful only if usepaw==1.
928 : : !! mpi_enreg<MPI_type>=Information about MPI parallelization
929 : : !! comm=MPI communicator. See notes
930 : : !! [check_hdr] <type(hdr_type)>=Optional. Used to compare with the hdr read from disk file
931 : : !! The routine will abort if restart cannot be performed.
932 : : !! [allow_interp]=If True, the density read from file will be interpolated if the mesh differs from the one
933 : : !! expected by the caller. This option is usually used in **self-consistent** calculations.
934 : : !! If False (default), the code stops if the two meshes are different.
935 : : !! [varname]=If present, check whether file contains varname
936 : : !!
937 : : !! OUTPUT
938 : : !! orhor(cplex*nfft,nspden)=The density on the real space mesh.
939 : : !! ohdr=Abinit header read from file.
940 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data. only
941 : : !! if pawread==1. The arrays is supposed to be already allocated in the caller and its
942 : : !! size must be consistent with the MPI communicator comm.
943 : : !!
944 : : !! NOTES
945 : : !! if xmpi_comm_size(comm) == 1, nfft shall be equal to nfftot, and len(pawrhoij) == natom
946 : : !! This means that one can call this routine with
947 : : !!
948 : : !! if (xmpi_comm_rank(comm) == 0) call read_rhor(...., comm=xmpi_comm_self)
949 : : !!
950 : : !! to get the full array and pawrhoij(natom) on the master node.
951 : : !!
952 : : !! if xmpi_comm_size(comm) > 1, nfft represents the number of FFT points treated by this processor,
953 : : !! and pawrhoij is dimensioned with my_natom
954 : : !! All the processors inside comm and comm_atom should call this routine.
955 : : !!
956 : : !! SOURCE
957 : :
958 [ - + - + : 2186 : subroutine read_rhor(fname, cplex, nspden, nfft, ngfft, pawread, mpi_enreg, orhor, ohdr, pawrhoij, comm, &
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - -
- ]
959 : : check_hdr, allow_interp, want_varname) ! Optional
960 : :
961 : : !Arguments ------------------------------------
962 : : !scalars
963 : : integer,intent(in) :: cplex,nfft,nspden,pawread,comm
964 : : character(len=*),intent(in) :: fname
965 : : type(MPI_type),intent(in) :: mpi_enreg
966 : : type(hdr_type),intent(out) :: ohdr
967 : : type(hdr_type),optional,intent(in) :: check_hdr
968 : : logical,optional,intent(in) :: allow_interp
969 : : character(len=*),optional,intent(in) :: want_varname
970 : : !arrays
971 : : integer,intent(in) :: ngfft(18)
972 : : real(dp),intent(out) :: orhor(cplex*nfft,nspden)
973 : : type(pawrhoij_type),intent(inout) :: pawrhoij(:)
974 : :
975 : : !Local variables-------------------------------
976 : : !scalars
977 : : integer,parameter :: master=0
978 : : integer :: unt,fform,iomode,my_rank,mybase,globase,cplex_file
979 : : integer :: ispden,ifft,nfftot_file,nprocs,ierr,i1,i2,i3,i3_local,n1,n2,n3
980 : : integer,parameter :: fform_den=52
981 : : integer :: restart, restartpaw
982 : : integer :: ncerr
983 : : real(dp) :: ratio,ucvol
984 : : real(dp) :: cputime,walltime,gflops
985 : : logical :: need_interp,have_mpifft,allow_interp__
986 : : character(len=500) :: msg,errmsg
987 : : character(len=fnlen) :: my_fname
988 : : character(len=nctk_slen) :: varname
989 : : !arrays
990 : 2186 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
991 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2)
992 : 2186 : real(dp),allocatable :: rhor_file(:,:),rhor_tmp(:,:)
993 : 2186 : type(pawrhoij_type),allocatable :: pawrhoij_file(:)
994 : : ! *************************************************************************
995 : :
996 : 2186 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
997 [ + + ]: 8744 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); have_mpifft = (nfft /= product(ngfft(1:3)))
998 [ + + ]: 2186 : allow_interp__ = .False.; if (present(allow_interp)) allow_interp__ = allow_interp
999 : :
1000 : 2186 : call timab(1280,1,tsec)
1001 : 2186 : call wrtout(std_out, sjoin(" About to read data(r) from:", fname), do_flush=.True.)
1002 : 2186 : call cwtime(cputime, walltime, gflops, "start")
1003 : :
1004 : : ! Master node opens the file, read the header and the FFT data
1005 : : ! This approach facilitates the interpolation of the density if in_ngfft(1:3) /= file_ngfft(1:3)
1006 [ + + ]: 2186 : if (my_rank == master) then
1007 [ + + ]: 1990 : my_fname = fname
1008 [ - + ]: 1990 : if (nctk_try_fort_or_ncfile(my_fname, msg) /= 0 ) then
1009 : 0 : ABI_ERROR(msg)
1010 : : end if
1011 : :
1012 : 1990 : iomode = iomode_from_fname(my_fname)
1013 : 3 : select case (iomode)
1014 : :
1015 : : case (IO_MODE_FORTRAN, IO_MODE_MPI)
1016 [ - + ]: 3 : if (open_file(my_fname, msg, newunit=unt, form='unformatted', status='old', action="read") /= 0) then
1017 : 0 : ABI_ERROR(msg)
1018 : : end if
1019 : :
1020 : 3 : call ohdr%fort_read(unt, fform)
1021 : :
1022 : : ! Check important dimensions.
1023 [ - + ]: 3 : ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1024 [ - + ]: 3 : if (present(want_varname)) then
1025 [ # # ]: 0 : ABI_CHECK(fform_contains(fform, want_varname, msg), msg)
1026 : : end if
1027 : :
1028 : 3 : cplex_file = 1
1029 [ + - ]: 3 : if (ohdr%pertcase /= 0) then
1030 [ - - ]: 3 : cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1
1031 : : end if
1032 [ - + ]: 3 : ABI_CHECK(cplex_file == cplex, "cplex_file != cplex")
1033 : :
1034 : : ! Read FFT array (full box)
1035 [ + + ]: 12 : nfftot_file = product(ohdr%ngfft(:3))
1036 [ - + + - : 12 : ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
- + + - +
- + - - +
- + - + ]
1037 [ + + ]: 6 : do ispden=1,ohdr%nspden
1038 [ + - ]: 6 : read(unt, err=10, iomsg=errmsg) (rhor_file(ifft,ispden), ifft=1,cplex*nfftot_file)
1039 : : end do
1040 : 3 : close(unt)
1041 : :
1042 : : case (IO_MODE_ETSF)
1043 [ - + ]: 1987 : NCF_CHECK(nctk_open_read(unt, my_fname, xmpi_comm_self))
1044 : 1987 : call ohdr%ncread(unt, fform)
1045 : :
1046 : : ! Check important dimensions.
1047 [ - + ]: 1987 : ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", my_fname))
1048 : : !if (fform /= fform_den) then
1049 : : ! write(msg, "(2a, 2(a, i0))")' File: ',trim(my_fname),' is not a density file: fform= ',fform,", expecting:", fform_den
1050 : : ! ABI_WARNING(msg)
1051 : : !end if
1052 : :
1053 : 1987 : cplex_file = 1
1054 [ + + ]: 1987 : if (ohdr%pertcase /= 0) then
1055 [ + + ]: 1987 : cplex_file = 2; if (ohdr%qptn(1)**2 + ohdr%qptn(2)**2 + ohdr%qptn(3)**2 <1.d-14) cplex_file= 1
1056 : : end if
1057 [ - + ]: 1987 : ABI_CHECK(cplex_file == cplex, "cplex_file != cplex")
1058 : :
1059 : : ! Read FFT array (full box)
1060 [ + + ]: 7948 : nfftot_file = product(ohdr%ngfft(:3))
1061 [ - + + - : 7948 : ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
- + + - +
- + - - +
- + - + ]
1062 : :
1063 : 1987 : varname = varname_from_fname(my_fname)
1064 : : ncerr= nf90_get_var(unt, nctk_idname(unt, varname), rhor_file, &
1065 [ + + ]: 11922 : count=[cplex, ohdr%ngfft(1), ohdr%ngfft(2), ohdr%ngfft(3), ohdr%nspden])
1066 [ - + ]: 1987 : NCF_CHECK(ncerr)
1067 [ - + ]: 1987 : NCF_CHECK(nf90_close(unt))
1068 : :
1069 : : case default
1070 [ + + - ]: 1990 : ABI_ERROR(sjoin("Wrong iomode:", itoa(iomode)))
1071 : : end select
1072 : :
1073 [ + + + + ]: 7909 : need_interp = any(ohdr%ngfft(1:3) /= ngfft(1:3))
1074 [ + + ]: 1990 : if (need_interp) then
1075 : : msg = sjoin("Different FFT meshes. Caller expects:", ltoa(ngfft(1:3)), &
1076 [ + - ]: 17 : ". File: ", ltoa(ohdr%ngfft(1:3)), ". Need to perform interpolation.")
1077 : 17 : ABI_COMMENT(msg)
1078 [ - + ]: 17 : if (.not. allow_interp__) then
1079 : : write(msg, "(5a)") &
1080 : 0 : " Cannot continue as allow_interp = .False. ", ch10, &
1081 [ # # ]: 0 : " Please set ngfft to: ", trim(ltoa(ohdr%ngfft(1:3))), " in the input file"
1082 : 0 : ABI_ERROR(msg)
1083 : : end if
1084 : :
1085 [ + + - + : 119 : ABI_MALLOC(rhor_tmp, (cplex*product(ngfft(1:3)), ohdr%nspden))
+ - - + +
- + - + -
- + - + ]
1086 : 17 : call timab(1281,1,tsec)
1087 : 17 : call interpolate_denpot(cplex, ohdr%ngfft(1:3), ohdr%nspden, rhor_file, ngfft(1:3), rhor_tmp)
1088 : 17 : call timab(1281,2,tsec)
1089 : :
1090 [ + + ]: 68 : ohdr%ngfft(1:3) = ngfft(1:3)
1091 [ + + ]: 68 : nfftot_file = product(ohdr%ngfft(:3))
1092 [ - + ]: 17 : ABI_FREE(rhor_file)
1093 [ - + + - : 68 : ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
- + + - +
- + - - +
- + ]
1094 [ + - + - : 223947 : rhor_file = rhor_tmp
- + - - +
+ + + ]
1095 [ - + ]: 17 : ABI_FREE(rhor_tmp)
1096 : :
1097 : : ! Renormalize charge to avoid errors due to the interpolation.
1098 : : ! Do this only for NC since for PAW we should add the onsite contribution.
1099 : : ! This is left to the caller.
1100 : : !if (ohdr%usepaw == 0) then
1101 [ + + + + ]: 17 : if (ohdr%usepaw == 0 .and. fform == fform_den) then
1102 : 4 : call metric(gmet, gprimd, -1, rmet, ohdr%rprimd, ucvol)
1103 [ + + + + ]: 46808 : ratio = ohdr%nelect / (sum(rhor_file(:,1))*ucvol/ product(ngfft(1:3)))
1104 [ + + + + ]: 46800 : rhor_file = rhor_file * ratio
1105 : 4 : write(msg,'(a,f8.2,a,f8.4)')' Expected nelect: ',ohdr%nelect,' renormalization ratio: ',ratio
1106 : 4 : call wrtout(std_out,msg)
1107 : : end if
1108 : : end if ! need_interp
1109 : :
1110 : : ! Read PAW Rhoij
1111 [ + + ]: 1990 : if (ohdr%usepaw == 1) then
1112 [ - + + - : 2149 : ABI_MALLOC(pawrhoij_file, (ohdr%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ ]
1113 : 427 : call pawrhoij_nullify(pawrhoij_file)
1114 : : call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex_rhoij, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, &
1115 : 427 : ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size, qphase=ohdr%pawrhoij(1)%qphase)
1116 : 427 : call pawrhoij_copy(ohdr%pawrhoij, pawrhoij_file)
1117 : : end if
1118 : :
1119 : : ! Check that restart is possible !
1120 : : ! This check must be done here because we may have changed hdr% if need_interp
1121 [ + + ]: 1990 : if (present(check_hdr)) then
1122 : : ! FIXME: Temporary hack: fform_den to make hdr_check happy!
1123 : 1958 : call hdr_check(fform_den, fform_den, check_hdr, ohdr, "COLL", restart, restartpaw)
1124 : : !call hdr_check(fform_den, fform, check_hdr, ohdr, "COLL", restart, restartpaw)
1125 : : end if
1126 : :
1127 : : end if ! master
1128 : :
1129 [ + + ]: 2186 : if (nprocs == 1) then
1130 [ + - ]: 1898 : if (ohdr%nspden == nspden) then
1131 [ + + + + ]: 21141320 : orhor = rhor_file
1132 : : else
1133 : 0 : call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform)
1134 : : end if
1135 [ + + ]: 1898 : if (pawread == 1) call pawrhoij_copy(pawrhoij_file, pawrhoij, keep_nspden=.true.)
1136 : :
1137 : : else
1138 : 288 : call ohdr%bcast(master, my_rank, comm)
1139 : 288 : call xmpi_bcast(fform, master, comm, ierr)
1140 : :
1141 : : ! Eventually copy (or distribute) PAW data
1142 [ + + + + ]: 288 : if (ohdr%usepaw == 1 .and. pawread == 1) then
1143 [ + + ]: 20 : if (my_rank /= master) then
1144 [ - + + - : 63 : ABI_MALLOC(pawrhoij_file, (ohdr%natom))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ ]
1145 : 15 : call pawrhoij_nullify(pawrhoij_file)
1146 : : call pawrhoij_alloc(pawrhoij_file, ohdr%pawrhoij(1)%cplex_rhoij, ohdr%pawrhoij(1)%nspden, ohdr%pawrhoij(1)%nspinor, &
1147 : 15 : ohdr%pawrhoij(1)%nsppol, ohdr%typat, lmnsize=ohdr%lmn_size, qphase=ohdr%pawrhoij(1)%qphase)
1148 : : end if
1149 [ + + ]: 20 : if (size(ohdr%pawrhoij) /= size(pawrhoij)) then
1150 : : call pawrhoij_copy(ohdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
1151 : 4 : keep_nspden=.true.)
1152 : : else
1153 : 16 : call pawrhoij_copy(ohdr%pawrhoij,pawrhoij, keep_nspden=.true.)
1154 : : end if
1155 : : end if
1156 : :
1157 [ + + ]: 288 : if (my_rank /= master) then
1158 [ + + ]: 784 : nfftot_file = product(ohdr%ngfft(1:3))
1159 [ - + + - : 784 : ABI_MALLOC(rhor_file, (cplex*nfftot_file, ohdr%nspden))
- + + - +
- + - - +
- + - + ]
1160 : : end if
1161 : 288 : call xmpi_bcast(rhor_file, master, comm,ierr)
1162 : :
1163 [ + + ]: 288 : if (have_mpifft) then
1164 : : ! Extract slice treated by this MPI-FFT process.
1165 : 24 : call ptabs_fourdp(mpi_enreg, ngfft(2), ngfft(3), fftn2_distrib, ffti2_local, fftn3_distrib, ffti3_local)
1166 [ + - ]: 24 : if (ohdr%nspden==nspden) then
1167 [ + + ]: 48 : do ispden=1,nspden
1168 [ + + ]: 744 : do i3=1,n3
1169 [ + + ]: 696 : if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle
1170 : 268 : i3_local = ffti3_local(i3)
1171 [ + + ]: 8524 : do i2=1,n2
1172 : 8232 : mybase = cplex * (n1 * (i2-1 + n2*(i3_local-1)))
1173 : 8232 : globase = cplex * (n1 * (i2-1 + n2*(i3-1)))
1174 [ + + ]: 293104 : do i1=1,n1*cplex
1175 : 292408 : orhor(i1+mybase,ispden) = rhor_file(i1+globase,ispden)
1176 : : end do
1177 : : end do
1178 : : end do
1179 : : end do
1180 : : else
1181 [ # # ]: 0 : do i3=1,n3
1182 [ # # ]: 0 : if (fftn3_distrib(i3) /= mpi_enreg%me_fft) cycle
1183 : 0 : i3_local = ffti3_local(i3)
1184 [ # # ]: 0 : do i2=1,n2
1185 : 0 : mybase = 1 + cplex * (n1 * (i2-1 + n2*(i3_local-1)))
1186 : 0 : globase = 1 + cplex * (n1 * (i2-1 + n2*(i3-1)))
1187 : : call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform,&
1188 : 0 : istart_in=globase,istart_out=mybase,nelem=n1*cplex)
1189 : : end do
1190 : : end do
1191 : : end if
1192 : : else
1193 [ + - ]: 264 : if (ohdr%nspden==nspden) then
1194 [ + + + + ]: 2998846 : orhor = rhor_file
1195 : : else
1196 : 0 : call denpot_spin_convert(rhor_file,ohdr%nspden,orhor,nspden,fform)
1197 : : end if
1198 : : end if
1199 : : end if ! nprocs > 1
1200 : :
1201 [ - + ]: 2186 : ABI_FREE(rhor_file)
1202 : :
1203 [ + + ]: 2186 : if (allocated(pawrhoij_file)) then
1204 : 442 : call pawrhoij_free(pawrhoij_file)
1205 [ + - + + : 1328 : ABI_FREE(pawrhoij_file)
- + - + -
+ - + - +
- + - + ]
1206 : : end if
1207 : :
1208 : : ! Non-collinear magnetism: avoid zero magnetization, because it produces numerical instabilities
1209 : : ! Add a small real to the magnetization
1210 [ + + + + ]: 16330 : if (nspden==4) orhor(:,4)=orhor(:,4)+tol14
1211 [ + + + + ]: 2186 : if (ohdr%usepaw==1.and.size(pawrhoij)>0) then
1212 [ + + ]: 421 : if (pawrhoij(1)%nspden==4) then
1213 [ + + ]: 5 : do i1=1,size(pawrhoij)
1214 [ + + ]: 711 : pawrhoij(i1)%rhoijp(:,4)=pawrhoij(i1)%rhoijp(:,4)+tol10
1215 : : end do
1216 : : end if
1217 : : end if
1218 : :
1219 : 2186 : call timab(1280,2,tsec)
1220 : 2186 : call cwtime_report(" read_rhor", cputime, walltime, gflops)
1221 : 2186 : return
1222 : :
1223 : : ! Handle Fortran IO error
1224 : : 10 continue
1225 : 0 : ABI_ERROR(errmsg)
1226 : :
1227 [ + + - + : 4382 : end subroutine read_rhor
- + - + -
- - - - -
- - - - -
- - - ]
1228 : : !!***
1229 : :
1230 : : !----------------------------------------------------------------------
1231 : :
1232 : : !!****f* m_ioarr/fort_denpot_skip
1233 : : !! NAME
1234 : : !! fort_denpot_skip
1235 : : !!
1236 : : !! FUNCTION
1237 : : !! Skip the header and the DEN/POT records. Mainly used to append data to a pre-existent file.
1238 : : !! Return exit code.
1239 : : !!
1240 : : !! INPUTS
1241 : : !! unit=Fortran unit number (already opened in the caller).
1242 : : !! msg=Error message if ierr /= 0
1243 : : !!
1244 : : !! SOURCE
1245 : :
1246 : 0 : integer function fort_denpot_skip(unit, msg) result(ierr)
1247 : :
1248 : : !Arguments ------------------------------------
1249 : : integer,intent(in) :: unit
1250 : : character(len=*),intent(out) :: msg
1251 : :
1252 : : !Local variables-------------------------------
1253 : : integer :: ii,fform,nspden
1254 : 0 : type(hdr_type) :: hdr
1255 : :
1256 : : ! *********************************************************************
1257 : :
1258 : 0 : ierr = 1
1259 : 0 : call hdr%fort_read(unit, fform)
1260 [ # # ]: 0 : if (fform == 0) then
1261 [ # # # # ]: 0 : msg = "hdr_fort_read returned fform == 0"; return
1262 : : end if
1263 : :
1264 : 0 : nspden = hdr%nspden
1265 : 0 : call hdr%free()
1266 : :
1267 : : ! Skip the records with v1.
1268 [ # # ]: 0 : do ii=1,nspden
1269 : 0 : read(unit, iostat=ierr, iomsg=msg)
1270 [ # # ]: 0 : if (ierr /= 0) return
1271 : : end do
1272 : :
1273 : 0 : ierr = 0
1274 : :
1275 [ # # # # : 0 : end function fort_denpot_skip
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # ]
1276 : : !!***
1277 : :
1278 : : !----------------------------------------------------------------------
1279 : :
1280 : : !!****f* m_ioarr/denpot_spin_convert
1281 : : !! NAME
1282 : : !! denpot_spin_convert
1283 : : !!
1284 : : !! FUNCTION
1285 : : !! Convert a density/potential from a spin representation to another
1286 : : !!
1287 : : !! INPUTS
1288 : : !! denpot_in(:,nspden_in)=input density//potential
1289 : : !! nspden_in=number of spin-component of the input density/potential
1290 : : !! fform=file format (density or potential)
1291 : : !! [istart_in]= --optional-- starting index in the denpot_in array; default is 1
1292 : : !! [istart_out]= --optional-- starting index in the denpot_out array; default is 1
1293 : : !! [nelem]= --optional-- number of elements to copy from denpot_in to denpot_out; default is all
1294 : : !!
1295 : : !! OUTPUT
1296 : : !! denpot_out(:,nspden_out)=output density//potential
1297 : : !! nspden_out=number of spin-component of the output density/potential
1298 : : !!
1299 : : !! NOTES
1300 : : !! More explicitely:
1301 : : !! We copy denpot_in(istar_in+1:istart_in+nelem,:)
1302 : : !! into denpot_out(istart_out+1:istart_out+nelem,:)
1303 : : !!
1304 : : !! SOURCE
1305 : :
1306 [ # # # # ]: 0 : subroutine denpot_spin_convert(denpot_in,nspden_in,denpot_out,nspden_out,fform,&
1307 : : & istart_in,istart_out,nelem) ! optional arguments
1308 : :
1309 : : !Arguments ------------------------------------
1310 : : !scalars
1311 : : integer,intent(in) :: nspden_in,nspden_out,fform
1312 : : integer,intent(in),optional :: istart_in,istart_out,nelem
1313 : : !arrays
1314 : : real(dp),intent(in) :: denpot_in(:,:)
1315 : : real(dp),intent(out) :: denpot_out(:,:)
1316 : :
1317 : : !Local variables-------------------------------
1318 : : integer :: iend_in,iend_out,ispden,my_istart_in,my_istart_out,my_nelem
1319 : : character(len=500) :: msg
1320 : :
1321 : : ! *********************************************************************
1322 : :
1323 : : !Optional arguments
1324 [ # # ]: 0 : my_istart_in=1;if (present(istart_in)) my_istart_in=istart_in
1325 [ # # ]: 0 : my_istart_out=1;if (present(istart_out)) my_istart_out=istart_out
1326 : 0 : iend_in=size(denpot_in,1) ; iend_out=size(denpot_out,1)
1327 : 0 : my_nelem=min(iend_in-my_istart_in+1,iend_out-my_istart_out+1)
1328 [ # # ]: 0 : if (present(nelem)) my_nelem=nelem
1329 : :
1330 : : !Checks
1331 [ # # ]: 0 : if (size(denpot_in,2)/=nspden_in) then
1332 : 0 : msg='size(denpot_in,2)/=nspden_in!'
1333 : 0 : ABI_BUG(msg)
1334 : : end if
1335 [ # # ]: 0 : if (size(denpot_out,2)/=nspden_out) then
1336 : 0 : msg='size(denpot_out,2)/=nspden_out!'
1337 : 0 : ABI_BUG(msg)
1338 : : end if
1339 [ # # ]: 0 : if (my_istart_in+my_nelem-1>size(denpot_in,1)) then
1340 : 0 : msg='istart_in+nelem>size(denpot_in,1)!'
1341 : 0 : ABI_BUG(msg)
1342 : : end if
1343 [ # # ]: 0 : if (my_istart_out+my_nelem-1>size(denpot_out,1)) then
1344 : 0 : msg='istart_out+nelem>size(denpot_out,1)!'
1345 : 0 : ABI_BUG(msg)
1346 : : end if
1347 : :
1348 : : !Simple copy if the number of spin-components is unchanged...
1349 [ # # ]: 0 : if (nspden_in==nspden_out) then
1350 [ # # ]: 0 : do ispden=1,nspden_in
1351 : : denpot_out(my_istart_out:my_istart_out+my_nelem-1,ispden)= &
1352 [ # # ]: 0 : & denpot_in(my_istart_in:my_istart_in+my_nelem-1,ispden)
1353 : : end do
1354 : 0 : return
1355 : : end if
1356 : :
1357 : : !...otherwise, we need to convert.
1358 [ # # ]: 0 : if ((fform-1)/2==25) then
1359 : :
1360 : : ! First case: DENSITY
1361 : :
1362 [ # # # # ]: 0 : if (nspden_in==1.and.nspden_out==2) then
1363 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1364 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half
1365 [ # # # # ]: 0 : else if (nspden_in==1.and.nspden_out==4) then
1366 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1367 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero
1368 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1369 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1370 [ # # # # ]: 0 : else if (nspden_in==2.and.nspden_out==1) then
1371 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1372 [ # # # # ]: 0 : else if (nspden_in==2.and.nspden_out==4) then
1373 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1374 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=zero
1375 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1376 : : denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*two &
1377 [ # # ]: 0 : & -denpot_in(my_istart_in:my_istart_in+my_nelem,1)
1378 [ # # # # ]: 0 : else if (nspden_in==4.and.nspden_out==1) then
1379 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1380 [ # # # # ]: 0 : else if (nspden_in==4.and.nspden_out==2) then
1381 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1382 : : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half &
1383 [ # # ]: 0 : & +denpot_in(my_istart_in:my_istart_in+my_nelem-1,4)*half
1384 : : end if
1385 : :
1386 : : else
1387 : :
1388 : : ! Second case: POTENTIAL
1389 : :
1390 [ # # # # ]: 0 : if (nspden_in==1.and.nspden_out==2) then
1391 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1392 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1393 [ # # # # ]: 0 : else if (nspden_in==1.and.nspden_out==4) then
1394 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1395 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1396 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1397 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1398 [ # # # # ]: 0 : else if (nspden_in==2.and.nspden_out==1) then
1399 : : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half &
1400 [ # # ]: 0 : & +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half
1401 [ # # # # ]: 0 : else if (nspden_in==2.and.nspden_out==4) then
1402 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1403 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)
1404 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,3)=zero
1405 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,4)=zero
1406 [ # # # # ]: 0 : else if (nspden_in==4.and.nspden_out==1) then
1407 : : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)*half &
1408 [ # # ]: 0 : & +denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)*half
1409 [ # # # # ]: 0 : else if (nspden_in==4.and.nspden_out==2) then
1410 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,1)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,1)
1411 [ # # ]: 0 : denpot_out(my_istart_out:my_istart_out+my_nelem-1,2)=denpot_in(my_istart_in:my_istart_in+my_nelem-1,2)
1412 : : end if
1413 : :
1414 : : end if
1415 : :
1416 : : end subroutine denpot_spin_convert
1417 : : !!***
1418 : :
1419 : : end module m_ioarr
1420 : : !!***
|
