Line data Source code
1 : !!****m* ABINIT/m_inwffil
2 : !! NAME
3 : !! m_inwffil
4 : !!
5 : !! FUNCTION
6 : !! Initialization of wavefunctions.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, AR, MB, MVer, ZL, MB, TD, 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_inwffil
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_wffile
28 : use m_wfk
29 : use m_errors
30 : use m_xomp
31 : use m_xmpi
32 : use m_nctk
33 : use m_hdr
34 : use m_dtset
35 : USE_MPI
36 :
37 : use defs_abitypes, only : MPI_type
38 : use m_fstrings, only : sjoin, itoa
39 : use m_time, only : timab, cwtime, cwtime_report
40 : use m_io_tools, only : file_exists, get_unit
41 : use m_geometry, only : getspinrot
42 : use m_pptools, only : prmat
43 : use m_matrix, only : matr3inv, mati3inv
44 : use m_cgtools, only : cg_envlop, pw_orthon
45 : use m_fftcore, only : kpgsph, sphere, sphereboundary
46 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_io
47 : use m_mpinfo, only : destroy_mpi_enreg, copy_mpi_enreg, proc_distrb_cycle
48 : use m_kg, only : kpgio, ph1d3d, getph
49 : use m_kpts, only : listkk
50 : use m_rwwf, only : rwwf, WffReadSkipK
51 : use m_wvl_wfsinp, only : wvl_wfsinp_disk, wvl_wfsinp_scratch
52 :
53 : implicit none
54 :
55 : #if defined HAVE_MPI1
56 : include 'mpif.h'
57 : #endif
58 :
59 : private
60 : !!***
61 :
62 : public :: inwffil ! Do initialization of wavefunctions.
63 : public :: cg_from_atoms !Initialize wave functions using Bloch sums of atomic orbitals.
64 : !!***
65 :
66 : contains
67 : !!***
68 :
69 : !!****f* m_inwffil/inwffil
70 : !! NAME
71 : !! inwffil
72 : !!
73 : !! FUNCTION
74 : !! Do initialization of wavefunctions.
75 : !! Also call other relevant routines for this initialisation
76 : !! (initialization of wavefunctions from scratch or from file, translations of wavefunctions, ...)
77 : !!
78 : !! INPUTS
79 : !! ask_accurate= if 1, the wavefunctions and eigenvalues must be
80 : !! accurate, that is, they must come from a k point that is
81 : !! symmetric of the needed k point, with a very small tolerance,
82 : !! the disk file contained sufficient bands to initialize all of them,
83 : !! the spinor and spin-polarisation characteristics must be identical
84 : !! dtset <type(dataset_type)>=all input variables for this dataset
85 : !! ecut=effective kinetic energy planewave cutoff (hartree), beyond
86 : !! which the coefficients of plane waves are zero
87 : !! ecut_eff=effective kinetic energy planewave cutoff (hartree), needed
88 : !! to generate the sphere of plane wave
89 : !! exchn2n3d=if 1, n2 and n3 are exchanged
90 : !! formeig=explained above
91 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
92 : !! ireadwf=option parameter described above for wf initialization
93 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs to be initialized here.
94 : !! kg(3,mpw*my_nkpt)=dimensionless coords of G vecs in basis sphere at k point
95 : !! kptns(3,nkpt)=reduced coords of k points
96 : !! localrdwf=(for parallel case) if 1, the wffnm file is local to each machine
97 : !! mband=maximum number of bands
98 : !! mband_mem=maximum number of bands for this cpu
99 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband_mem*mkmem*nsppol
100 : !! mkmem=number of k-points in core memory
101 : !! mpi_enreg=information about MPI parallelization
102 : !! mpw=maximum number of planewaves as dimensioned in calling routine
103 : !! nband(nkpt*nsppol)=number of bands at each k point
104 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
105 : !! nkpt=number of k points
106 : !! npwarr(nkpt)=array holding npw for each k point.
107 : !! nsppol=1 for unpolarized, 2 for spin-polarized
108 : !! nsym=number of symmetry elements in space group
109 : !! occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
110 : !! optorth= 1 if the WFS have to be orthogonalized; 0 otherwise
111 : !! prtvol=control print volume and debugging
112 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
113 : !! symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
114 : !! tnons(3,nsym)=nonsymmorphic translations for symmetry operations
115 : !! unkg=unit number for storage of basis sphere data: stores indirect
116 : !! indexing array and integer coordinates for all planewaves in basis
117 : !! sphere for each k point being considered
118 : !! unwff1,unwfnow= unit numbers for files wffnm and wft1nm.
119 : !! wffnm=name (character data) of file for input wavefunctions.
120 : !!
121 : !! OUTPUT
122 : !! wff1 = structure information for files wffnm .
123 : !! wffnow= structure information for wf file wft1nm
124 : !! if ground state format (formeig=0):
125 : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
126 : !! if respfn format (formeig=1):
127 : !! eigen(2*mband*mband*nkpt*nsppol)=matrix of eigenvalues
128 : !! (input or init to large number), (Ha)
129 : !! Conditional output (returned if mkmem/=0):
130 : !! cg(2,mcg)=complex wf array
131 : !! be careful : an array of size cg(2,npw*nspinor), as used
132 : !! in the response function code, is not enough !
133 : !! wvl <type(wvl_data)>=all wavelets data.
134 : !!
135 : !! NOTES
136 : !! Detailed description:
137 : !! Initialize unit wff1%unwff for input of wf data if ireadwf=1
138 : !! Opens file on unit wffnow%unwff
139 : !! if the storage on disk is needed (mkmem==0)
140 : !! Initializes wf data on wffnow%unwff, by calling the appropriate routine.
141 : !!
142 : !! formeig option (format of the eigenvalues and occupations) :
143 : !! 0 => ground-state format (initialisation of
144 : !! eigenvectors with random numbers, vector of eigenvalues,
145 : !! occupations are present)
146 : !! 1 => respfn format (initialisation of
147 : !! eigenvectors with 0 s, hermitian matrix of eigenvalues)
148 : !!
149 : !! ireadwf options:
150 : !! 0 => initialize with random numbers or 0 s
151 : !! 1 => read from disk file wff1, initializing higher bands
152 : !! with random numbers or 0 s if not provided in disk file
153 : !!
154 : !! The wavefunctions after this initialisation are stored in unit wffnow%unwff
155 : !!
156 : !! WARNINGS
157 : !!
158 : !! * The symmetry operations are used to translate the data from one
159 : !! k point to another, symmetric, k point.
160 : !! They can be completely different from the symmetry operations
161 : !! contained on the disk file. No check is performed between the two sets.
162 : !!
163 : !! * Occupations will not be modified nor output, in the present status of this routine.
164 : !!
165 : !! * If ground state format (formeig=0) occ(mband*nkpt*nsppol) was output.
166 : !! NOT OUTPUT NOW!
167 : !!
168 : !! SOURCE
169 :
170 6619 : subroutine inwffil(ask_accurate,cg,dtset,ecut,ecut_eff,eigen,exchn2n3d,&
171 6619 : & formeig,hdr,ireadwf,istwfk,kg,kptns,localrdwf,mband,&
172 6619 : & mcg,mkmem,mpi_enreg,mpw,nband,ngfft,nkpt,npwarr,&
173 6619 : & nsppol,nsym,occ,optorth,symafm,symrel,tnons,unkg,wff1,&
174 : & wffnow,unwff1,wffnm,wvl)
175 :
176 : !Arguments ------------------------------------
177 : integer,intent(in) :: ask_accurate,exchn2n3d,formeig,ireadwf,localrdwf,mband,mcg,mkmem,mpw
178 : integer,intent(in) :: nkpt,nsppol,nsym,optorth,unkg,unwff1
179 : real(dp),intent(in) :: ecut,ecut_eff
180 : character(len=*),intent(in) :: wffnm
181 : type(MPI_type),intent(inout),target :: mpi_enreg
182 : type(dataset_type),intent(in) :: dtset
183 : type(hdr_type),intent(inout) :: hdr
184 : type(wffile_type),intent(inout) :: wff1
185 : type(wffile_type),intent(inout) :: wffnow
186 : type(wvl_data),intent(inout) :: wvl
187 : integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),ngfft(18)
188 : integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym)
189 : integer,intent(in),target :: nband(nkpt*nsppol)
190 : real(dp),intent(inout),target :: cg(2,mcg),eigen((2*mband)**formeig*mband*nkpt*nsppol)
191 : real(dp),intent(in) :: kptns(3,nkpt),tnons(3,nsym)
192 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
193 :
194 : !Local variables-------------------------------
195 : integer,parameter :: master=0
196 : integer :: iomode,accurate,ceksp,debug,doorth,fform,fform_dum,fill
197 : integer :: headform0,iband,ibg,ibg0,icg,icg0,icgsft,ieigsft,ierr,ii
198 : integer :: ikassoc,ikpt,ikpt0,ikptsp,ikptsp0,imax,increase_nkassoc,isppol,isppol0
199 : integer :: mband0,mband0_rd,mband_eff,mcg_disk,me,me0,mkmem0,mpw0
200 : integer :: my_nkpt,my_nspinor,my_nspinor0,nband_k,nband0_k
201 : integer :: nkassoc,nkpt0,npw,npw0,nspinor0,nspinor_eff,nsppol0,nsppol_eff,nsppol2nspinor
202 : integer :: rdwr,randalg,restart,restartpaw,spaceComm,spaceComm_io,sppoldbl,sppoldbl_eff,squeeze
203 : logical :: out_of_core
204 : real(dp) :: dksqmax,ecut0
205 : character(len=500) :: msg
206 6619 : type(hdr_type) :: hdr0
207 : integer :: ngfft0(18)
208 6619 : integer,allocatable :: indkk0(:,:),indx(:),istwfk0(:),kg0(:,:)
209 6619 : integer,allocatable :: nband0_rd(:),npwarr0(:),npwi(:),npwtot0(:)
210 6619 : integer,allocatable,target :: indkk(:,:),nband0(:)
211 6619 : integer, pointer :: indkk_eff(:,:),nband_eff(:)
212 6619 : logical,allocatable :: my_kpt(:)
213 : real(dp) :: gmet(3,3),gmet0(3,3),gprim0(3,3),rprim0(3,3),tsec(2)
214 6619 : real(dp),allocatable :: cg_disk(:,:),kptns0(:,:)
215 6619 : real(dp),pointer :: cg_eff(:,:),eigen_eff(:)
216 : type(MPI_type),pointer :: mpi_enreg0
217 :
218 : ! *************************************************************************
219 :
220 : DBG_ENTER("COLL")
221 :
222 : !Keep track of total time spent in inwffil
223 6619 : call timab(710,1,tsec)
224 6619 : call timab(711,1,tsec)
225 :
226 : !Check the validity of formeig
227 6619 : if (formeig/=0.and.formeig/=1) then
228 0 : write(msg,'(a,i0,a)')' formeig = ',formeig,', but the only allowed values are 0 or 1.'
229 0 : ABI_BUG(msg)
230 : end if
231 :
232 : !Init mpi_comm
233 6619 : spaceComm=mpi_enreg%comm_cell
234 6619 : spaceComm_io=xmpi_comm_self
235 6619 : if (mpi_enreg%paral_kgb==1) spaceComm_io= mpi_enreg%comm_bandspinorfft
236 6619 : if (mpi_enreg%paral_hf ==1) spaceComm_io= mpi_enreg%comm_hf
237 6619 : me=xmpi_comm_rank(spaceComm)
238 :
239 : !Determine number of k points processed by current node
240 131325 : my_nkpt=nkpt;if (size(mpi_enreg%my_kpttab)>0) my_nkpt=maxval(mpi_enreg%my_kpttab)
241 6619 : out_of_core=(mkmem==0.and.my_nkpt/=0)
242 :
243 6619 : ngfft0(:)=ngfft(:)
244 6619 : headform0=0 !Default value for headform0 (will be needed later, to read wf blocks)
245 :
246 : !Chebyshev is more sensitive to the quality of input random numbers, so use a new algorithm
247 6619 : if(dtset%wfoptalg == 1 .or. dtset%wfoptalg == 2 .or. dtset%wfoptalg == 111 .or. dtset%wfoptalg == 112) then
248 208 : randalg = 1
249 : else
250 : ! Otherwise, use compatibility mode
251 6411 : randalg = 0
252 : end if
253 :
254 : !If the input data are on disk, determine the kind of restart
255 6619 : wff1%fname = wffnm
256 :
257 : !Checking the existence of data file
258 6619 : if (ireadwf==1 .and. .not.file_exists(wff1%fname)) then
259 : ! Trick needed to run Abinit test suite in netcdf mode.
260 2248 : if (file_exists(nctk_ncify(wff1%fname))) then
261 2248 : write(std_out,"(3a)")"- File: ",trim(wff1%fname)," does not exist but found netcdf file with similar name."
262 4496 : wff1%fname = nctk_ncify(wff1%fname)
263 : end if
264 2248 : if (localrdwf/=0 .and. .not. file_exists(wff1%fname)) then
265 0 : ABI_ERROR('Missing data file: '//TRIM(wff1%fname))
266 : end if
267 : end if
268 :
269 : !Compute reciprocal space metric gmet
270 6619 : call matr3inv(hdr%rprimd,gprim0) ! gprim0 is used as temporary storage
271 264760 : gmet=matmul(transpose(gprim0),gprim0)
272 :
273 6619 : if (ireadwf==1)then
274 :
275 2252 : iomode=dtset%iomode
276 2252 : if (localrdwf==0) then
277 : ! This is in case the wff file must be read by only the master proc
278 7 : if (iomode /= IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER
279 : !iomode=IO_MODE_FORTRAN_MASTER
280 : end if
281 :
282 2252 : call WffOpen(iomode,spaceComm,wff1%fname,ierr,wff1,master,me,unwff1,spaceComm_io)
283 :
284 : ! Initialize hdr0 (sent to all procs), thanks to reading of wff1
285 2252 : rdwr=1
286 2252 : if ( ANY(wff1%iomode == (/IO_MODE_FORTRAN_MASTER, IO_MODE_FORTRAN, IO_MODE_MPI/) )) then
287 4 : call hdr_io(fform_dum,hdr0,rdwr,wff1)
288 2248 : else if (wff1%iomode == IO_MODE_ETSF) then
289 2248 : call hdr0%ncread(wff1%unwff, fform_dum)
290 : end if
291 :
292 : ! Handle IO Error.
293 2252 : if (fform_dum == 0) then
294 : write(msg,"(4a)")&
295 0 : "hdr_io returned fform == 0 while trying to read the wavefunctions from file: ",trim(wff1%fname),ch10,&
296 0 : "This usually means that the file does not exist or that you don't have enough privileges to read it"
297 0 : ABI_ERROR(msg)
298 : end if
299 :
300 2252 : call wrtout(std_out,' inwffil: examining the header of disk file: '//trim(wff1%fname),'COLL')
301 :
302 : ! Check hdr0 versus hdr (and from now on ignore header consistency and write new info to header for each file)
303 2252 : if (dtset%usewvl == 0) then
304 : ! wait for plane waves.
305 2252 : fform=2
306 : else
307 : ! wait for wavelets.
308 0 : fform = 200
309 : end if
310 2252 : call hdr_check(fform,fform_dum,hdr,hdr0,'PERS',restart,restartpaw)
311 :
312 2252 : nkpt0=hdr0%nkpt
313 2252 : nsppol0=hdr0%nsppol
314 2252 : headform0=hdr0%headform
315 :
316 : ! Extended plane waves energy shift is required in order to restart from wave function.
317 2252 : hdr%extfpmd_eshift=hdr0%extfpmd_eshift
318 :
319 2252 : write(msg,'(2a)')'-inwffil : will read wavefunctions from disk file ',trim(wff1%fname)
320 2252 : call wrtout(std_out,msg,'COLL')
321 2252 : call wrtout(ab_out,msg,'COLL')
322 :
323 : else
324 4367 : restart=1; restartpaw=0
325 :
326 : ! Fill some data concerning an hypothetical file to be read
327 : ! This is to allow the safe use of same routines than with ireadwf==1.
328 4367 : nkpt0=nkpt ; nsppol0=nsppol
329 : end if ! end ireadwf
330 :
331 6619 : sppoldbl=1
332 137558 : if(minval(symafm(:))==-1)then
333 121 : if(nsppol0==1 .and. nsppol==2)sppoldbl=2
334 : end if
335 :
336 19857 : ABI_MALLOC(indkk,(nkpt*sppoldbl,6))
337 19857 : ABI_MALLOC(istwfk0,(nkpt0))
338 19857 : ABI_MALLOC(kptns0,(3,nkpt0))
339 19857 : ABI_MALLOC(nband0,(nkpt0*nsppol0))
340 13238 : ABI_MALLOC(npwarr0,(nkpt0))
341 :
342 6619 : if(restart==2)then ! restart with translations
343 :
344 517 : ecut0=hdr0%ecut_eff
345 5748 : istwfk0(1:nkpt0)=hdr0%istwfk(1:nkpt0)
346 21441 : kptns0(1:3,1:nkpt0)=hdr0%kptns(1:3,1:nkpt0)
347 5803 : nband0(1:nkpt0*nsppol0)=hdr0%nband(1:nkpt0*nsppol0)
348 2068 : ngfft0(1:3)=hdr0%ngfft(1:3)
349 5748 : npwarr0(1:nkpt0)=hdr0%npwarr(1:nkpt0)
350 517 : nspinor0=hdr0%nspinor
351 6721 : rprim0(:,:)=hdr0%rprimd(:,:)
352 5748 : mpw0=maxval(npwarr0(:))
353 :
354 : ! Compute reciprocal space metric gmet for unit cell of disk wf
355 517 : call matr3inv(rprim0,gprim0)
356 20680 : gmet0=matmul(transpose(gprim0),gprim0)
357 :
358 517 : if ((mpi_enreg%paral_kgb==1).or.(mpi_enreg%paral_hf==1)) then
359 0 : ABI_MALLOC(mpi_enreg0,)
360 0 : call copy_mpi_enreg(mpi_enreg,mpi_enreg0)
361 0 : ABI_MALLOC(kg0,(3,mpw0*nkpt0))
362 0 : ABI_MALLOC(npwtot0,(nkpt0))
363 0 : msg="tmpfil"
364 : call kpgio(ecut0,dtset%exchn2n3d,gmet0,istwfk0,kg0, &
365 : & kptns0,nkpt0,nband0,nkpt0,'PERS',mpi_enreg0,&
366 0 : & mpw0,npwarr0,npwtot0,nsppol0)
367 :
368 0 : ABI_FREE(kg0)
369 0 : ABI_FREE(npwtot0)
370 : else
371 : mpi_enreg0 => mpi_enreg
372 : end if
373 :
374 : ! At this stage, the header of the file wff1i%unwff is read, and
375 : ! the pointer is ready to read the first wavefunction block.
376 :
377 : ! Compute k points from input file closest to the output file
378 517 : call listkk(dksqmax,gmet0,indkk,kptns0,kptns,nkpt0,nkpt,nsym,sppoldbl,symafm,symrel,1,spaceComm)
379 :
380 6102 : else if (restart==1) then ! direct restart
381 :
382 : ! Fill variables that must be the same, as determined by hdr_check.f
383 : ! This is to allow the safe use of the same routines than with restart==2.
384 6102 : nspinor0=dtset%nspinor
385 6102 : ecut0=ecut_eff
386 6102 : gmet0(:,:)=gmet(:,:)
387 117180 : istwfk0(:)=istwfk(:)
388 450414 : kptns0(:,:)=kptns(:,:)
389 117180 : npwarr0(:)=npwarr(:)
390 6102 : mpw0=mpw
391 :
392 12204 : do isppol=1,sppoldbl
393 123282 : do ikpt=1,nkpt
394 111078 : indkk(ikpt+(isppol-1)*nkpt,1)=ikpt
395 672570 : indkk(ikpt+(isppol-1)*nkpt,2:6)=0
396 : end do
397 : end do
398 6102 : dksqmax=0.0_dp
399 :
400 : ! The treatment of nband0 asks for some care
401 6102 : if(ireadwf==0)then
402 59644 : nband0(:)=0
403 : else
404 62695 : nband0(1:nkpt0*nsppol0)=hdr0%nband(1:nkpt0*nsppol0)
405 : end if
406 :
407 : mpi_enreg0 => mpi_enreg
408 :
409 : else
410 : mpi_enreg0 => mpi_enreg
411 : end if
412 :
413 6619 : if(mpi_enreg0%paral_pert == 1.and.mpi_enreg0%me_pert/=-1) then
414 0 : me0 = mpi_enreg0%me_pert
415 : else
416 6619 : me0 = mpi_enreg0%me_cell
417 : end if
418 :
419 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420 : !Before hdr_free:
421 : !If restartpaw==1, store hdr0%pawrhoij in hdr%pawrhoij; else if restartpaw==0,
422 : !hdr%pawrhoij(:)has been initialized in hdr_init.
423 6619 : if(restartpaw==1) then
424 541 : call pawrhoij_copy(hdr0%pawrhoij,hdr%pawrhoij,keep_itypat=.true.)
425 : end if
426 :
427 6619 : call timab(711,2,tsec)
428 6619 : call timab(712,1,tsec)
429 :
430 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
431 : !At this stage, all the relevant information from the header of the disk file,
432 : !has been exploited, and stored in variables, on all processors.
433 : !It is also contained in hdr0
434 : !(on all processors, except if restart=1 and localrdwf=0,
435 : !in which case it is only on the master)
436 : !These information might be changed later, while processing the
437 : !wavefunction data, and converting it. The variable hdr0 might be kept
438 : !for further checking, or reference, or debugging, but at present,
439 : !it is simpler to close it. The other header, hdr, will be used for the new file, if any.
440 :
441 6619 : if(ask_accurate==1)then
442 :
443 : ! Check whether the accuracy requirements might be fulfilled
444 1038 : if(ireadwf==0)then
445 : write(msg,'(9a)')&
446 0 : & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
447 0 : & 'present calculation. It was asked that the wavefunctions be accurate,',ch10,&
448 0 : & 'but they were not even read.',ch10,&
449 0 : & 'Action: use a wf file, with ireadwf/=0.'
450 0 : ABI_ERROR(msg)
451 : end if
452 1038 : if(dksqmax>tol12)then
453 : write(msg, '(9a,es16.6,4a)' )&
454 0 : & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
455 0 : & 'present calculation. It was asked that the wavefunctions be accurate, but',ch10,&
456 0 : & 'at least one of the k points could not be generated from a symmetrical one.',ch10,&
457 0 : & 'dksqmax=',dksqmax,ch10,&
458 0 : & 'Action: check your wf file and k point input variables',ch10,&
459 0 : & ' (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.'
460 0 : ABI_ERROR(msg)
461 : end if
462 1038 : if(dtset%nspinor/=nspinor0)then
463 : write(msg,'(a,a, a,a,a,a,a, a,a,2i5,a,a)')&
464 0 : & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
465 0 : & 'present calculation. It was asked that the wavefunctions be accurate, but',ch10,&
466 0 : & 'nspinor differs in the file from the actual nspinor.',ch10,&
467 0 : & 'nspinor,nspinor0=',dtset%nspinor,nspinor0,ch10,&
468 0 : & 'Action: check your wf file, and nspinor input variables.'
469 0 : ABI_ERROR(msg)
470 : end if
471 1038 : if((nsppol>nsppol0 .and. sppoldbl==1) .or. nsppol<nsppol0 ) then
472 : write(msg,'(a,a, a,a,a,a,a, a,a,3i5,a,a)')&
473 0 : & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
474 0 : & 'present calculation. It was asked that the wavefunctions be accurate, but',ch10,&
475 0 : & 'the nsppol variables do not match in the file and in the actual calculation',ch10,&
476 0 : & 'nsppol,nsppol,sppoldbl=',dtset%nspinor,nspinor0,sppoldbl,ch10,&
477 0 : & 'Action: check your wf file, and nsppol input variables.'
478 0 : ABI_ERROR(msg)
479 : end if
480 :
481 : ! Now, check the number of bands
482 : accurate=1
483 2456 : do isppol=1,nsppol
484 60789 : do ikpt=1,nkpt
485 58333 : ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1)
486 58333 : ikptsp =ikpt +(isppol-1)*nkpt
487 58333 : ikptsp0=ikpt0+(isppol-1)*(2-sppoldbl)*nkpt0
488 59751 : if(nband0(ikptsp0)<nband(ikptsp))accurate=0
489 : end do
490 : end do
491 1038 : if(accurate==0)then
492 : write(msg,'(a,a, a,a,a,a,a, a,a)')&
493 0 : & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,&
494 0 : & 'present calculation. It was asked that the wavefunctions be accurate,',ch10,&
495 0 : & 'but the number of bands differ in the file and in the actual calculation.',ch10,&
496 0 : & 'Action: use a wf file with the correct characteristics.'
497 0 : ABI_ERROR(msg)
498 : end if
499 :
500 : end if
501 :
502 : !Flag: do we need to translate WF to (from) spinors ?
503 6619 : nsppol2nspinor=0
504 6619 : if (nsppol0==2.and.dtset%nspinor==2) nsppol2nspinor=+1
505 6619 : if (nspinor0==2.and.nsppol==2) nsppol2nspinor=-1
506 :
507 : !Take into account parallism over spinors
508 6619 : my_nspinor =max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
509 6619 : my_nspinor0=max(1,nspinor0/mpi_enreg0%nproc_spinor)
510 :
511 : !Not all bands might be read, if not needed to fill the wavefunctions
512 128142 : mband0=maxval(nband0(1:nkpt0*nsppol0))
513 6619 : mband0_rd=min(mband0,(mband/dtset%nspinor)*nspinor0)
514 :
515 : !****************************************************************************
516 : !If needed, transfer the input wf from disk to core memory
517 : !(in the parallel case, it allows to change localrdwf=0 in localrdwf=1)
518 :
519 6619 : mkmem0=0
520 :
521 : if(xmpi_paral == 1 .or. mpi_enreg%paral_kgb == 1 .or. mpi_enreg%paral_hf == 1) then
522 6619 : if(localrdwf==0 .and. out_of_core)then
523 0 : ABI_BUG('localrdwf==0 and mkmem==0 (out-of-core solution) are not allowed together (yet)')
524 : end if
525 : end if
526 :
527 6619 : call timab(712,2,tsec)
528 :
529 : !Here, treat reading wavefunctions with mkmem/=0, first step
530 6619 : if(ireadwf==1 .and. (.not.out_of_core))then
531 :
532 2252 : call timab(713,1,tsec)
533 :
534 : ! if(restart==1 .and. ireadwf==1 .and. mkmem/=0)then
535 :
536 : ! Compute table of k point associations. Make a trial choice for nkassoc.
537 2252 : nkassoc=(nkpt/nkpt0+1)*2
538 9008 : ABI_MALLOC(indkk0,(nkpt0,nkassoc))
539 : ! Infinite loops are allowed in F90
540 : do
541 306812 : indkk0(:,:)=0
542 2292 : increase_nkassoc=0
543 76056 : do ikpt=1,nkpt*sppoldbl
544 73804 : ikpt0=indkk(ikpt,1)
545 147220 : do ikassoc=1,nkassoc
546 147180 : if(indkk0(ikpt0,ikassoc)==0)then
547 73764 : indkk0(ikpt0,ikassoc)=ikpt
548 73764 : exit
549 : end if
550 73456 : if(nkassoc==ikassoc)increase_nkassoc=1
551 : end do
552 76056 : if(increase_nkassoc==1)then
553 40 : ABI_FREE(indkk0)
554 40 : nkassoc=2*nkassoc
555 160 : ABI_MALLOC(indkk0,(nkpt0,nkassoc))
556 40 : exit
557 : end if
558 : end do
559 2292 : if(increase_nkassoc==0)exit
560 : end do
561 :
562 : ! DEBUG
563 : ! write(std_out,*)' inwffil: indkk0, nkassoc=',nkassoc
564 : ! do ikpt0=1,nkpt0
565 : ! write(std_out,*)' ikpt0,indkk0(ikpt0,1)=',ikpt0,indkk0(ikpt0,1)
566 : ! end do
567 : ! ENDDEBUG
568 :
569 : ! DEBUG
570 : ! write(std_out,*)' inwffil : indkk(:,1)=',indkk(:,1)
571 : ! write(std_out,*)' inwffil : sppoldbl=',sppoldbl
572 : ! ENDDEBUG
573 :
574 : ! To treat the case (nsppol0=2,nspinor0=1)<->(nsppol=1,nspinor=2),
575 : ! apply the following trick:
576 : ! 1- We call wfsinp with fake arguments (nsppol_eff and nspinor_eff)
577 : ! 2- We transform collinear polarized WF into spinors
578 : ! or spinors into collinear polarized WF
579 : if (nsppol2nspinor/=0.and.out_of_core.and.dtset%usewvl==0) then
580 : write(msg, '(7a)')&
581 : & 'When mkmem=0 (out-of-core), the wavefunction translator is unable',ch10,&
582 : & 'to interchange spin-polarized wfs and spinor wfs.',ch10,&
583 : & 'Action: use a non-spin-polarized wf to start a spinor wf,',ch10,&
584 : & ' and a non-spinor wf to start a spin-polarized wf.'
585 : ABI_ERROR(msg)
586 : end if
587 :
588 : ! === Fake arguments definition for wfsinp
589 2252 : if (nsppol2nspinor==0.or.dtset%usewvl/=0) then
590 2251 : indkk_eff => indkk
591 2251 : nband_eff => nband
592 2251 : eigen_eff => eigen
593 2251 : cg_eff => cg
594 2251 : nspinor_eff=dtset%nspinor;nsppol_eff=nsppol;sppoldbl_eff=sppoldbl
595 76992 : mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff))
596 1 : else if (nsppol2nspinor==1.and.(.not.out_of_core)) then
597 1 : nsppol_eff=2;nspinor_eff=1;sppoldbl_eff=1
598 3 : ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6))
599 3 : ABI_MALLOC(nband_eff,(nkpt*nsppol_eff))
600 109 : indkk_eff(1:nkpt,1:6) =indkk(1:nkpt,1:6)
601 17 : nband_eff(1:nkpt) =nband(1:nkpt)/2
602 17 : nband_eff(1+nkpt:2*nkpt)=nband(1:nkpt)/2
603 17 : mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff))
604 1 : eigen_eff => eigen
605 1 : cg_eff => cg
606 : else if (nsppol2nspinor==-1.and.(.not.out_of_core)) then
607 : ! WARNING: MT 07072011 -> this is memory consuming
608 : ! A copy a spinorial WF (and eigenvalues) is temporary kept in memory;
609 : ! But the case (nspinor=2 => nsppol=2) might be rare
610 : ! and only useful for testing purposes.
611 : ! => print a warning for the user
612 : ! NOTE: in that case (nsppol=2), parallelization over spinors is not activated
613 :
614 : write(msg,'(5a)')&
615 0 : & 'In the case of spinor WF read from disk and converted into',ch10,&
616 0 : & 'spin-polarized non-spinor WF, the WF translator is memory',ch10,&
617 0 : & 'consuming (a copy of the spinor WF is temporarily stored in memory).'
618 0 : ABI_WARNING(msg)
619 :
620 0 : nsppol_eff=1;nspinor_eff=2;sppoldbl_eff=1
621 0 : ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6))
622 0 : ABI_MALLOC(nband_eff,(nkpt*nsppol_eff))
623 0 : indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6)
624 0 : nband_eff(1:nkpt) =2*nband(1:nkpt)
625 0 : mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff))
626 0 : ABI_MALLOC(eigen_eff,((2*mband_eff)**formeig*mband_eff*nkpt*nsppol_eff))
627 0 : ABI_MALLOC(cg_eff,(2,mpw0*nspinor_eff*mband_eff*mkmem*nsppol_eff))
628 : end if
629 :
630 : ! === nband0 argument definition for wfsinp
631 2252 : squeeze=0
632 2252 : ABI_MALLOC(cg_disk,(0,0))
633 : if(.not.out_of_core)then
634 6756 : ABI_MALLOC(nband0_rd,(nkpt0*nsppol0))
635 68498 : nband0_rd(:)=0
636 5000 : do isppol=1,nsppol_eff
637 79757 : do ikpt=1,nkpt
638 74757 : ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1)
639 74757 : isppol0=min(isppol,nsppol0)
640 74757 : ikptsp =ikpt +(isppol -1)*nkpt
641 74757 : ikptsp0=ikpt0+(isppol0-1)*(2-sppoldbl)*nkpt0
642 74757 : nband0_k=min(nband0(ikptsp0),(nband_eff(ikptsp)/nspinor_eff)*nspinor0)
643 74757 : nband0_rd(ikptsp0)=max(nband0_rd(ikptsp0),nband0_k)
644 74757 : npw0=npwarr0(ikpt0)
645 74757 : npw =npwarr (ikpt)
646 77505 : if(npw0*nspinor0*nband0_k > npw*nspinor_eff*nband_eff(ikptsp))squeeze=1
647 : end do
648 : end do
649 2252 : if(squeeze==1)then
650 86 : mcg_disk=mpw0*my_nspinor0*mband0_rd
651 86 : ABI_FREE(cg_disk)
652 258 : ABI_MALLOC(cg_disk,(2,mcg_disk))
653 : else
654 : if(xmpi_paral == 1 .or. mpi_enreg0%paral_kgb == 1 .or. mpi_enreg0%paral_hf == 1)then
655 2166 : if(localrdwf==0)then
656 7 : mcg_disk=mpw0*my_nspinor0*mband0_rd
657 7 : ABI_FREE(cg_disk)
658 21 : ABI_MALLOC(cg_disk,(2,mcg_disk))
659 : end if
660 : end if
661 : end if
662 : end if
663 :
664 2252 : call timab(713,2,tsec)
665 2252 : call timab(714,1,tsec)
666 :
667 : ! === call to wfsinp
668 2252 : if (dtset%usewvl == 0) then
669 : call wfsinp(cg_eff,cg_disk,ecut,ecut0,ecut_eff,eigen,&
670 : & exchn2n3d,formeig,gmet,gmet0,headform0,&
671 : & indkk_eff,indkk0,istwfk,istwfk0,kptns,kptns0,localrdwf,&
672 : & mband_eff,mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,&
673 : & nband_eff,nband0_rd,ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor_eff,nspinor0,&
674 : & nsppol_eff,nsppol0,nsym,occ,optorth,dtset%prtvol,randalg,restart,hdr%rprimd,sppoldbl_eff,squeeze,&
675 2252 : & symrel,tnons,wff1)
676 2252 : if (nsppol2nspinor/=0) then
677 1 : ABI_FREE(indkk_eff)
678 1 : ABI_FREE(nband_eff)
679 : end if
680 : else
681 : ! Read wavefunctions from file.
682 : call wvl_wfsinp_disk(dtset, hdr0, hdr, mpi_enreg, occ, 1, &
683 0 : & hdr%rprimd, wff1, wvl%wfs, wvl%descr, hdr%xred)
684 : end if
685 :
686 2252 : call timab(714,2,tsec)
687 2252 : call timab(715,1,tsec)
688 :
689 : ! Now, update xyz0 variables, for use in newkpt
690 68498 : nband0(:)=nband0_rd(:)
691 :
692 : ! If squeeze, the conversion was done in wfsinp, so no conversion left.
693 2252 : if(squeeze==1)then
694 86 : ecut0=ecut_eff
695 86 : gmet0(:,:)=gmet(:,:)
696 86 : ABI_FREE(kptns0)
697 86 : ABI_FREE(istwfk0)
698 86 : ABI_FREE(nband0)
699 86 : ABI_FREE(npwarr0)
700 258 : ABI_MALLOC(kptns0,(3,nkpt))
701 258 : ABI_MALLOC(istwfk0,(nkpt))
702 258 : ABI_MALLOC(nband0,(nkpt*nsppol))
703 172 : ABI_MALLOC(npwarr0,(nkpt))
704 19126 : kptns0(:,:)=kptns(:,:)
705 4846 : istwfk0(:)=istwfk(:)
706 4846 : npwarr0(:)=npwarr(:)
707 4846 : nband0(:)=0
708 172 : do isppol=1,nsppol
709 4932 : do ikpt=1,nkpt
710 4760 : ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1)
711 4760 : isppol0=min(isppol,nsppol0)
712 4760 : ikptsp =ikpt +(isppol -1)*nkpt
713 4760 : ikptsp0=ikpt0+(isppol0-1)*(sppoldbl-1)*nkpt0
714 4846 : nband0(ikptsp)=(nband0_rd(ikptsp0)/nspinor0)*dtset%nspinor
715 : end do
716 : end do
717 4846 : do ikpt=1,nkpt
718 4760 : indkk(ikpt,1)=ikpt
719 28646 : indkk(ikpt,2:6)=0
720 : end do
721 : ! This transfer must come after the nband0 transfer
722 86 : nspinor0=dtset%nspinor
723 86 : nkpt0=nkpt
724 86 : nsppol0=nsppol
725 : end if ! end squeeze == 1
726 :
727 : ! The input wavefunctions have been transferred from disk to core memory
728 2252 : mkmem0=mkmem
729 :
730 2252 : ABI_FREE(indkk0)
731 2252 : ABI_FREE(nband0_rd)
732 2252 : ABI_FREE(cg_disk)
733 :
734 2252 : call timab(715,2,tsec)
735 :
736 : else !ireadwf == 0
737 4367 : if (dtset%usewvl == 1) then
738 :
739 0 : call timab(714,1,tsec)
740 : ! Compute wavefunctions from input guess.
741 0 : call wvl_wfsinp_scratch(dtset, mpi_enreg, occ, hdr%rprimd, wvl, hdr%xred)
742 0 : call timab(714,2,tsec)
743 : end if
744 : end if
745 :
746 6619 : call timab(716,1,tsec)
747 :
748 : !=== Eventual conversion of WF into (from) spinors
749 6619 : if (dtset%usewvl==0) then
750 :
751 : ! ***** No conversion (standard case) ****
752 6619 : if (nsppol2nspinor==0) then
753 6618 : nspinor_eff=nspinor0;nsppol_eff=nsppol0;sppoldbl_eff=sppoldbl
754 6618 : indkk_eff => indkk
755 6618 : nband_eff => nband0
756 :
757 : ! ***** Conversion from collinear to spinorial WF ****
758 1 : else if (nsppol2nspinor==1.and.(.not.out_of_core)) then
759 : ! Translate the WF and eigenvalues from nsppol=2 to nspinor=2
760 : ! This is tricky (because we do not want to create a temporary array for cg)
761 1 : nsppol_eff=1;nspinor_eff=2;sppoldbl_eff=1
762 3 : ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6))
763 3 : ABI_MALLOC(nband_eff,(nkpt0*nsppol_eff))
764 109 : indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6)
765 17 : nband_eff(1:nkpt0)=2*nband0(1:nkpt0)
766 : ! Compute some shifts from isspol0=1 to isppol0=2
767 1 : imax=0;icgsft=0;ieigsft=0
768 3 : ABI_MALLOC(my_kpt,(nkpt0))
769 9 : do ikpt0=1,nkpt0
770 8 : nband0_k=nband0(ikpt0);nband_k=nband(ikpt0)
771 8 : my_kpt(ikpt0)=(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me0)))
772 8 : ieigsft=ieigsft+(2*nband0_k)**formeig*nband0_k
773 9 : if(my_kpt(ikpt0)) then
774 8 : imax=imax+nband0_k;icgsft=icgsft+nband0_k*npwarr0(ikpt0)
775 : end if
776 : end do
777 : ! --- First version: no parallelization over spinors
778 1 : if (mpi_enreg0%paral_spinor==0) then
779 : ! Compute some useful indexes
780 3 : ABI_MALLOC(indx,(2*imax))
781 3 : ABI_MALLOC(npwi,(imax))
782 1 : ii=0;icg=0
783 9 : do ikpt0=1,nkpt0
784 9 : if(my_kpt(ikpt0)) then
785 8 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
786 104 : do iband=1,nband0_k
787 96 : ii=ii+1;npwi(ii)=npw0
788 96 : indx(2*ii-1)=icg+mpw0;indx(2*ii)=icg+2*mpw0
789 104 : icg=icg+4*mpw0
790 : end do
791 : end if
792 : end do
793 : ! Expand WF in cg (try to use the whole array)
794 1 : ii=nsppol0*imax;icg0=nsppol0*icgsft
795 3 : do isppol=nsppol0,1,-1
796 19 : do ikpt0=nkpt0,1,-1
797 18 : if(my_kpt(ikpt0)) then
798 16 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
799 208 : do iband=nband0_k,1,-1
800 192 : icg0=icg0-npw0
801 192 : if (indx(ii)<icg0) then
802 0 : ABI_BUG("Unable to read WF!")
803 : end if
804 98400 : cg(:,indx(ii)+1:indx(ii)+npw0)=cg(:,icg0+1:icg0+npw0)
805 208 : ii=ii-1
806 : end do
807 : end if
808 : end do
809 : end do
810 : ! Convert polarized WF into spinors
811 1 : ii=1
812 9 : do ikpt0=1,nkpt0
813 9 : if(my_kpt(ikpt0)) then
814 8 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
815 104 : do iband=1,nband0_k
816 96 : npw0=npwi(ii)
817 49200 : cg(:,indx(2*ii-1)-mpw0+1:indx(2*ii-1)-mpw0+npw0)=cg(:,indx(ii)+1:indx(ii)+npw0)
818 49200 : cg(:,indx(2*ii )+mpw0+1:indx(2*ii )+mpw0+npw0)=cg(:,indx(ii+imax)+1:indx(ii+imax)+npw0)
819 104 : ii=ii+1
820 : end do
821 : end if
822 : end do
823 : ! Compress new cg array (from mpw to npw) and cancel zero-components
824 1 : icg0=0;icg=0
825 9 : do ikpt0=1,nkpt0
826 9 : if(my_kpt(ikpt0)) then
827 8 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
828 104 : do iband=1,nband0_k
829 49200 : cg(:,icg0 +1:icg0+ npw0)=cg(:,icg+1:icg+npw0)
830 24648 : cg(:,icg0+ npw0+1:icg0+2*npw0)=zero
831 24648 : cg(:,icg0+2*npw0+1:icg0+3*npw0)=zero
832 49200 : cg(:,icg0+3*npw0+1:icg0+4*npw0)=cg(:,icg+3*mpw0+1:icg+3*mpw0+npw0)
833 104 : icg0=icg0+4*npw0;icg=icg+4*mpw0
834 : end do
835 : end if
836 : end do
837 : ! --- Second version: parallelization over spinors
838 : else
839 : ! Compute some useful indexes
840 0 : ABI_MALLOC(indx,(imax))
841 0 : ABI_MALLOC(npwi,(imax))
842 0 : ii=0;icg=0
843 0 : do ikpt0=1,nkpt0
844 0 : if(my_kpt(ikpt0)) then
845 0 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
846 0 : do iband=1,nband0_k
847 0 : ii=ii+1;npwi(ii)=npw0
848 0 : indx(ii)=icg+mpi_enreg0%me_spinor*mpw0
849 0 : icg=icg+2*mpw0
850 : end do
851 : end if
852 : end do
853 : ! Expand WF in cg
854 0 : ii=(mpi_enreg0%me_spinor+1)*imax;icg0=(mpi_enreg0%me_spinor+1)*icgsft
855 0 : do ikpt0=nkpt0,1,-1
856 0 : if(my_kpt(ikpt0)) then
857 0 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
858 0 : do iband=nband0_k,1,-1
859 0 : icg0=icg0-npw0
860 0 : if (indx(ii)<icg0) then
861 0 : ABI_BUG("Unable to read WF!")
862 : end if
863 0 : cg(:,indx(ii)+1:indx(ii)+npw0)=cg(:,icg0+1:icg0+npw0)
864 0 : ii=ii-1
865 : end do
866 : end if
867 : end do
868 : ! Compress new cg array (from mpw to npw) and cancel zero-components
869 0 : icg0=0;icg=0
870 0 : do ikpt0=1,nkpt0
871 0 : if(my_kpt(ikpt0)) then
872 0 : nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0)
873 0 : do iband=1,nband0_k
874 0 : if (mpi_enreg0%me_spinor==0) then
875 0 : cg(:,icg0 +1:icg0+ npw0)=cg(:,icg+1:icg+npw0)
876 0 : cg(:,icg0+npw0+1:icg0+2*npw0)=zero
877 : else
878 0 : cg(:,icg0 +1:icg0+ npw0)=zero
879 0 : cg(:,icg0+npw0+1:icg0+2*npw0)=cg(:,icg+mpw0+1:icg+mpw0+npw0)
880 : end if
881 0 : icg0=icg0+2*npw0;icg=icg+2*mpw0
882 : end do
883 : end if
884 : end do
885 : end if
886 : ! Translate eigenvalues
887 1 : ibg0=2*ieigsft;ibg=2*ieigsft
888 9 : do ikpt0=nkpt0,1,-1
889 8 : nband0_k=nband0(ikpt0)
890 8 : ibg0=ibg0- nband0_k*(2*nband0_k)**formeig
891 8 : ibg =ibg -2*nband0_k*(2*nband0_k)**formeig
892 9 : if(my_kpt(ikpt0)) then
893 104 : do iband=nband0_k*(2*nband0_k)**formeig,1,-1
894 96 : eigen(2*iband-1+ibg)=eigen(iband+ibg0-ieigsft)
895 104 : eigen(2*iband +ibg)=eigen(iband+ibg0)
896 : end do
897 : end if
898 : end do
899 1 : ABI_FREE(indx)
900 1 : ABI_FREE(npwi)
901 1 : ABI_FREE(my_kpt)
902 :
903 : ! ***** Conversion from spinorial to collinear WF ****
904 0 : else if (nsppol2nspinor==-1.and.(.not.out_of_core)) then
905 : ! In that case parallelization over spinors is never activated
906 0 : nsppol_eff=2;nspinor_eff=1;sppoldbl_eff=1
907 0 : ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6))
908 0 : ABI_MALLOC(nband_eff,(nkpt0*nsppol_eff))
909 0 : indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6)
910 0 : nband_eff(1:nkpt0) =nband0(1:nkpt0)/2
911 0 : nband_eff(1+nkpt0:2*nkpt0)=nband0(1:nkpt0)/2
912 : ! Compute shifts from isspol0=1 to isppol0=2
913 0 : icgsft=0;ieigsft=0
914 0 : do ikpt0=1,nkpt0
915 0 : nband0_k=nband0(ikpt0);nband_k=nband(ikpt0)
916 0 : ieigsft=ieigsft+(nband0_k/2)*(nband0_k)**formeig
917 0 : if(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me))) &
918 0 : & icgsft=icgsft+(nband0_k/2)*npwarr0(ikpt0)
919 : end do
920 : ! Translate the WF and eigenvalues from nspinor=2 to nsppol=2
921 0 : icg0=0;icg=0;ibg=0
922 0 : do ikpt0=1,nkpt0
923 0 : nband0_k=nband0(ikpt0);nband_k=nband(ikpt0);npw0=npwarr0(ikpt0)
924 0 : if(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me))) then
925 0 : do iband=1,nband0_k/2
926 0 : do ii=1,npw0
927 0 : cg(:,ii+icg) =cg_eff(:,ii+icg0)
928 0 : cg(:,ii+icg+icgsft)=cg_eff(:,ii+icg0+3*npw0)
929 : end do
930 0 : icg0=icg0+4*npw0;icg=icg+npw0
931 : end do
932 0 : do iband=(nband0_k/2)*(nband0_k)**formeig,1,-1
933 0 : eigen(iband+ibg) =eigen_eff(2*iband-1+2*ibg)
934 0 : eigen(iband+ibg+ieigsft)=eigen_eff(2*iband +2*ibg)
935 : ! occ(iband+ibg) =occ_eff(2*iband-1+2*ibg)
936 : ! occ(iband+ibg+ieigsft)=occ_eff(2*iband +2*ibg)
937 : end do
938 : end if
939 0 : ibg=ibg+(nband0_k/2)*(nband0_k)**formeig
940 : end do
941 0 : ABI_FREE(cg_eff)
942 0 : ABI_FREE(eigen_eff)
943 :
944 : else
945 0 : ABI_BUG('unable to interchange nsppol and nspinor when mkmem=0')
946 : end if
947 : end if
948 :
949 : !Clean hdr0
950 6619 : call hdr0%free()
951 :
952 6619 : call timab(716,2,tsec)
953 6619 : call timab(717,1,tsec)
954 :
955 :
956 : !****************************************************************************
957 : !Now, treat translation of wavefunctions if wavefunctions are planewaves
958 :
959 6619 : ceksp=0; debug=0; doorth=1; fill=1
960 6619 : if (dtset%usewvl == 0) then
961 :
962 : call newkpt(ceksp,cg,debug,ecut0,ecut,ecut_eff,eigen,exchn2n3d,&
963 : & fill,formeig,gmet0,gmet,headform0,indkk_eff,&
964 : & ab_out,ireadwf,istwfk0,istwfk,kg,kptns0,kptns,&
965 : & mband,mcg,mkmem0,mkmem,mpi_enreg0,mpi_enreg,&
966 : & mpw0,mpw,my_nkpt,nband_eff,nband,ngfft0,ngfft,nkpt0,nkpt,npwarr0,npwarr,&
967 : & nspinor_eff,dtset%nspinor,nsppol_eff,nsppol,nsym,occ,optorth,&
968 6619 : & dtset%prtvol,randalg,restart,hdr%rprimd,sppoldbl_eff,symrel,tnons,unkg,wff1,wffnow)
969 :
970 6619 : if (nsppol2nspinor/=0) then
971 1 : ABI_FREE(indkk_eff)
972 1 : ABI_FREE(nband_eff)
973 : end if
974 :
975 : end if ! dtset%usewvl == 0
976 :
977 : !****************************************************************************
978 :
979 6619 : ABI_FREE(indkk)
980 6619 : ABI_FREE(istwfk0)
981 6619 : ABI_FREE(kptns0)
982 6619 : ABI_FREE(nband0)
983 6619 : ABI_FREE(npwarr0)
984 6619 : if (restart==2 .and.(mpi_enreg0%paral_kgb==1 .or. mpi_enreg0%paral_hf == 1)) then
985 0 : call destroy_mpi_enreg(mpi_enreg0)
986 0 : ABI_FREE(mpi_enreg0)
987 : else
988 6619 : nullify(mpi_enreg0)
989 : end if
990 :
991 6619 : call timab(717,2,tsec)
992 6619 : call timab(710,2,tsec)
993 :
994 : DBG_EXIT("COLL")
995 :
996 13238 : end subroutine inwffil
997 : !!***
998 :
999 : !!****f* m_inwffil/wfsinp
1000 : !! NAME
1001 : !! wfsinp
1002 : !!
1003 : !! FUNCTION
1004 : !! Do initialization of wavefunction files.
1005 : !! Also call other relevant routines for this initialisation.
1006 : !! Detailed description :
1007 : !! - Initialize unit wff1 for input of wf data
1008 : !!
1009 : !! formeig option (format of the eigenvalues and occupations) :
1010 : !! 0 => ground-state format (initialisation of
1011 : !! eigenvectors with random numbers, vector of eigenvalues,
1012 : !! occupations are present)
1013 : !! 1 => respfn format (initialisation of
1014 : !! eigenvectors with 0 s, hermitian matrix of eigenvalues)
1015 : !!
1016 : !! INPUTS
1017 : !! ecut0=kinetic energy cutoffs for basis sphere 0 (hartree) (if squeeze=1)
1018 : !! ecut=kinetic energy cutoffs beyond which the coefficients of cg vanish (Ha)
1019 : !! (needed only if squeeze=1)
1020 : !! ecut_eff=effective kinetic energy planewave cutoff (hartree), needed
1021 : !! to generate the sphere of plane wave
1022 : !! exchn2n3d=if 1, n2 and n3 are exchanged
1023 : !! formeig=explained above
1024 : !! gmet(3,3), gmet0(3,3)=reciprocal space metrics (bohr^-2)
1025 : !! headform0=header format (might be needed to read the block of wfs)
1026 : !! indkk(nkpt*sppoldbl,6)=describe k point number of kptns0 that allows to
1027 : !! generate wavefunctions closest to given kpt
1028 : !! indkk(:,1)=k point number of kptns0
1029 : !! indkk(:,2)=symmetry operation to be applied to kpt0, to give kpt0a
1030 : !! (if 0, means no symmetry operation, equivalent to identity )
1031 : !! indkk(:,3:5)=shift in reciprocal space to be given to kpt0a,
1032 : !! to give kpt0b, that is the closest to kpt.
1033 : !! indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
1034 : !! indkk0(nkpt0,nkassoc)=list of k points that will be generated by k point number ikpt0
1035 : !! istwfk(nkpt)=input parameter that describes the storage of wfs
1036 : !! istwfk0(nkpt0)=input parameter that describes the storage of wfs in set0
1037 : !! kptns(3,nkpt),kptns0(3,nkpt0)=k point sets (reduced coordinates)
1038 : !! localrdwf=(for parallel case) if 1, the wff1%unwff file is local to each machine
1039 : !! mband=maximum number of bands
1040 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*nsppol
1041 : !! mpi_enreg=information about MPI parallelization
1042 : !! mpi_enreg0=information about MPI parallelization in set0
1043 : !! mpw=maximum number of planewaves as dimensioned in calling routine
1044 : !! mpw0=maximum number of planewaves on disk file
1045 : !! nban_dp_rd(nkpt0*nsppol0)=number of bands to be read at each k point
1046 : !! nband(nkpt*nsppol)=number of bands at each k point
1047 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1048 : !! nkassoc=dimension of indkk0 array
1049 : !! nkpt=number of k points expected
1050 : !! nkpt0=number of k points on disk
1051 : !! npwarr(nkpt)=array holding npw for each k point.
1052 : !! npwarr0(nkpt0)=array holding npw for each k point, disk format.
1053 : !! nspinor=number of spinorial components of the wavefunctions
1054 : !! nspinor0=number of spinorial components of the wavefunctions on disk
1055 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1056 : !! nsppol0=1 for unpolarized, 2 for spin-polarized, when on disk
1057 : !! nsym=number of symmetry elements in space group
1058 : !! optorth= 1 if the WFS have to be orthogonalized; 0 otherwise
1059 : !! prtvol=control print volume and debugging
1060 : !! randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility
1061 : !! restart= if 2, want conversion between wavefunctions
1062 : !! if 1, direct restart is allowed (see hdr_check.f)
1063 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1064 : !! sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn
1065 : !! if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries
1066 : !! squeeze=1 if cg_disk is to be used, even with mkmem/=0
1067 : !! symrel(3,3,nsym)=symmetry operations in real space in terms
1068 : !! of primitive translations
1069 : !! tnons(3,nsym)=nonsymmorphic translations for symmetry operations
1070 : !! wff1, structure information for input and output files
1071 : !! dtset <type(dataset_type)>=all input variables for this dataset
1072 : !!
1073 : !! OUTPUT
1074 : !! if ground state format (formeig=0):
1075 : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
1076 : !! if respfn format (formeig=1):
1077 : !! eigen(2*mband*mband*nkpt*nsppol)=
1078 : !! matrix of eigenvalues (input or init to large number), (Ha)
1079 : !! Conditional output:
1080 : !! cg_disk(2,mpw*nspinor*mband*mkmem*nsppol)=complex wf array
1081 : !! be careful : an array of size cg(2,npw*nspinor), as used
1082 : !! in the response function code, is not enough !
1083 : !!
1084 : !! SIDE EFFECTS
1085 : !! if ground state format (formeig=0):
1086 : !! occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
1087 : !! NOT OUTPUT NOW !
1088 : !!
1089 : !! NOTES
1090 : !! occ will not be modified nor output, in the present status of this routine.
1091 : !!
1092 : !! WARNINGS
1093 : !! For parallelism : no distinction yet between nban_dp_rd and nband
1094 : !!
1095 : !! TODO
1096 : !! THE DESCRIPTION IS TO BE COMPLETELY REVISED, AS THIS ONE COMES FROM inwffil.f
1097 : !!
1098 : !! SOURCE
1099 :
1100 2252 : subroutine wfsinp(cg,cg_disk,ecut,ecut0,ecut_eff,eigen,exchn2n3d,&
1101 2252 : & formeig,gmet,gmet0,headform0,indkk,indkk0,istwfk,&
1102 2252 : & istwfk0,kptns,kptns0,localrdwf,mband,&
1103 2252 : & mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,nband,nban_dp_rd,&
1104 2252 : & ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor,&
1105 2252 : & nspinor0,nsppol,nsppol0,nsym,occ,optorth,prtvol,randalg,restart,rprimd,&
1106 2252 : & sppoldbl,squeeze,symrel,tnons,wff1)
1107 :
1108 : !Arguments ------------------------------------
1109 : integer, intent(in) :: exchn2n3d,formeig,headform0,localrdwf,mband,mcg,mcg_disk
1110 : integer, intent(in) :: mpw,mpw0,nkassoc,nkpt,nkpt0,nspinor,nspinor0,nsppol,nsppol0,nsym
1111 : integer, intent(in) :: optorth,prtvol,randalg,restart,sppoldbl,squeeze
1112 : real(dp), intent(in) :: ecut,ecut0,ecut_eff
1113 : type(MPI_type), intent(inout) :: mpi_enreg,mpi_enreg0
1114 : type(wffile_type), intent(inout) :: wff1
1115 : integer, intent(in) :: indkk(nkpt*sppoldbl,6),indkk0(nkpt0,nkassoc),istwfk(nkpt)
1116 : integer, intent(in) :: istwfk0(nkpt0),nband(nkpt*nsppol),nban_dp_rd(nkpt0*nsppol0)
1117 : integer, intent(in) :: ngfft(18),npwarr(nkpt),npwarr0(nkpt0),symrel(3,3,nsym)
1118 : real(dp), intent(in) :: gmet(3,3),gmet0(3,3),kptns(3,nkpt),kptns0(3,nkpt0),rprimd(3,3)
1119 : real(dp), intent(in) :: tnons(3,nsym)
1120 : real(dp), intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol)
1121 : real(dp), intent(inout) :: cg(2,mcg),cg_disk(2,mcg_disk) !vz_i pw_ortho
1122 : real(dp), intent(inout) :: occ(mband*nkpt*nsppol)
1123 :
1124 : !Local variables-------------------------------
1125 : integer :: band_index,band_index_trial,ceksp,debug,dim_eig_k,iband,icg
1126 : integer :: icg_disk,icg_trial,idum,ierr,ii,ikassoc,ikassoc_trial,ikpt,ikpt0
1127 : integer :: ikpt10,ikpt_trial,ikptsp,ikptsp_old,inplace,isp,isp_max,isppol,isppol0
1128 : ! integer :: ipw ! commented out below
1129 : integer :: isppol_trial,me,mgfft,my_nspinor,my_nspinor0
1130 : integer :: nban_dp_k,nban_dp_rdk,nband_k,nband_rdk,nband_trial,nbd,nbd_max
1131 : integer :: ncopy,nkpt_eff,nproc_max,npw0_k,npw_k,npw_ktrial
1132 : integer :: read_cg,read_cg_disk,sender,spaceComm
1133 : character(len=500) :: msg
1134 2252 : integer,allocatable :: band_index_k(:,:),icg_k(:,:),kg0_k(:,:),kg_k(:,:)
1135 : real(dp) :: tsec(2)
1136 2252 : real(dp),allocatable :: eig0_k(:),eig_k(:),occ0_k(:),occ_k(:)
1137 : #if defined HAVE_MPI
1138 : integer :: iproc,my_ikpt
1139 : integer :: tag,test_cycle
1140 : integer :: statux(MPI_STATUS_SIZE)
1141 2252 : integer,allocatable :: ikassoc_me(:),ikpt_me(:),isppol_me(:),nband_k_me(:)
1142 : #endif
1143 : integer :: nkpt_max=50
1144 :
1145 : ! *************************************************************************
1146 :
1147 : !DEBUG
1148 : !write(std_out,*)' wfsinp : enter'
1149 : !write(std_out,*)' wfsinp : nband=',nband(:)
1150 : !write(std_out,*)' wfsinp : nban_dp_rd=',nban_dp_rd(:)
1151 : !write(std_out,*)' wfsinp : localrdwf=',localrdwf
1152 : !write(std_out,*)' wfsinp : paralbd,formeig=',mpi_enreg%paralbd,formeig
1153 : !write(std_out,*)' wfsinp : indkk0(:,1)=',indkk0(:,1)
1154 : !ENDDEBUG
1155 :
1156 2252 : call timab(720,1,tsec)
1157 2252 : call timab(721,3,tsec)
1158 :
1159 2252 : nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1
1160 2252 : nbd_max=size(mpi_enreg%proc_distrb,2)
1161 2252 : isp_max=size(mpi_enreg%proc_distrb,3)
1162 :
1163 : !Init mpi_comm
1164 2252 : spaceComm=mpi_enreg%comm_cell
1165 2252 : nproc_max=xmpi_comm_size(spaceComm)
1166 2252 : me=mpi_enreg%me_kpt
1167 2252 : sender = 0
1168 :
1169 : #if defined HAVE_MPI
1170 2252 : if(localrdwf==0)then
1171 21 : ABI_MALLOC(ikpt_me,(nproc_max))
1172 14 : ABI_MALLOC(nband_k_me,(nproc_max))
1173 14 : ABI_MALLOC(ikassoc_me,(nproc_max))
1174 14 : ABI_MALLOC(isppol_me,(nproc_max))
1175 : end if
1176 : #endif
1177 :
1178 : !Check the validity of formeig
1179 2252 : if(formeig/=0.and.formeig/=1)then
1180 0 : write(msg, '(a,i0,a)' )' formeig=',formeig,' , but the only allowed values are 0 or 1.'
1181 0 : ABI_BUG(msg)
1182 : end if
1183 :
1184 2252 : my_nspinor =max(1,nspinor /mpi_enreg%nproc_spinor)
1185 2252 : my_nspinor0=max(1,nspinor0/mpi_enreg%nproc_spinor)
1186 :
1187 2252 : nkpt_eff=max(nkpt0,nkpt)
1188 2252 : if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max
1189 :
1190 9008 : ABI_MALLOC(icg_k,(nkpt,nsppol))
1191 6756 : ABI_MALLOC(band_index_k,(nkpt,nsppol))
1192 :
1193 : !write(std_out,*)' wfsinp : me,isppol,ikpt,icg_k,band_index_k'
1194 :
1195 : !Compute the locations of the blocks in cg, eig and occ
1196 2252 : icg=0
1197 2252 : band_index=0
1198 2252 : ikpt10=0
1199 :
1200 5000 : do isppol=1,nsppol
1201 79757 : do ikpt=1,nkpt
1202 :
1203 74757 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1204 74757 : band_index_k(ikpt,isppol)=band_index
1205 :
1206 : #if defined HAVE_MPI
1207 74757 : test_cycle=0;nbd=min(nband_k,nbd_max);isp=min(isppol,isp_max)
1208 74757 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbd,isp,me))test_cycle=1
1209 : if(test_cycle==1)then
1210 2866 : band_index=band_index+nband_k*(2*nband_k)**formeig
1211 : ! In the case this k point does not belong to me, cycle
1212 2866 : cycle
1213 : end if
1214 : #endif
1215 :
1216 71891 : npw_k=npwarr(ikpt)
1217 71891 : icg_k(ikpt,isppol)=icg
1218 71891 : icg=icg+npw_k*my_nspinor*nband_k
1219 :
1220 74639 : band_index=band_index+nband_k*(2*nband_k)**formeig
1221 : ! write(std_out,'(5i8)' )me,isppol,ikpt,icg_k(ikpt,isppol),band_index_k(ikpt,isppol)
1222 : end do ! End k point loop
1223 : end do ! End spin loop
1224 :
1225 2252 : band_index=0
1226 2252 : ikptsp_old=0
1227 :
1228 : !DEBUG
1229 : !write(std_out,*)' wfsinp: before loop'
1230 : !write(std_out,*)' nsppol0,nsppol,nkpt0',nsppol0,nsppol,nkpt0
1231 : !write(std_out,*)' mpw,mgfft,mpw,mpw0',mpw,mgfft,mpw,mpw0
1232 : !ENDDEBUG
1233 :
1234 : mgfft=maxval(ngfft(1:3))
1235 2252 : if(squeeze==1)then
1236 258 : ABI_MALLOC(kg_k,(3,mpw))
1237 258 : ABI_MALLOC(kg0_k,(3,mpw0))
1238 : end if
1239 :
1240 5031008 : eigen(:)=0.0_dp
1241 : !occ(:)=0.0_dp
1242 :
1243 2252 : call timab(721,2,tsec)
1244 :
1245 : !Loop over spins
1246 : !For the time being, do not allow nsppol=2 to nspinor=2 conversion
1247 : !MT 20110707: this can be done by a fake call to the routine: see inwffil
1248 4999 : do isppol0=1,min(nsppol0,nsppol)
1249 :
1250 : ! Loop on k points : get the cg then eventually write on unwfnow
1251 71240 : do ikpt0=1,nkpt0
1252 :
1253 66241 : call timab(722,1,tsec)
1254 :
1255 66241 : nban_dp_rdk=nban_dp_rd(ikpt0+(isppol0-1)*nkpt0)
1256 :
1257 : ! DEBUG
1258 : ! write(std_out,*)' wfsinp: ikpt0,isppol0,nkpt0=',ikpt0,isppol0,nkpt0
1259 : ! write(std_out,*)' nban_dp_rdk=',nban_dp_rdk
1260 : ! ENDDEBUG
1261 :
1262 66241 : npw0_k=npwarr0(ikpt0)
1263 66241 : if(ikpt0<=nkpt_eff)then
1264 404 : write(msg,'(a,a,2i4)')ch10,' wfsinp: inside loop, init ikpt0,isppol0=',ikpt0,isppol0
1265 404 : call wrtout(std_out,msg)
1266 : end if
1267 :
1268 : ! Must know whether this k point is needed, and in which
1269 : ! block (ikpt, isppol), the wavefunction is to be placed.
1270 : ! Select the one for which the number of bands is the biggest.
1271 66241 : ikpt=0
1272 66241 : isppol=0
1273 66241 : ikassoc=0
1274 66241 : nband_k=0
1275 : #if defined HAVE_MPI
1276 66241 : if(localrdwf==0)then
1277 56 : nband_k_me(:)=0
1278 56 : ikpt_me(:)=0
1279 56 : isppol_me(:)=0
1280 56 : ikassoc_me(:)=0
1281 56 : nband_k_me(:)=0
1282 : end if
1283 : #endif
1284 :
1285 139716 : do isppol_trial=1,nsppol
1286 :
1287 73475 : if(nsppol==2 .and. nsppol0==2 .and. isppol0/=isppol_trial)cycle
1288 :
1289 426460 : do ikassoc_trial=1,nkassoc
1290 :
1291 293976 : ikpt_trial=indkk0(ikpt0,ikassoc_trial)
1292 293976 : if(sppoldbl==2)then
1293 0 : if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle
1294 0 : if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle
1295 0 : if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt
1296 : end if
1297 :
1298 : #if defined HAVE_MPI
1299 293976 : if(localrdwf==1)then
1300 293724 : if(ikpt_trial/=0)then
1301 74645 : nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1302 74645 : nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
1303 74645 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me))ikpt_trial=0
1304 : end if
1305 : end if
1306 : #endif
1307 :
1308 367451 : if(ikpt_trial/=0)then
1309 71955 : nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1310 71955 : if(nband_k<nband_trial)then
1311 63814 : nband_k=nband_trial ; ikpt=ikpt_trial ; isppol=isppol_trial
1312 63814 : ikassoc=ikassoc_trial
1313 : end if
1314 :
1315 : #if defined HAVE_MPI
1316 71955 : if(localrdwf==0)then
1317 448 : do iproc=1,nproc_max
1318 336 : my_ikpt=1
1319 336 : nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1320 336 : nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
1321 336 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,(iproc-1))) my_ikpt=0
1322 112 : if(my_ikpt/=0)then
1323 : nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1324 112 : if(nband_k_me(iproc)<nband_trial)then
1325 42 : nband_k_me(iproc)=nband_trial
1326 42 : ikpt_me(iproc)=ikpt_trial
1327 42 : isppol_me(iproc)=isppol_trial
1328 42 : ikassoc_me(iproc)=ikassoc_trial
1329 : end if
1330 : end if
1331 : end do
1332 : end if
1333 : #endif
1334 :
1335 : end if
1336 :
1337 : end do ! ikassoc_trial
1338 : end do ! isppol_trial
1339 :
1340 : ! DEBUG
1341 : ! write(std_out,*)' wfsinp : me,select isppol,ikpt=',me,isppol,ikpt
1342 : #if defined HAVE_MPI
1343 : ! write(std_out,*)' wfsinp : me,ikpt_me(:)=',me,ikpt_me(:)
1344 : #endif
1345 : ! write(std_out,*)' wfsinp : me,isppol_me(:)=',me,isppol_me(:)
1346 : ! stop
1347 : ! ENDDEBUG
1348 :
1349 66241 : call timab(722,2,tsec)
1350 :
1351 : ! If the wavefunction block to be read is interesting ...
1352 68988 : if (ikpt/=0)then
1353 :
1354 63814 : call timab(723,3,tsec)
1355 63814 : sender = me
1356 63814 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1357 63814 : npw_k=npwarr(ikpt)
1358 :
1359 : #if defined HAVE_MPI
1360 63814 : if (localrdwf==1.or.(localrdwf==0.and.me==0)) then
1361 :
1362 63806 : if(ikpt<=nkpt_eff)then
1363 : write(msg,'(a,i6,a,i8,a,i4,a,i4)') &
1364 299 : & ' wfsinp: treating ',nband_k,' bands with npw=',npw_k,' for ikpt=',ikpt,' by node ',me
1365 299 : call wrtout(std_out,msg)
1366 63507 : else if(ikpt==nkpt_eff+1)then
1367 0 : call wrtout(std_out,' wfsinp: prtvol=0 or 1, do not print more k-points.')
1368 : end if
1369 :
1370 : end if
1371 : #endif
1372 :
1373 63814 : nband_rdk=nban_dp_rdk
1374 63814 : if(squeeze==1)nband_rdk=(nban_dp_rdk/nspinor0)*nspinor
1375 63814 : if(formeig==0)then
1376 28815 : ABI_MALLOC(eig_k,(nband_rdk))
1377 19210 : ABI_MALLOC(occ_k,(nband_rdk))
1378 28815 : ABI_MALLOC(eig0_k,(nban_dp_rdk))
1379 19210 : ABI_MALLOC(occ0_k,(nban_dp_rdk))
1380 63814 : dim_eig_k=nband_rdk
1381 54209 : else if(formeig==1)then
1382 162627 : ABI_MALLOC(eig_k,(2*nband_rdk*nband_rdk))
1383 162627 : ABI_MALLOC(eig0_k,(2*nban_dp_rdk*nban_dp_rdk))
1384 54209 : dim_eig_k=2*nband_rdk*nband_rdk
1385 54209 : ABI_MALLOC(occ0_k,(0))
1386 54209 : ABI_MALLOC(occ_k,(0))
1387 : else
1388 0 : ABI_MALLOC(occ0_k,(0))
1389 0 : ABI_MALLOC(eig0_k,(0))
1390 0 : ABI_MALLOC(occ_k,(0))
1391 63814 : ABI_MALLOC(eig_k,(0))
1392 : end if
1393 4993414 : eig_k(:)=0.0_dp
1394 4993414 : eig0_k(:)=0.0_dp
1395 :
1396 : ! Generate or read the cg for this k point
1397 : ! Either read into cg, or read into cg_disk
1398 63814 : read_cg=1 ; read_cg_disk=0
1399 63814 : if(squeeze==1)then
1400 999 : read_cg=0 ; read_cg_disk=1
1401 : end if
1402 : #if defined HAVE_MPI
1403 63814 : if(localrdwf==0)then
1404 14 : read_cg=0
1405 14 : read_cg_disk=0
1406 : ! XG20040106 The following condition is correct
1407 14 : if(me==0)read_cg_disk=1
1408 : end if
1409 : #endif
1410 :
1411 14 : icg=0
1412 63800 : if(read_cg==1)icg=icg_k(ikpt,isppol)
1413 :
1414 : ! DEBUG
1415 : ! write(std_out,*)' wfsinp: before initwf',wff1%offwff
1416 : ! write(std_out,*)' wfsinp: me,read_cg,read_cg_disk=',me,read_cg,read_cg_disk
1417 : ! write(std_out,*)' wfsinp: nban_dp_rdk=',nban_dp_rdk
1418 : ! ENDDEBUG
1419 63814 : call timab(723,2,tsec)
1420 63814 : call timab(724,3,tsec)
1421 :
1422 63814 : if(read_cg_disk==1)then
1423 : call initwf (cg_disk,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
1424 : & isppol0,mcg_disk,mpi_enreg0, &
1425 1005 : & nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1)
1426 : end if
1427 :
1428 63814 : if(read_cg==1)then
1429 : call initwf (cg,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,&
1430 : & isppol0,mcg,mpi_enreg0,&
1431 62801 : & nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1)
1432 : end if
1433 :
1434 63814 : call timab(724,2,tsec)
1435 63814 : call timab(725,3,tsec)
1436 :
1437 63814 : nban_dp_k=min(nban_dp_rdk,(nband_k/nspinor)*nspinor0)
1438 : ! This band_index is defined BEFORE the eventual redefinition
1439 : ! of ikpt and isppol, needed when localrdwf==0 in parallel
1440 63814 : band_index=band_index_k(ikpt,isppol)
1441 : ! DEBUG
1442 : ! write(std_out,*)' wfsinp: me,cg_disk(:,1)=',me,cg_disk(:,1)
1443 : ! ENDDEBUG
1444 :
1445 : ! DEBUG
1446 : ! if(me==0 .and. ikpt0==1)then
1447 : ! write(std_out,*)' wfsinp : cg array, before trial, ikpt0=',ikpt0
1448 : ! do ipw=1,15
1449 : ! write(std_out,'(i4,2es20.10)' )ipw,cg(:,ipw)
1450 : ! end do
1451 : ! end if
1452 : ! ENDDEBUG
1453 :
1454 :
1455 : #if defined HAVE_MPI
1456 63814 : if(localrdwf==0)then
1457 : ! Warning: In that case , not yet // on nspinors
1458 : ! Transmit to each of the other processors, when needed
1459 14 : if(nproc_max>=2)then
1460 40 : do iproc=2,nproc_max
1461 : ! Only me=0 and me=iproc-1 are concerned by this
1462 40 : if(me==0 .or. me==iproc-1)then
1463 :
1464 16 : ikpt=ikpt_me(iproc)
1465 16 : isppol=isppol_me(iproc)
1466 :
1467 16 : if(ikpt/=0)then
1468 : ! In this case, processor iproc-1 needs the data
1469 : ! Generate a common tag
1470 16 : tag=256*(ikpt-1)+iproc+1
1471 16 : if(isppol==2)tag=-tag
1472 16 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1473 16 : npw_k=npwarr(ikpt)
1474 : ! SEND
1475 16 : if(me==0)then
1476 8 : write(std_out,*)'SENDWFSINP ',me
1477 : call MPI_SEND(cg_disk,2*npw_k*my_nspinor*nband_k,&
1478 8 : & MPI_DOUBLE_PRECISION,iproc-1,tag,spaceComm,ierr)
1479 : end if
1480 : ! RECEIVE
1481 16 : if(me==iproc-1)then
1482 : call MPI_RECV(cg_disk,2*npw_k*my_nspinor*nband_k,&
1483 8 : & MPI_DOUBLE_PRECISION,0,tag,spaceComm,statux,ierr)
1484 8 : icg=icg_k(ikpt,isppol)
1485 8 : if(squeeze==0)then
1486 : cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)=&
1487 1448 : & cg_disk(:,1:npw_k*my_nspinor*nband_k)
1488 : end if
1489 8 : ikassoc=ikassoc_me(iproc)
1490 : end if
1491 : end if
1492 :
1493 : end if
1494 : end do ! iproc
1495 : end if
1496 :
1497 : ! DEBUG
1498 : ! write(std_out,*)' wfsinp: me, iproc loop finished',me
1499 : ! ENDDEBUG
1500 :
1501 : ! Take care of me=0 needing the data
1502 14 : if (me==0) then
1503 6 : ikpt=ikpt_me(me+1)
1504 6 : isppol=isppol_me(me+1)
1505 6 : if(ikpt/=0 )then
1506 6 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1507 6 : npw_k=npwarr(ikpt)
1508 : ! I am the master node, and I might need my own data
1509 6 : icg=icg_k(ikpt,isppol)
1510 6 : if(squeeze==0)then
1511 : ! Copy from cg_disk to cg
1512 : cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)= &
1513 1086 : & cg_disk(:,1:npw_k*my_nspinor*nband_k)
1514 : end if
1515 6 : ikassoc=ikassoc_me(me+1)
1516 : end if
1517 : end if
1518 : ! For the eigenvalues and occ, the transmission is much easier to write !
1519 : call MPI_BCAST(eig0_k,nban_dp_rdk*(2*nban_dp_rdk)**formeig ,&
1520 14 : & MPI_DOUBLE_PRECISION,0,spaceComm,ierr)
1521 : end if
1522 : #endif
1523 :
1524 63814 : if(formeig==0)then
1525 : ! The transfer from eig0_k to eig_k uses nban_dp_rdk, which contains
1526 : ! the maximal information.
1527 9605 : if(nspinor0==nspinor .or. squeeze==0)then
1528 82677 : eig_k(1:nban_dp_rdk)=eig0_k(1:nban_dp_rdk)
1529 82677 : occ_k(1:nban_dp_rdk)=occ0_k(1:nban_dp_rdk)
1530 0 : else if(nspinor0==1 .and. nspinor==2)then
1531 0 : do iband=1,nban_dp_rdk
1532 0 : eig_k(2*iband )=eig0_k(iband)
1533 0 : eig_k(2*iband-1)=eig0_k(iband)
1534 0 : occ_k(2*iband )=occ0_k(iband)*0.5_dp
1535 0 : occ_k(2*iband-1)=occ0_k(iband)*0.5_dp
1536 : end do
1537 0 : else if(nspinor0==2 .and. nspinor==1)then
1538 0 : do iband=1,nban_dp_rdk
1539 0 : eig_k(iband)=eig0_k(2*iband-1)
1540 0 : occ_k(iband)=occ0_k(2*iband-1)*2.0_dp
1541 : end do
1542 : end if
1543 :
1544 : ! DEBUG
1545 : ! write(std_out,*)' wfsinp: me,band_index,ikpt,isppol',me,band_index,ikpt,isppol
1546 : ! ENDDEBUG
1547 :
1548 : ! The transfer to eigen uses nban_dp_k, that is bound by the number
1549 : ! of bands for this k point.
1550 9605 : ncopy=min(dim_eig_k,(nban_dp_k/nspinor0)*nspinor)
1551 82677 : eigen(1+band_index:ncopy+band_index)=eig_k(1:ncopy)
1552 : ! The transfer of occ should be done.
1553 :
1554 : #if defined HAVE_MPI
1555 9605 : if(localrdwf==0 .and. ikpt/=0)then
1556 : ! Do not forget : ikpt,isppol were redefined ...
1557 14 : band_index=band_index_k(ikpt,isppol)
1558 70 : eigen(1+band_index:(nban_dp_k/nspinor0)*nspinor+band_index) = eig_k(1:(nban_dp_k/nspinor0)*nspinor)
1559 : ! The transfer of occ should be done.
1560 : end if
1561 : #endif
1562 :
1563 54209 : else if(formeig==1)then
1564 54209 : call wrtout(std_out,'wfsinp: transfer of first-order eigs not yet coded!',"COLL")
1565 : end if
1566 :
1567 : ! DEBUG
1568 : ! write(std_out,*)' wfsinp : me,transferred eig_k',me
1569 : ! write(std_out,*)' me,mkmem,nsppol,nsppol0,isppol',me,mkmem,nsppol,nsppol0,isppol
1570 : ! write(std_out,*)' me,nkassoc,ikassoc',me,nkassoc,ikassoc
1571 : ! ENDDEBUG
1572 :
1573 63814 : call timab(725,2,tsec)
1574 :
1575 : ! Write to disk if appropriate
1576 : ! The coding has to be done ... here, only fragments ...
1577 63814 : call timab(727,3,tsec)
1578 :
1579 134830 : do isppol_trial=1,nsppol
1580 71016 : if(nsppol==2 .and. nsppol0==2 .and. isppol_trial/=isppol)cycle
1581 404198 : do ikassoc_trial=1,nkassoc
1582 :
1583 : ! DEBUG
1584 : ! write(std_out,*)' wfsinp: me, for ikassoc,isppol',&
1585 : ! & me,ikassoc,isppol
1586 : ! write(std_out,*)' wfsinp: me, try ikassoc_trial,isppol_trial,nband_k',&
1587 : ! & me,ikassoc_trial,isppol_trial,nband_k
1588 : ! ENDDEBUG
1589 :
1590 : ! No conversion is to be done : it will be converted in newkpt
1591 : ! If squeeze==0, the block with the ikpt corresponding to ikassoc,
1592 : ! and with isppol, contains the wavefunction already
1593 : if( ikassoc_trial/=ikassoc .or. &
1594 276568 : & (isppol_trial/=isppol .and. sppoldbl==1 ).or. &
1595 71016 : & squeeze==1 )then
1596 :
1597 213753 : ikpt_trial=indkk0(ikpt0,ikassoc_trial)
1598 213753 : if(sppoldbl==2)then
1599 0 : if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle
1600 0 : if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle
1601 0 : if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt
1602 : end if
1603 :
1604 :
1605 : #if defined HAVE_MPI
1606 213753 : if(ikpt_trial/=0)then
1607 9729 : nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1608 9729 : nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max)
1609 9729 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me)) ikpt_trial=0
1610 : end if
1611 : #endif
1612 :
1613 213753 : if(ikpt_trial/=0 .and. ikpt_trial<=nkpt_eff)then
1614 171 : write(msg,'(2a,2i5)')ch10,' wfsinp: transfer to ikpt_trial,isppol_trial=',ikpt_trial,isppol_trial
1615 171 : call wrtout(std_out,msg)
1616 : end if
1617 :
1618 213753 : if(ikpt_trial/=0)then
1619 9076 : icg_trial=icg_k(ikpt_trial,isppol_trial)
1620 9076 : band_index_trial=band_index_k(ikpt_trial,isppol_trial)
1621 9076 : nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt)
1622 9076 : nban_dp_k=min(nban_dp_rdk,(nband_trial/nspinor)*nspinor0)
1623 :
1624 9076 : if(squeeze==0)then
1625 : ! GMR: modified to avoid compiler bug
1626 : ! cg(:,1+icg_trial:npw0_k*my_nspinor0*nband_trial+icg_trial)=&
1627 : ! & cg(:,1+icg:npw0_k*my_nspinor0*nband_trial+icg)
1628 3692046 : do ii=1,npw0_k*my_nspinor0*nband_trial
1629 11065714 : cg(:,ii+icg_trial)=cg(:,ii+icg)
1630 : end do
1631 : ! GMR
1632 5212 : if(formeig==0)then
1633 39575 : eigen(1+band_index_trial:nban_dp_k+band_index_trial)=eig_k(1:nban_dp_k)
1634 : ! occ(1+band_index_trial:nban_dp_k+band_index_trial)=&
1635 : ! & occ_k(1:nban_dp_k)
1636 : end if
1637 : ! RF transfer of eigenvalues still to be coded
1638 3864 : else if(squeeze==1)then
1639 3864 : npw_ktrial=npwarr(ikpt_trial)
1640 3864 : nband_k=(nban_dp_k/nspinor0)*nspinor
1641 : ! Conversion to be done
1642 3864 : ceksp=0 ; debug=0 ; icg_disk=0 ; idum=0 ; inplace=0
1643 : ! Note that this routine also convert eig and occ
1644 : ! even if the conversion had already been done
1645 :
1646 :
1647 : call wfconv(ceksp,cg_disk,cg,debug,ecut0,ecut,ecut_eff,&
1648 : & eig0_k,eig_k,exchn2n3d,formeig,gmet0,gmet,&
1649 : & icg_disk,icg_trial,ikpt0,ikpt10,ikpt_trial,indkk,&
1650 : & inplace,isppol_trial,istwfk0,istwfk,&
1651 : & kg0_k,kg_k,kptns0,kptns,nban_dp_rdk,nband_rdk,&
1652 : & mcg_disk,mcg,mpi_enreg0,mpi_enreg,mpw0,mpw,&
1653 : & nban_dp_rdk,nband_trial,ngfft,ngfft,nkpt0,nkpt,&
1654 : & npw0_k,npw_ktrial,nspinor0,nspinor,nsym,&
1655 : & occ0_k,occ_k,optorth,randalg,restart,rprimd,&
1656 3864 : & sppoldbl,symrel,tnons)
1657 :
1658 : ! DEBUG
1659 : ! write(std_out,*)' wfsinp: ikpt_trial=',ikpt_trial
1660 : ! write(std_out,*)' nband_k,band_index_trial',nband_k,band_index_trial
1661 : ! write(std_out,*)eig0_k(1:nban_dp_k)
1662 : ! write(std_out,*)eig_k(1:nband_k)
1663 : ! ENDDEBUG
1664 37278 : eigen(1+band_index_trial:nband_k+band_index_trial)=eig_k(1:nband_k)
1665 : ! occ(1+band_index_trial:nband_k+band_index_trial)=&
1666 : ! & occ_k(1:nband_k)
1667 : ! RF transfer of eigenvalues still to be coded
1668 :
1669 : ! Endif squeeze==1
1670 : end if
1671 :
1672 : ! DEBUG
1673 : ! if(ikpt_trial==2)then
1674 : ! write(std_out,*)' wfsinp: iband,ipw,cg for ikpt_trial=2'
1675 : ! write(std_out,*)' nband_trial,npw_ktrial=',nband_trial,npw_ktrial
1676 : ! do iband=1,nband_trial
1677 : ! do ipw=1,npw_ktrial
1678 : ! write(std_out,'(2i5,2es16.6)' )&
1679 : ! & iband,ipw,cg(:,ipw+(iband-1)*npw_ktrial+icg_trial)
1680 : ! end do
1681 : ! end do
1682 : ! end if
1683 : ! ENDDEBUG
1684 :
1685 : ! End if ikpt_trial/=0
1686 : end if
1687 :
1688 : ! End if ikpt_trial already initialized
1689 : end if
1690 :
1691 : end do ! ikassoc_trial
1692 : end do ! isppol_trial
1693 :
1694 63814 : call timab(727,2,tsec)
1695 :
1696 63814 : ABI_FREE(eig_k)
1697 63814 : ABI_FREE(eig0_k)
1698 : !if(formeig==0) then
1699 63814 : ABI_FREE(occ_k)
1700 63814 : ABI_FREE(occ0_k)
1701 : !end if
1702 :
1703 : end if ! End the condition of need of this k point
1704 : end do ! End of the k loop
1705 : end do ! End of spin loop
1706 :
1707 : #if defined HAVE_MPI
1708 2252 : call timab(67,1,tsec)
1709 :
1710 : !Still need to skip last k points to read history (MS)
1711 : !WARNING : not yet for formeig=1, but is it needed ?
1712 2252 : if(formeig==0)then
1713 2048 : do ikptsp=ikptsp_old+1,nkpt0*nsppol0
1714 725 : isppol=1 ; if(ikptsp>nkpt0)isppol=2
1715 725 : ikpt=ikptsp-nkpt0*(isppol-1)
1716 2048 : call WffReadSkipK(formeig,headform0,ikpt,isppol,mpi_enreg,wff1)
1717 : end do
1718 : end if
1719 :
1720 : !Transmit eigenvalues. This routine works in both localrdwf=0 or 1 cases.
1721 2252 : call pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,1)
1722 :
1723 :
1724 2252 : if(localrdwf==0)then
1725 7 : ABI_FREE(ikpt_me)
1726 7 : ABI_FREE(nband_k_me)
1727 7 : ABI_FREE(ikassoc_me)
1728 7 : ABI_FREE(isppol_me)
1729 : end if
1730 :
1731 2252 : call timab(67,2,tsec)
1732 : #endif
1733 :
1734 : !****************************************************************************
1735 :
1736 2252 : if(squeeze==1)then
1737 86 : ABI_FREE(kg_k)
1738 86 : ABI_FREE(kg0_k)
1739 : end if
1740 :
1741 :
1742 : !DEBUG
1743 : !if(me==0)then
1744 : !write(std_out,*)' wfsinp : cg array='
1745 : !icg=0
1746 : !do isppol=1,nsppol
1747 : !do ikpt=1,1
1748 : !nband_k=nband(ikpt+(isppol-1)*nkpt)
1749 : !npw_k=npwarr(ikpt)
1750 : !do iband=1,nband_k
1751 : !write(std_out,*)' new band, icg=',icg
1752 : !do ipw=1,npw_k
1753 : !write(std_out,'(4i4,2es20.10)' )isppol,ikpt,iband,ipw,cg(:,icg+ipw)
1754 : !end do
1755 : !icg=icg+npw_k
1756 : !end do
1757 : !end do
1758 : !end do
1759 : !end if
1760 : !if(ireadwf==1)stop
1761 : !write(std_out,*)' wfsinp : eigen array='
1762 : !do ikpt=1,nkpt
1763 : !do iband=1,mband
1764 : !write(std_out,*)'ikpt,iband,eigen',ikpt,iband,eigen(iband+(ikpt-1)*mband)
1765 : !end do
1766 : !end do
1767 : !ENDDEBUG
1768 :
1769 2252 : ABI_FREE(icg_k)
1770 2252 : ABI_FREE(band_index_k)
1771 :
1772 2252 : call timab(720,2,tsec)
1773 :
1774 2252 : end subroutine wfsinp
1775 : !!***
1776 :
1777 : !!****f* m_inwffil/initwf
1778 : !!
1779 : !! NAME
1780 : !! initwf
1781 : !!
1782 : !! FUNCTION
1783 : !! Initialization of wavefunctions.
1784 : !! If formeig==1, and partially filled case, I am not sure that the eig_k are initialized properly ...
1785 : !! formeig option (format of the eigenvalues and eigenvector) :
1786 : !! 0 => ground-state format (initialisation of
1787 : !! eigenvectors with random numbers, vector of eigenvalues)
1788 : !! 1 => respfn format (initialisation of
1789 : !! eigenvectors with 0 s, hermitian matrix of eigenvalues)
1790 : !!
1791 : !! INPUTS
1792 : !! formeig=see above
1793 : !! headform=header format (might be needed to read the block of wfs)
1794 : !! icg=shift to be given to the location of the data in the array cg
1795 : !! ikpt= number of the k point of which the wf is initialised
1796 : !! spin=spin index
1797 : !! mcg=dimension of the cg array
1798 : !! mpi_enreg=information about MPI parallelization
1799 : !! nband_k=number of bands at this particular k point
1800 : !! nkpt=number of k points
1801 : !! npw=number of plane waves
1802 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
1803 : !! wff1=structure info for file containing wavefunctions (when needed)
1804 : !!
1805 : !! OUTPUT
1806 : !! cg(2,mcg)=complex wf array
1807 : !! if ground state format (formeig=0):
1808 : !! eig_k(nband_k)=list of eigenvalues (input or init to large number), hartree
1809 : !! if respfn format (formeig=1):
1810 : !! eig_k(2*nband_k*nband_k)= matrix of eigenvalues (input or init to large number), hartree
1811 : !!
1812 : !! SIDE EFFECTS
1813 : !! Input/output:
1814 : !! occ_k(nband_k)=list of occupations (input or left to their initial value)
1815 : !! ikptsp_old=number of the previous spin-k point, or 0 if first call of present file
1816 : !!
1817 : !! SOURCE
1818 :
1819 63806 : subroutine initwf(cg,eig_k,formeig,headform,icg,ikpt,ikptsp_old,&
1820 63806 : & spin,mcg,mpi_enreg,nband_k,nkpt,npw,nspinor,occ_k,wff1)
1821 :
1822 : !Arguments ------------------------------------
1823 : !scalars
1824 : integer,intent(in) :: formeig,headform,icg,ikpt,spin,mcg,nband_k,nkpt,npw,nspinor
1825 : integer,intent(inout) :: ikptsp_old
1826 : type(MPI_type),intent(in) :: mpi_enreg
1827 : type(wffile_type),intent(inout) :: wff1
1828 : !arrays
1829 : real(dp),intent(inout) :: occ_k(nband_k)
1830 : real(dp),intent(inout) :: cg(2,mcg),eig_k((2*nband_k)**formeig*nband_k) !vz_i
1831 :
1832 : !Local variables-------------------------------
1833 : !scalars
1834 : integer,parameter :: nkpt_max=50
1835 : integer :: ikpt0,nband_disk,tim_rwwf
1836 : character(len=500) :: msg
1837 : !arrays
1838 63806 : integer,allocatable :: kg_dum(:,:)
1839 : real(dp) :: tsec(2)
1840 : #if 0
1841 : integer :: iomode,comm,funt,ierr
1842 : type(wfk_t) :: Wfk
1843 : #endif
1844 :
1845 : ! *************************************************************************
1846 :
1847 : !write(std_out,*)' initwf : enter, ikptsp_old,ikpt,spin,nkpt= ',ikptsp_old,ikpt,spin,nkpt
1848 : !stop
1849 :
1850 : #if 0
1851 : ABI_WARNING("Entering new IO section")
1852 : !call WffClose(wff1,ierr)
1853 : comm = MPI_enreg%comm_cell
1854 : iomode = iomode_from_fname(wff1%fname)
1855 : call wfk_open_read(Wfk,wff1%fname,formeig,iomode,get_unit(),comm)
1856 : call wfk_read_band_block(Wfk,(/1,nband_k/),ikpt,spin,xmpio_at,cg_k=cg(1:,icg:),eig_k=eig_k,occ_k=occ_k)
1857 : call wfk_close(Wfk)
1858 : !call clsopn(wff1)
1859 : RETURN
1860 : #endif
1861 :
1862 63806 : call timab(770,1,tsec)
1863 63806 : call timab(771,1,tsec)
1864 :
1865 63806 : ABI_MALLOC(kg_dum,(3,0))
1866 :
1867 : !Skip wavefunctions for k-points not treated by this proc.
1868 : !(from ikptsp_old+1 to ikpt+(spin-1)*nkpt-1)
1869 63806 : if (ikptsp_old<ikpt+(spin-1)*nkpt-1) then
1870 1626 : do ikpt0=ikptsp_old+1,ikpt+(spin-1)*nkpt-1
1871 1626 : call WffReadSkipK(formeig,headform,ikpt0,spin,mpi_enreg,wff1)
1872 : end do
1873 : end if
1874 :
1875 : !write(std_out,*)' initwf : before rwwf'
1876 : !write(std_out,*)' formeig,icg,ikpt,spin=',formeig,icg,ikpt,spin
1877 : !write(std_out,*)' nband_k,nband_disk,npw,nspinor=',nband_k,nband_disk,npw,nspinor
1878 : !write(std_out,*)' unwff1=',unwff1
1879 : !stop
1880 :
1881 63806 : if(mpi_enreg%paralbd==0)tim_rwwf=2
1882 63806 : if(mpi_enreg%paralbd==1)tim_rwwf=20
1883 :
1884 63806 : call timab(771,2,tsec)
1885 :
1886 : call rwwf(cg,eig_k,formeig,headform,icg,ikpt,spin,kg_dum,nband_k,mcg,mpi_enreg,nband_k,nband_disk,&
1887 63806 : & npw,nspinor,occ_k,1,0,tim_rwwf,wff1)
1888 :
1889 63806 : call timab(772,1,tsec)
1890 :
1891 63806 : if (ikpt<=nkpt_max) then
1892 29698 : write(msg,'(3(a,i0))')' initwf: disk file gives npw= ',npw,' nband= ',nband_disk,' for kpt number= ',ikpt
1893 29698 : call wrtout(std_out,msg)
1894 34108 : else if (ikpt==nkpt_max+1) then
1895 280 : call wrtout(std_out,' initwf: the number of similar msg is sufficient... stop printing them')
1896 : end if
1897 :
1898 : ! Check the number of bands on disk file against desired number. These are not required to agree)
1899 63806 : if (nband_disk /= nband_k .and. ikpt<=nkpt_max) then
1900 : write(msg,'(3(a,i0),3a)')&
1901 72 : 'For kpt number: ',ikpt,' disk file has: ',nband_disk,' bands but input file gave nband: ',nband_k,'.',ch10,&
1902 144 : 'This is not fatal. Bands are skipped or filled with random numbers.'
1903 72 : ABI_COMMENT(msg)
1904 : end if
1905 :
1906 63806 : if (ikpt<=nkpt_max) then
1907 29698 : write(msg,'(a,i0,a)')' initwf: ',nband_disk,' bands have been initialized from disk'
1908 29698 : call wrtout(std_out,msg)
1909 : end if
1910 :
1911 63806 : ikptsp_old=ikpt+(spin-1)*nkpt
1912 :
1913 63806 : ABI_FREE(kg_dum)
1914 :
1915 63806 : call timab(772,2,tsec)
1916 63806 : call timab(770,2,tsec)
1917 :
1918 63806 : end subroutine initwf
1919 : !!***
1920 :
1921 : !!****f* m_inwffil/newkpt
1922 : !! NAME
1923 : !! newkpt
1924 : !!
1925 : !! FUNCTION
1926 : !! This subroutine writes a starting guess for wave function (set 2)
1927 : !! It performs a "zero order" interpolation, ie simply
1928 : !! searches the nearest available k-point.
1929 : !! The data (set 1) associated with this point is either
1930 : !! read from a disk file (with a random access reading routine),
1931 : !! or input as argument.
1932 : !!
1933 : !! INPUTS
1934 : !! ceksp2=if 1, center the sphere of pw on Gamma; if 0, on each k-point.
1935 : !! doorth=1 to do orthogonalization
1936 : !! debug=>0 for debugging output
1937 : !! ecut1=kinetic energy cutoffs for basis sphere 1 (hartree)
1938 : !! ecut2=kinetic energy cutoffs beyond which the coefficients of wf2 vanish (Ha)
1939 : !! ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree)
1940 : !! exchn2n3d=if 1, n2 and n3 are exchanged
1941 : !! fill=if 1, fill the supplementary bands ; if 0, reduce the number of bands
1942 : !! Note : must have fill/=0 in the parallel execution
1943 : !! formeig=if 0, GS format for wfs, eig and occ ; if 1, RF format.
1944 : !! gmet1(3,3), gmet2(3,3)=reciprocal space metrics (bohr^-2)
1945 : !! headform1=header format (might be needed to read the block of wfs)
1946 : !! indkk(nkpt2*sppoldbl,6)=describe k point number of kptns1 that allows to
1947 : !! generate wavefunctions closest to given kpt2 (and possibly isppol2=2)
1948 : !! indkk(:,1)=k point number of kpt1
1949 : !! indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
1950 : !! (if 0, means no symmetry operation, equivalent to identity )
1951 : !! indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
1952 : !! to give kpt1b, that is the closest to ikpt2.
1953 : !! indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
1954 : !! iout=unit number for output file
1955 : !! ireadwf=if 0, no reading of disk wavefunction file (random or 0.0 initialisation)
1956 : !! istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1
1957 : !! istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2
1958 : !! kg2(3,mpw2*mkmem2)=dimensionless coords of G vecs in basis sphere at k point
1959 : !! kptns1(3,nkpt1), kptns2(3,nkpt2)=k point sets (reduced coordinates)
1960 : !! mband2= maximum number of bands of the output wavefunctions
1961 : !! mcg=dimension of the cg array
1962 : !! In case mkmem2/=0, all the output data must find their place in cg,
1963 : !! so that mcg must be at least Sum(ikpt,isppol) [npw*nspinor*nband](ikpt,isppol)
1964 : !! where these data are related to the output parameters
1965 : !! In case mkmem1/=0, the same is true, for the input parameters,
1966 : !! however, the maximum number of bands that will be read
1967 : !! will be at most (mband2/nspinor2)*nspinor1
1968 : !! In case mkmem1==0 and mkmem2==0, one must have at least mpw*nspinor*mband
1969 : !! for BOTH the input and output parameters, taking into account the
1970 : !! maximal number of band to be read, described above.
1971 : !! In case mkmem1/=0 and mkmem2/=0, it is expected that the input cg array
1972 : !! is organised using the output parameters nkpt2, nband2 ...
1973 : !! This is needed, in order to use the same pointer.
1974 : !! mkmem1= if 0, the input wf, eig, occ are available from disk
1975 : !! mkmem2= if 0, the output wf, eig, occ must be written onto disk
1976 : !! mpi_enreg1=information about MPI parallelization, for the input wf file
1977 : !! mpi_enreg2=information about MPI parallelization, for the output wf file
1978 : !! mpw1=maximum allowed number of planewaves at any k, for the input wf file
1979 : !! mpw2=maximum allowed number of planewaves at any k, for the output wf file
1980 : !! my_nkpt2= number of k points for the output wf file, handled by current processus
1981 : !! nband1(nkpt1*nsppol1)=number of bands, at each k point, on disk
1982 : !! nband2(nkpt2*nsppol2)=desired number of bands at each k point
1983 : !! ngfft1(18)=all needed information about 3D FFT, for the input wf file
1984 : !! ngfft2(18)=all needed information about 3D FFT, for the output wf file
1985 : !! see ~abinit/doc/variables/vargs.htm#ngfft
1986 : !! nkpt1, nkpt2=number of k points in each set
1987 : !! npwarr1(nkpt1)=array holding npw for each k point (input wf file).
1988 : !! npwarr2(nkpt2)=array holding npw for each k point (output wf file).
1989 : !! nspinor1,nspinor2=number of spinorial components of the wavefunctions
1990 : !! for each wf file (input or output)
1991 : !! nsppol1=1 for unpolarized, 2 for spin-polarized, input wf file
1992 : !! nsppol2=1 for unpolarized, 2 for spin-polarized, output wf file
1993 : !! nsym=number of symmetry elements in space group
1994 : !! optorth= 1 if the WFS have to be orthogonalized; 0 otherwise
1995 : !! prtvol=control print volume and debugging
1996 : !! randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility
1997 : !! restart= if 2, conversion between wavefunctions
1998 : !! if 1, direct restart is allowed (see hdr_check.f)
1999 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
2000 : !! sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn
2001 : !! if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries
2002 : !! symrel(3,3,nsym)=symmetry operations in real space in terms
2003 : !! of primitive translations
2004 : !! tnons(3,nsym)=nonsymmorphic translations for symmetry operations
2005 : !! unkg2=unit number for storage of basis sphere data: stores indirect
2006 : !! indexing array and integer coordinates for all planewaves in basis
2007 : !! sphere for each k point being considered (kptns2 set)
2008 : !! wffinp=structure info of input wf file unit number
2009 : !! wffout=structure info of output wf file unit number
2010 : !! dtset <type(dataset_type)>=all input variables for this dataset
2011 : !!
2012 : !! OUTPUT
2013 : !! (see side effects)
2014 : !!
2015 : !! SIDE EFFECTS
2016 : !! The following arrays are input if mkmem1/=0, otherwise their input
2017 : !! values are taken from disk, and are output if mkmem2/=0, otherwise
2018 : !! their output values are written on disk.
2019 : !! The location of the block for a given spin-k point at input MUST
2020 : !! be the same as the location of the corresponding spin-k point at output.
2021 : !! cg(2,mcg)=complex wf array
2022 : !! eigen(mband2*(2*mband2)**formeig *nkpt2*nsppol2)=
2023 : !! eigenvalues (input or init to large number for GS or init to 0.0 for RF), (Ha)
2024 : !! occ(mband2*nkpt2*nsppol2)=occupation (input or init to 0.0) NOT USED NOW
2025 : !!
2026 : !! NOTES
2027 : !! * When reading from disk, it is expected that the next record of
2028 : !! the wffinp%unwff disk unit is the first record of the first wavefunction block.
2029 : !!
2030 : !! * When the data is input as argument, it is assumed that the
2031 : !! data for each spin- k wavefunction block is located at the proper
2032 : !! corresponding location of the output array (this is to be described).
2033 : !!
2034 : !! * The information is pumped onto an fft box for the conversion.
2035 : !! This allows for changing the number of plane waves.
2036 : !!
2037 : !! * In the present status of this routine, occ is not output.
2038 : !!
2039 : !! SOURCE
2040 :
2041 6619 : subroutine newkpt(ceksp2,cg,debug,ecut1,ecut2,ecut2_eff,eigen,exchn2n3d,fill,&
2042 6619 : & formeig,gmet1,gmet2,headform1,indkk,iout,ireadwf,&
2043 6619 : & istwfk1,istwfk2,kg2,kptns1,kptns2,mband2,mcg,mkmem1,mkmem2,&
2044 6619 : & mpi_enreg1,mpi_enreg2,mpw1,mpw2,my_nkpt2,nband1,nband2,&
2045 6619 : & ngfft1,ngfft2,nkpt1,nkpt2,npwarr1,npwarr2,nspinor1,nspinor2,&
2046 6619 : & nsppol1,nsppol2,nsym,occ,optorth,prtvol,randalg,restart,rprimd,&
2047 6619 : & sppoldbl,symrel,tnons,unkg2,wffinp,wffout)
2048 :
2049 : !Arguments ------------------------------------
2050 : !scalars
2051 : integer,intent(in) :: ceksp2,debug,exchn2n3d,fill,formeig,headform1,iout
2052 : integer,intent(in) :: ireadwf,mband2,mcg,mkmem1,mkmem2,mpw1,mpw2,my_nkpt2,nkpt1,nkpt2
2053 : integer,intent(in) :: nspinor1,nspinor2,nsppol1,nsppol2,nsym,optorth,prtvol,restart
2054 : integer,intent(in) :: randalg,sppoldbl,unkg2
2055 : real(dp),intent(in) :: ecut1,ecut2,ecut2_eff
2056 : type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2
2057 : type(wffile_type),intent(inout) :: wffinp,wffout
2058 : !arrays
2059 : integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2)
2060 : integer,intent(in) :: kg2(3,mpw2*mkmem2),nband1(nkpt1*nsppol1)
2061 : integer,intent(in) :: nband2(nkpt2*nsppol2),ngfft1(18),ngfft2(18)
2062 : integer,intent(in) :: npwarr1(nkpt1),npwarr2(nkpt2),symrel(3,3,nsym)
2063 : real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
2064 : real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
2065 : real(dp),intent(inout) :: cg(2,mcg) !vz_i pw_orthon vecnm
2066 : real(dp),intent(inout) :: eigen(mband2*(2*mband2)**formeig*nkpt2*nsppol2)!vz_i newocc
2067 : real(dp),intent(inout) :: occ(mband2*nkpt2*nsppol2) !vz_i
2068 :
2069 : !Local variables-------------------------------
2070 : !scalars
2071 : integer,parameter :: init_random=-5,nkpt_max=50,tobox=1,tosph=-1,wr=2
2072 : integer :: aux_stor,band_index,iband,icg,icg_aux,idum
2073 : integer :: ii,ikg2,ikpt1,ikpt10,ikpt2,ikptsp_prev,inplace,iproc
2074 : integer :: isppol1,isppol2,istwf10_k,localrdwf
2075 : integer :: mband1,mband_rd,mband_rw,mcg_aux,me1,me2,mgfft1,mgfft2
2076 : integer :: my_nspinor1,my_nspinor2
2077 : integer :: nb_band,nbd1,nbd1_rd,nbd2,nkpt_eff,nproc2,npw1,npw2,nsp
2078 : integer :: test_cycle,tim_rwwf
2079 : logical :: out_of_core2
2080 : character(len=500) :: msg
2081 : !arrays
2082 6619 : integer,allocatable :: kg1(:,:),kg2_k(:,:),kg_dum(:,:)
2083 : real(dp) :: kpoint(3),tsec(2)
2084 6619 : real(dp),allocatable :: cg_aux(:,:),eig_k(:),occ_k(:)
2085 :
2086 : ! *************************************************************************
2087 :
2088 6619 : call timab(780,1,tsec)
2089 6619 : call timab(781,1,tsec)
2090 :
2091 6619 : icg=0
2092 :
2093 : !Init MPI data
2094 6619 : me1=mpi_enreg1%me_kpt
2095 6619 : me2=mpi_enreg2%me_kpt
2096 6619 : nproc2 = mpi_enreg2%nproc_cell
2097 6619 : out_of_core2=(my_nkpt2/=0.and.mkmem2==0)
2098 :
2099 :
2100 6619 : if((nsppol1==2.and.nspinor2==2).or.(nspinor1==2.and. nsppol2==2))then
2101 : ! This is not yet possible. See later for a msg about where to make the needed modifs.
2102 : ! EDIT MT 20110707: these modifs are no more needed as they are now done in inwffil
2103 : write(msg, '(5a,i2,a,i2,2a,i2,a,i2,4a)' ) &
2104 0 : & 'The wavefunction translator is (still) unable to interchange',ch10,&
2105 0 : & 'spin-polarized wfs and spinor wfs. However,',ch10,&
2106 0 : & 'the input variables are nsppol1=',nsppol1,', and nspinor1=',nspinor1,ch10,&
2107 0 : & 'the output variables are nsppol2=',nsppol2,', and nspinor2=',nspinor2,ch10,&
2108 0 : & 'Action: use a non-spin-polarized wf to start a spinor wf,',ch10,&
2109 0 : & ' and a non-spinor wf to start a spin-polarized wf.'
2110 0 : ABI_ERROR(msg)
2111 : end if
2112 :
2113 6619 : my_nspinor1=max(1,nspinor1/mpi_enreg1%nproc_spinor)
2114 6619 : my_nspinor2=max(1,nspinor2/mpi_enreg2%nproc_spinor)
2115 131176 : mband1=maxval(nband1(1:nkpt1*nsppol1))
2116 :
2117 6619 : if(mkmem1==0 .and. out_of_core2)then
2118 0 : mband_rd=min(mband1,(mband2/nspinor2)*nspinor1)
2119 0 : if(mcg<mpw1*my_nspinor1*mband_rd)then
2120 0 : write(msg,'(2(a,i0))')' The dimension mcg= ',mcg,', should be larger than mband_rd= ',mband_rd
2121 0 : ABI_BUG(msg)
2122 : end if
2123 0 : if(mcg<mband2*mpw2*my_nspinor2)then
2124 : write(msg,'(a,i0,a,a,a,i0,a,i0,a,i2)' )&
2125 0 : & 'The dimension mcg= ',mcg,', should be larger than',ch10,&
2126 0 : & 'the product of mband2= ',mband2,', mpw2= ',mpw2,', and nspinor2= ',my_nspinor2
2127 0 : ABI_BUG(msg)
2128 : end if
2129 : end if
2130 :
2131 6619 : idum=init_random
2132 6619 : ikpt10 = 0
2133 6619 : istwf10_k=0
2134 6619 : band_index=0
2135 : icg=0
2136 :
2137 6619 : nkpt_eff=nkpt2
2138 6619 : if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max
2139 :
2140 : mgfft1=maxval(ngfft1(1:3))
2141 : mgfft2=maxval(ngfft2(1:3))
2142 19857 : ABI_MALLOC(kg1,(3,mpw1))
2143 19857 : ABI_MALLOC(kg2_k,(3,mpw2))
2144 6619 : ABI_MALLOC(kg_dum,(3,0))
2145 :
2146 6619 : if (debug>0) then
2147 0 : if (me1==0) then
2148 0 : write(std_out,'(a)' ) ' newkpt: kptns1'
2149 0 : call prmat (kptns1, 3, nkpt1, 3)
2150 : end if
2151 0 : if (me2==0) then
2152 0 : write(std_out,'(a)' ) ' newkpt: kptns2'
2153 0 : call prmat (kptns2, 3, nkpt2, 3)
2154 : end if
2155 : end if
2156 :
2157 6619 : ikptsp_prev=0
2158 :
2159 6619 : call timab(781,2,tsec)
2160 :
2161 : !Do outer loop over spins
2162 14117 : do isppol2=1,nsppol2
2163 :
2164 7498 : if (nsppol2==2 .and. me2==0) then
2165 1662 : write(std_out,'(a,i5)' ) ' newkpt: spin channel isppol2 = ',isppol2
2166 : end if
2167 :
2168 7498 : if (restart==1 .and. out_of_core2) rewind (unkg2)
2169 7498 : ikg2=0
2170 :
2171 : ! Do loop over new k point set
2172 144143 : do ikpt2=1,nkpt2
2173 :
2174 130026 : call timab(782,1,tsec)
2175 :
2176 130026 : nbd2=nband2(ikpt2+(isppol2-1)*nkpt2)
2177 130026 : npw2=npwarr2(ikpt2)
2178 :
2179 130026 : if(restart==1)then
2180 :
2181 : ! Announce the treatment of k point ikpt
2182 116234 : if(ikpt2<=nkpt_eff)then
2183 : ! This msg might be overwritten in parallel
2184 73159 : write(msg, '(a,i6,a,i8,a,i4)' )'P newkpt: treating ',nbd2,' bands with npw=',npw2,' for ikpt=',ikpt2
2185 : ! This msg might be overwritten in parallel
2186 73159 : if(mpi_enreg2%paralbd==1)then
2187 664 : do iproc=0,nproc2-1
2188 332 : nb_band=0
2189 1660 : do iband=1,nbd2
2190 1660 : if(mpi_enreg2%proc_distrb(ikpt2,iband,isppol2) == iproc)nb_band=nb_band+1
2191 : end do
2192 664 : if(nb_band/=0)then
2193 : write(msg, '(a,i6,a,i8,a,i4,a,i4)' ) &
2194 332 : & 'P newkpt: treating ',nb_band,' bands with npw=',npw2,' for ikpt=',ikpt2,' by node ',iproc
2195 : end if
2196 : end do
2197 : end if
2198 73159 : if(mpi_enreg2%paralbd==0) then
2199 : write(msg, '(a,i6,a,i8,a,i4,a,i4)' )&
2200 72827 : & 'P newkpt: treating ',nbd2,' bands with npw=',npw2,&
2201 145654 : & ' for ikpt=',ikpt2,' by node ',mpi_enreg2%proc_distrb(ikpt2,1,isppol2)
2202 : end if
2203 73159 : if(prtvol>0)then
2204 2726 : call wrtout(iout,msg,'COLL')
2205 : end if
2206 : end if
2207 :
2208 : ! Cut the writing if the limit is reached
2209 116234 : if(ikpt2==nkpt_eff+1)then
2210 375 : if(prtvol>0)then
2211 0 : call wrtout(iout,' newkpt: prtvol=0 or 1, do not print more k-points.','COLL')
2212 : end if
2213 : end if
2214 :
2215 : ! End of restart==1
2216 : end if
2217 :
2218 130026 : test_cycle=0
2219 130026 : if(proc_distrb_cycle(mpi_enreg2%proc_distrb,ikpt2,1,nbd2,isppol2,me2)) test_cycle=1
2220 : if(test_cycle==1)then
2221 19772 : if(formeig==0)then
2222 100316 : eigen(1+band_index : nbd2+band_index) = zero
2223 : ! occ(1+band_index : nbd2+band_index) = zero
2224 : band_index=band_index+nbd2
2225 : else
2226 21120 : eigen(1+band_index : 2*nbd2**2+band_index) = 0.0_dp
2227 : band_index=band_index+2*nbd2**2
2228 : end if
2229 : ! In the case this k point does not belong to me, cycle
2230 19772 : if (my_nkpt2==0) cycle
2231 19699 : if ((mkmem1==0) .and. (ireadwf==1) .and. (mpi_enreg2%paralbd==1))then
2232 0 : call WffReadSkipK(formeig,headform1,ikpt2,isppol2,mpi_enreg2,wffinp)
2233 0 : ikptsp_prev=ikptsp_prev+1
2234 : end if
2235 : cycle
2236 : end if
2237 :
2238 110254 : if(restart==1)then
2239 :
2240 97838 : if(mkmem2/=0)then
2241 58098626 : kg2_k(:,1:npw2)=kg2(:,1+ikg2:npw2+ikg2)
2242 : else if(mkmem2==0)then
2243 : ! Read the first line of a block and performs some checks on the unkg file.
2244 0 : ABI_ERROR("mkmem2 == 0 and rdnpw are not supported anymore.")
2245 0 : nsp=nspinor2
2246 : !call rdnpw(ikpt2,isppol2,nbd2,npw2,nsp,0,unkg2)
2247 : ! Read k+g data
2248 0 : read (unkg2) kg2_k(1:3,1:npw2)
2249 : end if
2250 :
2251 : end if
2252 :
2253 : ! Get ikpt1, the closest k from original set, from indkk
2254 110254 : ikpt1=indkk(ikpt2,1)
2255 110254 : if(sppoldbl==2 .and. isppol2==2)ikpt1=indkk(ikpt2+nkpt2,1)
2256 :
2257 110254 : npw1=npwarr1(ikpt1)
2258 : kpoint(:)=kptns1(:,ikpt1)
2259 :
2260 : ! Determine the spin polarization of the input data
2261 110254 : isppol1=isppol2
2262 110254 : if(nsppol2==2 .and. nsppol1==1)isppol1=1
2263 :
2264 110254 : if(restart==2)then
2265 12416 : if(ikpt2<=nkpt_eff)then
2266 6717 : write(msg,'(a,i4,i8,a,i4,i8)')'- newkpt: read input wf with ikpt,npw=',ikpt1,npw1,', make ikpt,npw=',ikpt2,npw2
2267 6717 : call wrtout(std_out,msg)
2268 6717 : if(iout/=6 .and. me2==0 .and. prtvol>0)then
2269 320 : call wrtout(iout,msg)
2270 : end if
2271 5699 : else if(ikpt2==nkpt_eff+1)then
2272 57 : call wrtout(std_out, '- newkpt: prtvol=0 or 1, do not print more k-points.')
2273 57 : if(iout/=6 .and. me2==0 .and. prtvol>0)then
2274 1 : call wrtout(iout,msg)
2275 : end if
2276 : end if
2277 : end if
2278 :
2279 : ! Set up the number of bands to be read
2280 110254 : nbd1=nband1(ikpt1+(isppol1-1)*nkpt1)
2281 110254 : nbd1_rd=min(nbd1,(nbd2/nspinor2)*nspinor1)
2282 :
2283 : ! Check that number of bands is not being increased if fill==0 --if so
2284 : ! print warning and reset new wf file nband2 to only allowed number
2285 110254 : if ( nbd2/nspinor2 > nbd1/nspinor1 .and. fill==0) then
2286 0 : if(ikpt2<=nkpt_eff)then
2287 0 : write(msg, '(a,i8,a,i8,a,i8)' )' newkpt: nband2=',nbd2,' < nband1=',nbd1,' => reset nband2 to ',nbd1
2288 0 : call wrtout(std_out,msg)
2289 : end if
2290 0 : nbd2=nbd1
2291 : end if
2292 :
2293 : ! Prepare the reading of the wavefunctions: the correct record is selected
2294 : ! WARNING : works only for GS - for RF the number of record differs
2295 110254 : if(restart==2 .and. mkmem1==0)then
2296 0 : ABI_ERROR("mkmem1 == 0 has been removed.")
2297 :
2298 0 : if(debug>0)then
2299 0 : write(msg, '(a,a,a,a,i5,a,i5,a,a,i5,a,i5)' ) ch10,&
2300 0 : ' newkpt: about to call randac',ch10,&
2301 0 : ' for ikpt1=',ikpt1,', ikpt2=',ikpt2,ch10,&
2302 0 : ' and isppol1=',isppol1,', isppol2=',isppol2
2303 0 : call wrtout(std_out,msg)
2304 : end if
2305 :
2306 : !call randac(debug,headform1,ikptsp_prev,ikpt1,isppol1,nband1,nkpt1,nsppol1,wffinp)
2307 : end if
2308 :
2309 : ! Read the data for nbd2 bands at this k point
2310 : ! Must decide whether an auxiliary storage is needed
2311 : ! When mkmem1==0 and mkmem2==0 , the cg array should be large enough ...
2312 : ! When mkmem1==0 and mkmem2/=0 , each k-point block in cg might not be large enough
2313 : ! however, will read at most (nbd2/nspinor2)*nspinor1 bands from disk
2314 : ! When mkmem1/=0 , it is supposed that each input k-point block is smaller
2315 : ! than the corresponding output k-point block, so that the input data
2316 : ! have been placed already in cg, at the k-point location where they are needed
2317 110254 : aux_stor=0
2318 110254 : if(mkmem2/=0 .and. mkmem1==0)then
2319 38371 : mcg_aux=npw1*my_nspinor1*nbd1
2320 38371 : if(nbd1_rd<nbd1)mcg_aux=npw1*my_nspinor1*nbd1_rd
2321 38371 : if( mcg_aux > npw2*my_nspinor2*nbd2 )then
2322 0 : aux_stor=1 ; icg_aux=0
2323 0 : ABI_MALLOC(cg_aux,(2,mcg_aux))
2324 : end if
2325 : end if
2326 :
2327 110254 : mband_rw=max(nbd1_rd,nbd2)
2328 330762 : ABI_MALLOC(eig_k,(mband_rw*(2*mband_rw)**formeig))
2329 110254 : if(formeig==0) then
2330 168135 : ABI_MALLOC(occ_k,(mband_rw))
2331 : else
2332 54209 : ABI_MALLOC(occ_k,(0))
2333 : end if
2334 :
2335 110254 : if(mkmem1/=0 .and. ireadwf==1)then
2336 : ! Checks that nbd1 and nbd1_rd are equal if eig and occ are input
2337 71883 : if(nbd1/=nbd1_rd)then
2338 : write(msg,'(a,a,a,i6,a,i6)')&
2339 0 : & 'When mkmem1/=0, one must have nbd1=nbd1_rd, while',ch10,&
2340 0 : & 'nbd1 = ',nbd1,', and nbd1_rd = ',nbd1_rd
2341 0 : ABI_BUG(msg)
2342 : end if
2343 : ! Need to put eigenvalues in eig_k, same for occ
2344 : ! Note use of band_index, since it is assumed that eigen and occ
2345 : ! already have spin-k point location structure than output.
2346 71883 : if(formeig==0)then
2347 149520 : eig_k(1:nbd1_rd)=eigen(1+band_index : nbd1_rd+band_index)
2348 : ! occ_k(1:nbd1_rd)=occ(1+band_index : nbd1_rd+band_index)
2349 54209 : else if(formeig==1)then
2350 : ! The matrix of eigenvalues has size nbd1 , that must be equal
2351 : ! to nbd1_rd in the case mkmem1/=0)
2352 4910737 : eig_k(1:2*nbd1_rd**2)=eigen(1+band_index : 2*nbd1_rd**2+band_index)
2353 : end if
2354 : end if
2355 :
2356 110254 : call timab(782,2,tsec)
2357 :
2358 : ! Must read the wavefunctions if they are not yet in place
2359 110254 : if(mkmem1==0 .and. ireadwf==1)then
2360 :
2361 0 : if (debug>0 .and. restart==2) then
2362 : write(msg,'(a,i5,a,a,i5,a,i5,a)' ) &
2363 0 : & ' newkpt: about to call rwwf with ikpt1=',ikpt1,ch10,&
2364 0 : & ' and nband(ikpt1)=',nband1(ikpt1),' nbd2=',nbd2,'.'
2365 0 : call wrtout(std_out,msg)
2366 : end if
2367 :
2368 0 : if(mpi_enreg1%paralbd==0)tim_rwwf=21
2369 0 : if(mpi_enreg1%paralbd==1)tim_rwwf=22
2370 :
2371 0 : if(aux_stor==0)then
2372 : call rwwf(cg,eig_k,formeig,headform1,icg,ikpt1,isppol1,kg_dum,mband_rw,mcg,mpi_enreg1,&
2373 0 : & nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp)
2374 : else
2375 0 : icg_aux=0
2376 : call rwwf(cg_aux,eig_k,formeig,headform1,icg_aux,ikpt1,isppol1,kg_dum,mband_rw,mcg_aux,&
2377 0 : & mpi_enreg1,nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp)
2378 : end if
2379 : end if
2380 :
2381 110254 : call timab(783,1,tsec)
2382 :
2383 110254 : if(formeig==1 .and. nbd2/=nbd1_rd .and. ireadwf==1)then
2384 : ! Change the storage of eig_k
2385 0 : if(nbd1_rd<nbd2)then
2386 0 : do iband=nbd1_rd,1,-1
2387 : ! The factor of two is for complex eigenvalues
2388 0 : do ii=2*nbd2,2*nbd1_rd+1,-1
2389 0 : eig_k(ii+(iband-1)*2*nbd2)=huge(zero)/10.0_dp
2390 : end do
2391 0 : do ii=2*nbd1_rd,1,-1
2392 0 : eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd)
2393 : end do
2394 : end do
2395 : else if(nbd1_rd>nbd2)then
2396 0 : do iband=1,nbd2
2397 : ! The factor of two is for complex eigenvalues
2398 0 : do ii=1,2*nbd2
2399 0 : eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd)
2400 : end do
2401 : end do
2402 : end if
2403 : end if
2404 :
2405 : ! If change nsppol, must adapt the occupation numbers
2406 : ! if(nsppol1/=nsppol2)then
2407 : ! occ_k(1:nbd2)=occ_k(1:nbd2)*nsppol1/dbl(nsppol2)
2408 : ! then
2409 :
2410 : ! In case nsppol1=2 and nspinor2=2, one should read
2411 : ! the other spin-component, and form a spinor wf here, before calling
2412 : ! wfconv. One should treat eig_k and occ_k as well.
2413 : ! A similar operation is to be performed when nspino1=2 and nsppol2=2
2414 : ! EDIT - MT 20110707: the building of the spinor wf is now done in wfffil
2415 : ! no need to make it here....
2416 :
2417 : ! DEBUG
2418 : ! write(std_out,*)' newkpt: before wfconv'
2419 : ! write(std_out,*)' newkpt: mkmem2=',mkmem2
2420 : ! stop
2421 : ! ENDDEBUG
2422 :
2423 110254 : call timab(783,2,tsec)
2424 110254 : call timab(784,1,tsec)
2425 :
2426 : ! Note the use of mband2, while mband is used inside
2427 : ! write(std_out,*) 'in newkpt,before wfconv,npw1,npw2',npw1,npw2
2428 110254 : inplace=1
2429 110254 : if(aux_stor==0)then
2430 : call wfconv(ceksp2,cg,cg,debug,ecut1,ecut2,ecut2_eff,&
2431 : & eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg,icg,&
2432 : & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
2433 : & kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,&
2434 : & mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,&
2435 : & ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,&
2436 110254 : & occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons)
2437 : else
2438 : call wfconv(ceksp2,cg_aux,cg_aux,debug,ecut1,ecut2,ecut2_eff,&
2439 : & eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg_aux,icg_aux,&
2440 : & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
2441 : & kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,&
2442 : & mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,&
2443 : & ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,&
2444 0 : & occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons)
2445 : end if
2446 :
2447 110254 : call timab(784,2,tsec)
2448 :
2449 : ! Finally write new wf to disk file or save in permanent file
2450 110254 : if(mkmem2==0)then
2451 :
2452 : ! Note that in this case, we are sure aux_stor==0
2453 0 : if(mpi_enreg2%paralbd==0)tim_rwwf=21
2454 0 : if(mpi_enreg2%paralbd==1)tim_rwwf=22
2455 : call rwwf(cg,eig_k,formeig,0,0,ikpt2,isppol2,kg2_k,nbd2,mcg,mpi_enreg2,&
2456 0 : & nbd2,nbd2,npw2,my_nspinor2,occ_k,wr,1,tim_rwwf,wffout)
2457 :
2458 : end if
2459 :
2460 110254 : call timab(785,1,tsec)
2461 :
2462 110254 : if(mkmem2/=0)then
2463 110254 : if(aux_stor==1)then
2464 0 : cg(:,1+icg:npw2*nbd2*my_nspinor2+icg)=cg_aux(:,1:npw2*nbd2*my_nspinor2)
2465 0 : ABI_FREE(cg_aux)
2466 : end if
2467 :
2468 110254 : icg=icg+npw2*nbd2*my_nspinor2
2469 110254 : ikg2=ikg2+npw2
2470 : end if
2471 :
2472 5351209 : eigen(1+band_index:nbd2*(2*nbd2)**formeig+band_index) = eig_k(1:nbd2*(2*nbd2)**formeig)
2473 : ! occ(1+band_index:nbd2+band_index)=occ_k(1:nbd2)
2474 :
2475 110254 : if(formeig==0)then
2476 56045 : band_index=band_index+nbd2
2477 54209 : else if(formeig==1)then
2478 54209 : band_index=band_index+2*nbd2**2
2479 : end if
2480 :
2481 110254 : ABI_FREE(eig_k)
2482 110254 : ABI_FREE(occ_k)
2483 :
2484 117825 : call timab(785,2,tsec)
2485 :
2486 : end do ! ikpt2
2487 : end do ! isppol2
2488 :
2489 6619 : call timab(786,1,tsec)
2490 :
2491 : if(xmpi_paral==1)then
2492 : ! Transmit eigenvalues (not yet occupation numbers)
2493 : ! newkpt.F90 is not yet suited for RF format
2494 : ! This routine works in both localrdwf=0 or 1 cases.
2495 : ! However, in the present routine, localrdwf is to be considered
2496 : ! as 1 always, since the transfer has been made in wfsinp .
2497 6619 : localrdwf=1
2498 6619 : call pareigocc(eigen,formeig,localrdwf,mpi_enreg2,mband2,nband2,nkpt2,nsppol2,occ,1)
2499 : end if
2500 :
2501 6619 : ABI_FREE(kg1)
2502 6619 : ABI_FREE(kg2_k)
2503 6619 : ABI_FREE(kg_dum)
2504 :
2505 6619 : call timab(786,2,tsec)
2506 6619 : call timab(780,2,tsec)
2507 :
2508 6619 : end subroutine newkpt
2509 : !!***
2510 :
2511 : !!****f* m_inwffil/wfconv
2512 : !! NAME
2513 : !! wfconv
2514 : !!
2515 : !! FUNCTION
2516 : !! This subroutine treats the wavefunctions for one k point,
2517 : !! and converts them to other parameters.
2518 : !!
2519 : !! INPUTS
2520 : !! ceksp2=if 1, center the output sphere of pw on Gamma; if 0, on each k-point (usual).
2521 : !! cg1(2,mcg1)=wavefunction array
2522 : !! debug= if 1, print some messages ; otherwise, 0.
2523 : !! ecut1=kinetic energy cutoffs for basis sphere 1 (hartree)
2524 : !! ecut2=kinetic energy cutoff beyond which the coefficients of wf2 vanish (Ha)
2525 : !! ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree)
2526 : !! eig_k1(mband1*(2*mband1)**formeig)=eigenvalues
2527 : !! exchn2n3d=if 1, n2 and n3 are exchanged
2528 : !! formeig option (format of the eigenvalues and eigenvector) :
2529 : !! 0 => ground-state format (initialisation of
2530 : !! eigenvectors with random numbers, vector of eigenvalues)
2531 : !! 1 => respfn format (initialisation of
2532 : !! eigenvectors with 0 s, hermitian matrix of eigenvalues)
2533 : !! gmet1(3,3)=reciprocal space metric (bohr^-2) for input wf
2534 : !! gmet2(3,3)=reciprocal space metric (bohr^-2) for output wf
2535 : !! icg1=shift to be given to the location of the data in the array cg1
2536 : !! icg2=shift to be given to the location of the data in the array cg2
2537 : !! ikpt1=number of the k point actually treated (input wf numbering)
2538 : !! ikpt10=number of the k point previously treated (input wf numbering)
2539 : !! ikpt2=number of the k point actually treated (output numbering)
2540 : !! indkk(nkpt2*sppoldbl,6)=describe k point number of kptns1 that allows to
2541 : !! generate wavefunctions closest to given kpt2 (and possibly isppol2=2)
2542 : !! indkk(:,1)=k point number of kpt1
2543 : !! indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
2544 : !! (if 0, means no symmetry operation, equivalent to identity )
2545 : !! indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
2546 : !! to give kpt1b, that is the closest to kpt2.
2547 : !! indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise
2548 : !! inplace= if 0, cg1 and cg2 are different in the calling routine,
2549 : !! if 1, cg1 and cg2 are identical (they have the same memory location)
2550 : !! This is also true for the pairs (eig_k1,eig_k2) and (occ_k1,occ_k2)
2551 : !! isppol2=spin variable for output wavefunctions
2552 : !! istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1
2553 : !! istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2
2554 : !! kg1(3,mpw1)=dimensionless coords of G vecs in basis sphere at k point (input wf)
2555 : !! kg2(3,mpw2)=dimensionless coords of G vecs in basis sphere at k point (output wf)
2556 : !! kptns1(3,nkpt1)=k point set for input wavefunctions
2557 : !! kptns2(3,nkpt2)=k point set for output wavefunctions
2558 : !! mband1=dimension of eig_k1 and occ_k1 arrays
2559 : !! mband2=dimension of eig_k2 and occ_k2 arrays
2560 : !! mcg1=dimension of cg1 array (at least npw1*nspinor1*nbd1)
2561 : !! mcg2=dimension of cg2 array (at least npw2*nspinor2*nbd2)
2562 : !! mpi_enreg1=information about MPI parallelization for set 1
2563 : !! mpi_enreg2=information about MPI parallelization for set 2
2564 : !! mpw1=dimension of kg1, can be set to 0 if not needed
2565 : !! mpw2=dimension of kg2, can be set to 0 if not needed
2566 : !! nbd1=number of bands contained in cg1,eig_k1,occ_k1 at this k-point - spin (at input)
2567 : !! nbd2=number of bands contained in cg2,eig_k2,occ_k2 at this k-point - spin (at output)
2568 : !! ngfft1(18)=all needed information about 3D FFT, for input wavefunctions
2569 : !! ngfft2(18)=all needed information about 3D FFT, for output wavefunctions
2570 : !! see ~abinit/doc/variables/vargs.htm#ngfft
2571 : !! nkpt1=number of k points for input wavefunctions
2572 : !! nkpt2=number of k points for output wavefunctions
2573 : !! npw1=number of planewaves for input wavefunctions
2574 : !! npw2=number of planewaves for output wavefunctions
2575 : !! nspinor1=number of spinors for input wavefunctions
2576 : !! nspinor2=number of spinors for output wavefunctions
2577 : !! nsym=number of symmetry elements in space group
2578 : !! occ_k1(mband1)=occupation numbers
2579 : !! optorth=1 if the WFs are orthogonalized before leaving the routine
2580 : !! randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility
2581 : !! restart=if 2, conversion between wavefunctions
2582 : !! if 1, direct restart is allowed (see hdr_check.f)
2583 : !! rprimd2(3,3)=dimensional primitive translations for real space (bohr)
2584 : !! needed only for the spinor rotation
2585 : !! sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn
2586 : !! if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries
2587 : !! symrel(3,3,nsym)=symmetry operations in real space in terms
2588 : !! of primitive translations
2589 : !! tnons(3,nsym)=nonsymmorphic translations for symmetry operations
2590 : !! dtset <type(dataset_type)>=all input variables for this dataset
2591 : !!
2592 : !! OUTPUT
2593 : !! cg2(2,mcg2)=wavefunction array
2594 : !! eig_k2(mband2*(2*mband2)**formeig)=eigenvalues
2595 : !! occ_k2(mband2)=occupation (completed with zeros)
2596 : !!
2597 : !! SIDE EFFECTS
2598 : !! Input/Output:
2599 : !! ikpt10=at input, number of the k point previously treated (input wf numbering)
2600 : !! (if this is the first call for the present k point set, ikpt10 should be 0)
2601 : !! at output, number of the k point just treated (input wf numbering)
2602 : !! kg1, kg2, npw1 and npw2 should not be modified by kpgsph (TD).
2603 : !!
2604 : !! NOTES
2605 : !! Note that this routine can make an in-place conversion
2606 : !! (see the input variable "inplace"),
2607 : !! if cg1 and cg2 are equal, as well as the pairs (icg1,icg2),
2608 : !! (eig_k1,eig_k2),(occ_k1,occ_k2) and (mband1,mband2)
2609 : !!
2610 : !! It can also be used to fill or to initialize wavefunctions
2611 : !! at one k point
2612 : !! (filling with random numbers or 0''s, according to the value
2613 : !! of formeig), if the input number of bands (nbd1) is 0.
2614 : !! In the latter case, one should use the same values of input
2615 : !! wavefunction parameters
2616 : !! than for output wavefunction parameters, except nbd1.
2617 : !!
2618 : !! The input parameters are indexed with 1, the output parameters
2619 : !! are indexed with 2.
2620 : !!
2621 : !! Some of the arguments are arrays dimensioned with nkpt1 or nkpt2.
2622 : !! Note that for these, only the elements for ikpt1 or ikpt2 will be used.
2623 : !!
2624 : !! The number of input bands must already be minimal at the input.
2625 : !! This means, when input and output nspinor are equal : nbd1<nbd2
2626 : !! When the two nspinor differ, one must have nbd1/nspinor1<nbd2/nspinor2
2627 : !!
2628 : !! SOURCE
2629 :
2630 114118 : subroutine wfconv(ceksp2,cg1,cg2,debug,ecut1,ecut2,ecut2_eff,&
2631 114118 : & eig_k1,eig_k2,exchn2n3d,formeig,gmet1,gmet2,icg1,icg2,&
2632 114118 : & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,&
2633 114118 : & kg1,kg2,kptns1,kptns2,mband1,mband2,mcg1,mcg2,mpi_enreg1,mpi_enreg2,&
2634 : & mpw1,mpw2,nbd1,nbd2,ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,&
2635 114118 : & nsym,occ_k1,occ_k2,optorth,randalg,restart,rprimd2,sppoldbl,symrel,tnons)
2636 :
2637 : !Arguments ------------------------------------
2638 : !scalars
2639 : integer,intent(in) :: ceksp2,debug,exchn2n3d,formeig,icg1,icg2,ikpt1
2640 : integer,intent(in) :: ikpt2,inplace,isppol2,mband1,mband2,mcg1,mcg2,mpw1,mpw2
2641 : integer,intent(in) :: nbd1,nbd2,nkpt1,nkpt2,nspinor1,nspinor2,nsym
2642 : integer,intent(in) :: optorth,randalg,restart,sppoldbl
2643 : integer,intent(inout) :: ikpt10,npw1,npw2
2644 : real(dp),intent(in) :: ecut1,ecut2,ecut2_eff
2645 : type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2
2646 : !arrays
2647 : integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2)
2648 : integer,intent(in) :: ngfft1(18),ngfft2(18),symrel(3,3,nsym)
2649 : integer,intent(inout) :: kg1(3,mpw1),kg2(3,mpw2)
2650 : real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
2651 : real(dp),intent(in) :: rprimd2(3,3),tnons(3,nsym)
2652 : real(dp),intent(inout) :: cg1(2,mcg1),cg2(2,mcg2)
2653 : real(dp),intent(inout) :: eig_k1(mband1*(2*mband1)**formeig)
2654 : real(dp),intent(inout) :: eig_k2(mband2*(2*mband2)**formeig),occ_k1(mband1)
2655 : real(dp),intent(inout) :: occ_k2(mband2)
2656 :
2657 : !Local variables ------------------------------
2658 : !scalars
2659 : integer,parameter :: nkpt_max=50,tobox=1,tosph=-1
2660 : integer :: conv_tnons,convert,fftalg,fold1,fold2,foldim,foldre,i1,i2,iband
2661 : integer :: iband_first,iband_last,icgmod,ierr,index,ipw
2662 : integer :: ispinor,ispinor1,ispinor2,ispinor_first,ispinor_last
2663 : integer :: istwf10_k,istwf1_k,istwf2_k,isym,itimrev
2664 : integer :: mgfft1,mgfft2,n1,n2,n3,n4,n5,n6
2665 : integer :: nbremn,npwtot,nspinor_index,nspinor1_this_proc,nspinor2_this_proc
2666 : integer :: order,ortalgo
2667 : real(dp) :: ai,ar,arg,bi,br,eig_tmp,spinrots,spinrotx,spinroty,spinrotz
2668 : character(len=500) :: msg
2669 : integer, parameter :: int64 = selected_int_kind(18)
2670 : integer(KIND=int64) :: seed
2671 : !arrays
2672 : integer :: atindx(1),identity(3,3),ngfft_now(18),no_shift(3),shiftg(3)
2673 : integer :: symm(3,3),symrel_conv(3,3)
2674 114118 : integer,allocatable :: gbound1(:,:),gbound2(:,:)
2675 : real(dp) :: kpoint1(3),kpoint2_sph(3),phktnons(2,1),spinrot(4),tnons_conv(3),tsec(2)
2676 114118 : real(dp),allocatable :: cfft(:,:,:,:),dum(:,:),phase1d(:,:),phase3d(:,:)
2677 114118 : real(dp),allocatable :: wavef1(:,:),wavef2(:,:),wavefspinor(:,:)
2678 :
2679 : ! *************************************************************************
2680 :
2681 456472 : mgfft1=maxval(ngfft1(1:3))
2682 456472 : mgfft2=maxval(ngfft2(1:3))
2683 : if(.false.)write(std_out,*)occ_k1 ! just to keep occ_k1 as an argument before resolving the issue of its transfer
2684 :
2685 114118 : if(nspinor1/=1 .and. nspinor1/=2)then
2686 0 : write(msg,'(a,i0)')'The argument nspinor1 must be 1 or 2, while it is nspinor1 = ',nspinor1
2687 0 : ABI_BUG(msg)
2688 : end if
2689 :
2690 114118 : if(nspinor2/=1 .and. nspinor2/=2)then
2691 0 : write(msg,'(a,i0)')' The argument nspinor2 must be 1 or 2, while it is nspinor2=',nspinor2
2692 0 : ABI_BUG(msg)
2693 : end if
2694 :
2695 114118 : if(nspinor1==2 .and. mod(nbd1,2)/=0)then
2696 0 : write(msg,'(a,i0)')' When nspinor1 is 2, nbd1 must be even, while it is nbd1 = ',nbd1
2697 0 : ABI_BUG(msg)
2698 : end if
2699 :
2700 114118 : if(nspinor2==2 .and. mod(nbd2,2)/=0)then
2701 0 : write(msg,'(a,i0)')' When nspinor2 is 2, nbd2 must be even, while it is nbd2=',nbd2
2702 0 : ABI_BUG(msg)
2703 : end if
2704 :
2705 114118 : if(nbd1/nspinor1>nbd2/nspinor2)then
2706 : write(msg, '(3a,2i6,3a,2i6,a)' )&
2707 0 : & 'In wfconv, the nbd/nspinor ratio cannot decrease. However,',ch10,&
2708 0 : & 'the initial quantities are nbd1,nspinor1=',nbd1,nspinor1,', and',ch10,&
2709 0 : & 'the requested final quantities are nbd2,nspinor2=',nbd2,nspinor2,'.'
2710 0 : ABI_BUG(msg)
2711 : end if
2712 :
2713 456472 : ngfft_now(1:3)=ngfft1(1:3)
2714 1369416 : ngfft_now(8:18)=ngfft1(8:18)
2715 : !This line is the reason why ngfft_now has to be introduced
2716 114118 : ngfft_now(7)=101
2717 456472 : ngfft_now(4:6)=ngfft_now(1:3)
2718 114118 : n1=ngfft_now(1) ; n2=ngfft_now(2) ; n3=ngfft_now(3)
2719 114118 : n4=ngfft_now(4) ; n5=ngfft_now(5) ; n6=ngfft_now(6)
2720 114118 : fftalg=ngfft_now(7)
2721 :
2722 : !Parallelization over spinors management
2723 114118 : nspinor1_this_proc=max(1,nspinor1/mpi_enreg1%nproc_spinor)
2724 114118 : nspinor2_this_proc=max(1,nspinor2/mpi_enreg2%nproc_spinor)
2725 :
2726 : !In order to generate IN PLACE new wfs from old wfs, the loop
2727 : !over bands and spinors must be done in one direction or the other,
2728 : !depending on npw1 and npw2, nspinor1 and nspinor2.
2729 : !If nspinor1=1 and nspinor2=2 , note that one will generate
2730 : !from nbd1 states of npw1 coefficients,
2731 : !2*nbd1 states of 2*npw2 coefficients. nbd1 cancels in comparing
2732 : !these expressions, but nspinor2 appears squared.
2733 : !The same line of thought works for the case nspinor1=2 and nspinor2=1
2734 114118 : order=1
2735 114118 : iband_first=1 ; iband_last=nbd1
2736 114118 : ispinor_first=1 ; ispinor_last=nspinor1
2737 114118 : if(nspinor1==2 .and. nspinor2==1)then
2738 0 : order=2 ; iband_last=nbd1-1 ; ispinor_last=1
2739 : end if
2740 : !Here, reverse the order if needed
2741 114118 : if( npw2*nspinor2**2 > npw1*nspinor1**2 )then
2742 2477 : order=-order
2743 2477 : iband_first=iband_last ; iband_last=1
2744 2477 : ispinor_first=ispinor_last ; ispinor_last=1
2745 : end if
2746 :
2747 456472 : kpoint1(:)=kptns1(:,ikpt1)
2748 114118 : istwf1_k=istwfk1(ikpt1)
2749 :
2750 114118 : kpoint2_sph(:)=0.0_dp
2751 456472 : if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2)
2752 114118 : istwf2_k=istwfk2(ikpt2)
2753 :
2754 : !DEBUG
2755 : !write(std_out,*)'ecut1,ecut2_eff=',ecut1,ecut2_eff
2756 : !write(std_out,*)'gmet1,gmet2=',gmet1,gmet2
2757 : !write(std_out,*)'kpoint1,kpoint2_sph=',kpoint1,kpoint2_sph
2758 : !write(std_out,*)'nspinor1,nspinor2',nspinor1,nspinor2
2759 : !write(std_out,*)'istwf1_k,istwf2_k=',istwf1_k,istwf2_k
2760 : !write(std_out,*)'nbd1,tol8=',nbd1,tol8
2761 : !ENDDEBUG
2762 :
2763 : !Determine whether it will be needed to convert the existing
2764 : !wavefunctions, or simply to complete them.
2765 :
2766 114118 : convert=0
2767 114118 : if(nbd1/=0)then
2768 75747 : if(abs(ecut2_eff-ecut1)>tol8)convert=convert+1
2769 984711 : if(sum(abs(gmet2(:,:)-gmet1(:,:)))>tol8)convert=convert+2
2770 302988 : if(sum(abs(kpoint2_sph(:)-kpoint1(:)))>tol8)convert=convert+4
2771 75747 : if(nspinor2/=nspinor1)convert=convert+8
2772 75747 : if(istwf2_k/=istwf1_k)convert=convert+16
2773 : end if
2774 :
2775 : !This is a supplementary check
2776 114118 : if(restart==1 .and. convert/=0)then
2777 0 : ABI_BUG('Restart==1 and convert/=0 are exclusive')
2778 : end if
2779 :
2780 : !Determine whether symmetries must be used
2781 114118 : conv_tnons=0
2782 114118 : no_shift(:)=0
2783 114118 : identity(:,:)=0
2784 114118 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2785 114118 : isym=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,2)
2786 : !write(std_out,*)' wfconv : isym=',isym
2787 114118 : itimrev=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,6)
2788 114118 : if(isym/=0)then
2789 161408 : symrel_conv(:,:)=symrel(:,:,isym)
2790 12416 : call mati3inv(symrel_conv,symm)
2791 49664 : shiftg(:)=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,3:5)
2792 49664 : tnons_conv(:)=tnons(:,isym)
2793 49664 : if(sum(tnons_conv(:)**2)>tol8)then
2794 : ! Need to compute phase factors associated with nonsymmorphic translations.
2795 443 : conv_tnons=1
2796 1329 : ABI_MALLOC(phase3d,(2,npw1))
2797 1329 : ABI_MALLOC(phase1d,(2,(2*n1+1)+(2*n2+1)+(2*n3+1)))
2798 : ! Although the routine getph is originally written for
2799 : ! atomic phase factors, it does precisely what we want
2800 443 : atindx(1)=1
2801 443 : call getph(atindx,1,n1,n2,n3,phase1d,tnons_conv)
2802 : end if
2803 12416 : if(nspinor1==2 .and. nspinor2==2)then
2804 : ! Compute rotation in spinor space
2805 115 : call getspinrot(rprimd2,spinrot,symrel_conv)
2806 : end if
2807 : else
2808 101702 : shiftg(:)=0
2809 101702 : symm(:,:)=identity(:,:)
2810 101702 : spinrot(:)=zero
2811 101702 : spinrot(1)=one
2812 : end if
2813 114118 : if(itimrev/=0)then
2814 34619 : symm(:,:)=-symm(:,:)
2815 : end if
2816 :
2817 : !DEBUG
2818 : !write(std_out,'(a,i3,2x,3i3,2x,9i3)')' wfconv : isym,shiftg,symm=',isym,shiftg,symm
2819 : !write(std_out,*)' wfconv : ecut2_eff,ecut1=',ecut2_eff,ecut1
2820 : !write(std_out,*)' wfconv : istwf1_k,istwf2_k=',istwf1_k,istwf2_k
2821 : !write(std_out,*)' wfconv : kpoint1(:),kpoint2_sph(:)=',&
2822 : !& kpoint1(:),kpoint2_sph(:)
2823 : !write(std_out,*)' wfconv : nspinor1,nspinor2=',nspinor1,nspinor2
2824 : !ENDDEBUG
2825 :
2826 : !if (mpi_enreg1%fft_option_lob==0) mpi_enreg1%fft_option_lob=1
2827 : !if (mpi_enreg2%fft_option_lob==0) mpi_enreg2%fft_option_lob=1
2828 :
2829 114118 : if (restart==2.and.(convert/=0.or.(nbd2/nspinor2>nbd1/nspinor1.and.formeig==0))) then
2830 : ! kg2 is needed both for FFT grid conversion and for envlop
2831 : ! Choose the center of the sphere : either gamma, or each k-point
2832 10520 : kpoint2_sph(:)=0.0_dp
2833 42080 : if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2)
2834 : istwf2_k=istwfk2(ikpt2)
2835 10520 : call kpgsph(ecut2_eff,exchn2n3d,gmet2,0,ikpt2,istwf2_k,kg2,kpoint2_sph,1,mpi_enreg2,mpw2,npw2)
2836 : end if
2837 :
2838 114118 : if(convert/=0)then
2839 10514 : istwf10_k=0
2840 10514 : if(ikpt10/=0)istwf10_k=istwfk1(ikpt10)
2841 :
2842 : ! Only need G sphere if different from last time
2843 10514 : if ( ikpt1/=ikpt10 .or. istwf1_k/=istwf10_k ) then
2844 :
2845 6625 : call kpgsph (ecut1,exchn2n3d,gmet1,0,ikpt1,istwf1_k,kg1,kpoint1,1,mpi_enreg1,mpw1,npw1)
2846 6625 : if (debug>0) then
2847 : write(msg, '(a,f8.3,a,a,3f8.5,a,a,i3,a,3(a,3es16.8,a),a,3i4,a,i5,a)' )&
2848 0 : & ' wfconv: called kpgsph with ecut1=',ecut1,ch10,&
2849 0 : & ' kpt1=',kptns1(1:3,ikpt1),ch10,&
2850 0 : & ' istwf1_k=',istwf1_k,ch10,&
2851 0 : & ' gmet1= ',gmet1(1:3,1),ch10,&
2852 0 : & ' ',gmet1(1:3,2),ch10,&
2853 0 : & ' ',gmet1(1:3,3),ch10,&
2854 0 : & ' ngfft=',ngfft_now(1:3),' giving npw1=',npw1,'.'
2855 0 : call wrtout(std_out,msg)
2856 : end if
2857 6625 : ikpt10 = ikpt1
2858 6625 : istwf10_k=istwf1_k
2859 : end if
2860 :
2861 10514 : if(conv_tnons==1)then
2862 443 : arg=two_pi*(kpoint1(1)*tnons_conv(1)+ kpoint1(2)*tnons_conv(2)+ kpoint1(3)*tnons_conv(3) )
2863 443 : phktnons(1,1)=cos(arg)
2864 443 : phktnons(2,1)=sin(arg)
2865 : ! Convert 1D phase factors to 3D phase factors exp(i 2 pi (k+G).tnons )
2866 443 : call ph1d3d(1,1,kg1,1,1,npw1,n1,n2,n3,phktnons,phase1d,phase3d)
2867 : end if
2868 :
2869 52570 : ABI_MALLOC(cfft,(2,n4,n5,n6))
2870 31542 : ABI_MALLOC(wavef1,(2,npw1))
2871 31542 : ABI_MALLOC(wavef2,(2,npw2))
2872 10514 : if(nspinor1==2 .and. nspinor2==2) then
2873 264 : ABI_MALLOC(wavefspinor,(2,2*npw2))
2874 : end if
2875 42056 : ABI_MALLOC(gbound1,(2*mgfft1+8,2))
2876 42056 : ABI_MALLOC(gbound2,(2*mgfft2+8,2))
2877 10514 : call sphereboundary(gbound1,istwf1_k,kg1,mgfft1,npw1)
2878 10514 : call sphereboundary(gbound2,istwf2_k,kg2,mgfft2,npw2)
2879 :
2880 : ! Take old wf from sphere->box, the new from box->sphere
2881 : ! One pays attention not to have a problem of erasing data when replacing
2882 : ! a small set of coefficient by a large set, or the reverse.
2883 : ! This is the reason of the use of order, _first and _last variables,
2884 : ! defined earlier.
2885 10514 : nspinor_index=mpi_enreg1%me_spinor+1
2886 84623 : do iband=iband_first,iband_last,order
2887 149960 : do ispinor1=ispinor_first,ispinor_last,order
2888 75851 : ispinor=ispinor1
2889 75851 : if (mpi_enreg1%paral_spinor==1) then
2890 0 : if (ispinor1==nspinor_index) then
2891 : ispinor=1
2892 : else
2893 0 : if (nspinor1==2.and.nspinor2==2) wavefspinor(:,(ispinor1-1)*npw2+1:ispinor1*npw2)=zero
2894 : cycle
2895 : end if
2896 : end if
2897 :
2898 : ! Copy input wf
2899 75851 : i1=(ispinor-1)*npw1+(iband-1)*nspinor1_this_proc*npw1+icg1
2900 67728491 : wavef1(:,1:npw1)=cg1(:,i1+1:i1+npw1)
2901 :
2902 : ! Make symmetry-induced conversion, if needed (translation part)
2903 75851 : if(conv_tnons==1)then
2904 : !$OMP PARALLEL DO PRIVATE(ai,ar)
2905 563587 : do ipw=1,npw1
2906 560472 : ar=phase3d(1,ipw)*wavef1(1,ipw)-phase3d(2,ipw)*wavef1(2,ipw)
2907 560472 : ai=phase3d(2,ipw)*wavef1(1,ipw)+phase3d(1,ipw)*wavef1(2,ipw)
2908 560472 : wavef1(1,ipw)=ar
2909 563587 : wavef1(2,ipw)=ai
2910 : end do
2911 : end if
2912 :
2913 : ! Take into account time-reversal symmetry, if needed, in the scalar case
2914 75851 : if(itimrev==1 .and. (nspinor1==1 .or. nspinor2==1))then
2915 : !$OMP PARALLEL DO
2916 4176504 : do ipw=1,npw1
2917 4235018 : wavef1(2,ipw)=-wavef1(2,ipw)
2918 : end do
2919 : end if
2920 :
2921 : ! DEBUG
2922 : ! write(std_out,*)' wfconv : before sphere, isym,ispinor=',isym,ispinor
2923 : ! write(std_out,*)' no_shift,identity=',no_shift,identity
2924 : ! write(std_out,*)' shiftg,symm=',shiftg,symm
2925 : ! stop
2926 : ! This debugging sequence is an attempt to rotate spinors,
2927 : ! and works indeed for test13, when symmetry 9 is used ...
2928 : ! if(isym==9 .and. ispinor==1)then
2929 : ! write(std_out,*)' wfconv : gives a 120 degree rotation to first component'
2930 : ! do ipw=1,npw1
2931 : ! ar=- half*wavef1(1,ipw)-sqrt(three)*half*wavef1(2,ipw)
2932 : ! ai= sqrt(three)*half*wavef1(1,ipw)- half*wavef1(2,ipw)
2933 : ! wavef1(1,ipw)=ar
2934 : ! wavef1(2,ipw)=ai
2935 : ! end do
2936 : ! end if
2937 : ! ENDDEBUG
2938 :
2939 : ! Convert wf, and also include the symmetry operation and shiftg.
2940 : call sphere(wavef1,1,npw1,cfft,n1,n2,n3,n4,n5,n6,kg1,istwf1_k,tobox,&
2941 75851 : & mpi_enreg1%me_g0,no_shift,identity,one)
2942 :
2943 : call sphere(wavef2,1,npw2,cfft,n1,n2,n3,n4,n5,n6,kg2,istwf2_k,tosph,&
2944 75851 : & mpi_enreg2%me_g0,shiftg,symm,one)
2945 :
2946 149960 : if(nspinor2==1 )then
2947 71910 : i2=(ispinor-1)*npw2+(iband-1)*nspinor2_this_proc*npw2+icg2
2948 67253655 : cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
2949 3941 : else if(nspinor1==2.and.nspinor2==2)then
2950 : ! Will treat this case outside of the ispinor loop
2951 3484 : i2=(ispinor1-1)*npw2
2952 1231492 : wavefspinor(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
2953 457 : else if(nspinor1==1 .and. nspinor2==2)then
2954 : ! The number of bands is doubled, and the number of coefficients
2955 : ! is doubled also
2956 457 : if (mpi_enreg2%paral_spinor==0) then
2957 457 : i2=(iband-1)*nspinor2_this_proc*nspinor2_this_proc*npw2+icg2
2958 481312 : cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
2959 481312 : cg2(:,i2+npw2+1:i2+2*npw2)=zero
2960 481312 : cg2(:,i2+2*npw2+1:i2+3*npw2)=zero
2961 481312 : cg2(:,i2+3*npw2+1:i2+4*npw2)=wavef2(:,1:npw2)
2962 : else
2963 0 : i2=(iband-1)*nspinor2_this_proc*npw2+icg2
2964 0 : if (nspinor_index==1) then
2965 0 : cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2)
2966 0 : cg2(:,i2+npw2+1:i2+2*npw2)=zero
2967 : else
2968 0 : cg2(:,i2+1:i2+npw2)=zero
2969 0 : cg2(:,i2+npw2+1:i2+2*npw2)=wavef2(:,1:npw2)
2970 : end if
2971 : end if
2972 : end if
2973 : end do ! ispinor=ispinor_first,ispinor_last,order
2974 :
2975 84623 : if(nspinor1==2.and.nspinor2==2)then
2976 : ! Take care of possible parallelization over spinors
2977 1742 : if (mpi_enreg2%paral_spinor==1) then
2978 0 : call xmpi_sum(wavefspinor,mpi_enreg2%comm_spinor,ierr)
2979 : end if
2980 : ! Take care of time-reversal symmetry, if needed
2981 1742 : if(itimrev==1)then
2982 : ! Exchange spin-up and spin-down
2983 : ! Make complex conjugate of one component,
2984 : ! and change sign of other component
2985 : !$OMP PARALLEL DO PRIVATE(ipw,ar,ai) SHARED(wavefspinor,npw2)
2986 2640 : do ipw=1,npw2
2987 : ! Here, change sign of real part
2988 2600 : ar=-wavefspinor(1,ipw)
2989 2600 : ai= wavefspinor(2,ipw)
2990 2600 : wavefspinor(1,ipw)= wavefspinor(1,npw2+ipw)
2991 : ! Here, change sign of imaginary part
2992 2600 : wavefspinor(2,ipw)=-wavefspinor(2,npw2+ipw)
2993 2600 : wavefspinor(1,npw2+ipw)=ar
2994 2640 : wavefspinor(2,npw2+ipw)=ai
2995 : end do
2996 : end if ! itimrev==1
2997 :
2998 : ! Rotation in spinor space
2999 : !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(npw2,spinrot,wavefspinor)
3000 1742 : spinrots=spinrot(1)
3001 1742 : spinrotx=spinrot(2)
3002 1742 : spinroty=spinrot(3)
3003 1742 : spinrotz=spinrot(4)
3004 : !$OMP DO
3005 206410 : do ipw=1,npw2
3006 204668 : ar=wavefspinor(1,ipw)
3007 204668 : ai=wavefspinor(2,ipw)
3008 204668 : br=wavefspinor(1,npw2+ipw)
3009 204668 : bi=wavefspinor(2,npw2+ipw)
3010 204668 : wavefspinor(1,ipw) = spinrots*ar - spinrotz*ai + spinroty*br - spinrotx*bi
3011 204668 : wavefspinor(2,ipw) = spinrots*ai + spinrotz*ar + spinroty*bi + spinrotx*br
3012 204668 : wavefspinor(1,npw2+ipw)= -spinroty*ar - spinrotx*ai + spinrots*br + spinrotz*bi
3013 206410 : wavefspinor(2,npw2+ipw)= -spinroty*ai + spinrotx*ar + spinrots*bi - spinrotz*br
3014 : end do
3015 : !$OMP END DO
3016 : !$OMP END PARALLEL
3017 :
3018 : ! Save wavefunction
3019 1742 : i2=(iband-1)*nspinor2_this_proc*npw2+icg2
3020 1742 : if (mpi_enreg2%paral_spinor==0) then
3021 615746 : cg2(:,i2 +1:i2+ npw2)=wavefspinor(:,1:npw2)
3022 615746 : cg2(:,i2+npw2+1:i2+2*npw2)=wavefspinor(:,npw2+1:2*npw2)
3023 : else
3024 0 : if (nspinor_index==1) then
3025 0 : cg2(:,i2+1:i2+npw2)=wavefspinor(:,1:npw2)
3026 : else
3027 0 : cg2(:,i2+1:i2+npw2)=wavefspinor(:,npw2+1:2*npw2)
3028 : end if
3029 : end if
3030 : end if ! nspinor1==2 .and. nspinor2==2
3031 :
3032 : end do
3033 :
3034 : ! Take care of copying eig and occ when nspinor increases or decreases
3035 10514 : if(nspinor1==1.and.nspinor2==2)then
3036 45 : if(formeig==0)then
3037 : ! Note the reverse order, needed in case inplace=1
3038 502 : do iband=nbd1,1,-1
3039 : ! use eig_tmp to avoid bug on ifort10.1 x86_64
3040 457 : eig_tmp=eig_k1(iband)
3041 1416 : eig_k2(2*iband-1:2*iband)=eig_tmp
3042 : ! occ_tmp=occ_k1(iband)*0.5_dp
3043 : ! occ_k2(2*iband-1:2*iband )=occ_tmp
3044 : end do
3045 : else
3046 0 : call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL")
3047 : end if
3048 : end if
3049 10514 : if(nspinor1==2 .and. nspinor2==1)then
3050 0 : if(formeig==0)then
3051 0 : do iband=1,nbd1
3052 : ! use eig_tmp to avoid bug on ifort10.1 x86_64
3053 0 : eig_tmp=eig_k1(2*iband-1)
3054 0 : eig_k2(iband)=eig_tmp
3055 : ! occ_tmp=occ_k1(2*iband-1)*2.0_dp
3056 : ! occ_k2(iband)=occ_tmp
3057 : end do
3058 : else
3059 0 : call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL")
3060 : end if
3061 : end if
3062 :
3063 10514 : ABI_FREE(cfft)
3064 10514 : ABI_FREE(gbound1)
3065 10514 : ABI_FREE(gbound2)
3066 10514 : ABI_FREE(wavef1)
3067 10514 : ABI_FREE(wavef2)
3068 10514 : if(nspinor1==2 .and. nspinor2==2) then
3069 88 : ABI_FREE(wavefspinor)
3070 : end if
3071 :
3072 : else if(convert==0)then
3073 :
3074 103604 : if(inplace==0)then
3075 : ! Must copy cg, eig and occ if not in-place while convert==0
3076 : ! Note that npw1=npw2, nspinor1=nspinor2
3077 : cg2(:,1+icg2:npw1*nspinor1_this_proc*nbd1+icg2)=&
3078 3896675 : & cg1(:,1+icg1:npw1*nspinor1_this_proc*nbd1+icg1)
3079 3890 : eig_k2(:)=eig_k1(:)
3080 : ! occ_k2(:)=occ_k1(:)
3081 : end if
3082 :
3083 : end if ! End of if convert/=0
3084 :
3085 114118 : if(conv_tnons==1) then
3086 443 : ABI_FREE(phase1d)
3087 443 : ABI_FREE(phase3d)
3088 : end if
3089 :
3090 :
3091 : !If not enough bands, complete with random numbers or zeros
3092 114118 : if(nbd2/nspinor2>nbd1/nspinor1)then
3093 38494 : if(formeig==0)then
3094 :
3095 : ! Ground state wf and eig case
3096 290618 : eig_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=huge(zero)/10.0_dp
3097 290618 : occ_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=0.0_dp
3098 38494 : index=(nbd1/nspinor1)*nspinor2*npw2*nspinor2_this_proc
3099 :
3100 : ! Initialisation of wavefunctions
3101 : ! One needs to initialize wfs in such a way to avoid symmetry traps,
3102 : ! and to avoid linear dependencies between wavefunctions
3103 : ! No need for a difference for different k points and/or spin-polarization
3104 :
3105 38494 : npwtot=npw2
3106 38494 : if (mpi_enreg1%paral_kgb == 1) then
3107 2559 : call timab(539,1,tsec)
3108 2559 : call xmpi_sum(npwtot, mpi_enreg1%comm_bandfft, ierr)
3109 2559 : call timab(539,2,tsec)
3110 : end if
3111 :
3112 290618 : do iband=(nbd1/nspinor1)*nspinor2+1,nbd2
3113 528808 : do ispinor2=1,nspinor2_this_proc
3114 276684 : ispinor=ispinor2;if (nspinor2_this_proc/=nspinor2) ispinor=mpi_enreg2%me_spinor+1
3115 :
3116 70935555 : do ipw=1,npw2
3117 70406747 : index=index+1
3118 : ! Different seed for different planewave and band
3119 : ! DEBUG seq==par
3120 : ! if(.false.) then
3121 : ! ENDDEBUG seq==par
3122 :
3123 70406747 : if ( mpi_enreg2%paral_kgb /= 1.or.mpi_enreg2%nproc_cell == 1) then
3124 65617151 : seed=(iband-1)*npw2*nspinor2 + (ispinor-1)*npw2 + ipw
3125 : else
3126 4789596 : seed=kg2(1,ipw)*npwtot*npwtot + kg2(2,ipw)*npwtot + kg2(3,ipw)
3127 4789596 : seed=(iband*nspinor2+ispinor-1)*seed
3128 : end if
3129 :
3130 70683431 : if(randalg == 0) then
3131 : ! For portability, use only integer numbers
3132 : ! The series of couples (fold1,fold2) is periodic with a period of
3133 : ! 3x5x7x11x13x17x19x23x29x31, that is, larger than 2**32, the largest integer*4
3134 : ! fold1 is between 0 and 34, fold2 is between 0 and 114. As sums of five
3135 : ! uniform random variables, their distribution is close to a gaussian
3136 67087555 : fold1=modulo(seed,3)+modulo(seed,5)+modulo(seed,7)+modulo(seed,11)+modulo(seed,13)
3137 67087555 : fold2=modulo(seed,17)+modulo(seed,19)+modulo(seed,23)+modulo(seed,29)+modulo(seed,31)
3138 :
3139 : ! The gaussian distributions are folded, in order to be back to a uniform distribution
3140 : ! foldre is between 0 and 20, foldim is between 0 and 18
3141 67087555 : foldre=mod(fold1+fold2,21)
3142 67087555 : foldim=mod(3*fold1+2*fold2,19)
3143 :
3144 67087555 : cg2(1,index+icg2)=dble(foldre)
3145 67087555 : cg2(2,index+icg2)=dble(foldim)
3146 : else
3147 : ! (Antoine Levitt) Simple linear congruential generator from
3148 : ! numerical recipes, modulo'ed and 64bit'ed to avoid
3149 : ! overflows (NAG doesn't like overflows, even though
3150 : ! they are perfectly legitimate here). Then, we get some
3151 : ! lowest order bits and sum them, as the previous
3152 : ! generator, to get quasi-normal numbers.
3153 : ! This is clearly suboptimal and might cause problems,
3154 : ! but at least it doesn't seem to create linear
3155 : ! dependencies and local minima like the previous one.
3156 : ! it's not trivial to generate good reproductible random
3157 : ! numbers in parallel. Patches welcome !
3158 : ! Note a fun fortran fact : MOD simply ignores 64 bits integer
3159 : ! and casts them into 32bits, so we use MODULO.
3160 3319192 : fold1 = modulo(1664525_int64 * seed + 1013904223_int64, 2147483648_int64)
3161 3319192 : fold2 = modulo(1664525_int64 * fold1 + 1013904223_int64, 2147483648_int64)
3162 3319192 : fold1=modulo(fold1,3)+modulo(fold1,5)+modulo(fold1,7)+modulo(fold1,11)+modulo(fold1,13)
3163 3319192 : fold2=modulo(fold2,3)+modulo(fold2,5)+modulo(fold2,7)+modulo(fold2,11)+modulo(fold2,13)
3164 :
3165 3319192 : cg2(1,index+icg2)=dble(fold1)/34-0.5
3166 3319192 : cg2(2,index+icg2)=dble(fold2)/34-0.5
3167 : end if
3168 : end do
3169 : end do
3170 :
3171 : ! XG030513: Time-reversal symmetry for k=gamma imposes zero imaginary part at G=0
3172 : ! XG: I do not know what happens for spin-orbit here.
3173 290618 : if (istwf2_k == 2 .and. mpi_enreg2%me_g0 == 1) then
3174 7484 : cg2(2,1+(iband-1)*npw2*nspinor2_this_proc+icg2)=zero
3175 : end if
3176 : end do ! iband
3177 :
3178 : ! Multiply with envelope function to reduce kinetic energy
3179 38494 : icgmod=icg2+npw2*nspinor2_this_proc*(nbd1/nspinor1)
3180 38494 : nbremn=nbd2-nbd1
3181 38494 : call cg_envlop(cg2,ecut2,gmet2,icgmod,kg2,kpoint2_sph,mcg2,nbremn,npw2,nspinor2_this_proc)
3182 :
3183 38494 : if(ikpt2<=nkpt_max)then
3184 30726 : write(msg,'(3(a,i6))')' wfconv:',nbremn,' bands initialized randomly with npw=',npw2,', for ikpt=',ikpt2
3185 30726 : call wrtout(std_out,msg)
3186 : end if
3187 :
3188 0 : else if(formeig==1)then
3189 :
3190 : ! For response function, put large numbers in the remaining of the
3191 : ! eigenvalue array (part of it was already filled in calling routine)
3192 : ! WARNING : Change of nspinor not yet coded
3193 0 : eig_k2(1+2*nbd1*nbd2 : 2*nbd2*nbd2)=huge(zero)/10.0_dp
3194 : ! Initialisation of wfs with 0 s
3195 0 : index=npw2*nbd1*nspinor2_this_proc
3196 0 : do iband=nbd1+1,nbd2
3197 0 : do ipw=1,npw2*nspinor2_this_proc
3198 0 : index=index+1
3199 0 : cg2(:,index+icg2)=zero
3200 : end do
3201 : end do
3202 :
3203 0 : if(ikpt2<=nkpt_max)then
3204 0 : nbremn=nbd2-nbd1
3205 0 : write(msg,'(3(a,i0))')' wfconv:',nbremn,' bands set=0 with npw=',npw2,', for ikpt=',ikpt2
3206 0 : call wrtout(std_out,msg)
3207 : end if
3208 :
3209 : end if ! End of initialisation to 0
3210 : end if
3211 :
3212 : !Orthogonalize GS wfs
3213 114118 : if (optorth==1.and.formeig==0.and.(mpi_enreg2%paral_kgb/=1.or.mpi_enreg2%nproc_cell==1)) then
3214 : !if (.True.) then
3215 51001 : ABI_MALLOC(dum,(2,0))
3216 51001 : ortalgo=0 !;ortalgo=3
3217 : call pw_orthon(icg2,0,istwf2_k,mcg2,0,npw2*nspinor2_this_proc,nbd2,ortalgo,dum,0,cg2,&
3218 51001 : & mpi_enreg2%me_g0,mpi_enreg2%comm_bandspinorfft)
3219 51001 : ABI_FREE(dum)
3220 : end if
3221 :
3222 114118 : end subroutine wfconv
3223 : !!***
3224 :
3225 : !!****f* m_inwffil/pareigocc
3226 : !! NAME
3227 : !! pareigocc
3228 : !!
3229 : !! FUNCTION
3230 : !! This subroutine transmit to all processors, using MPI:
3231 : !! - the eigenvalues and,
3232 : !! - if ground-state, the occupation numbers
3233 : !! (In fact, in the present status of the routine,
3234 : !! occupation numbers are NOT transmitted)
3235 : !! transmit_occ = 2 is used in case the occ should be transmitted.
3236 : !! Yet the code is not already written.
3237 : !!
3238 : !! INPUTS
3239 : !! formeig=format of eigenvalues (0 for GS, 1 for RF)
3240 : !! localrdwf=(for parallel case) if 1, the eig and occ initial values
3241 : !! are local to each machine, if 0, they are on proc me=0.
3242 : !! mband=maximum number of bands of the output wavefunctions
3243 : !! mpi_enreg=information about MPI parallelization
3244 : !! nband(nkpt*nsppol)=desired number of bands at each k point
3245 : !! nkpt=number of k points
3246 : !! nsppol=1 for unpolarized, 2 for spin-polarized, output wf file processors,
3247 : !! Warning : defined only when paralbd=1
3248 : !! transmit_occ/=2 transmit only eigenvalues, =2 for transmission of occ also
3249 : !! (yet transmit_occ=2 is not safe or finished at all)
3250 : !!
3251 : !! OUTPUT
3252 : !! (see side effects)
3253 : !!
3254 : !! SIDE EFFECTS
3255 : !! eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
3256 : !! occ(mband*nkpt*nsppol)=occupation (input or init to 0.0) NOT USED NOW
3257 : !!
3258 : !! NOTES
3259 : !! * The case paralbd=1 with formeig=0 is implemented, but not yet used.
3260 : !!
3261 : !! * The transmission of occ is not activated yet !
3262 : !!
3263 : !! * The routine takes the eigenvalues in the eigen array on one of the
3264 : !! processors that possess the wavefunctions, and transmit it to all procs.
3265 : !! If localrdwf==0, me=0 has the full array at start,
3266 : !! If localrdwf==1, the transfer might be more complex.
3267 : !!
3268 : !! * This routine should not be used for RF wavefunctions, since
3269 : !! it does not treat the eigenvalues as a matrix.
3270 : !!
3271 : !! SOURCE
3272 :
3273 8871 : subroutine pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,transmit_occ)
3274 :
3275 : !Arguments ------------------------------------
3276 : !scalars
3277 : integer,intent(in) :: formeig,localrdwf,mband,nkpt,nsppol,transmit_occ
3278 : type(MPI_type),intent(in) :: mpi_enreg
3279 : !arrays
3280 : integer,intent(in) :: nband(nkpt*nsppol)
3281 : real(dp),intent(inout) :: eigen(mband*(2*mband)**formeig*nkpt*nsppol)
3282 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
3283 :
3284 : !Local variables-------------------------------
3285 : !scalars
3286 : integer :: band_index,iband,ierr,ikpt,isppol,me,nbks,spaceComm
3287 : !character(len=500) :: msg
3288 : !arrays
3289 : real(dp) :: tsec(2)
3290 8871 : real(dp),allocatable :: buffer1(:),buffer2(:)
3291 :
3292 : ! *************************************************************************
3293 :
3294 : if(xmpi_paral==1)then
3295 :
3296 : ! Init mpi_comm
3297 8871 : spaceComm=mpi_enreg%comm_cell
3298 8871 : if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
3299 8871 : if(mpi_enreg%paral_hf==1) spaceComm=mpi_enreg%comm_kpt
3300 : ! Init me
3301 8871 : me=mpi_enreg%me_kpt
3302 :
3303 8871 : if(localrdwf==0)then
3304 7 : call xmpi_bcast(eigen,0,spaceComm,ierr)
3305 :
3306 8864 : else if(localrdwf==1)then
3307 :
3308 : ! Prepare transmission of eigen (and occ)
3309 26592 : ABI_MALLOC(buffer1,(2*mband**(formeig+1)*nkpt*nsppol))
3310 17728 : ABI_MALLOC(buffer2,(2*mband**(formeig+1)*nkpt*nsppol))
3311 10996874 : buffer1(:)=zero
3312 10996874 : buffer2(:)=zero
3313 :
3314 8864 : band_index=0
3315 19103 : do isppol=1,nsppol
3316 223774 : do ikpt=1,nkpt
3317 204671 : nbks=nband(ikpt+(isppol-1)*nkpt)
3318 :
3319 214910 : if(mpi_enreg%paralbd==0)then
3320 :
3321 202799 : if(formeig==0)then
3322 710246 : buffer1(2*band_index+1:2*band_index+nbks) = eigen(band_index+1:band_index+nbks)
3323 94701 : if(transmit_occ==2) then
3324 0 : buffer1(2*band_index+nbks+1:2*band_index+2*nbks) = occ(band_index+1:band_index+nbks)
3325 : end if
3326 94701 : band_index=band_index+nbks
3327 108098 : else if(formeig==1)then
3328 9810914 : buffer1(band_index+1:band_index+2*nbks**2) = eigen(band_index+1:band_index+2*nbks**2)
3329 108098 : band_index=band_index+2*nbks**2
3330 : end if
3331 :
3332 1872 : else if(mpi_enreg%paralbd==1)then
3333 :
3334 : ! Skip this k-point if not the proper processor
3335 1872 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbks,isppol,me)) then
3336 0 : if(formeig==0) then
3337 0 : band_index=band_index+nbks
3338 : else
3339 0 : band_index=band_index+2*nbks**2
3340 : end if
3341 : cycle
3342 : end if
3343 : ! Loop on bands
3344 9360 : do iband=1,nbks
3345 7488 : if(mpi_enreg%proc_distrb(ikpt, iband,isppol) /= me)cycle
3346 9360 : if(formeig==0)then
3347 1088 : buffer1(2*band_index+iband)=eigen(band_index+iband)
3348 : ! if(transmit_occ==2) buffer1(2*band_index+iband+nbdks)=occ(band_index+iband)
3349 6400 : else if (formeig==1)then
3350 : buffer1(band_index+(iband-1)*2*nbks+1:band_index+(iband-1)*2*nbks+2*nbks) = &
3351 57600 : & eigen(band_index+(iband-1)*2*nbks+1:band_index+(iband-1)*2*nbks+2*nbks)
3352 : end if
3353 : end do
3354 1872 : if(formeig==0)then
3355 272 : band_index=band_index+nbks
3356 : else
3357 1600 : band_index=band_index+2*nbks**2
3358 : end if
3359 : end if
3360 :
3361 : end do
3362 : end do
3363 :
3364 : ! Build sum of everything
3365 8864 : call timab(48,1,tsec)
3366 8864 : if(formeig==0)band_index=band_index*2
3367 8864 : call xmpi_sum(buffer1,buffer2,band_index,spaceComm,ierr)
3368 8864 : call timab(48,2,tsec)
3369 :
3370 8864 : band_index=0
3371 19103 : do isppol=1,nsppol
3372 223774 : do ikpt=1,nkpt
3373 204671 : nbks=nband(ikpt+(isppol-1)*nkpt)
3374 214910 : if(formeig==0)then
3375 711606 : eigen(band_index+1:band_index+nbks) = buffer2(2*band_index+1:2*band_index+nbks)
3376 94973 : if(transmit_occ==2) then
3377 0 : occ(band_index+1:band_index+nbks) = buffer2(2*band_index+nbks+1:2*band_index+2*nbks)
3378 : end if
3379 94973 : band_index=band_index+nbks
3380 109698 : else if(formeig==1)then
3381 9863714 : eigen(band_index+1:band_index+2*nbks**2) = buffer1(band_index+1:band_index+2*nbks**2)
3382 109698 : band_index=band_index+2*nbks**2
3383 : end if
3384 : end do
3385 : end do
3386 :
3387 8864 : ABI_FREE(buffer1)
3388 17728 : ABI_FREE(buffer2)
3389 : end if
3390 : end if
3391 :
3392 8871 : end subroutine pareigocc
3393 : !!***
3394 :
3395 : !!****f* m_inwffil/cg_from_atoms
3396 : !! NAME
3397 : !! cg_from_atoms
3398 : !!
3399 : !! FUNCTION
3400 : !! Initialize wave functions at a given (k-point, spin) using Bloch sums of atomic orbitals.
3401 : !!
3402 : !! INPUTS
3403 : !! ikpt,isppol=k-point index, spin index
3404 : !! rprimd(3,3)=Direct lattice vectors in Bohr.
3405 : !! xred(3,natom)=Atomic positions.
3406 : !! kg_k(3,npw_k)=reduced planewave coordinates.
3407 : !! dtset <type(dataset_type)>=all input variables for this dataset
3408 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
3409 : !! mpi_enreg=information about MPI parallelization
3410 : !! nband=number of bands at this k point for that spin polarization
3411 : !! npw=number of plane waves at this k point
3412 : !! my_nspinor=number of spinors treated by this MPI proc
3413 : !!
3414 : !! OUTPUT
3415 : !! eig(nband)=array for holding eigenvalues (hartree)
3416 : !!
3417 : !! SIDE EFFECTS
3418 : !! cg(2,*)=updated wavefunctions
3419 : !!
3420 : !! SOURCE
3421 :
3422 8 : subroutine cg_from_atoms(ikpt, isppol, rprimd, xred, kg_k, cg, dtset, psps, eig, gs_hamk, &
3423 : mpi_enreg, nband, npw, my_nspinor)
3424 :
3425 : use defs_datatypes, only : pseudopotential_type
3426 : use m_geometry, only : metric
3427 : use m_splines, only : splfit
3428 : use m_initylmg, only : initylmg_k
3429 : use m_hamiltonian, only : gs_hamiltonian_type
3430 : use m_getghc, only : getghc
3431 : use m_nonlop, only : nonlop
3432 : use m_pawcprj, only : pawcprj_type
3433 : use m_rmm_diis, only : subspace_rotation
3434 : use m_cgtools, only : cgpaw_cholesky, cgnc_cholesky
3435 :
3436 : !Arguments ------------------------------------
3437 : integer,intent(in) :: ikpt, isppol, nband, npw, my_nspinor
3438 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
3439 : type(dataset_type),intent(in) :: dtset
3440 : type(pseudopotential_type),intent(in) :: psps
3441 : type(mpi_type),intent(inout) :: mpi_enreg
3442 : real(dp),intent(in) :: rprimd(3,3), xred(3,dtset%natom)
3443 : integer, intent(in) :: kg_k(3,npw)
3444 : real(dp),intent(inout) :: cg(2,npw*my_nspinor,nband)
3445 : real(dp),intent(out) :: eig(nband)
3446 :
3447 : !Local variables-------------------------------
3448 : integer,parameter :: optder0 = 0, ider0 = 0, icg0 = 0
3449 : !integer,parameter :: idir0 = 0, !, type_calc0 = 0, option1 = 1, option2 = 2, tim_getghc = 0
3450 : integer :: npwsp !, ortalgo, ierr,
3451 : integer :: istwf_k, usepaw, mcg !, mgsc
3452 : integer :: me_g0, me_cell !prev_mixprec,
3453 : integer :: comm_bsf, savemem, ll
3454 : integer :: iatom, itypat, iln, ig, iband, ilmn, ilm, im, lnmax
3455 : real(dp) :: kpg1, kpg2, kpg3, kpgc1, kpgc2, kpgc3
3456 : complex(dp) :: cfact
3457 : logical :: supported !, use_fft_mixprec
3458 : real(dp) :: ucvol, arg ! cpu, wall, gflops,
3459 : !character(len=500) :: msg
3460 : !arrays
3461 : real(dp) :: gmet(3,3), gprimd(3,3), rmet(3,3), kpt(3), phase_l(2), ri(2)
3462 8 : real(dp) :: enlx(nband) !, tsec(2)
3463 4 : real(dp),allocatable :: ghc(:,:), gvnlxc(:,:)
3464 4 : real(dp),allocatable :: kpg_k(:,:), tphiq(:,:,:), sf(:,:) !,ph3d(:,:,:)
3465 4 : real(dp),allocatable :: ylm(:,:), ylm_gr(:,:,:), gsc(:,:)
3466 4 : real(dp),allocatable :: kpgnorm(:), wk_ffnl2(:)
3467 :
3468 : ! *************************************************************************
3469 :
3470 : ! Define useful vars.
3471 4 : usepaw = dtset%usepaw; istwf_k = gs_hamk%istwf_k
3472 4 : me_g0 = mpi_enreg%me_g0; comm_bsf = mpi_enreg%comm_bandspinorfft
3473 4 : npwsp = npw * my_nspinor; mcg = npwsp * nband !; mgsc = npwsp * nband * usepaw
3474 4 : me_cell = mpi_enreg%me_cell
3475 16 : kpt = dtset%kptns(:,ikpt)
3476 :
3477 4 : supported = .True.
3478 4 : if (dtset%usepaw == 0) then
3479 12 : lnmax = maxval(psps%nctab(:)%num_tphi)
3480 12 : supported = supported .and. minval(psps%nctab(:)%num_tphi) > 0
3481 : end if
3482 :
3483 : ! Test whether cg initialization from ps atomic orbitals is coded/supported.
3484 4 : if (dtset%nspinor == 2) supported = .False.
3485 4 : if (dtset%usepaw /= 0) supported = .False.
3486 4 : if (.not. supported) then
3487 0 : if (me_cell == 0 .and. ikpt == 1) then
3488 0 : call wrtout(std_out, " cg initialization from atomic orbitals not available. returning")
3489 : end if
3490 0 : return
3491 : end if
3492 :
3493 4 : if (me_cell == 0 .and. ikpt == 1) then
3494 1 : call wrtout(std_out, sjoin(" Initializing cg from atomic orbitals for ikpt:", itoa(ikpt), ", spin:", itoa(isppol)))
3495 : end if
3496 : !call cwtime(cpu, wall, gflops, "start")
3497 :
3498 4 : call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
3499 :
3500 16 : ABI_MALLOC(ylm, (npw, psps%mpsang**2))
3501 16 : ABI_MALLOC(ylm_gr, (npw, 3+6*(optder0/2), psps%mpsang**2))
3502 :
3503 4 : call initylmg_k(npw, psps%mpsang, optder0, rprimd, gprimd, kpt, kg_k, ylm, ylm_gr)
3504 4 : ABI_SFREE(ylm_gr)
3505 :
3506 : ! Compute nonlocal form factors at (k+G)
3507 : ! Note that we need to work with useylm = 1 to keep the m-dependency
3508 : ! even when Vnl is applied with Legendre polynomials (useylm = 0)
3509 :
3510 : ! Get |k+G|
3511 12 : ABI_MALLOC(kpgnorm, (npw))
3512 :
3513 : !$OMP PARALLEL DO PRIVATE(kpg1, kpg2, kpg3, kpgc1, kpgc2, kpgc3)
3514 9101 : do ig=1,npw
3515 9097 : kpg1=kpt(1)+dble(kg_k(1,ig))
3516 9097 : kpg2=kpt(2)+dble(kg_k(2,ig))
3517 9097 : kpg3=kpt(3)+dble(kg_k(3,ig))
3518 9097 : kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
3519 9097 : kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
3520 9097 : kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
3521 9101 : kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
3522 : end do
3523 :
3524 20 : ABI_MALLOC(tphiq, (npw, lnmax, psps%ntypat))
3525 8 : ABI_MALLOC(wk_ffnl2, (npw))
3526 :
3527 12 : do itypat=1,psps%ntypat
3528 28 : do iln=1,psps%nctab(itypat)%num_tphi
3529 : call splfit(psps%qgrid_ff, wk_ffnl2, psps%nctab(itypat)%tphi_qspl(:,:,iln), &
3530 24 : ider0, kpgnorm, tphiq(:,iln,itypat), psps%mqgrid_ff, npw)
3531 : end do
3532 : end do
3533 :
3534 4 : ABI_FREE(kpgnorm)
3535 4 : ABI_FREE(wk_ffnl2)
3536 :
3537 : !call getph(atindx, natom, n1, n2, n3, ph1d, xred)
3538 : !call ph1d3d(iatom, jatom, kg_k, matblk, natom, npw_k, n1, n2, n3, phkxred, ph1d, ph3d)
3539 :
3540 : ! Now init cg. Assuming cg has been already filled with random numbers previously
3541 : ! so we only need to init the first states. We don't take into account the occupancies
3542 : ! in the isolated atom. We just loop over all nlm states until we have filled max nband states.
3543 12 : ABI_MALLOC(sf, (2, npw))
3544 36 : iband = 0
3545 36 : iatom_loop: do iatom=1,dtset%natom
3546 36 : itypat = dtset%typat(iatom)
3547 :
3548 : ! Structure factor.
3549 81909 : do ig=1,npw
3550 327492 : arg = -two_pi * dot_product(xred(:,iatom), kpt + kg_k(:,ig))
3551 81873 : sf(1,ig) = cos(arg)
3552 81909 : sf(2,ig) = sin(arg)
3553 : end do
3554 :
3555 36 : ilmn = 0
3556 104 : do iln=1,psps%nctab(itypat)%num_tphi
3557 : !if (psps%nctab(itypat)%tphi_occ(iln) < zero) cycle
3558 68 : ll = psps%nctab(itypat)%tphi_l(iln)
3559 : !cfact = (j_dpc ** ll) * four_pi / sqrt(ucvol)
3560 68 : cfact = (-j_dpc ** ll) * four_pi / sqrt(ucvol)
3561 68 : phase_l(1) = dble(cfact)
3562 68 : phase_l(2) = aimag(cfact)
3563 228 : do im=1, 2*ll+1
3564 132 : ilmn = ilmn + 1
3565 132 : ilm = im + ll**2
3566 132 : iband = iband + 1
3567 : ! Another good reason why nband should be > nbocc.
3568 132 : if (iband > nband) exit iatom_loop
3569 :
3570 : ! NB: Assuming nspinor == 1
3571 291232 : do ig=1,npw ! *my_nspinor
3572 291104 : ri(1) = phase_l(1) * sf(1, ig) - phase_l(2) * sf(2, ig)
3573 291104 : ri(2) = phase_l(1) * sf(2, ig) + phase_l(2) * sf(1, ig)
3574 291104 : if (dtset%wfinit == 1) call randomize(ri)
3575 873440 : cg(:, ig, iband) = ri(:) * ylm(ig, ilm) * tphiq(ig, iln, itypat)
3576 : !wfcatom (ig, 1, n_starting_wfc) = phase_l * sf(1, ig) * ylm(ig, ilm) * chiq(ig, nb, nt)
3577 : end do ! ig
3578 :
3579 : ! XG030513: Time-reversal symmetry for k=gamma imposes zero imaginary part at G=0
3580 : ! XG: I do not know what happens for spin-orbit here.
3581 192 : if (istwf_k == 2 .and. mpi_enreg%me_g0 == 1) cg(2, 1, iband) = zero
3582 : end do ! im
3583 : end do ! iln
3584 : end do iatom_loop
3585 :
3586 : !call cg_envlop(cg, dtset%ecut, gmet, icg0, kg_k, kpt, mcg, nband, npw, my_nspinor)
3587 :
3588 4 : ABI_FREE(sf)
3589 4 : ABI_FREE(tphiq)
3590 : ABI_SFREE(kpg_k)
3591 4 : ABI_FREE(ylm)
3592 :
3593 : ! Use mixed precisions if requested by the user but only for low accuracy_level
3594 : !use_fft_mixprec = dtset%mixprec == 1 .and. accuracy_level < 2
3595 : !if (use_fft_mixprec) prev_mixprec = fftcore_set_mixprec(1)
3596 :
3597 : ! =========================
3598 : ! === Subspace rotation ===
3599 : ! =========================
3600 4 : savemem = 1
3601 12 : ABI_MALLOC(gsc, (2, npw*my_nspinor*nband*dtset%usepaw))
3602 : call subspace_rotation(gs_hamk, dtset%prtvol, mpi_enreg, nband, npw, my_nspinor, savemem, &
3603 4 : enlx, eig, cg, gsc, ghc, gvnlxc)
3604 :
3605 4 : ABI_SFREE(ghc)
3606 4 : ABI_SFREE(gvnlxc)
3607 :
3608 : ! Revert mixprec to previous status before returning.
3609 : !if (use_fft_mixprec) prev_mixprec = fftcore_set_mixprec(prev_mixprec)
3610 :
3611 : ! Ortoghonalization is in principle not needed but it seems to improve a bit.
3612 4 : if (dtset%wfinit < 0) then
3613 : !ortalgo = 3 !; ortalgo = mpi_enreg%paral_kgb
3614 : !call pw_orthon(0, 0, istwf_k, mcg, mgsc, npwsp, nband, ortalgo, gsc, usepaw, cg, me_g0, comm_bsf)
3615 :
3616 : ! TODO: Merge the two routines.
3617 : if (usepaw == 1) then
3618 : call cgpaw_cholesky(npwsp, nband, cg, gsc, istwf_k, me_g0, comm_bsf)
3619 : else
3620 0 : call cgnc_cholesky(npwsp, nband, cg, istwf_k, me_g0, comm_bsf, use_gemm=.False.)
3621 : end if
3622 : end if
3623 :
3624 4 : ABI_FREE(gsc)
3625 : !call cwtime_report(" cg_from_atoms:", cpu, wall, gflops)
3626 :
3627 : contains
3628 0 : subroutine randomize(ri)
3629 : real(dp),intent(inout) :: ri(2)
3630 :
3631 : !Local variables-------------------------------
3632 : real(dp) :: arg, rr
3633 : complex(dp) :: c2, c1
3634 : ! *************************************************************************
3635 :
3636 0 : call random_number(arg)
3637 0 : arg = two_pi * arg
3638 0 : call random_number(rr)
3639 0 : c1 = cmplx(ri(1), ri(2), kind=dp)
3640 0 : c2 = one + 0.05_dp * cmplx(rr*cos(arg), rr*sin(arg), kind=dp)
3641 0 : c2 = c1 * c2
3642 0 : ri(1) = real(c2)
3643 0 : ri(2) = aimag(c2)
3644 :
3645 4 : end subroutine randomize
3646 :
3647 : end subroutine cg_from_atoms
3648 : !!***
3649 :
3650 : end module m_inwffil
3651 : !!***
|
