Line data Source code
1 : !!****m* ABINIT/m_iowf
2 : !! NAME
3 : !! m_iowf
4 : !!
5 : !! FUNCTION
6 : !! Procedures for the IO of the WFK file.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, AR, MB, MVer, MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : MODULE m_iowf
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_errors
28 : use m_dtset
29 : use m_xmpi
30 : use m_wffile
31 : use m_abi_etsf
32 : use m_nctk
33 : use m_wfk
34 : use netcdf
35 : use m_hdr
36 : use m_ebands
37 :
38 : use m_time, only : cwtime, cwtime_report, timab
39 : use m_io_tools, only : get_unit, flush_unit, iomode2str
40 : use m_fstrings, only : endswith, sjoin
41 : use m_numeric_tools, only : mask2blocks
42 : use defs_abitypes, only : MPI_type
43 : use defs_datatypes, only : pseudopotential_type
44 : use m_cgtools, only : cg_zcopy
45 : use m_crystal, only : crystal_t
46 : use m_rwwf, only : rwwf
47 : use m_mpinfo, only : proc_distrb_cycle
48 : use m_vkbr, only : calc_vkb
49 : use m_wvl_rwwf, only : wvl_write
50 :
51 : implicit none
52 :
53 : private
54 :
55 : public :: outwf
56 : public :: outresid
57 : public :: prtkbff ! Write KB form factors to WFK in netcdf format.
58 :
59 : !!***
60 :
61 : CONTAINS !====================================================================================================
62 : !!***
63 :
64 : !!****f* m_iowf/outresid
65 : !! NAME
66 : !! outresid
67 : !!
68 : !! FUNCTION
69 : !! - Compute the maximal residual and eventually print it
70 : !!
71 : !! INPUTS
72 : !! dtset <type(dataset_type)>=all input variables for this dataset
73 : !! kptns(3,nkpt)=k points in terms of recip primitive translations
74 : !! mband=maximum number of bands
75 : !! nband=number of bands
76 : !! nkpt=number of k points
77 : !! nsppol=1 for unpolarized, 2 for spin-polarized
78 : !! resid(mband*nkpt*nsppol)=squared residuals for each band and k point
79 : !! where resid(n,k)=|<C(n,k)|(H-e(n,k))|C(n,k)>|^2 for the ground state
80 : !!
81 : !! OUTPUT
82 : !! (only writing)
83 : !!
84 : !! SOURCE
85 :
86 7984 : subroutine outresid(dtset,kptns,mband,nband,nkpt,nsppol,resid)
87 :
88 : !Arguments ------------------------------------
89 : !scalars
90 : integer,intent(in) :: mband,nkpt,nsppol
91 : type(dataset_type),intent(in) :: dtset
92 : !arrays
93 : integer, intent(in) :: nband(nkpt*nsppol)
94 : real(dp), intent(in) :: kptns(3,nkpt)
95 : real(dp), intent(in) :: resid(mband*nkpt*nsppol)
96 :
97 : !Local variables-------------------------------
98 : integer,parameter :: nkpt_max=50
99 : integer :: band_index,spin,ikpt,ibdkpt,nkpt_eff,nband_k,nband_eff
100 : integer :: iband, ii
101 : real(dp) :: resims,residm, residk
102 : character(len=500) :: msg
103 :
104 : !Compute mean square and maximum residual over all bands and k points and spins
105 : !(disregard k point weights and occupation numbers here)
106 :
107 : !Find largest residual over bands, k points, and spins, except for nbdbuf highest bands
108 : !Already AVAILABLE in hdr ?!
109 7984 : ibdkpt=0
110 7984 : residm=zero
111 7984 : resims=zero
112 7984 : band_index=0
113 16676 : do spin=1,nsppol
114 239605 : do ikpt=1,nkpt
115 222929 : nband_k=nband(ikpt+(spin-1)*nkpt)
116 222929 : if (dtset%nbdbuf>0) then
117 70931 : nband_eff=max(1,nband_k-dtset%nbdbuf)
118 : else
119 : nband_eff=nband_k
120 : end if
121 1574851 : residm=max(residm,maxval(resid(ibdkpt+1:ibdkpt+nband_eff)))
122 1351922 : resims=resims + sum(resid(ibdkpt+1:ibdkpt+nband_eff))
123 222929 : ibdkpt=ibdkpt+nband_k
124 231621 : band_index=band_index + nband_eff
125 : end do
126 : end do
127 7984 : resims=resims/dble(band_index)
128 :
129 7984 : write(msg,'(a,2p,e12.4,a,e12.4)')' Mean square residual over all n,k,spin= ',resims,'; max=',residm
130 23952 : call wrtout([std_out, ab_out], msg)
131 :
132 7984 : band_index=0
133 7984 : nkpt_eff=nkpt
134 7984 : if( (dtset%prtvol==0 .or. dtset%prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max
135 :
136 : !Loop over spin again
137 16676 : do spin=1,nsppol
138 : ! Give (squared) residuals for all bands at each k
139 239605 : do ikpt=1,nkpt
140 222929 : nband_k=nband(ikpt+(spin-1)*nkpt)
141 : ! Will not print all residuals when prtvol=0 or 1
142 222929 : if(ikpt<=nkpt_eff)then
143 : ! Find largest residual over all bands for given k point
144 1059729 : residk=maxval(resid(1+band_index:nband_k+band_index))
145 128430 : write(msg,'(1x,3f8.4,3x,i2,1p,e13.5,a)')kptns(1:3,ikpt),spin,residk,' kpt; spin; max resid(k); each band:'
146 128430 : if(dtset%prtvol>=2) call wrtout(ab_out, msg)
147 128430 : call wrtout(std_out, msg)
148 295846 : do ii=0,(nband_k-1)/8
149 970285 : write(msg,'(1x,1p,8e9.2)')(resid(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8))
150 167416 : if(dtset%prtvol>=2) call wrtout(ab_out, msg)
151 295846 : call wrtout(std_out, msg)
152 : end do
153 94499 : else if(ikpt==nkpt_eff+1)then
154 758 : write(msg,'(2a)')' outresid : prtvol=0 or 1, do not print more k-points.',ch10
155 758 : if(dtset%prtvol>=2) call wrtout(ab_out, msg)
156 758 : call wrtout(std_out, msg)
157 : end if
158 231621 : band_index=band_index+nband_k
159 : end do
160 : end do
161 :
162 7984 : end subroutine outresid
163 : !!***
164 :
165 : !!****f* m_iowf/outwf
166 : !! NAME
167 : !! outwf
168 : !!
169 : !! FUNCTION
170 : !! Conduct output of a "wave-functions" file.
171 : !! - Compute the maximal residual
172 : !! - Then open a permanent file wff2 for final output of wf data
173 : !! - Create a new header for the file.
174 : !! - Write wave-functions (and energies)
175 : !!
176 : !! INPUTS
177 : !! cg(2,mcg)=wavefunction array (storage if nkpt>1)
178 : !! dtset <type(dataset_type)>=all input variables for this dataset
179 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
180 : !! eigen( (2*mband)**response *mband*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
181 : !! filnam= character string giving the root to form the name of the
182 : !! output WFK or WFQ file if response==0, otherwise it is the filename.
183 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
184 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
185 : !! kptns(3,nkpt)=k points in terms of recip primitive translations
186 : !! mband=maximum number of bands
187 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
188 : !! mkmem=Number of k-points treated by this node.
189 : !! mpi_enreg=information about MPI parallelization
190 : !! mpw=maximum number of plane waves
191 : !! natom=number of atoms in unit cell
192 : !! nband=number of bands
193 : !! nkpt=number of k points
194 : !! npwarr(nkpt)=number of planewaves in basis and on boundary for each k
195 : !! nsppol=1 for unpolarized, 2 for spin-polarized
196 : !! nstep=desired number of electron iteration steps
197 : !! occ(mband*nkpt*nsppol)=occupations for all bands at each k point
198 : !! response: if == 0, GS wavefunctions , if == 1, RF wavefunctions
199 : !! unwff2=unit for output of wavefunction
200 : !! wfs <type(wvl_projector_type)>=wavefunctions information for wavelets.
201 : !!
202 : !! OUTPUT
203 : !! (only writing)
204 : !!
205 : !! NOTES
206 : !! * The name of the file wff2 might be the same as that of the file wff1.
207 : !!
208 : !! SOURCE
209 :
210 5691 : subroutine outwf(cg,dtset,psps,eigen,filnam,hdr,kg,kptns,mband,mcg,mkmem,&
211 5691 : mpi_enreg,mpw,natom,nband,nkpt,npwarr,&
212 5691 : nsppol,occ,response,unwff2,&
213 : wfs,wvl,force_write) ! optional
214 :
215 : !Arguments ------------------------------------
216 : !scalars
217 : integer,intent(in) :: mband,mcg,mkmem,mpw,natom,nkpt,nsppol,response,unwff2
218 : !integer,intent(in) :: nstep
219 : character(len=*),intent(in) :: filnam
220 : type(MPI_type),intent(in) :: mpi_enreg
221 : type(dataset_type),intent(in) :: dtset
222 : type(pseudopotential_type),intent(in) :: psps
223 : type(hdr_type), intent(inout) :: hdr
224 : type(wvl_wf_type),intent(in) :: wfs
225 : type(wvl_internal_type), intent(in) :: wvl
226 : logical, intent(in), optional :: force_write
227 : !arrays
228 : integer, intent(in) :: kg(3,mpw*mkmem),nband(nkpt*nsppol),npwarr(nkpt)
229 : real(dp), intent(inout) :: cg(2,mcg)
230 : real(dp), intent(in) :: eigen((2*mband)**response*mband*nkpt*nsppol),kptns(3,nkpt)
231 : real(dp), intent(in) :: occ(mband*nkpt*nsppol)
232 :
233 : !Local variables-------------------------------
234 : integer :: iomode,action,band_index,fform,formeig,iband,icg !,iat,iproj
235 : integer :: ierr,ikg,ikpt,spin,master,mcg_disk,me,me0,mtag,my_nspinor
236 : integer :: nband_k,nmaster,npw_k,option,rdwr,sender,source !npwtot_k,
237 : integer :: spaceComm,spaceComm_io,spacecomsender,spaceWorld,sread,sskip,tim_rwwf,xfdim2
238 : #ifdef HAVE_MPI
239 : integer :: ipwnbd
240 : #endif
241 : real(dp) :: cpu,wall,gflops
242 : logical :: ihave_data,iwrite,iam_master,done,prtwf
243 : character(len=500) :: msg
244 : type(wffile_type) :: wff2
245 : !character(len=fnlen) :: path
246 : !arrays
247 5691 : integer,allocatable :: kg_disk(:,:)
248 : real(dp) :: tsec(2)
249 5691 : real(dp),allocatable :: cg_disk(:,:),eig_k(:),occ_k(:)
250 :
251 : ! *************************************************************************
252 : !For readability of the source file, define a "me" variable also in the sequential case
253 :
254 : DBG_ENTER("COLL")
255 :
256 : ABI_UNUSED(kptns(1,1))
257 :
258 5691 : xfdim2 = natom+4
259 : !Init mpi_comm
260 5691 : spaceWorld= mpi_enreg%comm_cell
261 5691 : spaceComm=spaceWorld
262 5691 : spaceComm_io=xmpi_comm_self
263 :
264 5691 : if (mpi_enreg%paral_kgb==1 ) spaceComm_io= mpi_enreg%comm_bandspinorfft
265 : if (mpi_enreg%paral_kgb==1 ) spaceComm= mpi_enreg%comm_cell
266 :
267 : !Paral_kgb=1 and Fortran-I/O is not supported (only for testing purpose)
268 5691 : if (mpi_enreg%paral_kgb==1.and.dtset%iomode==IO_MODE_FORTRAN) then
269 132 : spaceWorld=mpi_enreg%comm_kpt
270 : write(msg,'(7a)') &
271 132 : 'WF file is written using standard Fortran I/O',ch10,&
272 132 : 'and Kpt-band-FFT parallelization is active !',ch10,&
273 132 : 'This is only allowed for testing purposes.',ch10,&
274 264 : 'The produced WF file will be incomplete and not useable.'
275 132 : ABI_WARNING(msg)
276 : end if
277 :
278 : !If parallel HF calculation
279 5691 : if (mpi_enreg%paral_hf==1 ) spaceComm_io= mpi_enreg%comm_hf
280 5691 : if (mpi_enreg%paral_hf==1 ) spaceComm= mpi_enreg%comm_cell
281 :
282 : !Paral_hf=1 and Fortran-I/O is not supported (copy from paral_kgb... not tested)
283 5691 : if (mpi_enreg%paral_hf==1.and.dtset%iomode==IO_MODE_FORTRAN) then
284 0 : spaceWorld=mpi_enreg%comm_kpt
285 : write(msg,'(7a)') &
286 0 : 'WF file is written using standard Fortran I/O',ch10,&
287 0 : 'and HF parallelization is active !',ch10,&
288 0 : 'This is only allowed for testing purposes.',ch10,&
289 0 : 'The produced WF file will be incomplete and not useable.'
290 0 : ABI_WARNING(msg)
291 : end if
292 :
293 : ! check consistency between dimensions and input hdr.
294 : !ABI_CHECK(mband == maxval(hdr%nband), "hdr:mband")
295 : !ABI_CHECK(nkpt == hdr%nkpt, "hdr:nkpt")
296 : !ABI_CHECK(nsppol == hdr%nsppol, "hdr:nsppol")
297 : !ABI_CHECK(all(hdr%npwarr == npwarr), "hdr:npwarr")
298 : !ABI_CHECK(all(hdr%nband == nband), "hdr:nband")
299 : !ABI_CHECK(maxval(hdr%npwarr) == mpw, "hdr:nband")
300 :
301 : !Init me
302 5691 : me=mpi_enreg%me_kpt
303 5691 : me0=me
304 : !Define master
305 5691 : master=0
306 :
307 5691 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
308 5691 : tim_rwwf =0
309 5691 : source = master
310 5691 : sread = master
311 5691 : iam_master=(master==me)
312 5691 : iwrite=iam_master
313 5691 : sender=-1
314 :
315 : !Will write the wavefunction file only when nstep>0
316 : !MT 07 2015: writing reactivated when nstep=0
317 : !if (nstep>0 .and. dtset%prtwf/=0) then
318 : !FB 03/2022: Added an option to force writing (used in RT-TDDFT)
319 5691 : prtwf = dtset%prtwf/=0
320 5691 : if (present(force_write)) then
321 49 : if (force_write) prtwf = .true.
322 : end if
323 5642 : if (prtwf) then
324 :
325 : ! Only the master write the file, except if MPI I/O, but the
326 : ! full wff dataset should be provided to WffOpen in this case
327 3218 : iomode=IO_MODE_FORTRAN_MASTER
328 3218 : if (dtset%iomode==IO_MODE_MPI) iomode = IO_MODE_MPI
329 3218 : if (dtset%iomode==IO_MODE_ETSF) iomode = IO_MODE_ETSF
330 :
331 3218 : if (dtset%iomode == IO_MODE_ETSF .and. dtset%usewvl == 0) then
332 : call cg_ncwrite(filnam,hdr,dtset,response,mpw,mband,nband,nkpt,nsppol,&
333 3058 : dtset%nspinor,mcg,mkmem,eigen,occ,cg,npwarr,kg,mpi_enreg,done)
334 : #ifdef HAVE_NETCDF_DEFAULT
335 3058 : ABI_CHECK(done, "cg_ncwrite must handle the output of the WFK file.")
336 : #endif
337 :
338 : ! Write KB form factors. Only master works. G-vectors are read from file to avoid
339 : ! having to deal with paral_kgb distribution.
340 3058 : if (me == master .and. dtset%prtkbff == 1 .and. dtset%iomode == IO_MODE_ETSF .and. dtset%usepaw == 0) then
341 1 : ABI_CHECK(done, "cg_ncwrite was not able to generate WFK.nc in parallel. Perhaps hdf5 is not working")
342 1 : call prtkbff(filnam, hdr, psps, dtset%prtvol)
343 : end if
344 :
345 3058 : if (done) return
346 : ! If cg_ncwrite cannot handle the IO because HDF5 + MPI-IO support is missing, we fallback to Fortran + MPI-IO.
347 0 : msg = "Could not produce a netcdf file in parallel (MPI-IO support is missing). Will fallback to MPI-IO with Fortran"
348 0 : ABI_WARNING(msg)
349 0 : iomode=IO_MODE_MPI
350 : end if
351 :
352 160 : call cwtime(cpu, wall, gflops, "start")
353 160 : call wrtout(std_out, sjoin(ch10,' outwf: writing wavefunctions to:', trim(filnam), "with iomode:", iomode2str(iomode)))
354 :
355 : ! Create an ETSF file for the wavefunctions
356 160 : if (iomode == IO_MODE_ETSF) then
357 0 : ABI_CHECK(xmpi_comm_size(spaceComm) == 1, "Legacy etsf-io code does not support nprocs > 1")
358 0 : ABI_ERROR("ETSF_IO is not activated")
359 : ABI_UNUSED(psps%ntypat)
360 : end if
361 :
362 160 : call WffOpen(iomode,spaceComm,filnam,ierr,wff2,master,me0,unwff2,spaceComm_io)
363 : ! Conduct wavefunction output to wff2
364 :
365 480 : ABI_MALLOC(kg_disk,(3,mpw))
366 :
367 160 : mcg_disk=mpw*my_nspinor*mband
368 160 : formeig=0; if (response==1) formeig=1
369 :
370 480 : ABI_MALLOC(eig_k,( (2*mband)**formeig * mband))
371 480 : ABI_MALLOC(occ_k,(mband))
372 :
373 : #ifdef HAVE_MPI
374 160 : call xmpi_barrier(spaceComm)
375 : ! Compute mband and mpw
376 480 : ABI_MALLOC_OR_DIE(cg_disk,(2,mcg_disk), ierr)
377 : #endif
378 :
379 160 : band_index=0
380 160 : icg=0
381 160 : if(mpi_enreg%paralbd==0) tim_rwwf=6
382 160 : if(mpi_enreg%paralbd==1) tim_rwwf=12
383 :
384 : ! Write header info for new wf file
385 160 : rdwr=2
386 160 : if (dtset%usewvl==0) then
387 160 : fform=2
388 : else
389 0 : fform = 200 ! Use 200 as radical for naming file format used by wavelets.
390 : end if
391 :
392 160 : if (wff2%iomode < 2) then
393 160 : call hdr_io(fform,hdr,rdwr,wff2)
394 160 : call WffKg(wff2,1)
395 0 : else if (wff2%iomode==IO_MODE_ETSF .and. iam_master) then
396 0 : NCF_CHECK(hdr%ncwrite(wff2%unwff, fform, nc_define=.True.))
397 : end if
398 :
399 340 : do spin=1,nsppol
400 180 : ikg=0
401 :
402 932 : do ikpt=1,nkpt
403 592 : nband_k=nband(ikpt+(spin-1)*nkpt)
404 592 : npw_k=npwarr(ikpt)
405 :
406 : #ifdef HAVE_MPI
407 592 : if (dtset%usewvl == 0) then
408 592 : mtag=ikpt+(spin-1)*nkpt
409 592 : call xmpi_barrier(spaceWorld)
410 :
411 : ! Must transfer the wavefunctions to the master processor
412 : ! Separate sections for paralbd=1 or other values ; might be merged
413 592 : if(mpi_enreg%paralbd==0)then
414 592 : nmaster=0
415 9648 : source=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,spin))
416 592 : ihave_data=.false.
417 592 : if(source==me)ihave_data=.true.
418 416 : action=0
419 : ! I am the master node, and I have the data in cg or cg_disk
420 416 : if((iam_master).and.(ihave_data))action=1
421 : ! I am not the master, and I have the data => send to master
422 592 : if((.not.iam_master).and.(ihave_data))action=2
423 : ! I am the master, and I receive the data
424 592 : if((iam_master).and.(.not.ihave_data))action=3
425 :
426 : ! I have the data in cg or cg_disk ( MPI_IO case)
427 592 : if (iomode==IO_MODE_MPI) then
428 60 : action = 0
429 60 : sender=-1
430 60 : iwrite=.false.
431 60 : if (ihave_data)then
432 52 : action=1
433 52 : iwrite=.true.
434 52 : sender=me
435 : end if
436 : end if
437 :
438 : ! I am the master node, and I have the data in cg or cg_disk
439 : ! I have the data in cg or cg_disk ( MPI_IO case)
440 532 : if(action==1)then
441 : ! Copy from kg to kg_disk
442 182300 : kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
443 : ! Copy from cg to cg_disk
444 874956 : do ipwnbd=1,nband_k*npw_k*my_nspinor
445 874604 : cg_disk(1,ipwnbd)=cg(1,ipwnbd+icg)
446 874956 : cg_disk(2,ipwnbd)=cg(2,ipwnbd+icg)
447 : end do
448 : end if
449 :
450 : ! I am not the master, and I have the data => send to master
451 : ! I am the master, and I receive the data
452 592 : if ( action==2.or.action==3) then
453 : !write(std_out,*)npw_k,nband_k
454 128 : call timab(48,1,tsec)
455 128 : if(action==2)then
456 64 : call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_disk,nmaster,spaceWorld,2*mtag+1,ierr)
457 : call xmpi_exch(cg(:,icg+1:icg+nband_k*npw_k*my_nspinor),2*nband_k*npw_k*my_nspinor, &
458 64 : & source,cg_disk,nmaster,spaceWorld,2*mtag+2,ierr)
459 : else
460 64 : call xmpi_exch(kg_disk,3*npw_k,source,kg_disk,nmaster,spaceWorld,2*mtag+1,ierr)
461 64 : call xmpi_exch(cg_disk,2*nband_k*npw_k*my_nspinor,source,cg_disk,nmaster,spaceWorld,2*mtag+2,ierr)
462 : end if
463 128 : call timab(48,2,tsec)
464 : end if
465 :
466 :
467 0 : else if(mpi_enreg%paralbd==1)then
468 0 : nmaster=0
469 : #ifdef HAVE_MPI_IO
470 0 : sender=IO_MODE_FORTRAN_MASTER
471 0 : if( iomode==IO_MODE_MPI) then
472 0 : nmaster=mpi_enreg%proc_distrb(ikpt,1,spin)
473 0 : sender=nmaster
474 : end if
475 : #endif
476 :
477 : ! Note the loop over bands
478 0 : do iband=1,nband_k
479 :
480 : ! The message passing related to kg is counted as one band
481 0 : action=0
482 :
483 : ! I am the master node, and I have the data in cg or cg_disk
484 0 : if( mpi_enreg%proc_distrb(ikpt,iband,spin)==nmaster .and. me==nmaster) then
485 : action=1
486 : ! I am not the master, and I have the data => send to master
487 0 : elseif( mpi_enreg%proc_distrb(ikpt,iband,spin)==me .and. me/=nmaster ) then
488 : action = 2
489 : ! I am the master, and I receive the data
490 0 : elseif( mpi_enreg%proc_distrb(ikpt,iband,spin)/=me .and. me==nmaster ) then
491 : action=3
492 : end if
493 :
494 : if(action==1) then
495 : ! I am the master node, and I have the data in cg or cg_disk
496 : ! Copy from kg to kg_disk
497 0 : if(iband==1)kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
498 : ! Copy from cg to cg_disk
499 0 : do ipwnbd=1,npw_k*my_nspinor
500 0 : cg_disk(1,(iband-1)*npw_k*my_nspinor+ipwnbd) = cg(1,(iband-1)*npw_k*my_nspinor+ipwnbd+icg)
501 0 : cg_disk(2,(iband-1)*npw_k*my_nspinor+ipwnbd) = cg(2,(iband-1)*npw_k*my_nspinor+ipwnbd+icg)
502 : end do
503 : end if ! action=1
504 :
505 0 : if ( action==2.or.action==3) then
506 : ! action=2 : I am not the master, and I have the data => send to master
507 : ! action=3 : I am the master, and I receive the data
508 0 : call timab(48,1,tsec)
509 0 : if ( iband == 1 ) then
510 0 : if (action==2) then
511 : call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,spin), &
512 0 : & kg_disk,nmaster,spaceWorld,iband*(mtag-1)+1,ierr)
513 : else
514 : call xmpi_exch(kg_disk,3*npw_k,mpi_enreg%proc_distrb(ikpt,iband,spin), &
515 0 : & kg_disk,nmaster,spaceWorld,iband*(mtag-1)+1,ierr)
516 : end if
517 : end if ! iband =1
518 0 : ipwnbd=(iband-1)*npw_k*my_nspinor
519 0 : if (action==2) then
520 : call xmpi_exch( cg(:,ipwnbd+icg+1:ipwnbd+icg+npw_k*my_nspinor),2*npw_k*my_nspinor &
521 : & ,mpi_enreg%proc_distrb(ikpt,iband,spin) &
522 0 : & ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),nmaster,spaceWorld,iband*(mtag-1)+2,ierr)
523 : else
524 : call xmpi_exch( cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),2*npw_k*my_nspinor &
525 : & ,mpi_enreg%proc_distrb(ikpt,iband,spin) &
526 0 : & ,cg_disk(:,ipwnbd+1:ipwnbd+npw_k*my_nspinor),nmaster,spaceWorld,iband*(mtag-1)+2,ierr)
527 : end if
528 :
529 0 : call timab(48,2,tsec)
530 : end if ! action=2 or action=3
531 :
532 0 : if(iomode==IO_MODE_MPI) then
533 : ! I have the data in cg or cg_disk
534 0 : iwrite=.false.
535 0 : if (nmaster == me) iwrite=.true.
536 : end if
537 :
538 : end do ! End of loop over bands
539 : end if ! End of paralbd=1
540 : end if
541 : #endif
542 :
543 : ! Only the master will write to disk the final output wf file.
544 : ! in MPI_IO case only iwrite will write to disk the final output wf file.
545 592 : if(iwrite) then
546 : ! write(std_out,*) 'outwf : I am master and will write wf file'
547 416 : if(formeig==0)then
548 7056 : eig_k(1:nband_k)=eigen(1+band_index:nband_k+band_index)
549 7056 : occ_k(1:nband_k)=occ(1+band_index:nband_k+band_index)
550 : else
551 0 : eig_k(1:2*nband_k*nband_k)=eigen(1+band_index:2*nband_k*nband_k+band_index)
552 : end if
553 416 : option=2
554 416 : if(dtset%prtwf==3)option=5
555 : ! if (dtset%prtwf == 2 .and. mkmem/=0) option=4
556 :
557 416 : if (dtset%usewvl == 0) then
558 : #ifdef HAVE_MPI
559 : call rwwf(cg_disk,eig_k,formeig,0,0,ikpt,spin,kg_disk,mband,mcg_disk,mpi_enreg, &
560 416 : & nband_k, nband_k,npw_k,my_nspinor,occ_k,option,1,tim_rwwf,wff2)
561 :
562 : #else
563 : kg_disk(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
564 : call rwwf(cg,eig_k,formeig,0,icg,ikpt,spin,kg_disk,mband,mcg,mpi_enreg,nband_k, &
565 : & nband_k, npw_k,my_nspinor,occ_k,option,1,tim_rwwf,wff2)
566 : #endif
567 : else
568 0 : call wvl_write(dtset,eigen,mpi_enreg,option,hdr%rprimd,wff2,wfs,wvl,hdr%xred)
569 : end if
570 : end if
571 :
572 : ! The wavefunctions for the present k point and spin are written
573 592 : if(response==0)band_index=band_index+nband_k
574 592 : if(response==1)band_index=band_index+2*nband_k*nband_k
575 :
576 592 : sskip=1
577 : #ifdef HAVE_MPI
578 592 : if (dtset%usewvl == 0) then
579 592 : sskip=0
580 592 : if(.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me)))sskip=1
581 : end if
582 : #endif
583 : if(sskip==1)then
584 416 : icg=icg+npw_k*my_nspinor*nband_k
585 416 : ikg=ikg+npw_k
586 : end if
587 :
588 :
589 : #ifdef HAVE_MPI_IO
590 592 : spacecomsender=spaceComm
591 592 : if (mpi_enreg%paral_kgb==1) spacecomsender =mpi_enreg%comm_kpt
592 592 : if (mpi_enreg%paral_hf==1) spacecomsender =mpi_enreg%comm_kpt
593 772 : call WffOffset(wff2,sender,spacecomsender,ierr)
594 : #endif
595 :
596 : end do ! ikpt
597 : end do ! spin
598 :
599 160 : ABI_FREE(kg_disk)
600 : #ifdef HAVE_MPI
601 160 : ABI_FREE(cg_disk)
602 : #endif
603 :
604 160 : ABI_FREE(eig_k)
605 160 : ABI_FREE(occ_k)
606 :
607 : ! Close the wavefunction file (and do NOT delete it !)
608 : !if (wff2%iomode /= IO_MODE_NETCDF) then
609 160 : call WffClose(wff2,ierr)
610 : !end if
611 :
612 480 : call cwtime_report(" WFK output", cpu, wall, gflops)
613 : end if ! End condition of nstep>0
614 :
615 : ! Block here because we might need to read the WFK file in the caller.
616 2633 : call xmpi_barrier(mpi_enreg%comm_cell)
617 :
618 : DBG_EXIT("COLL")
619 :
620 2633 : end subroutine outwf
621 : !!***
622 :
623 : !----------------------------------------------------------------------
624 :
625 : !!****f* m_iowf/cg_ncwrite
626 : !! NAME
627 : !! cg_ncwrite
628 : !!
629 : !! FUNCTION
630 : !! Conduct output of a "wave-functions" file with netcdf
631 : !!
632 : !! INPUTS
633 : !! fname= character string giving the root to form the name of the
634 : !! output WFK or WFQ file if response==0, otherwise it is the filename.
635 : !! dtset <type(dataset_type)>=all input variables for this dataset
636 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
637 : !! response: if == 0, GS wavefunctions , if == 1, RF wavefunctions
638 : !! mpw=maximum number of plane waves
639 : !! mband=maximum number of bands
640 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
641 : !! mkmem=maximum number of k-points treated by this node
642 : !! nkpt=number of k points
643 : !! nsppol=1 for unpolarized, 2 for spin-polarized
644 : !! eigen((2*mband)**response *mband*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
645 : !! occ(mband*nkpt*nsppol)=occupations for all bands at each k point
646 : !! cg(2,mcg)=wavefunction array (storage if nkpt>1)
647 : !! npwarr(nkpt)=number of planewaves in basis and on boundary for each k
648 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
649 : !! mpi_enreg=information about MPI parallelization
650 : !!
651 : !! OUTPUT
652 : !! done=.True if cg_ncwrite can handle the output of the WFK file in parallel.
653 : !!
654 : !! SOURCE
655 :
656 3058 : subroutine cg_ncwrite(fname,hdr,dtset,response,mpw,mband,nband,nkpt,nsppol,nspinor,mcg,&
657 3058 : mkmem,eigen,occ,cg,npwarr,kg,mpi_enreg,done)
658 :
659 : !Arguments ------------------------------------
660 : !scalars
661 : integer,intent(in) :: response,mband,mcg,mkmem,mpw,nkpt,nsppol,nspinor
662 : character(len=*),intent(in) :: fname
663 : logical,intent(out) :: done
664 : type(dataset_type),intent(in) :: dtset
665 : type(MPI_type),intent(in) :: mpi_enreg
666 : type(hdr_type),intent(in) :: hdr
667 : !arrays
668 : integer, intent(in) :: nband(nkpt*nsppol),kg(3,mpw*mkmem),npwarr(nkpt)
669 : real(dp),intent(in) :: cg(2,mcg)
670 : real(dp),intent(in) :: eigen((2*mband)**response*mband*nkpt*nsppol),occ(mband*nkpt*nsppol)
671 :
672 : !Local variables-------------------------------
673 : !scalars
674 : integer,parameter :: master=0,fform2=2
675 : integer :: ii,iomode,icg,iband,ikg,ikpt,spin,me_cell,me_kpt,me_band,me_spinor,my_nspinor,nband_k,npw_k
676 : integer :: comm_cell,comm_fft,comm_bandfft,mtag,formeig
677 : integer :: cnt,min_cnt,max_cnt,ierr,action,source,ncid,ncerr,cg_varid,kg_varid !,eig_varid,
678 : integer :: paral_kgb,npwtot_k !,start_pwblock !,start_cgblock !count_pwblock,
679 : integer :: ipw,ispinor_index,npwso,npwsotot,npwtot,nspinortot,ikpt_this_proc,ispinor
680 : integer :: bandpp,nproc_band,nproc_fft,nproc_spinor,me_fft,nproc_cell,nwrites
681 : integer :: comm_mpiio,nranks,bstart,bcount !nbdblock,nblocks,
682 : !integer :: band_blocksize,band
683 : real(dp) :: cpu,wall,gflops
684 : logical :: ihave_data,iam_master,single_writer,same_layout,use_collective
685 : character(len=500) :: msg
686 : character(len=fnlen) :: path
687 12232 : type(wfk_t) :: wfk
688 : type(crystal_t) :: crystal
689 3058 : type(ebands_t) :: gs_ebands
690 : !arrays
691 3058 : integer,allocatable :: kg_k(:,:),iter2kscgkg(:,:),ind_cg_mpi_to_seq(:),rank_has_cg(:),ranks_io(:)!,gblock(:,:)
692 : real(dp) :: tsec(2)
693 3058 : real(dp),allocatable :: eigen3d(:,:,:),occ3d(:,:,:),cg_k(:,:),my_cgblock(:,:,:)
694 :
695 : ! *************************************************************************
696 :
697 : DBG_ENTER("COLL")
698 3058 : done = .False.
699 :
700 3058 : path = nctk_ncify(fname)
701 3058 : call wrtout(std_out, sjoin(" In cg_ncwrite with path:", path))
702 :
703 : ! communicators and ranks
704 3058 : comm_cell = mpi_enreg%comm_cell; me_cell = mpi_enreg%me_cell; nproc_cell = mpi_enreg%nproc_cell
705 3058 : comm_fft = mpi_enreg%comm_fft; me_fft = mpi_enreg%me_fft; nproc_fft = mpi_enreg%nproc_fft
706 3058 : comm_bandfft = mpi_enreg%comm_bandfft
707 3058 : me_kpt = mpi_enreg%me_kpt; me_band = mpi_enreg%me_band; me_spinor = mpi_enreg%me_spinor
708 3058 : iam_master = (me_kpt == master)
709 :
710 3058 : paral_kgb = dtset%paral_kgb
711 3058 : nproc_band = mpi_enreg%nproc_band
712 3058 : bandpp = mpi_enreg%bandpp
713 3058 : nproc_spinor = mpi_enreg%nproc_spinor
714 :
715 : ! FIXME
716 3058 : my_nspinor = max(1, nspinor/nproc_spinor)
717 3058 : if (nspinor == 2 .and. my_nspinor == 1) then
718 0 : ABI_ERROR("Spinor parallelization not coded yet")
719 : end if
720 :
721 : ! TODO: Be careful with response == 1 in parallel because the distribution of the cg
722 : ! can be **VERY** different from cg if nprocs > nkpt * nsppol
723 : !ABI_CHECK(response==0, "response == 1 not coded")
724 3058 : if (size(hdr%nband) == size(nband)) then
725 44696 : ABI_CHECK(all(Hdr%nband == nband),"nband")
726 : else
727 0 : ABI_ERROR("hdr%nband and nband have different size!")
728 : end if
729 :
730 3058 : if (xmpi_comm_size(comm_cell) == 1) then
731 32660 : ABI_CHECK(all(npwarr == hdr%npwarr), "npwarr != hdr%npwarr")
732 : end if
733 :
734 3058 : crystal = hdr%get_crystal()
735 :
736 : ! TODO
737 : ! Be careful with response == 1.
738 : ! gs_ebands contains the GS eigenvalues and occupation and will be written if this is a
739 : ! GS wfk. If we have a DFPT file, eigen stores the GKK matrix element, in this case
740 : ! we don't write gs_ebands but we define new variables in the netcdf file to store the GKK
741 15290 : ABI_MALLOC(occ3d, (mband,nkpt,nsppol))
742 3058 : call unpack_eneocc(nkpt,nsppol,mband,nband,occ,occ3d)
743 :
744 3058 : if (response == 0) then
745 3045 : formeig = 0
746 12180 : ABI_MALLOC(eigen3d, (mband,nkpt,nsppol))
747 3045 : call unpack_eneocc(nkpt,nsppol,mband,nband,eigen,eigen3d)
748 388462 : call gs_ebands%from_hdr(hdr, mband, eigen3d); gs_ebands%occ = occ3d
749 3045 : ABI_FREE(eigen3d)
750 : else
751 13 : formeig = 1
752 : end if
753 :
754 : ! same_layout is set to True if the internal representation of the cgs
755 : ! is compatible with the representation on file.
756 3058 : iomode = IO_MODE_ETSF
757 3058 : if (response == 0) then
758 5878 : same_layout = (paral_kgb == 0 .or. (paral_kgb == 1 .and. all([nproc_fft, nproc_band, nproc_spinor] == 1)))
759 : else
760 : ! For the time being, these cases are not implemented in the DFPT part.
761 13 : ABI_CHECK(nproc_fft==1, "nproc_fft != 1 not coded")
762 13 : ABI_CHECK(nproc_band==1, "nproc_band != 1 not coded")
763 13 : ABI_CHECK(nproc_spinor==1, "nproc_spinor != 1 not coded")
764 :
765 : ! Note: It would be possible to use collective IO if the cg1 are block-distributed
766 13 : same_layout = .True.
767 26 : spin_loop: do spin=1,nsppol
768 1722 : do ikpt=1,nkpt
769 1696 : nband_k = nband(ikpt + (spin-1)*nkpt)
770 1709 : if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then
771 8736 : if (any(mpi_enreg%proc_distrb(ikpt,:nband_k,spin) /= me_kpt)) then
772 0 : same_layout = .False.; exit spin_loop
773 : end if
774 : end if
775 : end do
776 : end do spin_loop
777 13 : call xmpi_land(same_layout, comm_cell)
778 : end if
779 :
780 3058 : call cwtime(cpu, wall, gflops, "start")
781 :
782 3058 : if (same_layout) then
783 2846 : single_writer = .True.
784 2846 : if (nctk_has_mpiio) single_writer = .False.
785 : !single_writer = .True.
786 :
787 : write(msg,'(5a,l1)')&
788 2846 : " same layout --> writing WFK file: ",trim(path),", with iomode: ",trim(iomode2str(iomode)),", single writer: ",single_writer
789 2846 : call wrtout(std_out,msg, 'PERS', do_flush=.True.)
790 :
791 2846 : if (.not. single_writer) then
792 :
793 : ! master opens the file and write the metadata.
794 2846 : if (xmpi_comm_rank(comm_cell) == master) then
795 : ncerr = nf90_einval
796 : #ifdef HAVE_NETCDF_MPI
797 : ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
798 2630 : comm=xmpi_comm_self, info=xmpio_info, ncid=ncid)
799 : #endif
800 2630 : NCF_CHECK_MSG(ncerr, sjoin("create_par: ", path))
801 :
802 2630 : call wfk_ncdef_dims_vars(ncid, hdr, fform2, write_hdr=.True.)
803 2630 : NCF_CHECK(crystal%ncwrite(ncid))
804 :
805 2630 : if (response == 0) then
806 : ! Write Gs bands
807 2617 : NCF_CHECK(gs_ebands%ncwrite(ncid))
808 : else
809 : ! Write H1 matrix elements and occupancies.
810 : ! Note that GS eigens are not written here.
811 13 : call ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d)
812 : end if
813 :
814 2630 : NCF_CHECK(nf90_close(ncid))
815 : end if
816 :
817 : ! Compute the table for collective IO.
818 207668 : ABI_CALLOC(iter2kscgkg, (4, nkpt*nsppol))
819 2846 : cnt = 0; icg = 0
820 5965 : do spin=1,nsppol
821 3119 : ikg = 0
822 45791 : do ikpt=1,nkpt
823 39826 : nband_k = nband(ikpt + (spin-1)*nkpt)
824 39826 : npw_k = npwarr(ikpt)
825 42945 : if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then
826 : ! FIXME v2[01], v3[6], with 4 procs
827 : ! v6[7] and v7[68], v7[69] fail but due to a extra line with nprocs
828 : ! v7[96] with np=4 seems to be more serious but it crashes also in trunk.
829 : ! v2[88] is slow likely due to outscfv
830 318040 : ABI_CHECK(all(mpi_enreg%proc_distrb(ikpt,:nband_k,spin) == me_kpt), "bands are distributed")
831 34044 : cnt = cnt + 1
832 170220 : iter2kscgkg(:,cnt) = [ikpt, spin, icg, ikg]
833 34044 : icg = icg + npw_k*my_nspinor*nband_k
834 34044 : ikg = ikg + npw_k
835 : end if
836 : end do
837 : end do
838 2846 : if (cnt == 0) then
839 36 : write(std_out,*)"cnt == 0 for me_cell, me, me_kpt",mpi_enreg%me_cell, mpi_enreg%me, mpi_enreg%me_kpt
840 : !ABI_CHECK(cnt > 0, "This processor does not have wavefunctions!")
841 : end if
842 :
843 2846 : call xmpi_min(cnt, min_cnt, comm_cell, ierr)
844 2846 : call xmpi_max(cnt, max_cnt, comm_cell, ierr)
845 :
846 : ! Handle idle procs, i.e. processors that do not have wavefunctions
847 : ! This happens if paral_kgb == 0 and nprocs > nkpt * nsppol (Abinit does not stop anymore!)
848 2846 : comm_mpiio = comm_cell
849 :
850 2846 : if (min_cnt <= 0) then
851 72 : ABI_COMMENT("Will create subcommunicator to exclude idle processors from MPI-IO collective calls")
852 72 : ABI_CHECK(paral_kgb == 0, "paral_kgb == 1 with idle processors should never happen")
853 :
854 : ! Prepare the call to xmpi_subcomm that will replace comm_mpiio.
855 496 : ABI_CALLOC(rank_has_cg, (0:nproc_cell-1))
856 160 : do spin=1,nsppol
857 258 : do ikpt=1,nkpt
858 134 : nband_k = nband(ikpt + (spin-1)*nkpt)
859 186 : if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) then
860 36 : rank_has_cg(me_kpt) = 1
861 36 : exit
862 : end if
863 : end do
864 : end do
865 :
866 72 : call xmpi_sum(rank_has_cg,comm_cell,ierr)
867 352 : nranks = count(rank_has_cg == 1)
868 216 : ABI_MALLOC(ranks_io, (nranks))
869 72 : cnt = 0
870 352 : do ii=0,nproc_cell-1
871 352 : if (rank_has_cg(ii) == 1) then
872 140 : cnt = cnt + 1
873 140 : ranks_io(cnt) = ii
874 : end if
875 : end do
876 : !write(std_out,*)"nranks, ranks_io:", nranks, ranks_io
877 72 : comm_mpiio = xmpi_subcomm(comm_cell, nranks, ranks_io)
878 72 : if (.not. rank_has_cg(me_kpt) == 1) then
879 36 : comm_mpiio = xmpi_comm_null
880 36 : ABI_CHECK(rank_has_cg(me_kpt) == 0, "rank_has_cg must be 0 or 1")
881 : end if
882 72 : ABI_FREE(ranks_io)
883 144 : ABI_FREE(rank_has_cg)
884 : end if
885 :
886 : ! Open the file in parallel inside comm_mpiio.
887 2846 : call xmpi_barrier(comm_cell)
888 2846 : if (comm_mpiio == xmpi_comm_null) goto 100
889 :
890 : ncerr = nf90_einval
891 : #ifdef HAVE_NETCDF_MPI
892 : ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write),&
893 2810 : comm=comm_mpiio, info=xmpio_info, ncid=ncid)
894 : #endif
895 2810 : NCF_CHECK_MSG(ncerr, sjoin("open_par: ", path))
896 :
897 : ! Use individual IO (default) for the G-vectors [3, mpw, nkpt]
898 2810 : NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
899 2810 : spin = 1; ikg=0
900 41441 : do ikpt=1,nkpt
901 38631 : npw_k = npwarr(ikpt)
902 38631 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,spin,me_kpt)) cycle
903 230636 : ncerr = nf90_put_var(ncid, kg_varid, kg(:, 1+ikg:), start=[1,1,ikpt], count=[3,npw_k,1])
904 32948 : NCF_CHECK_MSG(ncerr, "put_kg_k")
905 41441 : ikg = ikg + npw_k
906 : end do
907 :
908 2810 : NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", cg_varid))
909 :
910 2810 : use_collective = (response == 0) ! or (response == 1 .and. nproc_cell == 1)
911 :
912 2810 : if (use_collective) then
913 2797 : call wrtout(std_out,"Using collective IO for the CGs")
914 : ! Use collective IO for the CGs
915 : ncerr = nf90_einval
916 : #ifdef HAVE_NETCDF_MPI
917 2797 : ncerr = nf90_var_par_access(ncid, cg_varid, nf90_collective)
918 : #endif
919 2797 : NCF_CHECK(ncerr)
920 :
921 35175 : do cnt=1,max_cnt
922 32378 : ikpt = iter2kscgkg(1,cnt)
923 32378 : spin = iter2kscgkg(2,cnt)
924 32378 : icg = iter2kscgkg(3,cnt)
925 32378 : ikg = iter2kscgkg(4,cnt)
926 :
927 : ! The array on file has shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
928 35175 : if (ikpt /= 0) then
929 32348 : nband_k = nband(ikpt + (spin-1)*nkpt)
930 32348 : npw_k = npwarr(ikpt)
931 :
932 : !write(std_out,*)"About to write ikpt, spin, npw_k, icg: ",ikpt, spin, npw_k, icg
933 : ncerr = nf90_put_var(ncid, cg_varid, cg(:, 1+icg:), start=[1,1,1,1,ikpt,spin], &
934 420524 : count=[2,npw_k,nspinor,nband_k,1,1])
935 32348 : NCF_CHECK_MSG(ncerr, "writing cg")
936 : else
937 : ! This happens when nkpt * nsppol // nprocs != 0
938 : ! Note that we are using collective MPI-IO hence all processors must call put_var
939 : ! Here we re-write the ug(0) of the first (k-point, spin) treated by the node.
940 30 : ikpt = iter2kscgkg(1,1)
941 30 : spin = iter2kscgkg(2,1)
942 210 : ncerr = nf90_put_var(ncid, cg_varid, cg, start=[1,1,1,1,ikpt,spin], count=[1,1,1,1,1,1])
943 30 : NCF_CHECK_MSG(ncerr, "re-writing cg")
944 : end if
945 : end do
946 :
947 : else
948 13 : call wrtout(std_out, "Using individual IO for the CGs")
949 : ! Individual IO of the CGs (for debugging purposes)
950 13 : icg = 0
951 26 : do spin=1,nsppol
952 1722 : do ikpt=1,nkpt
953 1696 : nband_k = nband(ikpt + (spin-1)*nkpt)
954 1696 : npw_k = npwarr(ikpt)
955 1696 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) cycle
956 22048 : ncerr = nf90_put_var(ncid, cg_varid, cg(:, 1+icg:), start=[1,1,1,1,ikpt,spin], count=[2,npw_k,nspinor,nband_k,1,1])
957 1696 : NCF_CHECK_MSG(ncerr, "put_var")
958 1709 : icg = icg+npw_k*my_nspinor*nband_k
959 : end do
960 : end do
961 : end if
962 :
963 2810 : NCF_CHECK(nf90_close(ncid))
964 :
965 2846 : 100 call xmpi_barrier(comm_cell)
966 2846 : ABI_FREE(iter2kscgkg)
967 :
968 2846 : done = .True.
969 5692 : call cwtime_report(" collective ncwrite", cpu, wall, gflops)
970 :
971 : else ! single_writer
972 0 : if (nproc_cell > 1) then
973 0 : ABI_WARNING("Slow version without MPI-IO support. Processors send data to master...")
974 : else
975 0 : call wrtout(std_out, "Using netcdf library without MPI-IO support")
976 : end if
977 :
978 0 : ABI_MALLOC(kg_k,(3,mpw))
979 0 : ABI_MALLOC_OR_DIE(cg_k,(2,mpw*my_nspinor*mband), ierr)
980 :
981 0 : if (iam_master) then
982 0 : call wfk%open_write(hdr,path,formeig,iomode,get_unit(),xmpi_comm_self,write_hdr=.True.)
983 :
984 0 : NCF_CHECK(crystal%ncwrite(wfk%fh))
985 : !write(std_out,*)"after crystal_ncwrite"
986 :
987 : ! Write eigenvalues and occupations (these arrays are not MPI-distributed)
988 0 : if (response == 0) then
989 0 : NCF_CHECK(gs_ebands%ncwrite(wfk%fh))
990 : else
991 0 : call ncwrite_eigen1_occ(wfk%fh, nband, mband, nkpt, nsppol, eigen, occ3d)
992 : end if
993 : end if
994 :
995 0 : icg = 0
996 0 : do spin=1,nsppol
997 0 : ikg = 0
998 0 : do ikpt=1,nkpt
999 0 : nband_k = nband(ikpt + (spin-1)*nkpt)
1000 0 : npw_k = npwarr(ikpt)
1001 0 : mtag = ikpt+(spin-1)*nkpt
1002 :
1003 0 : call xmpi_barrier(comm_cell)
1004 :
1005 : ! Transfer the wavefunctions and the g-vectors to the master processor
1006 0 : source = minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,spin))
1007 0 : ihave_data = (source==me_kpt)
1008 :
1009 0 : action=0
1010 0 : if (iam_master .and. ihave_data) action=1 ! I am the master node, and I have the data in cg
1011 0 : if (.not.iam_master.and.ihave_data) action=2 ! I am not the master, and I have the data => send to master
1012 0 : if (iam_master.and..not.ihave_data) action=3 ! I am the master, and I receive the data
1013 :
1014 0 : if (action==1) then ! Copy from kg and cg
1015 0 : kg_k(:,1:npw_k) = kg(:,ikg+1:ikg+npw_k)
1016 0 : call cg_zcopy(npw_k*my_nspinor*nband_k, cg(1,icg+1), cg_k)
1017 : end if
1018 :
1019 : ! Exchange data
1020 0 : if (action==2.or.action==3) then
1021 0 : call timab(48,1,tsec)
1022 0 : if (action==2) then
1023 0 : call xmpi_exch(kg(:,1+ikg:npw_k+ikg),3*npw_k,source,kg_k,master,comm_cell,2*mtag+1,ierr)
1024 : call xmpi_exch(cg(:,icg+1:icg+nband_k*npw_k*my_nspinor),2*nband_k*npw_k*my_nspinor,&
1025 0 : source,cg_k,master,comm_cell,2*mtag+2,ierr)
1026 : else
1027 0 : call xmpi_exch(kg_k,3*npw_k,source,kg_k,master,comm_cell,2*mtag+1,ierr)
1028 0 : call xmpi_exch(cg_k,2*nband_k*npw_k*my_nspinor,source,cg_k,master,comm_cell,2*mtag+2,ierr)
1029 : end if
1030 0 : call timab(48,2,tsec)
1031 : end if
1032 :
1033 : ! Master writes this block of bands.
1034 0 : if (iam_master) then
1035 0 : if (response == 0) then
1036 : call wfk%write_band_block([1,nband_k],ikpt,spin,xmpio_single,kg_k=kg_k,cg_k=cg_k,&
1037 0 : eig_k=gs_ebands%eig(:,ikpt,spin),occ_k=gs_ebands%occ(:,ikpt,spin))
1038 : !write(std_out,*)"cg_k",cg_k(:,1:2)
1039 : else
1040 0 : call wfk%write_band_block([1,nband_k],ikpt,spin,xmpio_single,kg_k=kg_k,cg_k=cg_k)
1041 : end if
1042 : end if
1043 :
1044 0 : if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt))) then
1045 0 : icg = icg+npw_k*my_nspinor*nband_k
1046 0 : ikg = ikg+npw_k
1047 : end if
1048 :
1049 : end do !ikpt
1050 : end do !spin
1051 :
1052 0 : ABI_FREE(kg_k)
1053 0 : ABI_FREE(cg_k)
1054 :
1055 0 : if (iam_master) call wfk%close()
1056 0 : call xmpi_barrier(comm_cell)
1057 :
1058 0 : done = .True.
1059 0 : call cwtime_report(" individual ncwrite", cpu, wall, gflops)
1060 : end if
1061 :
1062 : else ! not same_layout
1063 :
1064 212 : if (nctk_has_mpiio) then
1065 : call wrtout(std_out, &
1066 212 : sjoin("scattered data. writing WFK file",trim(path),", with iomode: ",iomode2str(iomode)), 'PERS', do_flush=.True.)
1067 :
1068 : ! master write the metadata.
1069 212 : if (xmpi_comm_rank(comm_cell) == master) then
1070 : ncerr = nf90_einval
1071 : #ifdef HAVE_NETCDF_MPI
1072 : ncerr = nf90_create(path, cmode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
1073 47 : comm=xmpi_comm_self, info=xmpio_info, ncid=ncid)
1074 : #endif
1075 47 : NCF_CHECK_MSG(ncerr, sjoin("create_par:", path))
1076 :
1077 47 : call wfk_ncdef_dims_vars(ncid, hdr, fform2, write_hdr=.True.)
1078 47 : NCF_CHECK(crystal%ncwrite(ncid))
1079 :
1080 : ! Write eigenvalues and occupations (these arrays are not MPI-distributed)
1081 47 : if (response == 0) then
1082 47 : NCF_CHECK(gs_ebands%ncwrite(ncid))
1083 : else
1084 0 : call ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d)
1085 : end if
1086 :
1087 47 : NCF_CHECK(nf90_close(ncid))
1088 : end if
1089 :
1090 : ! Reopen the file inside comm_cell
1091 212 : call xmpi_barrier(comm_cell)
1092 : ncerr = nf90_einval
1093 : #ifdef HAVE_NETCDF_MPI
1094 : ncerr = nf90_open(path, mode=ior(ior(nf90_netcdf4, nf90_mpiio), nf90_write), &
1095 212 : comm=comm_cell, info=xmpio_info, ncid=ncid)
1096 : #endif
1097 212 : NCF_CHECK_MSG(ncerr, sjoin("create_par:", path))
1098 :
1099 : ! Get var ids
1100 212 : NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
1101 212 : NCF_CHECK(nf90_inq_varid(ncid, "coefficients_of_wavefunctions", cg_varid))
1102 :
1103 : ! Write the G-vectors
1104 212 : ikg = 0
1105 2016 : do ikpt=1,nkpt
1106 1804 : npw_k = npwarr(ikpt)
1107 1804 : npwtot_k = hdr%npwarr(ikpt)
1108 1804 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,1,me_kpt)) cycle
1109 1014 : if (me_spinor /= 0) cycle
1110 : !write(std_out,*)"In G-vector loop ",ikpt,", with me_cell me_kpt me_band, me_spinor ",me_cell,me_kpt,me_band,me_spinor
1111 : !write(std_out,*)" with ik, npw_k, npwtot_k: ",ikpt, npw_k, npwtot_k
1112 :
1113 3042 : ABI_MALLOC(ind_cg_mpi_to_seq, (npw_k))
1114 1014 : if (allocated(mpi_enreg%my_kgtab)) then
1115 1014 : ikpt_this_proc = mpi_enreg%my_kpttab(ikpt)
1116 152418 : ind_cg_mpi_to_seq = mpi_enreg%my_kgtab(1:npw_k,ikpt_this_proc)
1117 : else
1118 0 : ABI_CHECK(nproc_fft==1, "nproc_fft !=1 and my_kgtab not allocated")
1119 0 : ind_cg_mpi_to_seq(1:npw_k) = [(ipw, ipw=1,npw_k)]
1120 : end if
1121 :
1122 1996218 : ABI_CALLOC(kg_k, (3, npwtot_k))
1123 151404 : do ipw=1,npw_k
1124 602574 : kg_k(:, ind_cg_mpi_to_seq(ipw)) = kg(:,ikg+ipw)
1125 : end do
1126 1014 : call xmpi_sum_master(kg_k,master,comm_bandfft,ierr)
1127 1014 : if (xmpi_comm_rank(comm_bandfft) == master) then
1128 2401 : ncerr = nf90_put_var(ncid, kg_varid, kg_k, start=[1,1,ikpt], count=[3,npwtot_k,1])
1129 343 : NCF_CHECK_MSG(ncerr, "putting kg_k")
1130 : else
1131 2684 : ncerr = nf90_put_var(ncid, kg_varid, kg_k, start=[1,1,ikpt], count=[0,0,0])
1132 671 : NCF_CHECK_MSG(ncerr, "skiping kg_k")
1133 : end if
1134 1014 : ABI_FREE(kg_k)
1135 :
1136 : ! gblock contains block-distributed G-vectors inside comm_fft
1137 : !call kg2seqblocks(npwtot_k,npw_k,kg(:,ikg+1:),ind_cg_mpi_to_seq,comm_fft,start_pwblock,count_pwblock,gblock)
1138 : !write(std_out,*)"gblock(:,2)",gblock(:,2)
1139 : !ncerr = nf90_put_var(ncid, kg_varid, gblock, start=[1,start_pwblock,ikpt], count=[3,count_pwblock,1])
1140 : !NCF_CHECK_MSG(ncerr, "putting kg_k")
1141 : !ABI_FREE(gblock)
1142 :
1143 1014 : ABI_FREE(ind_cg_mpi_to_seq)
1144 :
1145 2016 : ikg = ikg+npw_k
1146 : end do
1147 :
1148 : ! Write wavefunctions
1149 212 : if (response == 1) then
1150 : ! The cg1 is allocated with size mcg1=mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol (see dfpt_looppert)
1151 : ! hence bands are not MPI distributed, this is the reason why we have to use the loop over bands
1152 : ! and the check on nwrites
1153 0 : icg = 0
1154 0 : do spin=1,nsppol
1155 0 : do ikpt=1,nkpt
1156 0 : nband_k = nband(ikpt + (spin-1)*nkpt)
1157 0 : npw_k = npwarr(ikpt)
1158 0 : npwtot_k = hdr%npwarr(ikpt)
1159 0 : nwrites = 0
1160 0 : do iband=1,nband_k
1161 0 : if (mpi_enreg%proc_distrb(ikpt,iband,spin)/=me_kpt) cycle
1162 : ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
1163 0 : nwrites = nwrites + 1
1164 0 : ii = 1 + (iband-1)*npw_k*my_nspinor + icg
1165 0 : ncerr = nf90_put_var(ncid, cg_varid, cg(1:,ii:), start=[1,1,1,iband,ikpt,spin], count=[2,npwtot_k,nspinor,1,1,1])
1166 0 : NCF_CHECK_MSG(ncerr, "put_var cg")
1167 : end do ! iband
1168 0 : if (nwrites /= 0) icg = icg + npw_k*my_nspinor*nband_k
1169 : end do !ikpt
1170 : end do !spin
1171 :
1172 : else
1173 212 : icg = 0
1174 432 : do spin=1,nsppol
1175 2244 : do ikpt=1,nkpt
1176 1812 : nband_k = nband(ikpt + (spin-1)*nkpt)
1177 1812 : npw_k = npwarr(ikpt)
1178 1812 : npwtot_k = hdr%npwarr(ikpt)
1179 :
1180 1812 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,spin,me_kpt)) cycle
1181 : !write(std_out,*)"In u(g)-vector loop ",ikpt,", with me_cell me_kpt me_band, me_spinor ",me_cell,me_kpt,me_band,me_spinor
1182 : !write(std_out,*)"nband_k, npw_k, npwtot_k: ",nband_k, npw_k, npwtot_k
1183 :
1184 : ! Part taken from writewf
1185 : ! Note that ind_cg_mpi_to_seq is wrong in nspinor > 1
1186 1020 : npwtot=npw_k; npwso=npw_k*nspinor
1187 1020 : npwsotot=npwso
1188 1020 : nspinortot = min(2,(1+mpi_enreg%paral_spinor)*nspinor)
1189 :
1190 3060 : ABI_MALLOC(ind_cg_mpi_to_seq, (npwso))
1191 1020 : if (allocated(mpi_enreg%my_kgtab)) then
1192 1020 : ikpt_this_proc=mpi_enreg%my_kpttab(ikpt)
1193 2040 : do ispinor=1,nspinor
1194 1020 : ispinor_index=ispinor
1195 1020 : if (nproc_spinor > 1) ispinor_index = mpi_enreg%me_spinor + 1
1196 : ind_cg_mpi_to_seq(1+npw_k*(ispinor-1):npw_k*ispinor)=npwtot*(ispinor_index-1) &
1197 152819 : & + mpi_enreg%my_kgtab(1:npw_k,ikpt_this_proc)
1198 : end do
1199 : else
1200 0 : ABI_CHECK(nproc_fft==1, "nproc_fft !=1 and my_kgtab not allocated")
1201 0 : ind_cg_mpi_to_seq(1:npwso) = [(ipw, ipw=1,npwso)]
1202 : end if
1203 :
1204 : ! TODO: Blocking + collective IO in comm_cell
1205 : !band_blocksize = nband_k/nbdblock
1206 : !ibandmin = 1
1207 : !step=min(ii,MAXBAND, nband_disk)
1208 : !do iblock=1,nband_disk/step+1
1209 : ! ibandmax = min(ibandmin+step-1, nband_disk)
1210 : ! nband_block = ibandmax-ibandmin+1
1211 : !gbase = icg + npw_k*my_nspinor*nband_k
1212 :
1213 : ! Each MPI proc write bcount bands with npwtot_k G-vectors starting from bstart.
1214 1020 : call cg2seqblocks(npwtot_k,npw_k,nband_k,cg(:,icg+1:),ind_cg_mpi_to_seq,comm_bandfft,bstart,bcount,my_cgblock)
1215 :
1216 : ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
1217 13260 : ncerr = nf90_put_var(ncid, cg_varid, my_cgblock, start=[1,1,1,bstart,ikpt,spin], count=[2,npwtot_k,nspinor,bcount,1,1])
1218 1020 : NCF_CHECK_MSG(ncerr, "put_var cg")
1219 :
1220 1020 : ABI_FREE(my_cgblock)
1221 1020 : ABI_FREE(ind_cg_mpi_to_seq)
1222 :
1223 2032 : icg = icg+npw_k*my_nspinor*nband_k
1224 : end do !ikpt
1225 : end do !spin
1226 : end if
1227 :
1228 212 : NCF_CHECK(nf90_close(ncid))
1229 212 : done = .True.
1230 :
1231 212 : call cwtime_report(" scattered ncwrite", cpu, wall, gflops)
1232 : end if !nctk_has_mpiio
1233 : end if
1234 :
1235 3058 : call crystal%free()
1236 3058 : if (response == 0) call gs_ebands%free()
1237 :
1238 3058 : ABI_FREE(occ3d)
1239 :
1240 : DBG_EXIT("COLL")
1241 :
1242 6116 : end subroutine cg_ncwrite
1243 : !!***
1244 :
1245 : !!****f* m_iowf/prtkbff
1246 : !! NAME
1247 : !! prtkbff
1248 : !!
1249 : !! FUNCTION
1250 : !! Write KB form factors to WFK in netcdf format.
1251 : !! Only master works. G-vectors are read from file to avoid
1252 : !! having to deal with paral_kgb distribution.
1253 : !!
1254 : !! INPUTS
1255 : !!
1256 : !! OUTPUT
1257 : !!
1258 : !! SOURCE
1259 :
1260 2 : subroutine prtkbff(filnam, hdr, psps, prtvol)
1261 :
1262 : !Arguments ------------------------------------
1263 : character(len=*),intent(in) :: filnam
1264 : type(hdr_type),intent(in) :: hdr
1265 : type(pseudopotential_type),intent(in) :: psps
1266 : integer,intent(in) :: prtvol
1267 :
1268 : !Local variables-------------------------------
1269 : !scalars
1270 : character(len=fnlen) :: path
1271 : character(len=500) :: msg
1272 : integer :: ncid, ncerr, kg_varid, mpw_disk, npwk_disk, ikpt, iat, iproj
1273 2 : integer,allocatable :: kg_disk(:,:)
1274 2 : real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:),vkbsign(:,:)
1275 104 : type(crystal_t) :: crystal
1276 :
1277 : ! *************************************************************************
1278 :
1279 2 : path = nctk_ncify(filnam)
1280 2 : call wrtout(std_out, sjoin("Writing KB form factors to:", path))
1281 2 : NCF_CHECK(nctk_open_modify(ncid, path, xmpi_comm_self))
1282 2 : NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
1283 268 : mpw_disk = maxval(hdr%npwarr)
1284 :
1285 : ! Dimensions needed by client code to allocate memory when reading.
1286 : ncerr = nctk_def_dims(ncid, [ &
1287 : nctkdim_t("mproj", psps%mproj), &
1288 : nctkdim_t("mpsang", psps%mpsang), &
1289 : nctkdim_t("mpssoang", psps%mpssoang), &
1290 : nctkdim_t("lnmax", psps%lnmax), &
1291 : nctkdim_t("lmnmax", psps%lmnmax) &
1292 12 : ])
1293 2 : NCF_CHECK(ncerr)
1294 :
1295 4 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "mpspso"])
1296 2 : NCF_CHECK(ncerr)
1297 :
1298 : ! Write indlmn table (needed to access vkb arrays)
1299 : ncerr = nctk_def_arrays(ncid, [ &
1300 : nctkarr_t("indlmn", "int", "six, lmnmax, number_of_atom_species"), &
1301 : nctkarr_t("vkbsign", "dp", "lnmax, number_of_atom_species"), &
1302 : nctkarr_t("vkb", "dp", "max_number_of_coefficients, lnmax, number_of_atom_species, number_of_kpoints"), &
1303 : nctkarr_t("vkbd", "dp", "max_number_of_coefficients, lnmax, number_of_atom_species, number_of_kpoints") &
1304 10 : ], defmode=.True.)
1305 2 : NCF_CHECK(ncerr)
1306 :
1307 : ! Switch to write mode.
1308 2 : NCF_CHECK(nctk_set_datamode(ncid))
1309 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "indlmn"), psps%indlmn))
1310 :
1311 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
1312 : "mpspso"], &
1313 6 : [psps%mpspso])
1314 2 : NCF_CHECK(ncerr)
1315 :
1316 : ! Calculate KB form factors and derivatives.
1317 : ! The arrays are allocated with lnmax to support pseudos with more than projector.
1318 : ! Note that lnmax takes into account lloc hence arrays are in packed form and should be
1319 : ! accessed with the indices provided by psps%indlmn.
1320 8 : ABI_MALLOC(vkbsign, (psps%lnmax, psps%ntypat))
1321 10 : ABI_MALLOC(vkb, (mpw_disk, psps%lnmax, psps%ntypat))
1322 8 : ABI_MALLOC(vkbd, (mpw_disk, psps%lnmax, psps%ntypat))
1323 6 : ABI_MALLOC(kg_disk, (3, mpw_disk))
1324 :
1325 2 : crystal = hdr%get_crystal()
1326 :
1327 : ! For each k-point: read full G-vector list from file, compute KB data and write to file.
1328 268 : do ikpt=1,hdr%nkpt
1329 266 : npwk_disk = hdr%npwarr(ikpt)
1330 1862 : NCF_CHECK(nf90_get_var(ncid, kg_varid, kg_disk, start=[1, 1, ikpt], count=[3, npwk_disk, 1]))
1331 3518682 : vkb = zero; vkbd = zero
1332 266 : call calc_vkb(crystal, psps, hdr%kptns(:, ikpt), npwk_disk, mpw_disk, kg_disk, vkbsign, vkb, vkbd)
1333 :
1334 266 : if (ikpt == 1) then
1335 : ! This for the automatic tests.
1336 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkbsign"), vkbsign))
1337 2 : if(prtvol >= 2) then
1338 1 : write(msg,'(a)') 'prtkbff: writing first and last G-components of the KB form factors'
1339 1 : call wrtout(ab_out, msg)
1340 3 : do iat=1,psps%ntypat
1341 2 : write(msg,'(a10,i5)') 'atype ',iat
1342 2 : call wrtout(ab_out, msg)
1343 5 : do iproj=1,psps%lnmax
1344 : write(msg,'(a10,i5,a,a10,e12.4,a,2(a10,2e12.4,a))') &
1345 2 : 'projector ', iproj,ch10, &
1346 2 : 'vkbsign ', vkbsign(iproj,iat), ch10, &
1347 2 : 'vkb ', vkb(1,iproj,iat), vkb(npwk_disk,iproj,iat), ch10, &
1348 4 : 'vkbd ', vkbd(1,iproj,iat), vkbd(npwk_disk,iproj,iat), ''
1349 4 : call wrtout(ab_out, msg)
1350 : end do
1351 : end do
1352 : end if
1353 : end if
1354 :
1355 1330 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkb"), vkb, start=[1, 1, 1, ikpt]))
1356 1332 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vkbd"), vkbd, start=[1, 1, 1, ikpt]))
1357 : end do
1358 2 : NCF_CHECK(nf90_close(ncid))
1359 :
1360 2 : ABI_FREE(kg_disk)
1361 2 : ABI_FREE(vkbsign)
1362 2 : ABI_FREE(vkb)
1363 2 : ABI_FREE(vkbd)
1364 2 : call crystal%free()
1365 :
1366 2 : end subroutine prtkbff
1367 : !!***
1368 :
1369 : !----------------------------------------------------------------------
1370 :
1371 : !!****f* m_iowf/ncwrite_eigen1_occ
1372 : !! NAME
1373 : !! ncwrite_eigen1_occ
1374 : !!
1375 : !! FUNCTION
1376 : !! Write the first order DFPT eigenvalues and the occupations.
1377 : !!
1378 : !! INPUTS
1379 : !! ncid=Netcdf file handler.
1380 : !! nband(nkpt*nsppol)=Number of bands.
1381 : !! mband=maximum number of bands
1382 : !! nkpt=Total number of k points
1383 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1384 : !! eigen((2*mband**2*nkpt*nsppol)= eigenvalues (hartree) for all bands at each k point
1385 : !! occ3d(mband*nkpt*nsppol)=occupations for all bands at each k point
1386 : !! Note that occ3d differes from the occ arrays used in the rest of the code.
1387 : !! occ3d is an array with constant stride `mband` whereas abinit (for reasons that are not clear to me)
1388 : !! packs the occupations in a 1d vector with a k-dependent separation (nband_k).
1389 : !! occ and occ3d differ only if nband is k-dependent but you should never assume this, hence
1390 : !! remember to *convert* occ into occ3d before calling this routine.
1391 : !!
1392 : !! SOURCE
1393 :
1394 13 : subroutine ncwrite_eigen1_occ(ncid, nband, mband, nkpt, nsppol, eigen, occ3d)
1395 :
1396 : !Arguments ------------------------------------
1397 : !scalars
1398 : integer,intent(in) :: ncid,mband,nkpt,nsppol
1399 : !arrays
1400 : integer, intent(in) :: nband(nkpt*nsppol)
1401 : real(dp),intent(in) :: eigen(2*mband**2*nkpt*nsppol),occ3d(mband,nkpt,nsppol)
1402 :
1403 : !Local variables-------------------------------
1404 : !scalars
1405 : integer :: idx,spin,ikpt,nband_k,ib2,ib1
1406 : integer :: ncerr,occ_varid,h1mat_varid,eigens_varid
1407 : !arrays
1408 13 : real(dp),allocatable :: h1mat(:,:,:,:,:), fake_eigens(:,:,:)
1409 :
1410 : ! *************************************************************************
1411 :
1412 : ! Declare h1 array with abinit conventions
1413 : ! Cannot use a 3D array since eigen1 are packed and one could have different number of bands
1414 : ! per k-points. (this is not an official etsf-io variable)!
1415 : ! elements in eigen are packed in the first positions.
1416 : ! Remember that eigen are not MPI-distributed so no communication is needed
1417 13 : idx=1
1418 99451 : ABI_CALLOC(h1mat, (2,mband,mband,nkpt,nsppol))
1419 26 : do spin=1,nsppol
1420 1722 : do ikpt=1,nkpt
1421 1696 : nband_k = nband(ikpt + (spin-1)*nkpt)
1422 8749 : do ib2=1,nband_k
1423 38944 : do ib1=1,nband_k
1424 90624 : h1mat(:,ib1,ib2,ikpt,spin) = eigen(idx:idx+1)
1425 37248 : idx=idx+2
1426 : end do
1427 : end do
1428 : end do
1429 : end do
1430 :
1431 13 : NCF_CHECK(nctk_set_defmode(ncid))
1432 :
1433 : ncerr = nctk_def_arrays(ncid, nctkarr_t('h1_matrix_elements', "dp", &
1434 13 : &"complex, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins"))
1435 13 : NCF_CHECK(ncerr)
1436 :
1437 13 : NCF_CHECK(nf90_inq_varid(ncid, "occupations", occ_varid))
1438 13 : NCF_CHECK(nf90_inq_varid(ncid, "h1_matrix_elements", h1mat_varid))
1439 13 : NCF_CHECK(nf90_inq_varid(ncid, "eigenvalues", eigens_varid))
1440 :
1441 : ! Write data
1442 13 : NCF_CHECK(nctk_set_datamode(ncid))
1443 13 : NCF_CHECK_MSG(nf90_put_var(ncid, occ_varid, occ3d), "putting occ3d")
1444 13 : NCF_CHECK_MSG(nf90_put_var(ncid, h1mat_varid, h1mat), "putting h1mat")
1445 :
1446 13 : ABI_FREE(h1mat)
1447 :
1448 : ! GS eigenvalues are set to zero.
1449 8814 : ABI_CALLOC(fake_eigens, (mband,nkpt,nsppol))
1450 13 : NCF_CHECK_MSG(nf90_put_var(ncid, eigens_varid, fake_eigens), "putting fake eigens")
1451 13 : ABI_FREE(fake_eigens)
1452 :
1453 13 : end subroutine ncwrite_eigen1_occ
1454 : !!***
1455 :
1456 : !----------------------------------------------------------------------
1457 :
1458 : !!****f* m_iowf/kg2seqblocks
1459 : !! NAME
1460 : !! kg2seqblocks
1461 : !!
1462 : !! FUNCTION
1463 : !!
1464 : !! INPUTS
1465 : !! npwtot_k
1466 : !! npw_k=
1467 : !! kg(3,npw_k)=reduced planewave coordinates.
1468 : !! gmpi2seq
1469 : !! comm_fft=FFT communicator
1470 : !!
1471 : !! OUTPUT
1472 : !! start_pwblock
1473 : !! count_pwblock
1474 : !! gblock(:,:)
1475 : !!
1476 : !! SOURCE
1477 :
1478 : subroutine kg2seqblocks(npwtot_k,npw_k,kg_k,gmpi2seq,comm_fft,start_pwblock,count_pwblock,gblock)
1479 :
1480 : !Arguments ------------------------------------
1481 : !scalars
1482 : integer,intent(in) :: npwtot_k,npw_k,comm_fft
1483 : integer,intent(out) :: start_pwblock,count_pwblock
1484 : !arrays
1485 : integer,intent(in) :: kg_k(3,npw_k),gmpi2seq(npw_k)
1486 : integer,allocatable,intent(out) :: gblock(:,:)
1487 :
1488 : !Local variables-------------------------------
1489 : !scalars
1490 : integer :: me_fft,nproc_fft,rank,ig,igseq,maxnpw,ierr
1491 : !arrays
1492 : integer,allocatable :: igstart_rank(:),gbuf(:,:)
1493 :
1494 : ! *************************************************************************
1495 :
1496 : me_fft = xmpi_comm_rank(comm_fft); nproc_fft = xmpi_comm_size(comm_fft)
1497 :
1498 : ! Handle sequential case.
1499 : if (nproc_fft == 1) then
1500 : start_pwblock = 1; count_pwblock = npw_k
1501 : ABI_MALLOC(gblock, (3, npw_k))
1502 : gblock(:,:) = kg_k
1503 : return ! DOH
1504 : end if
1505 :
1506 : ABI_MALLOC(igstart_rank, (0:nproc_fft))
1507 : igstart_rank = [(1 + (npwtot_k/nproc_fft)*rank, rank=0,nproc_fft-1), 1 + npwtot_k]
1508 :
1509 : ! Get max dimension for workspace array
1510 : ! Cannot use npwtot_k / nproc_fft because G-vectors are not equally distributed.
1511 : call xmpi_max(npw_k, maxnpw, comm_fft, ierr)
1512 : ABI_MALLOC(gbuf, (3, maxnpw))
1513 :
1514 : do rank=0,nproc_fft-1
1515 : start_pwblock = igstart_rank(rank)
1516 : count_pwblock = igstart_rank(rank+1) - igstart_rank(rank)
1517 :
1518 : gbuf = 0
1519 : do ig=1,npw_k
1520 : igseq = gmpi2seq(ig)
1521 : if (igseq >= igstart_rank(rank) .and. igseq < igstart_rank(rank+1)) then
1522 : igseq = igseq - start_pwblock + 1
1523 : !ABI_CHECK(igseq <= maxnpw, "boom")
1524 : gbuf(:,igseq) = kg_k(:,ig)
1525 : end if
1526 : end do
1527 : call xmpi_sum_master(gbuf,rank,comm_fft,ierr)
1528 :
1529 : if (me_fft == rank) then
1530 : ABI_MALLOC_OR_DIE(gblock, (3, count_pwblock), ierr)
1531 : gblock = gbuf(:, :count_pwblock)
1532 : end if
1533 : end do
1534 :
1535 : start_pwblock = igstart_rank(me_fft)
1536 : count_pwblock = igstart_rank(me_fft+1) - igstart_rank(me_fft)
1537 :
1538 : ABI_FREE(gbuf)
1539 : ABI_FREE(igstart_rank)
1540 :
1541 : end subroutine kg2seqblocks
1542 : !!***
1543 :
1544 : !----------------------------------------------------------------------
1545 :
1546 : !!****f* m_iowf/cg2seqblocks
1547 : !! NAME
1548 : !! cg2seqblocks
1549 : !!
1550 : !! FUNCTION
1551 : !!
1552 : !! INPUTS
1553 : !! npwtot_k
1554 : !! npw_k=
1555 : !! cg_k(2,npw_k*nband_k)
1556 : !! gmpi2seq
1557 : !! comm_fft=FFT communicator
1558 : !!
1559 : !! OUTPUT
1560 : !! bstart
1561 : !! bcount
1562 : !!
1563 : !! SOURCE
1564 :
1565 1020 : subroutine cg2seqblocks(npwtot_k,npw_k,nband,cg_k,gmpi2seq,comm_bandfft,bstart,bcount,my_cgblock)
1566 :
1567 : !Arguments ------------------------------------
1568 : !scalars
1569 : integer,intent(in) :: npwtot_k,npw_k,nband,comm_bandfft
1570 : integer,intent(out) :: bstart,bcount
1571 : !arrays
1572 : integer,intent(in) :: gmpi2seq(npw_k)
1573 : real(dp),intent(in) :: cg_k(2,npw_k,nband)
1574 : real(dp),allocatable,intent(out) :: my_cgblock(:,:,:)
1575 :
1576 : !Local variables-------------------------------
1577 : !scalars
1578 : integer :: me,nprocs,rank,ig,ib,igseq,ierr,nbb,nbmax,band
1579 : !arrays
1580 1020 : integer,allocatable :: bstart_rank(:)
1581 1020 : real(dp),allocatable :: cgbuf(:,:,:)
1582 :
1583 : ! *************************************************************************
1584 :
1585 1020 : me = xmpi_comm_rank(comm_bandfft); nprocs = xmpi_comm_size(comm_bandfft)
1586 :
1587 : ! Handle sequential case.
1588 1020 : if (nprocs == 1) then
1589 0 : bstart = 1; bcount = nband
1590 0 : ABI_MALLOC_OR_DIE(my_cgblock, (2, npwtot_k, nband), ierr)
1591 0 : my_cgblock = cg_k
1592 : return ! DOH
1593 : end if
1594 :
1595 3060 : ABI_MALLOC(bstart_rank, (0:nprocs))
1596 14160 : bstart_rank = [(1 + (nband/nprocs)*rank, rank=0,nprocs-1), 1 + nband]
1597 :
1598 : ! Allocate MPI buffer (same size on each MPI proc)
1599 1020 : nbmax = 0
1600 4380 : do rank=0,nprocs-1
1601 4380 : nbmax = max(nbmax, bstart_rank(rank+1) - bstart_rank(rank))
1602 : end do
1603 4080 : ABI_MALLOC_OR_DIE(cgbuf, (2, npwtot_k, nbmax), ierr)
1604 :
1605 1020 : nbb = bstart_rank(me+1) - bstart_rank(me)
1606 4080 : ABI_MALLOC_OR_DIE(my_cgblock, (2, npwtot_k, nbb), ierr)
1607 :
1608 : ! TODO: This should be replaced by gatherv but premature optimization....
1609 4380 : do rank=0,nprocs-1
1610 3360 : bstart = bstart_rank(rank)
1611 3360 : bcount = bstart_rank(rank+1) - bstart_rank(rank)
1612 :
1613 15840 : do band=bstart, bstart+bcount-1
1614 12480 : ib = band - bstart + 1
1615 28849152 : cgbuf(:,:,ib) = zero
1616 2784236 : do ig=1,npw_k
1617 2768396 : igseq = gmpi2seq(ig)
1618 8317668 : cgbuf(:,igseq,ib) = cg_k(:,ig,band)
1619 : end do
1620 : end do ! band
1621 :
1622 3360 : call xmpi_sum_master(cgbuf,rank,comm_bandfft,ierr)
1623 8314604 : if (me == rank) my_cgblock = cgbuf(:,:,:bcount)
1624 : end do ! rank
1625 :
1626 1020 : bstart = bstart_rank(me)
1627 1020 : bcount = bstart_rank(me+1) - bstart_rank(me)
1628 :
1629 1020 : ABI_FREE(cgbuf)
1630 1020 : ABI_FREE(bstart_rank)
1631 :
1632 : end subroutine cg2seqblocks
1633 : !!***
1634 :
1635 : !----------------------------------------------------------------------
1636 :
1637 : END MODULE m_iowf
1638 : !!***
|
