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 7834 : 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 7834 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:),fftn3_distrib(:),ffti3_local(:)
696 7834 : integer(XMPI_OFFSET_KIND) :: bsize_frecord(nspden)
697 : ! *************************************************************************
698 :
699 7834 : abifile = abifile_from_varname(varname)
700 7834 : 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 7834 : fform = abifile%fform
705 :
706 7834 : comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
707 7834 : 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 7834 : my_iomode = iomode
714 7834 : if (my_iomode /= IO_MODE_ETSF .and. nproc_fft == 1) my_iomode = IO_MODE_FORTRAN
715 7834 : if (nproc_fft > 1 .and. my_iomode == IO_MODE_FORTRAN) my_iomode = IO_MODE_MPI
716 :
717 7834 : call wrtout(std_out, sjoin(ch10, "fftdatar_write: About to write data to:", path, "with iomode:",iomode2str(my_iomode)))
718 7834 : call cwtime(cputime, walltime, gflops, "start")
719 :
720 : ! Get MPI-FFT tables from input ngfft
721 7834 : 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 7672 : file_etsf = nctk_ncify(path)
802 :
803 : ! Write datar.
804 : ncerr = nctk_write_datar(varname,file_etsf,ngfft,cplex,nfft,nspden, &
805 7672 : comm_fft,fftn3_distrib,ffti3_local,datar,action="create")
806 7672 : NCF_CHECK(ncerr)
807 7672 : call xmpi_barrier(comm_fft)
808 :
809 : ! Master writes the header.
810 7672 : if (xmpi_comm_rank(comm_fft) == master) then
811 7515 : NCF_CHECK(nctk_open_modify(ncid, file_etsf, xmpi_comm_self))
812 7515 : NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.))
813 : ! Add information on the crystalline structure.
814 7515 : NCF_CHECK(crystal%ncwrite(ncid))
815 7515 : if (present(ebands)) then
816 3673 : 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 7515 : NCF_CHECK(nf90_close(ncid))
825 : end if
826 :
827 : case default
828 7968 : ABI_ERROR(sjoin("Wrong iomode:",itoa(my_iomode)))
829 : end select
830 :
831 7834 : 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 7834 : 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 3877 : 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 3877 : type(ebands_t) :: ebands
882 : !arrays
883 3877 : real(dp),allocatable :: ene3d(:,:,:)
884 : ! *************************************************************************
885 :
886 3877 : crystal = hdr%get_crystal()
887 :
888 3877 : 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 3841 : call fftdatar_write(varname,path,iomode,hdr,crystal,ngfft,cplex,nfft,nspden,datar,mpi_enreg)
900 : end if
901 :
902 3877 : call crystal%free()
903 :
904 3877 : 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 4386 : 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 : !!***
|
