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