Branch data Line data Source code
1 : : !!****m* ABINIT/m_io_kss
2 : : !! NAME
3 : : !! m_io_kss
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains procedured dealing with the IO of the KSS file.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (MG, MT, VO, AR, LR, RWG, MM, XG, RShaltaf)
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 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : MODULE m_io_kss
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_xmpi
28 : : use m_errors
29 : : use m_nctk
30 : : use netcdf
31 : : use m_hdr
32 : : use m_wfk
33 : : use m_cgtools
34 : : use m_hamiltonian
35 : : use m_electronpositron
36 : : use m_pawtab
37 : : use m_paw_ij
38 : : use m_pawcprj
39 : : use m_pawfgr
40 : : use m_dtfil
41 : : use m_dtset
42 : :
43 : : use defs_datatypes, only : pseudopotential_type
44 : : use defs_abitypes, only : MPI_type
45 : : use m_time, only : timab
46 : : use m_io_tools, only : open_file
47 : : use m_fstrings, only : sjoin, itoa, strcat
48 : : use m_hide_lapack, only : xheevx_cplex, xhegvx_cplex
49 : : use m_geometry, only : metric, remove_inversion
50 : : use m_mpinfo, only : destroy_mpi_enreg, proc_distrb_cycle
51 : : use m_fftcore, only : get_kg, sphere
52 : : use m_fft, only : fftpac
53 : : use m_crystal , only : crystal_t
54 : : use m_gsphere, only : table_gbig2kg, merge_and_sort_kg
55 : : use m_kg, only : mkkin, mkkpg
56 : : use m_ksdiago, only : ksdiago, init_ddiago_ctl, ddiago_ctl_type
57 : : use m_mkffnl, only : mkffnl
58 : : use m_getghc, only : getghc
59 : : use m_fourier_interpol, only : transgrid
60 : :
61 : : implicit none
62 : :
63 : : private
64 : :
65 : : public :: write_kss_header ! Writes the header of the KSS file.
66 : : !private :: write_vkb ! Writes the KB form factors and derivates on file for a single k-point.
67 : : public :: write_kss_wfgk ! Write the Gamma-centered wavefunctions and energies on the KSS file for a single k-point.
68 : : public :: k2gamma_centered ! Convert a set of wavefunctions from the k-centered to the gamma-centered basis set.
69 : : public :: make_gvec_kss ! Build the list of G-vectors for the KSS file.
70 : : public :: outkss ! Generate KSS file
71 : :
72 : : CONTAINS !===========================================================
73 : : !!***
74 : :
75 : : !!****f* m_io_kss/write_kss_header
76 : : !! NAME
77 : : !! write_kss_header
78 : : !!
79 : : !! FUNCTION
80 : : !! Write the header of the KSS file either using plain Fortran-IO or netcdf with ETSF-IO format.
81 : : !! Returns the unit number to be used for further writing.
82 : : !! It should be executed by master node only.
83 : : !!
84 : : !! INPUTS
85 : : !! filekss(len=fnlen)=The name of the KSS file.
86 : : !! kss_npw=Number of planewaves used for the wavefunctions in the KSS files.
87 : : !! ishm=Max number of shells written on file
88 : : !! shlim(ishm)=The cumulative number of G"s in each shell.
89 : : !! nbandksseff=Number of bands to be written.
90 : : !! mband=The maximum number of bands treated by abinit.
91 : : !! nsym2=Number of symmetry operations to be written on the header.
92 : : !! symrel2(3,3,nsym2)=The symmetry operations in real space to be written.
93 : : !! tnons2(3,nsym2)=The fractional translations associateed to symrel2.
94 : : !! gbig(3,kss_npw)=The set of G-vectors for the KSS wavefunctions (Gamma-centered)
95 : : !! Hdr<hdr_type>=The abinit header.
96 : : !! Dtset <dataset_type>=all input variables for this dataset
97 : : !! Psps<pseudopotential_type>=Structure gathering info on the pseudopotentials.
98 : : !! iomode=Input variables specifying the fileformat. (0-->Fortran,3-->netcdf with ETSF-IO format).
99 : : !! occ(mband*nkpt*nsppol)=The occupation factors.
100 : : !!
101 : : !! OUTPUT
102 : : !! kss_unt=The unit number of the opened file.
103 : : !!
104 : : !! SIDE EFFECTS
105 : : !! The KSS Header is written on file.
106 : : !!
107 : : !! SOURCE
108 : :
109 : 3 : subroutine write_kss_header(filekss,kss_npw,ishm,nbandksseff,mband,nsym2,symrel2,tnons2,occ,gbig,shlim,&
110 : : crystal,Dtset,Hdr,Psps,iomode,kss_unt)
111 : :
112 : : !Arguments ------------------------------------
113 : : !scalars
114 : : integer,intent(in) :: iomode,kss_npw,nbandksseff,ishm,nsym2,mband
115 : : integer,intent(out) :: kss_unt
116 : : character(len=fnlen),intent(in) :: filekss
117 : : type(crystal_t),intent(in) :: crystal
118 : : type(pseudopotential_type),intent(in) :: Psps
119 : : type(Hdr_type),intent(in) :: Hdr
120 : : type(Dataset_type),intent(in) :: Dtset
121 : : !arrays
122 : : integer,intent(in) :: symrel2(3,3,nsym2)
123 : : integer,intent(in) :: gbig(3,kss_npw),shlim(ishm)
124 : : real(dp),intent(in) :: tnons2(3,nsym2)
125 : : real(dp),intent(in) :: occ(mband*Dtset%nkpt*Dtset%nsppol)
126 : :
127 : : !Local variables-------------------------------
128 : : !scalars
129 : : integer :: nspinor,nsppol,nkpt,itypat,ierr
130 : : integer :: nb,isppol,ik,fform,ii,jj,kk,ig
131 : : integer :: il,il0,ilmn,in,ind1,ind2
132 : : character(len=80) :: title
133 : : character(len=500) :: msg
134 : 3 : type(hdr_type) :: my_Hdr
135 [ + + + + : 384 : type(dataset_type) :: Dtset_cpy
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + +
+ ]
136 : : integer :: ncerr
137 : : !arrays
138 : 3 : integer,allocatable :: vkbsign_int(:,:,:)
139 : 3 : real(dp),allocatable :: vkbsign(:,:)
140 : :
141 : : ! *********************************************************************
142 : :
143 : : DBG_ENTER("COLL")
144 : :
145 : 3 : nsppol = Dtset%nsppol
146 : 3 : nkpt = Dtset%nkpt
147 : 3 : nspinor= Dtset%nspinor
148 : :
149 : 3 : write(msg,'(3a)')ch10,' Opening file for KS structure output: ',TRIM(filekss)
150 : 3 : call wrtout(std_out,msg,'COLL')
151 : :
152 : 3 : write(msg,'(a,i6)') ' number of Gamma centered plane waves ',kss_npw
153 : 3 : call wrtout(std_out,msg,'COLL')
154 : 3 : call wrtout(ab_out,msg,'COLL')
155 : 3 : write(msg,'(a,i6)') ' number of Gamma centered shells ',ishm
156 : 3 : call wrtout(std_out,msg,'COLL')
157 : 3 : call wrtout(ab_out,msg,'COLL')
158 : 3 : write(msg,'(a,i6)') ' number of bands ',nbandksseff
159 : 3 : call wrtout(std_out,msg,'COLL')
160 : 3 : call wrtout(ab_out,msg,'COLL')
161 : 3 : write(msg,'(a,i6)') ' maximum angular momentum components ',Psps%mpsang
162 : 3 : call wrtout(std_out,msg,'COLL')
163 : 3 : call wrtout(ab_out,msg,'COLL')
164 : 3 : write(msg,'(a,i2,a)')' number of symmetry operations ',nsym2,' (without inversion)'
165 : 3 : call wrtout(std_out,msg,'COLL')
166 : :
167 : : !Copy the header so that we can change some basic dimensions using the KSS values:
168 : : !(bantot, npwarr, nband) and the occupation factors
169 : :
170 : : !Note that nsym and symrel might have been changed this has to be fixed
171 : : !carefully in the next patch since in the new implementation symmorphy=0 should be dafault
172 : 3 : call hdr%copy(my_Hdr)
173 : :
174 [ + + ]: 34 : my_Hdr%npwarr =kss_npw
175 [ + + ]: 34 : my_Hdr%nband =nbandksseff
176 [ + + ]: 34 : my_hdr%mband = maxval(my_hdr%nband)
177 : 3 : my_Hdr%bantot =nbandksseff*nkpt*nsppol
178 : :
179 [ + + ]: 34 : my_Hdr%istwfk = 1 ! KSS file does not support istwfk/=1 even though the GS run
180 : : ! can take advantage of time-reversal symmetry.
181 : :
182 : : !Copy the occ number in the new header with correct dimensions
183 : : !fill with zero the rest since mband can be < nbandksseff
184 : : !write(std_out,*)associated(my_Hdr%occ)
185 [ - + ]: 3 : ABI_FREE(my_Hdr%occ)
186 [ - + + - : 9 : ABI_MALLOC(my_Hdr%occ,(my_Hdr%bantot))
+ - + - -
+ - + ]
187 : : !mband = MAXVAL(Hdr%nband)
188 : :
189 [ + + ]: 1323 : my_Hdr%occ=zero; nb=MIN(mband,nbandksseff)
190 [ + + ]: 6 : do isppol=1,nsppol
191 [ + + ]: 37 : do ik=1,nkpt
192 : 31 : ind1=1+(ik-1)*nbandksseff+(isppol-1)*nkpt*nbandksseff
193 : 31 : ind2=1+(ik-1)*mband +(isppol-1)*nkpt*mband
194 [ + + ]: 444 : my_Hdr%occ(ind1:ind1+nb-1) = occ(ind2:ind2+nb-1)
195 : : end do
196 : : end do
197 : :
198 : : !Change dimension in the local Dtset_cpy as well.
199 : 3 : dtset_cpy = Dtset%copy()
200 : 3 : Dtset_cpy%mpw = kss_npw
201 : 3 : Dtset_cpy%mband = nbandksseff
202 : :
203 : 3 : fform=502
204 : :
205 : 3 : SELECT CASE (iomode)
206 : :
207 : : CASE (IO_MODE_FORTRAN)
208 : :
209 [ # # ]: 0 : if (open_file(filekss, msg, newunit=kss_unt, form="unformatted") /= 0) then
210 : 0 : ABI_ERROR(msg)
211 : : end if
212 : :
213 : 0 : call my_hdr%fort_write(kss_unt, fform, ierr)
214 [ # # ]: 0 : ABI_CHECK(ierr == 0, "hdr_Fort_write returned ierr != 0")
215 : :
216 : 0 : title='Results from ABINIT code'; write(kss_unt) title(1:80)
217 : 0 : title='Ab-initio plane waves calculation'; write(kss_unt) title(1:80)
218 : :
219 : 0 : write(kss_unt) nsym2,nbandksseff,kss_npw,ishm,Psps%mpsang ! To be modified to deal with more than one projector
220 : 0 : write(kss_unt) (((symrel2(ii,jj,kk),ii=1,3),jj=1,3),kk=1,nsym2)
221 : 0 : write(kss_unt) ((tnons2(ii,kk),ii=1,3),kk=1,nsym2)
222 : 0 : write(kss_unt) ((gbig(ii,ig),ii=1,3),ig=1,kss_npw)
223 : 0 : write(kss_unt) (shlim(in),in=1,ishm)
224 : :
225 : : ! Write vkbsign for NC pseudos with Fortran IO
226 : : ! MG FIXME: only one projector in each angular channel is treated.
227 : : ! Moreover the allocation is done in the wrong order for dimensions...
228 : : ! but if I change this code, compatibility with external codes is broken.
229 [ # # ]: 0 : if (Psps%usepaw==0) then
230 [ # # # # : 0 : ABI_MALLOC(vkbsign,(Psps%ntypat,Psps%mpsang))
# # # # #
# # # # #
# # ]
231 [ # # # # ]: 0 : vkbsign(:,:)=zero
232 [ # # ]: 0 : do itypat=1,Psps%ntypat
233 : 0 : il0=0
234 [ # # ]: 0 : do ilmn=1,Psps%lmnmax
235 : 0 : il=1+Psps%indlmn(1,ilmn,itypat)
236 : 0 : in=Psps%indlmn(3,ilmn,itypat)
237 [ # # ]: 0 : if (il/=il0 .and. in==1) then
238 : 0 : il0=il
239 : 0 : vkbsign(itypat,il)=DSIGN(one,Psps%ekb(ilmn,itypat))
240 : : end if
241 : : end do
242 : : end do
243 [ # # # # ]: 0 : write(kss_unt) ((vkbsign(itypat,il),il=1,Psps%mpsang),itypat=1,Psps%ntypat)
244 [ # # ]: 0 : ABI_FREE(vkbsign)
245 : : end if
246 : :
247 : : CASE (IO_MODE_ETSF)
248 : :
249 : : ! Create file.
250 [ - + ]: 3 : NCF_CHECK(nctk_open_create(kss_unt, nctk_ncify(filekss), xmpi_comm_self))
251 : :
252 : : ! Add additional info from abinit header.
253 [ - + ]: 3 : NCF_CHECK(my_hdr%ncwrite(kss_unt, fform, nc_define=.True.))
254 : :
255 : : ! Add info on crystalline structure
256 : : ! FIXME: Check symmorphi trick and crystal%symrel!
257 : : ! We currently use the dataset symmetries, as defined in the Hdr structure
258 : : ! instead of the symmetries recomputed in outkss.
259 [ - + ]: 3 : NCF_CHECK(crystal%ncwrite(kss_unt))
260 : :
261 : : ! Defined G-vectors and wavefunctions.
262 : 3 : call wfk_ncdef_dims_vars(kss_unt, my_hdr, fform, iskss=.True.)
263 : : !call abi_etsf_init(Dtset_cpy, filekss, 4, .false., my_Hdr%lmn_size, Psps, Dummy_wfs)
264 : :
265 : : ! If NC pseudos, write vkbsign.
266 : : ! Here multi-projectors are supported, array is dimensioned according to etsf-io standard.
267 [ + + ]: 3 : if (psps%usepaw == 0) then
268 : :
269 : : ! Define dims and variables needed for KB matrix elements.
270 : : ncerr = nctk_def_dims(kss_unt, [ &
271 : : nctkdim_t("max_number_of_angular_momenta", psps%mpsang), &
272 : : nctkdim_t("max_number_of_projectors", psps%mproj) &
273 [ + + ]: 6 : ])
274 [ - + ]: 2 : NCF_CHECK(ncerr)
275 : :
276 : : ncerr = nctk_def_arrays(kss_unt, [ &
277 : : nctkarr_t("kb_formfactor_sign", "int", &
278 : : &"max_number_of_projectors, max_number_of_angular_momenta, number_of_atom_species"), &
279 : : nctkarr_t("kb_formfactors", "dp", &
280 : : &"max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,&
281 : : &max_number_of_angular_momenta, number_of_atom_species"), &
282 : : nctkarr_t("kb_formfactor_derivative", "dp", &
283 : : &"max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,&
284 : : &max_number_of_angular_momenta, number_of_atom_species") &
285 [ + + ]: 8 : ])
286 [ - + ]: 2 : NCF_CHECK(ncerr)
287 : :
288 [ - + + - : 10 : ABI_MALLOC(vkbsign_int, (psps%mproj, Psps%mpsang, Psps%ntypat))
- + + - -
+ + - + -
+ - - + -
+ ]
289 [ + + + + : 19 : vkbsign_int=0
+ + ]
290 [ + + ]: 5 : do itypat=1,Psps%ntypat
291 [ + + ]: 9 : do ilmn=1,Psps%lmnmax
292 : 4 : il=1+Psps%indlmn(1,ilmn,itypat)
293 : 4 : in=Psps%indlmn(3,ilmn,itypat)
294 : 7 : vkbsign_int(in,il,itypat)=NINT(DSIGN(one,Psps%ekb(ilmn,itypat)))
295 : : end do
296 : : end do
297 : :
298 [ - + ]: 2 : NCF_CHECK(nctk_set_datamode(kss_unt))
299 : :
300 : : ! Write KB sign here
301 [ - + ]: 2 : NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "kb_formfactor_sign"), vkbsign_int))
302 [ - + ]: 2 : ABI_FREE(vkbsign_int)
303 : : end if
304 : :
305 [ - + ]: 3 : NCF_CHECK(nctk_set_datamode(kss_unt))
306 : :
307 : : CASE DEFAULT
308 [ - + - ]: 3 : ABI_ERROR(sjoin("Unsupported value for iomode:", itoa(iomode)))
309 : : END SELECT
310 : :
311 : 3 : call Dtset_cpy%free()
312 : 3 : call my_Hdr%free()
313 : :
314 : : DBG_EXIT("COLL")
315 : :
316 [ - + - + : 3 : end subroutine write_kss_header
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
- - - - -
- - - - -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
317 : : !!***
318 : :
319 : : !----------------------------------------------------------------------
320 : :
321 : : !!****f* m_io_kss/write_vkb
322 : : !! NAME
323 : : !! write_vkb
324 : : !!
325 : : !! FUNCTION
326 : : !! Writes the KB form factors and derivates on file for a single k-point.
327 : : !! Supports plain Fortran IO and netcdf with ETSF-IO format
328 : : !!
329 : : !! INPUTS
330 : : !! kss_unt=The unit number of the file
331 : : !! ikpt=The index of the k-point
332 : : !! kpoint(3)=The k-point in reduced coordinates.
333 : : !! kss_npw=Number of planewaves used for the wavefunctions in the KSS files.
334 : : !! npw_k=Number of planewaves at this k-point in the k-centered basis set used in abinit (ecut).
335 : : !! ecut=cutoff energy used in abinit.
336 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr).
337 : : !! Psps<Pseudopotential_type>=Datatype gathering data on the Pseudopotentials.
338 : : !! iomode=Input variables specifying the fileformat. (0-->Fortran,3-->netcdf with ETSF-IO).
339 : : !! gbig(3,kss_npw)=Set of G-vectors used in the KSS file.
340 : : !!
341 : : !! OUTPUT
342 : : !! Only writing.
343 : : !!
344 : : !! SOURCE
345 : :
346 : 20 : subroutine write_vkb(kss_unt,ikpt,kpoint,kss_npw,gbig,rprimd,Psps,iomode)
347 : :
348 : : !Arguments ------------------------------------
349 : : !scalars
350 : : integer,intent(in) :: ikpt,iomode,kss_npw,kss_unt
351 : : type(Pseudopotential_type),intent(in) :: Psps
352 : : !arrays
353 : : integer,intent(in) :: gbig(3,kss_npw)
354 : : real(dp),intent(in) :: kpoint(3),rprimd(3,3)
355 : :
356 : : !Local variables-------------------------------
357 : : !scalars
358 : : integer :: itypat,il,ig,mpsang,ntypat
359 : : !array
360 : 20 : real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:)
361 : 20 : real(dp),allocatable :: dum_vkbsign(:,:)
362 : : integer :: ncerr,varid
363 : 20 : real(dp),allocatable,target :: vkb_tgt(:,:,:,:), vkbd_tgt(:,:,:,:)
364 : : ! *********************************************************************
365 : :
366 : 20 : mpsang = Psps%mpsang; ntypat = Psps%ntypat
367 : :
368 [ - + + - : 100 : ABI_MALLOC(vkb ,(kss_npw,ntypat,mpsang))
- + + - -
+ + - + -
+ - - + -
+ ]
369 [ - + + - : 80 : ABI_MALLOC(vkbd,(kss_npw,ntypat,mpsang))
- + + - -
+ + - - +
- + ]
370 [ - + + - : 80 : ABI_MALLOC(dum_vkbsign,(ntypat,mpsang))
- + + - +
- + - - +
- + ]
371 : :
372 : 20 : call kss_calc_vkb(Psps,kpoint,kss_npw,gbig,rprimd,dum_vkbsign,vkb,vkbd)
373 [ - + ]: 20 : ABI_FREE(dum_vkbsign)
374 : :
375 : 20 : SELECT CASE (iomode)
376 : :
377 : : CASE (IO_MODE_FORTRAN)
378 [ # # ]: 0 : do itypat=1,ntypat
379 [ # # ]: 0 : do il=1,mpsang
380 : 0 : write(kss_unt) (vkb (ig,itypat,il),ig=1,kss_npw)
381 : 0 : write(kss_unt) (vkbd(ig,itypat,il),ig=1,kss_npw)
382 : : end do
383 : : end do
384 : :
385 : : CASE (IO_MODE_ETSF)
386 [ - + + - : 100 : ABI_MALLOC(vkb_tgt ,(kss_npw,1,mpsang,ntypat))
- + + - -
+ + - + -
+ - - + -
+ ]
387 [ - + + - : 80 : ABI_MALLOC(vkbd_tgt,(kss_npw,1,mpsang,ntypat))
- + + - -
+ + - - +
- + ]
388 [ + + ]: 50 : do itypat=1,ntypat
389 [ + + ]: 120 : do il=1,mpsang
390 [ + + ]: 89990 : do ig=1,kss_npw
391 : 89890 : vkb_tgt (ig,1,il,itypat)=vkb (ig,itypat,il)
392 : 89960 : vkbd_tgt(ig,1,il,itypat)=vkbd(ig,itypat,il)
393 : : end do
394 : : end do
395 : : end do
396 : :
397 : : ! FIXME: Multiple projectors
398 : : !ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat))
399 : : !call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd)
400 : : !ABI_FREE(vkbsign)
401 : : !ABI_FREE(vkb)
402 : : !ABI_FREE(vkbd)
403 : :
404 : : ! The shape of the variable on disk is (Fortran API):
405 : : ! (max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,
406 : : ! max_number_of_angular_momenta, number_of_atom_species)
407 : 20 : varid = nctk_idname(kss_unt, "kb_formfactors")
408 [ + + + + ]: 220 : ncerr = nf90_put_var(kss_unt, varid, vkb_tgt, start=[1,ikpt,1,1,1], count=[kss_npw,1,1,mpsang,ntypat])
409 [ - + ]: 20 : NCF_CHECK(ncerr)
410 : :
411 : 20 : varid = nctk_idname(kss_unt, "kb_formfactor_derivative")
412 [ + + + + ]: 220 : ncerr = nf90_put_var(kss_unt, varid, vkbd_tgt, start=[1,ikpt,1,1,1], count=[kss_npw,1,1,mpsang,ntypat])
413 [ - + ]: 20 : NCF_CHECK(ncerr)
414 : :
415 [ - + ]: 20 : ABI_FREE(vkb_tgt)
416 [ - + ]: 20 : ABI_FREE(vkbd_tgt)
417 : :
418 : : CASE DEFAULT
419 [ - + - ]: 20 : ABI_ERROR(sjoin("Unsupported value for iomode:", itoa(iomode)))
420 : : END SELECT
421 : :
422 [ - + ]: 20 : ABI_FREE(vkb)
423 [ - + ]: 20 : ABI_FREE(vkbd)
424 : :
425 [ - + - + ]: 20 : end subroutine write_vkb
426 : : !!***
427 : :
428 : : !----------------------------------------------------------------------
429 : :
430 : : !!****f* m_io_kss/write_kss_wfgk
431 : : !! NAME
432 : : !! write_kss_wfgk
433 : : !!
434 : : !! FUNCTION
435 : : !! Write the Gamma-centered wavefunctions and energies on the KSS file for a single k-point.
436 : : !! (Only the master node should call this routine).
437 : : !!
438 : : !! INPUTS
439 : : !! kss_unt=The unit number of the file
440 : : !! ikpt=The index of the k-point
441 : : !! isppol=The spin index.
442 : : !! nspinor=number of spinorial components (on current proc)
443 : : !! kss_npw=Number of planewaves used for the wavefunctions in the KSS files.
444 : : !! npw_k=Number of plane-waves in the k-centered basis set.
445 : : !! nbandksseff=Number of bands to be written.
446 : : !! natom=Number of atoms.
447 : : !! Psps<Pseudopotential_type>=Structure gathering pseudopotential data.
448 : : !! kpoint(3)=The k-points in reduced coordinates.
449 : : !! ene_k(nbandksseff)=Energies at this k-point
450 : : !! occ_k(nbandksseff)=Occupation factors at this k-point.
451 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr).
452 : : !! gbig(3,kss_npw)=The set of G-vectors for the KSS wavefunctions (Gamma-centered)
453 : : !! wfg(2,kss_npw*nspinor,nbandksseff)=The wavefunction Fourier coefficients.
454 : : !! iomode=Input variables specifying the fileformat. (0-->Fortran,3--> netcdf with ETSF-IO format).
455 : : !!
456 : : !! OUTPUT
457 : : !! Only writing.
458 : : !!
459 : : !! SOURCE
460 : :
461 : 31 : subroutine write_kss_wfgk(kss_unt,ikpt,isppol,kpoint,nspinor,kss_npw,&
462 : 31 : & nbandksseff,natom,Psps,ene_k,occ_k,rprimd,gbig,wfg,Cprjnk_k,iomode)
463 : :
464 : : !Arguments ------------------------------------
465 : : !scalars
466 : : integer,intent(in) :: ikpt,isppol,iomode,kss_npw,nspinor,kss_unt,nbandksseff
467 : : integer,intent(in) :: natom
468 : : type(pseudopotential_type),intent(in) :: Psps
469 : : !arrays
470 : : integer,intent(in) :: gbig(3,kss_npw)
471 : : real(dp),intent(in) :: kpoint(3),rprimd(3,3)
472 : : real(dp),intent(in) :: ene_k(nbandksseff),occ_k(nbandksseff)
473 : : real(dp),intent(in) :: wfg(2,kss_npw*nspinor,nbandksseff)
474 : : type(pawcprj_type),intent(in) :: Cprjnk_k(natom,nspinor*nbandksseff*Psps%usepaw)
475 : :
476 : : !Local variables-------------------------------
477 : : !scalars
478 : : integer :: ib,ibsp,ig,ispinor,iatom,ii !,ierr
479 : : integer :: kg_varid,cg_varid,ncerr
480 : : character(len=nctk_slen) :: kdep
481 : : ! *********************************************************************
482 : :
483 : : ! Calculate and write KB form factors and derivative at this k-point.
484 [ + + ]: 31 : if (Psps%usepaw==0) call write_vkb(kss_unt,ikpt,kpoint,kss_npw,gbig,rprimd,Psps,iomode)
485 : :
486 : : ! ============================================================
487 : : ! ==== Write wavefunctions and PAW matrix elements on disk ====
488 : : ! ============================================================
489 : 31 : SELECT CASE (iomode)
490 : :
491 : : CASE (IO_MODE_FORTRAN)
492 : 0 : write(kss_unt) (ene_k(ib),ib=1,nbandksseff)
493 : :
494 : 0 : ibsp=0
495 [ # # ]: 0 : do ib=1,nbandksseff
496 [ # # # # ]: 0 : write(kss_unt) (wfg(:,ig,ib),ig=1,kss_npw*nspinor)
497 [ # # ]: 0 : if (Psps%usepaw==1) then ! Remember that cprj are unsorted.
498 [ # # ]: 0 : do ispinor=1,nspinor
499 : 0 : ibsp=ibsp+1
500 [ # # ]: 0 : do iatom=1,natom
501 : 0 : ii=Cprjnk_k(iatom,ibsp)%nlmn
502 [ # # # # ]: 0 : write(kss_unt) (Cprjnk_k(iatom,ibsp)%cp(:,1:ii))
503 : : end do
504 : : end do
505 : : end if
506 : : end do
507 : :
508 : : CASE (IO_MODE_ETSF)
509 [ + + ]: 31 : if (Psps%usepaw==1) then
510 : 11 : ABI_WARNING("PAW output with ETSF-IO netcdf: cprj won't be written")
511 : : end if
512 : :
513 : : ! Write G-vectors (gbig because it's not k-dependent)
514 [ - + ]: 31 : NCF_CHECK(nf90_inq_varid(kss_unt, "reduced_coordinates_of_plane_waves", kg_varid))
515 [ - + ]: 31 : NCF_CHECK(nf90_get_att(kss_unt, kg_varid, "k_dependent", kdep))
516 [ + - ]: 31 : if (kdep == "no") then
517 [ + + ]: 93 : ncerr = nf90_put_var(kss_unt, kg_varid, gbig, start=[1,1], count=[3,kss_npw])
518 : : else
519 [ # # # # ]: 0 : ncerr = nf90_put_var(kss_unt, kg_varid, gbig, start=[1,1,ikpt], count=[3,kss_npw,1])
520 : : end if
521 [ - + ]: 31 : NCF_CHECK_MSG(ncerr, "putting gibg")
522 : :
523 : : ! Write wavefunctions
524 : : ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
525 [ - + ]: 31 : NCF_CHECK(nf90_inq_varid(kss_unt, "coefficients_of_wavefunctions", cg_varid))
526 : : ncerr = nf90_put_var(kss_unt, cg_varid, wfg, start=[1,1,1,1,ikpt,isppol], &
527 [ + + + + ]: 403 : count=[2,kss_npw,nspinor,nbandksseff,1,1])
528 [ - + ]: 31 : NCF_CHECK_MSG(ncerr, "putting cg_k")
529 : :
530 : : ! Write eigenvalues and occupations
531 [ + + - + : 124 : NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "eigenvalues"), ene_k, start=[1,ikpt,isppol]))
- - ]
532 [ + + - + : 124 : NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "occupations"), occ_k, start=[1,ikpt,isppol]))
- - ]
533 : :
534 : : CASE DEFAULT
535 [ - + - ]: 31 : ABI_ERROR(sjoin("Unsupported iomode:", itoa(iomode)))
536 : : END SELECT
537 : :
538 : 31 : end subroutine write_kss_wfgk
539 : : !!***
540 : :
541 : : !----------------------------------------------------------------------
542 : :
543 : : !!****f* m_io_kss/k2gamma_centered
544 : : !! NAME
545 : : !! k2gamma_centered
546 : : !!
547 : : !! FUNCTION
548 : : !! Helper function to translate a set of wavefunctions from the k-centered G-sphere
549 : : !! to the Gamma-centered G-sphere used for GW calculations.
550 : : !!
551 : : !! INPUTS
552 : : !! npw_k=Number of planewaves in the k-centered basis set.
553 : : !! kss_npw=Number of planewaves in the Gamma-centered G-sphere.
554 : : !! nspinor=Number of spinorial component.
555 : : !! nbandksseff=Number of bands in input-output arrays.
556 : : !! [icg]=Shift to be used when accessing the cg array. 0 if not specified (usually k_index).
557 : : !! [eig_vec(2,npw_k*nspinor,nbandksseff)]=wavefunctions defined on the k-centered G-sphere.
558 : : !! [cg(2,ikg+1:ikg+npw_k*nspinor*nbandksseff)]=wavefunctions defined on the k-centered G-sphere.
559 : : !! ngfft(18)=Info on the FFT.
560 : : !! MPI_enreg<MPI_type>=Structure gathering info on the parallelization.
561 : : !! istwf_k
562 : : !! ecut
563 : : !! gbig(3,kss_npw)
564 : : !! kg_k(3,npw_k)
565 : : !! gmet(3,3)
566 : : !! kpoint(3)
567 : : !!
568 : : !! OUTPUT
569 : : !! wfg(2,kss_npw*nspinor,nbandksseff)=Wavefunctions in the Gamma-centered representation.
570 : : !!
571 : : !! NOTES
572 : : !! 1) icg is used only if cg is present.
573 : : !! 2) cg and eig_vec are mutually exclusive. One and only one can be passed to the routine.
574 : : !!
575 : : !! SOURCE
576 : :
577 : 62 : subroutine k2gamma_centered(kpoint,npw_k,istwf_k,ecut,kg_k,kss_npw,nspinor,nbandksseff,ngfft,gmet,&
578 [ + + + - : 52 : & MPI_enreg,gbig,ug,icg,cg,eig_vec)
- + ]
579 : :
580 : : !Arguments ------------------------------------
581 : : !scalars
582 : : integer,intent(in) :: nbandksseff,nspinor,kss_npw,npw_k,istwf_k
583 : : integer,optional,intent(in) :: icg
584 : : real(dp),intent(in) :: ecut
585 : : type(MPI_type),intent(inout) :: MPI_enreg
586 : : !arrays
587 : : integer,intent(in) :: gbig(3,kss_npw)
588 : : integer,intent(in) :: kg_k(3,npw_k)
589 : : integer,intent(in) :: ngfft(18)
590 : : real(dp),intent(in) :: gmet(3,3),kpoint(3)
591 : : real(dp),intent(out) :: ug(2,kss_npw*nspinor,nbandksseff)
592 : : real(dp),optional,intent(in) :: eig_vec(2,npw_k*nspinor,nbandksseff)
593 : : real(dp),optional,intent(in) :: cg(:,:)
594 : :
595 : : !Local variables-------------------------------
596 : : !scalars
597 : : integer,parameter :: tobox=1,tosph=-1
598 : : integer :: band,ispinor,spinor_shift2,spinor_shift1,ig,my_icg,ierr
599 : : integer :: n1,n2,n3,n4,n5,n6,ndat,full_npw_k,ii
600 : : character(len=500) :: msg
601 : : !arrays
602 : : integer :: identity(3,3)=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/))
603 : : integer :: no_shift(3)=(/0,0,0/)
604 : 31 : integer,allocatable :: trsl(:),full_kg_k(:,:)
605 : 31 : real(dp),allocatable :: cfft(:,:,:,:)
606 : 31 : real(dp),allocatable :: full_cg(:,:),tmp_cg(:,:)
607 : :
608 : : ! *********************************************************************
609 : :
610 [ - + ]: 31 : if (PRESENT(cg).and.PRESENT(eig_vec)) then
611 : 0 : ABI_ERROR("Both cg and eig_vec are present!")
612 : : end if
613 : :
614 : : ! Mapping between the gamma-centered basis set and the k-centered one.
615 : : ! trsl(ig)=npw_k+1 if vector ig is not inside the k-centered G-sphere.
616 [ - + + - : 93 : ABI_MALLOC(trsl,(kss_npw))
+ - + - -
+ - + ]
617 : :
618 : 31 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
619 : 31 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
620 : :
621 [ + - ]: 31 : if (istwf_k==1) then ! Full k-centered G-sphere.
622 : 31 : call table_gbig2kg(npw_k,kg_k,kss_npw,gbig,trsl,ierr)
623 [ - + - - ]: 31 : if (ierr/=0.and.(kss_npw>=npw_k)) then
624 : 0 : ABI_ERROR(' The set of G vectors is inconsistent')
625 : : end if
626 : :
627 : : else ! Calculate full kg with istwf_k=1 then do the mapping.
628 : 0 : call get_kg(kpoint,1,ecut,gmet,full_npw_k,full_kg_k)
629 : :
630 : 0 : call table_gbig2kg(full_npw_k,full_kg_k,kss_npw,gbig,trsl,ierr)
631 [ # # # # ]: 0 : if (ierr/=0.and.(kss_npw>=npw_k)) then
632 : 0 : ABI_ERROR(' The set of G vectors is inconsistent')
633 : : end if
634 : : end if
635 : : !
636 : : ! Branching, depending on optional arguments.
637 [ + + ]: 31 : if (PRESENT(cg)) then
638 [ + - ]: 10 : my_icg=0; if (PRESENT(icg)) my_icg=icg
639 : :
640 : : SELECT CASE (istwf_k)
641 : :
642 : : CASE (1)
643 [ + + ]: 110 : do band=1,nbandksseff
644 [ + + ]: 210 : do ispinor=1,nspinor
645 : 100 : spinor_shift1=(ispinor-1)*kss_npw
646 : 100 : spinor_shift2=(ispinor-1)*npw_k
647 [ + + ]: 188700 : do ig=1,kss_npw ! Retrieve the correct components
648 [ + + ]: 188600 : if (trsl(ig)<=npw_k) then
649 [ + + ]: 414660 : ug(:,ig+spinor_shift1,band)=cg(:,trsl(ig)+spinor_shift2+(band-1)*npw_k*nspinor+my_icg)
650 : : else
651 [ + + ]: 150840 : ug(:,ig+spinor_shift1,band)=zero
652 : : end if
653 : : end do
654 : : end do
655 : : end do
656 : :
657 : : CASE (2:9)
658 : :
659 [ # # ]: 0 : ABI_CHECK(nspinor==1,"nspinor/=1!")
660 : : !
661 : : ! Convert input wfs from reduced to full G-sphere.
662 : 0 : ndat=1
663 [ # # # # : 0 : ABI_MALLOC(cfft,(2,n4,n5,n6*ndat))
# # # # #
# # # # #
# # # # #
# ]
664 [ # # # # : 0 : ABI_MALLOC(full_cg,(2,full_npw_k*ndat))
# # # # #
# # # ]
665 [ # # # # : 0 : ABI_MALLOC(tmp_cg,(2,npw_k*ndat))
# # # # #
# # # ]
666 : :
667 : : !write(std_out,*)"npw_k, full_kg_k",npw_k,full_npw_k
668 : :
669 [ # # ]: 0 : do band=1,nbandksseff
670 : 0 : ii = (band-1)*npw_k
671 [ # # # # : 0 : tmp_cg = cg(:,my_icg+ii+1:my_icg+ii+npw_k)
# # # # #
# # # ]
672 : : !write(776,*)"band= ",band,tmp_cg !cg(1:,my_icg+1+ii:my_icg+ii+npw_k)
673 : :
674 : 0 : call sphere(tmp_cg,ndat,npw_k,cfft,n1,n2,n3,n4,n5,n6,kg_k,istwf_k,tobox,MPI_enreg%me_g0,no_shift,identity,one)
675 : :
676 : 0 : call sphere(full_cg,ndat,full_npw_k,cfft,n1,n2,n3,n4,n5,n6,full_kg_k,1,tosph,MPI_enreg%me_g0,no_shift,identity,one)
677 : : !write(777,*)"band= ",band,full_cg(:,:)
678 : :
679 [ # # ]: 0 : do ig=1,kss_npw ! Retrieve the correct components
680 [ # # ]: 0 : if (trsl(ig)<=full_npw_k) then
681 [ # # ]: 0 : ug(:,ig,band)=full_cg(:,trsl(ig))
682 : : else
683 [ # # ]: 0 : ug(:,ig,band)=zero
684 : : end if
685 : : end do
686 : : end do !band
687 : :
688 [ # # ]: 0 : ABI_FREE(cfft)
689 [ # # ]: 0 : ABI_FREE(tmp_cg)
690 [ # # ]: 0 : ABI_FREE(full_cg)
691 : :
692 : : CASE DEFAULT
693 [ + - - ]: 10 : ABI_BUG("Wrong istwf_k")
694 : : END SELECT
695 : :
696 [ + - ]: 21 : else if (PRESENT(eig_vec)) then
697 : :
698 : : SELECT CASE (istwf_k)
699 : :
700 : : CASE (1)
701 [ + + ]: 1241 : do band=1,nbandksseff
702 [ + + ]: 2461 : do ispinor=1,nspinor
703 : 1220 : spinor_shift1=(ispinor-1)*kss_npw
704 : 1220 : spinor_shift2=(ispinor-1)*npw_k
705 [ + + ]: 623380 : do ig=1,kss_npw ! Retrieve the correct components
706 [ + + ]: 622160 : if (trsl(ig)<=npw_k) then
707 [ + + ]: 1176480 : ug(:,ig+spinor_shift1,band)=eig_vec(:,trsl(ig)+spinor_shift2,band)
708 : : else
709 [ + + ]: 686340 : ug(:,ig+spinor_shift1,band)=zero
710 : : end if
711 : : end do
712 : : end do
713 : : end do
714 : :
715 : : CASE DEFAULT
716 : 0 : write(msg,'(a,i0)')" Unsupported value for istwf_k: ",istwf_k
717 [ + - ]: 21 : ABI_ERROR(msg)
718 : : END SELECT
719 : :
720 : : else
721 : 0 : ABI_ERROR("neither cg not eig_vec are in input")
722 : : end if
723 : :
724 [ - + ]: 31 : ABI_FREE(trsl)
725 [ - + ]: 31 : if (allocated(full_kg_k)) then
726 : 0 : ABI_FREE(full_kg_k)
727 : : end if
728 : :
729 [ + + - + ]: 62 : end subroutine k2gamma_centered
730 : : !!***
731 : :
732 : : !----------------------------------------------------------------------
733 : :
734 : : !!****f* m_io_kss/make_gvec_kss
735 : : !! NAME
736 : : !! make_gvec_kss
737 : : !!
738 : : !! FUNCTION
739 : : !! Build the list of G-vectors using the KSS convention.
740 : : !!
741 : : !! INPUTS
742 : : !! nkpt=Number of k-points.
743 : : !! nsym=Number of symmetries.
744 : : !! prtvol=Verbosity option.
745 : : !! symmorphi=
746 : : !! 0 : Old (Obsolete) implementation => Suppress inversion from symmetries list
747 : : !! 1 : Use input symrel, tnons.
748 : : !! ecut_eff=Effective cutoff
749 : : !! symrel(3,3,nsym)= Symmetry operation in real space.
750 : : !! tnons(3,nsym)=Fractional translations
751 : : !! kptns(3,nkpt)=K-points in reduced coordinates.
752 : : !!
753 : : !! OUTPUT
754 : : !! npwkss = Input: Initial guess for the number of G-vectors required. Use 0 to have the
755 : : !! full list of G-vectors that form a closed shell.
756 : : !! Output: Actual number of G-vectors that form a set of closed shells
757 : : !! gvec_kss(:,:) = Input: null pointer. Output: gvec_kss(3,npwkss), list of G-vectors (closed shells)
758 : : !! ierr=Status error
759 : : !!
760 : : !! SOURCE
761 : :
762 : 364 : subroutine make_gvec_kss(nkpt,kptns,ecut_eff,symmorphi,nsym,symrel,tnons,gprimd,prtvol,npwkss,gvec_kss,ierr)
763 : :
764 : : !Arguments ------------------------------------
765 : : !scalars
766 : : integer,intent(in) :: nkpt,nsym,prtvol,symmorphi
767 : : integer,intent(out) :: ierr
768 : : integer,intent(inout) :: npwkss
769 : : real(dp),intent(in) :: ecut_eff
770 : : !arrays
771 : : integer,intent(in) :: symrel(3,3,nsym)
772 : : integer,pointer :: gvec_kss(:,:)
773 : : real(dp),intent(in) :: tnons(3,nsym),kptns(3,nkpt)
774 : : real(dp),intent(in) :: gprimd(3,3)
775 : :
776 : : !Local variables-------------------------------
777 : : !scalars
778 : : integer :: ii,ishm,maxpw,nbase
779 : : integer :: nrst1,nrst2,nsym2,pinv
780 : 364 : integer,pointer :: gbig(:,:)
781 : : character(len=500) :: msg
782 : : !arrays
783 : 364 : integer,pointer :: symrel2(:,:,:),shlim(:)
784 : 364 : real(dp),pointer :: tnons2(:,:)
785 : : ! *********************************************************************
786 : :
787 : 364 : ierr = 0
788 : 364 : write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.'
789 : 364 : call wrtout(std_out,msg,'COLL')
790 : :
791 : : !ecut_eff = ecut * Dtset%dilatmx**2 ! Use ecut_eff instead of ecut_eff since otherwise
792 : : !
793 : : !============================================================
794 : : !=== Prepare set containing all G-vectors sorted by stars ===
795 : : !============================================================
796 : : !
797 : : !=== Analyze symmetry operations ===
798 [ + + ]: 364 : if (symmorphi==0) then ! Old (Obsolete) implementation: Suppress inversion from symmetries list:
799 : 15 : nullify(symrel2,tnons2)
800 : 15 : call remove_inversion(nsym,symrel,tnons,nsym2,symrel2,tnons2,pinv)
801 [ + + + + : 1279 : if (ANY(ABS(tnons2(:,1:nsym2))>tol8)) then
+ - - + ]
802 : : write(msg,'(3a)')&
803 : 0 : & ' Non-symmorphic operations still remain in the symmetries list ',ch10,&
804 : 0 : & ' Program does not stop but _KSS file will not be created...'
805 : 0 : ABI_WARNING(msg)
806 : 0 : ierr=ierr+1 ; RETURN
807 : : end if
808 [ + - ]: 349 : else if (symmorphi==1) then
809 : : ! If in the input file symmorphi==1 all the symmetry operations are retained:
810 : : ! both identity and inversion (if any) as well as non-symmorphic operations.
811 : 349 : nsym2=nsym ; pinv=1
812 [ - + + - : 1047 : ABI_MALLOC(symrel2,(3,3,nsym))
+ - + - -
+ - + ]
813 [ - + + - : 1047 : ABI_MALLOC(tnons2,(3,nsym))
+ - + - -
+ - + ]
814 [ + + + + : 187289 : symrel2(:,:,:)=symrel(:,:,1:nsym)
+ + ]
815 [ + + + + ]: 57869 : tnons2(:,:) =tnons(:,1:nsym)
816 : : else
817 : : write(msg,'(a,i4,3a)')&
818 : 0 : & ' symmorphi = ',symmorphi,' while it must be 0 or 1',ch10,&
819 : 0 : & ' Program does not stop but KSS file will not be created...'
820 : 0 : ABI_WARNING(msg)
821 : 0 : ierr=ierr+1 ; RETURN
822 : : end if
823 : : !
824 : : !===================================================================
825 : : !==== Merge the set of k-centered G-spheres into a big set gbig ====
826 : : !===================================================================
827 : : !* Vectors in gbig are ordered by shells
828 : : !
829 : 364 : nullify(gbig,shlim)
830 [ - + ]: 364 : call merge_and_sort_kg(nkpt,kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim)
831 : :
832 : 364 : nbase = SIZE(shlim) ! Number of independent G in the big sphere.
833 : 364 : maxpw = shlim(nbase) ! Total number of G"s in the big sphere.
834 : : !
835 : : ! * Determine optimal number of bands and G"s to be written.
836 : : !npwkss=Dtset%npwkss
837 [ + + + - ]: 364 : if ((npwkss==0).or.(npwkss>=maxpw)) then
838 : 364 : npwkss=maxpw
839 : : write(msg,'(5a)')&
840 : 364 : & ' Since the number of g''s to be written on file',ch10,&
841 : 364 : & ' was 0 or too large, it has been set to the max. value.,',ch10,&
842 : 728 : & ' computed from the union of the sets of G vectors for the different k-points.'
843 : 364 : call wrtout(std_out,msg,'COLL')
844 : : end if
845 : :
846 : 364 : ishm=0
847 [ + + ]: 20863 : do ii=1,nbase
848 [ + - ]: 20863 : if (shlim(ii)<=npwkss) then
849 : 20499 : ishm=ii
850 : : else
851 : : EXIT
852 : : end if
853 : : end do
854 : : !ishm=bisect(shlim,npwkss)
855 : :
856 [ - + ]: 364 : if (shlim(ishm)/=npwkss) then
857 : 0 : nrst1=shlim(ishm)
858 : 0 : nrst2=MIN0(shlim(MIN0(ishm+1,nbase)),maxpw)
859 [ # # ]: 0 : if (IABS(npwkss-nrst2)<IABS(npwkss-nrst1)) nrst1=nrst2
860 : 0 : npwkss=nrst1
861 : 0 : if (shlim(ishm)<npwkss) ishm=ishm+1
862 : : write(msg,'(3a)')&
863 : 0 : & ' The number of G''s to be written on file is not a whole number of stars ',ch10,&
864 : 0 : & ' the program set it to the nearest star limit.'
865 : 0 : call wrtout(std_out,msg,'COLL')
866 : : end if
867 : :
868 : 364 : write(msg,'(a,i5)')' Number of G-vectors is: ',npwkss
869 : 364 : call wrtout(std_out,msg,'COLL')
870 : :
871 [ - + + - : 1092 : ABI_MALLOC(gvec_kss,(3,npwkss))
+ - + - -
+ - + ]
872 [ + - + + : 2085804 : gvec_kss = gbig(:,1:npwkss)
+ + + + +
+ ]
873 : :
874 [ - + ]: 364 : ABI_FREE(gbig)
875 [ - + ]: 364 : ABI_FREE(symrel2)
876 [ - + ]: 364 : ABI_FREE(tnons2)
877 [ - + ]: 364 : ABI_FREE(shlim)
878 : :
879 : 364 : end subroutine make_gvec_kss
880 : : !!***
881 : :
882 : : !----------------------------------------------------------------------
883 : :
884 : : !!****f* m_io_kss/kss_calc_vkb
885 : : !! NAME
886 : : !! kss_calc_vkb
887 : : !!
888 : : !! FUNCTION
889 : : !! This routine calculates the Kleynman-Bylander form factors and its derivatives
890 : : !! needed for the evaluation of the matrix elements of the dipole operator <phi1|r|phi2>.
891 : : !!
892 : : !! INPUTS
893 : : !! npw_k=Number of plane waves for this k-point.
894 : : !! Psps<pseudopotential_type>=Structured datatype gathering information on the pseudopotentials.
895 : : !! kg_k(3,npw_k)=Reduced coordinates of the G-vectors.
896 : : !! kpoint(3)=The k-point in reduced coordinates.
897 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
898 : : !!
899 : : !! OUTPUT
900 : : !! vkb (npw_k,Psps%ntypat,Psps%mpsang)=KB form factors.
901 : : !! vkbd(npw_k,Psps%ntypat,Psps%mpsang)=KB form factor derivatives.
902 : : !! vkbsign(Psps%mpsang,Psps%ntypat) =KS dyadic sign.
903 : : !!
904 : : !! NOTES
905 : : !! This piece of code has been extracted from outkss.F90. The implementation is consistent
906 : : !! with the KSS file formata (Fortran version) but it presents two design flaws.
907 : : !!
908 : : !! 1) Pseudo with more that one projector per l-channel are not supported.
909 : : !! 2) Ordering of dimensions in vkb and vkbd is not optimal. We are not programming C!!!
910 : : !!
911 : : !! TODO
912 : : !! *) Spinorial case is not implemented.
913 : : !!
914 : : !! SOURCE
915 : :
916 : 20 : subroutine kss_calc_vkb(Psps,kpoint,npw_k,kg_k,rprimd,vkbsign,vkb,vkbd)
917 : :
918 : : !Arguments ------------------------------------
919 : : !scalars
920 : : integer,intent(in) :: npw_k
921 : : type(Pseudopotential_type),intent(in) :: Psps
922 : : !arrays
923 : : integer,intent(in) :: kg_k(3,npw_k)
924 : : real(dp),intent(in) :: kpoint(3),rprimd(3,3)
925 : : real(dp),intent(out) :: vkb (npw_k,Psps%ntypat,Psps%mpsang)
926 : : real(dp),intent(out) :: vkbd(npw_k,Psps%ntypat,Psps%mpsang)
927 : : real(dp),intent(out) :: vkbsign(Psps%mpsang,Psps%ntypat)
928 : :
929 : : !Local variables ------------------------------
930 : : !scalars
931 : : integer :: dimffnl,ider,idir,itypat,nkpg,il0,in
932 : : integer :: il,ilmn,ig,is
933 : : real(dp) :: ucvol,effmass_free,ecutsm,ecut
934 : : !arrays
935 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
936 : 20 : real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:), ylm(:,:),ylm_gr(:,:,:),ylm_k(:,:)
937 : : ! *************************************************************************
938 : :
939 : : DBG_ENTER("COLL")
940 : :
941 [ - + ]: 20 : ABI_CHECK(Psps%usepaw==0,"You should not be here!")
942 [ - + ]: 20 : ABI_CHECK(Psps%useylm==0,"useylm/=0 not considered!")
943 : :
944 : 20 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
945 : : !
946 : : ! === Save KB dyadic sign (integer-valued) ===
947 [ + + + + ]: 120 : vkbsign=zero
948 [ + + ]: 50 : do itypat=1,Psps%ntypat
949 : 30 : il0=0
950 [ + + ]: 90 : do ilmn=1,Psps%lmnmax
951 : 40 : il=1+Psps%indlmn(1,ilmn,itypat)
952 : 40 : in=Psps%indlmn(3,ilmn,itypat)
953 [ + - ]: 70 : if (il/=il0 .and. in==1) then
954 : 40 : il0=il
955 : 40 : vkbsign(il,itypat)=DSIGN(one,Psps%ekb(ilmn,itypat))
956 : : end if
957 : : end do
958 : : end do
959 : :
960 : : ! === Allocate KB form factor and derivative wrt k+G ===
961 : : ! * Here we do not use correct ordering for dimensions
962 : :
963 : 20 : ider=1; dimffnl=2 ! To retrieve the first derivative.
964 : 20 : idir=0; nkpg=0
965 : : !
966 : : ! Quantities used only if useylm==1
967 [ - + + - : 80 : ABI_MALLOC(ylm,(npw_k,Psps%mpsang**2*Psps%useylm))
+ - - - +
- - + - +
- + ]
968 [ - + + - : 100 : ABI_MALLOC(ylm_gr,(npw_k,3+6*(ider/2),Psps%mpsang**2*Psps%useylm))
+ - + - -
- + - - +
- + - + ]
969 [ - + + - : 60 : ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm))
+ - - - -
+ - + ]
970 [ - + + - : 40 : ABI_MALLOC(kpg_dum,(npw_k,nkpg))
- + - + ]
971 : :
972 [ - + + - : 120 : ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat))
+ - - + +
- - + + -
+ - + - -
+ - + ]
973 : :
974 : : call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,&
975 : : kg_k,kpg_dum,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,&
976 : 20 : Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)
977 : :
978 [ - + ]: 20 : ABI_FREE(kpg_dum)
979 [ - + ]: 20 : ABI_FREE(ylm)
980 [ - + ]: 20 : ABI_FREE(ylm_gr)
981 [ - + ]: 20 : ABI_FREE(ylm_k)
982 : :
983 [ - + + - : 60 : ABI_MALLOC(modkplusg,(npw_k))
+ - + - -
+ - + ]
984 : :
985 : 20 : effmass_free=one; ecutsm=zero; ecut=HUGE(one)
986 : : ! call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,modkplusg,kpoint,npw_k)
987 : 20 : call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,modkplusg,kpoint,npw_k,0,0)
988 [ + + ]: 23700 : modkplusg(:)=SQRT(half/pi**2*modkplusg(:))
989 [ + + ]: 23700 : modkplusg(:)=MAX(modkplusg(:),tol10)
990 : :
991 : : !do ig=1,npw_k
992 : : ! kpg(:)= kpoint(:)+kg_k(:,ig)
993 : : ! modkplusg(ig) = normv(kpg,gmet,"G")
994 : : !end do
995 : :
996 : : ! Calculate matrix elements.
997 [ + + + + : 180040 : vkb=zero; vkbd=zero
+ + + + +
+ + + ]
998 : :
999 [ + + ]: 50 : do is=1,Psps%ntypat
1000 : : il0=0
1001 [ + + ]: 90 : do ilmn=1,Psps%lmnmax
1002 : 40 : il=1+Psps%indlmn(1,ilmn,is)
1003 : 40 : in=Psps%indlmn(3,ilmn,is)
1004 [ + - ]: 70 : if ((il/=il0).and.(in==1)) then
1005 : 40 : il0=il
1006 [ + - ]: 40 : if (ABS(Psps%ekb(ilmn,is))>1.0d-10) then
1007 [ + + ]: 40 : if (il==1) then
1008 [ + + ]: 42560 : vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)
1009 [ + + ]: 42560 : vkbd(1:npw_k,is,il) = ffnl(:,2,ilmn,is)*modkplusg(:)/two_pi
1010 [ + - ]: 10 : else if (il==2) then
1011 [ + + ]: 4840 : vkb(1:npw_k,is,il) = ffnl(:,1,ilmn,is)*modkplusg(:)
1012 [ + + ]: 4840 : do ig=1,npw_k
1013 : : vkbd(ig,is,il) = ((ffnl(ig,2,ilmn,is)*modkplusg(ig)*modkplusg(ig))+&
1014 : 4840 : ffnl(ig,1,ilmn,is) )/two_pi
1015 : : end do
1016 [ # # ]: 0 : else if (il==3) then
1017 [ # # ]: 0 : vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)*modkplusg(:)**2
1018 : : vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**3+&
1019 [ # # ]: 0 : 2*ffnl(:,1,ilmn,is)*modkplusg(:) )/two_pi
1020 [ # # ]: 0 : else if (il==4) then
1021 [ # # ]: 0 : vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)*modkplusg(:)**3
1022 : : vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**4+&
1023 [ # # ]: 0 : 3*ffnl(:,1,ilmn,is)*modkplusg(:)**2 )/two_pi
1024 : : end if
1025 [ + + ]: 47400 : vkb (:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkb (:,is,il)
1026 [ + + ]: 47400 : vkbd(:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkbd(:,is,il)
1027 : : else
1028 [ # # ]: 0 : vkb (:,is,il)=zero
1029 [ # # ]: 0 : vkbd(:,is,il)=zero
1030 : : end if
1031 : : end if
1032 : : end do
1033 : : end do
1034 : :
1035 [ - + ]: 20 : ABI_FREE(ffnl)
1036 [ - + ]: 20 : ABI_FREE(modkplusg)
1037 : :
1038 : : DBG_EXIT("COLL")
1039 : :
1040 : 20 : end subroutine kss_calc_vkb
1041 : : !!***
1042 : :
1043 : : !!****f* m_io_kss/outkss
1044 : : !! NAME
1045 : : !! outkss
1046 : : !!
1047 : : !! FUNCTION
1048 : : !! This routine creates an output file containing the Kohn-Sham electronic Structure
1049 : : !! for a large number of eigenstates (energies and eigen-functions).
1050 : : !! The resulting file (_KSS) is needed for a GW post-treatment.
1051 : : !!
1052 : : !! The routine drives the following operations:
1053 : : !! - Re-ordering G-vectors according to stars (sets of Gs related by symmetry operations).
1054 : : !! A set of g for all k-points is created.
1055 : : !! - Creating and opening the output "_KSS'" file
1056 : : !! - Printing out output file header information...
1057 : : !! ... and, for each k-point:
1058 : : !! According to 'kssform', either
1059 : : !! - Re-computing <G|H|G_prim> matrix elements for all (G, G_prim).
1060 : : !! Diagonalizing H in the plane-wave basis.
1061 : : !! or - Taking eigenvalues/vectors from congugate-gradient ones.
1062 : : !! - Writing out eigenvalues and eigenvectors.
1063 : : !!
1064 : : !! INPUTS
1065 : : !! cg(2,mcg)=planewave coefficients of wavefunctions.
1066 : : !! usecprj=1 if cprj datastructure has been allocated (ONLY PAW)
1067 : : !! Cprj(natom,mcprj*usecprj) <type(pawcprj_type)>=
1068 : : !! projected input wave functions <Proj_i|Cnk> with all NL projectors (only for PAW)
1069 : : !! NOTE that Cprj are unsorted, see ctoprj.F90
1070 : : !! Dtfil <type(datafiles_type)>=variables related to files
1071 : : !! Dtset <type(dataset_type)>=all input variables for this dataset
1072 : : !! ecut=cut-off energy for plane wave basis sphere (Ha)
1073 : : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
1074 : : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
1075 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1076 : : !! Hdr <type(hdr_type)>=the header of wf, den and pot files
1077 : : !! kssform=govern the Kohn-Sham Structure file format
1078 : : !! mband=maximum number of bands
1079 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
1080 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
1081 : : !! mgfft=maximum size of 1D FFTs
1082 : : !! mkmem =number of k points treated by this node.
1083 : : !! MPI_enreg=information about MPI parallelization
1084 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
1085 : : !! mpw=maximum dimensioned size of npw.
1086 : : !! my_natom=number of atoms treated by current processor
1087 : : !! natom=number of atoms in cell.
1088 : : !! nfft=(effective) number of FFT grid points (for this processor)
1089 : : !! nkpt=number of k points.
1090 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
1091 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
1092 : : !! nspden=number of density components
1093 : : !! nsym=number of symmetries in space group
1094 : : !! ntypat=number of types of atoms in unit cell.
1095 : : !! occ(mband*nkpt*nsppol)=occupation number for each band (usually 2) for each k.
1096 : : !! Pawtab(Psps%ntypat*Psps%usepaw) <type(pawtab_type)>=paw tabulated starting data
1097 : : !! Pawfgr<pawfgr_type>=fine grid parameters and related data
1098 : : !! prtvol=control print volume and debugging output
1099 : : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
1100 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1101 : : !! vtrial(nfft,nspden)=the trial potential
1102 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
1103 : : !!
1104 : : !! OUTPUT
1105 : : !! Output is written on file.
1106 : : !! ierr=Status error.
1107 : : !!
1108 : : !! NOTES
1109 : : !!
1110 : : !! * This routine is maintained for legacy reasons. Abinit8 is not able to read KSS files
1111 : : !! anymore hence KSS files should be used only to interface Abinit with external codes
1112 : : !! that are still using the old KSS format.
1113 : : !!
1114 : : !! * The routine can be time consuming (in particular when computing
1115 : : !! <G|H|G_prim> elements for all (G, G_prim)) (kssform=1).
1116 : : !! So, it is recommended to call it once per run...
1117 : : !!
1118 : : !! * The IO code is not parallelized and this represents a serious bottleneck when np is large.
1119 : : !!
1120 : : !! * when kssform==1, the routine RE-computes all Hamiltonian terms.
1121 : : !! So it is equivalent to an additional electronic SC cycle.
1122 : : !! (This has no effect is convergence was reach...
1123 : : !! If not, eigenvalues/vectors may differs from the congugaste gradient ones)
1124 : : !!
1125 : : !! * The KB form factors and derivatives are not calculated correctly if there are
1126 : : !! pseudos with more than one projector in an angular momentum channel.
1127 : : !!
1128 : : !! * In the ETSF output format (Dtset%iomode == 3), the complete symmetry set
1129 : : !! is output. So, if reading programs need only the symmorphic symmetries, they
1130 : : !! will need to remove themselves the non-symmorphic ones.
1131 : : !!
1132 : : !! * There exists two file formats:
1133 : : !! kssform==1 diagonalized file _KSS in real(dp) is generated.
1134 : : !! kssform==3 same as kssform=1 but the wavefunctions are not diagonalized
1135 : : !! (they are taken from conjugate-gradient ones)
1136 : : !! Old kssform=0 and kssform=2 are obsolete and no longer available
1137 : : !!
1138 : : !! TESTS
1139 : : !! * ETSF_IO output is tested in tests/etsf_io/t02.
1140 : : !!
1141 : : !! SOURCE
1142 : :
1143 : 6 : subroutine outkss(crystal,Dtfil,Dtset,ecut,gmet,gprimd,Hdr,&
1144 : : & kssform,mband,mcg,mcprj,mgfft,mkmem,MPI_enreg,mpsang,mpw,my_natom,natom,&
1145 : 6 : & nfft,nkpt,npwarr,nspden,nsppol,nsym,ntypat,occ,Pawtab,Pawfgr,Paw_ij,&
1146 : 6 : & prtvol,Psps,rprimd,vtrial,xred,cg,usecprj,Cprj,eigen,ierr)
1147 : :
1148 : : use m_linalg_interfaces
1149 : :
1150 : : !Arguments ------------------------------------
1151 : : !scalars
1152 : : integer,intent(in) :: kssform,mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,usecprj
1153 : : integer,intent(in) :: nfft,nkpt,nsppol,nspden,nsym,ntypat,prtvol
1154 : : integer,intent(out) :: ierr
1155 : : real(dp),intent(in) :: ecut
1156 : : type(MPI_type),intent(inout) :: MPI_enreg
1157 : : type(Datafiles_type),intent(in) :: Dtfil
1158 : : type(Dataset_type),intent(in) :: Dtset
1159 : : type(Hdr_type),intent(inout) :: Hdr
1160 : : type(Pseudopotential_type),intent(in) :: Psps
1161 : : type(pawfgr_type), intent(in) :: Pawfgr
1162 : : type(crystal_t),intent(in) :: crystal
1163 : : !arrays
1164 : : integer,intent(in),target :: npwarr(nkpt)
1165 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol)
1166 : : real(dp),intent(in) :: rprimd(3,3)
1167 : : real(dp),intent(inout) :: vtrial(nfft,nspden)
1168 : : real(dp),intent(in) :: xred(3,natom)
1169 : : real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol)
1170 : : type(pawcprj_type),intent(in) :: Cprj(natom,mcprj*usecprj)
1171 : : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
1172 : : type(paw_ij_type),intent(inout),target :: Paw_ij(my_natom*Psps%usepaw)
1173 : :
1174 : : !Local variables-------------------------------
1175 : : !scalars
1176 : : integer,parameter :: tim_rwwf=0,bufnb=20
1177 : : integer :: untkss,onband_diago
1178 : : integer :: bdtot_index,i,ib,ibp,iomode
1179 : : integer :: ibsp,ibsp1,ibsp2,ibg,ig,ii,ikpt
1180 : : integer :: master,receiver,sender,spinor_shift1,shift
1181 : : integer :: ishm,ispinor,isppol,istwf_k,my_rank,j
1182 : : integer :: k_index,maxpw,mproj,mtag,n1,n2,n2dim,n3,n4,n5,n6,nband_k
1183 : : integer :: nbandkss_k,nbandksseff,nbase,nprocs,npw_k,onpw_k,npwkss
1184 : : integer :: nrst1,nrst2,nsym2,ntemp,pinv,sizepw,spaceComm,comm_self
1185 : : integer :: pad1,pad2
1186 : : integer :: bufrt,bufsz
1187 : : real(dp) :: cinf=1.0e24,csup=zero,einf=1.0e24,esup=zero
1188 : : real(dp) :: norm,cfact,ecut_eff
1189 : : logical :: do_diago,found,ltest,lhack
1190 : : logical,parameter :: skip_test_ortho=.FALSE.
1191 : : character(len=500) :: msg
1192 : : character(len=80) :: frmt1,frmt2
1193 : : character(len=10) :: stag(2)=(/' ',' '/)
1194 : : !arrays
1195 : 12 : integer :: nbandkssk(nkpt)
1196 : 6 : integer,pointer :: symrel2(:,:,:)
1197 : 6 : integer,pointer :: gbig(:,:)
1198 : 6 : integer,pointer :: shlim(:)
1199 : 6 : integer,allocatable :: kg_k(:,:)
1200 : 6 : integer,allocatable :: dimlmn(:)
1201 : : integer :: nattyp_dum(0)
1202 : : real(dp) :: ovlp(2),kpoint(3),tsec(2)
1203 : 6 : real(dp),pointer :: tnons2(:,:)
1204 : 6 : real(dp),allocatable :: ene(:)
1205 : 6 : real(dp),pointer :: eig_ene(:),eig_vec(:,:,:)
1206 : 6 : real(dp),allocatable :: occ_k(:)
1207 : 6 : real(dp),allocatable,target :: wfg(:,:,:)
1208 : 6 : real(dp),contiguous, pointer :: ug1(:,:),ug2(:,:)
1209 : 6 : type(pawcprj_type),allocatable :: Cprjnk_k(:,:)
1210 : 6 : type(pawcprj_type),pointer :: Cprj_diago_k(:,:)
1211 : : type(ddiago_ctl_type) :: Diago_ctl
1212 : 6 : type(paw_ij_type),pointer :: Paw_ij_all(:)
1213 : : ! *********************************************************************
1214 : :
1215 : : ABI_UNUSED(mkmem)
1216 : :
1217 : : DBG_ENTER("COLL")
1218 : :
1219 : 6 : call timab(933,1,tsec) ! outkss
1220 : 6 : call timab(934,1,tsec) ! outkss(Gsort+hd)
1221 : :
1222 : 6 : spaceComm=MPI_enreg%comm_cell
1223 : 6 : my_rank=xmpi_comm_rank(spaceComm)
1224 : 6 : nprocs=xmpi_comm_size(spaceComm)
1225 : 6 : master=0
1226 : :
1227 : 6 : iomode = Dtset%iomode
1228 : 6 : nullify(eig_ene)
1229 : 6 : nullify(eig_vec)
1230 : 6 : nullify(Cprj_diago_k)
1231 : :
1232 : : ! JB: Valgrind complains about non initialized value. Set to -1 so if an array
1233 : : ! should be allocated with this "unintialized value" it crashes
1234 : 6 : onband_diago = -1
1235 : :
1236 : : !MG: since in seq case MPI_enreg%proc_distrb is not defined
1237 : : !we hack a bit the data type in order to get rid of MPI preprocessing options.
1238 : : !The previous status of %proc_distrb is restored before exiting.
1239 : : !Note that in case of seq run MPI_enreg%proc_distrb is nullified at the very beginning of abinit.F90
1240 : : !
1241 : : !FIXME this is a design flaw that should be solved: proc_distrb should always
1242 : : !be allocated and filled with my_rank in case of sequential run otherwise checks like
1243 : : !if (nprocs>1.and.MPI_enreg%proc_distrb(ii)==me) leads to SIGFAULT under gfortran.
1244 : : !as the second array is not allocated.
1245 : 6 : lhack=.FALSE.
1246 [ + + ]: 6 : if (nprocs==1) then
1247 : 2 : ltest=allocated(MPI_enreg%proc_distrb)
1248 [ - + ]: 2 : if (.not.ltest) then
1249 [ # # # # : 0 : ABI_MALLOC(MPI_enreg%proc_distrb,(nkpt,mband,nsppol))
# # # # #
# # # # #
# # # # #
# # # ]
1250 [ # # # # : 0 : MPI_enreg%proc_distrb=my_rank
# # ]
1251 : : lhack=.TRUE.
1252 : : end if
1253 [ + + + + : 213 : ltest=ALL(MPI_enreg%proc_distrb==my_rank)
+ + + - ]
1254 [ - + ]: 2 : ABI_CHECK(ltest,'wrong values in %proc_distrb')
1255 : : end if
1256 : : !
1257 : : !============================
1258 : : !==== Perform some tests ====
1259 : : !============================
1260 : 6 : ierr=0
1261 : :
1262 [ + - ]: 6 : if (iomode==IO_MODE_ETSF) then
1263 : : write(msg,'(3a)')&
1264 : 6 : & 'when iomode==3 in outkss, support for netcdf ',ch10,&
1265 : 12 : & 'must be compiled. Use --enable-netcdf when configuring '
1266 : : end if
1267 : :
1268 [ + + ]: 6 : if (kssform==3) then
1269 : 5 : write(msg,'(a,70("="),4a)')ch10,ch10,&
1270 : 5 : & ' Calculating and writing out Kohn-Sham electronic Structure file',ch10, &
1271 : 10 : & ' Using conjugate gradient wavefunctions and energies (kssform=3)'
1272 [ + - ]: 1 : else if (kssform==1) then
1273 : 1 : write(msg,'(a,70("="),4a,i1,a)') ch10,ch10, &
1274 : 1 : & ' Calculating and writing out Kohn-Sham electronic Structure file',ch10, &
1275 : 2 : & ' Using diagonalized wavefunctions and energies (kssform=',kssform,')'
1276 : : else
1277 : : write(msg,'(a,i0,2a)')&
1278 : 0 : & " Unsupported value for kssform: ",kssform,ch10,&
1279 : 0 : & " Program does not stop but _KSS file will not be created..."
1280 : 0 : ierr=ierr+1
1281 : : end if
1282 : 6 : call wrtout(std_out,msg,'COLL')
1283 : 6 : call wrtout(ab_out,msg,'COLL')
1284 : : !
1285 : : !* Check whether nband is constant in metals
1286 [ + + + - : 70 : if ( (Dtset%occopt>=2.and.Dtset%occopt<=8) .and. (ANY(Dtset%nband(1:nkpt*nsppol)/=Dtset%nband(1))) ) then
+ + - + ]
1287 : : write(msg,'(3a,i4,a,i3,a,i4,3a)')&
1288 : 0 : & ' The number of bands must be the same for all k-points ',ch10,&
1289 : 0 : & ' but nband(1)=',Dtset%nband(1),' is different of nband(',&
1290 : 0 : & ikpt+(isppol-1)*nkpt,')=',Dtset%nband(ikpt+(isppol-1)*nkpt),'.',ch10,&
1291 : 0 : & ' Program does not stop but _KSS file will not be created...'
1292 : 0 : ABI_WARNING(msg)
1293 : 0 : ierr=ierr+1
1294 : : end if
1295 : : !* istwfk must be 1 for each k-point
1296 [ + + + - : 70 : if (ANY(Dtset%istwfk(1:nkpt)/=1).and.kssform/=3) then
- + - - ]
1297 : : write(msg,'(7a)')&
1298 : 0 : & ' istwfk/=1 not allowed when kssform/=3 :',ch10,&
1299 : 0 : & ' States output not programmed for time-reversal symmetry.',ch10,&
1300 : 0 : & ' Action : change istwfk in input file (put it to 1 for all kpt).',ch10,&
1301 : 0 : & ' Program does not stop but _KSS file will not be created...'
1302 : 0 : ABI_WARNING(msg)
1303 : 0 : ierr=ierr+1
1304 : : end if
1305 : : !* Check spin-orbit
1306 [ - + ]: 6 : if (Psps%mpssoang/=mpsang) then
1307 : : write(msg,'(3a)')&
1308 : 0 : & ' Variable mpspso should be 1 !',ch10,&
1309 : 0 : & ' Program does not stop but _KSS file will not be created...'
1310 : 0 : ABI_WARNING(msg)
1311 : 0 : ierr=ierr+1
1312 : : end if
1313 : : !* Check mproj
1314 [ + + + + ]: 89 : mproj=MAXVAL(Psps%indlmn(3,:,:))
1315 [ + + - + ]: 6 : if (mproj>1.and.Psps%usepaw==0) then ! TODO One has to derive the expression for [Vnl,r], in particular HGH and GTH psps
1316 : 0 : write(msg,'(8a)')ch10,&
1317 : 0 : & ' outkss : COMMENT - ',ch10,&
1318 : 0 : & ' At least one NC pseudopotential has more that one projector per angular channel',ch10,&
1319 : 0 : & ' Note that inclvkb==0 should be used in screening, since the evaluation of the commutator',ch10,&
1320 : 0 : & ' for this particular case is not implemented yet'
1321 : 0 : call wrtout(std_out,msg,'COLL')
1322 : 0 : call wrtout(ab_out,msg,'COLL')
1323 : : end if
1324 : : !* Check max angular momentum
1325 [ + + + + : 89 : if (MAXVAL(Psps%indlmn(1,:,:))+1 >= 5) then
- + ]
1326 : : write(msg,'(3a)')&
1327 : 0 : & ' Pseudopotentials with f-projectors not implemented',ch10,&
1328 : 0 : & ' Program does not stop but _KSS file will not be created...'
1329 : 0 : ABI_WARNING(msg)
1330 : 0 : ierr=ierr+1
1331 : : end if
1332 : : !* Check useylm
1333 [ + + - + ]: 6 : if (Psps%useylm/=0.and.Psps%usepaw==0) then
1334 : : write(msg,'(3a)')&
1335 : 0 : & ' The present version of outkss does not work with useylm/=0 !',ch10,&
1336 : 0 : & ' Program does not stop but _KSS file will not be created...'
1337 : 0 : ABI_WARNING(msg)
1338 : 0 : ierr=ierr+1
1339 : : end if
1340 : : !* Check PAW and kssform value
1341 [ + + ]: 6 : if (Psps%usepaw/=0) then
1342 [ + - - + ]: 4 : if (nprocs>1.and.kssform==1) then
1343 : : write(msg,'(3a)')&
1344 : 0 : & ' Parallel PAW with kssform=1, not yet allowed',ch10,&
1345 : 0 : & ' Program does not stop but _KSS file will not be created...'
1346 : 0 : ABI_WARNING(msg)
1347 : 0 : ierr=ierr+1
1348 : : end if
1349 [ + - - + ]: 4 : if (kssform==3.and.usecprj/=1) then
1350 : : write(msg,'(3a)')&
1351 : 0 : & ' If PAW and kssform=3, usecprj must be 1',ch10,&
1352 : 0 : & ' Program does not stop but _KSS file will not be created...'
1353 : 0 : ABI_WARNING(msg)
1354 : 0 : ierr=ierr+1
1355 : : end if
1356 : : end if
1357 : : !* Check parallelization
1358 [ - + ]: 6 : if (MPI_enreg%paralbd/=0) then
1359 : : write(msg,'(3a)')&
1360 : 0 : & ' outkss cannot be used with parallelization on bands (paralbd/=0) !',ch10,&
1361 : 0 : & ' Program does not stop but _KSS file will not be created...'
1362 : 0 : ABI_WARNING(msg)
1363 : 0 : ierr=ierr+1
1364 : : end if
1365 [ - + ]: 6 : if (MPI_enreg%paral_spinor/=0) then
1366 : : write(msg,'(3a)')&
1367 : 0 : & ' outkss cannot be used yet with parallelization on nspinors !',ch10,&
1368 : 0 : & ' Program does not stop but _KSS file will not be created...'
1369 : 0 : ABI_WARNING(msg)
1370 : 0 : ierr=ierr+1
1371 : :
1372 : : endif
1373 [ - + ]: 6 : if (ierr/=0) then
1374 : : write(msg,'(3a)')&
1375 : 0 : & ' outkss: Not allowed options found !',ch10,&
1376 : 0 : & ' Program does not stop but _KSS file will not be created...'
1377 : 0 : call wrtout(std_out,msg,'COLL')
1378 : 0 : call wrtout(ab_out,msg,'COLL')
1379 : 0 : write(msg,'(a)')' outkss: see the log file for more information.'
1380 : 0 : call wrtout(ab_out,msg,'COLL')
1381 : 0 : RETURN ! Houston we have a problem!
1382 : : end if
1383 : : !
1384 : : !Estimate required memory in case of diagonalization.
1385 : : !TODO to be modified to take into account the case nsppol=2
1386 [ + + ]: 6 : if (kssform/=3) then
1387 : 1 : call memkss(mband,mgfft,mproj,Psps%mpssoang,mpw,natom,Dtset%ngfft,nkpt,dtset%nspinor,nsym,ntypat)
1388 : : end if
1389 : : !
1390 : : !=== Initialize some variables ===
1391 [ - + - - ]: 6 : if (nsppol==2) stag(:)=(/'SPIN UP: ','SPIN DOWN:'/)
1392 : 6 : n1=Dtset%ngfft(1); n2=Dtset%ngfft(2); n3=Dtset%ngfft(3)
1393 : 6 : n4=Dtset%ngfft(4); n5=Dtset%ngfft(5); n6=Dtset%ngfft(6)
1394 : 6 : ecut_eff = ecut * Dtset%dilatmx**2 ! Use ecut_eff instead of ecut_eff since otherwise
1395 : : !one cannot restart from a previous density file
1396 : 6 : sizepw=2*mpw ; do_diago=(kssform/=3)
1397 [ + + + - : 18 : ABI_MALLOC(dimlmn,(natom*Psps%usepaw))
+ - + + -
+ - + ]
1398 [ + + ]: 6 : if (Psps%usepaw==1) then
1399 : 4 : call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,Dtset%typat,pawtab,'R')
1400 : : end if
1401 : : !
1402 : : !============================================================
1403 : : !=== Prepare set containing all G-vectors sorted by stars ===
1404 : : !============================================================
1405 : 6 : write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.'
1406 : 6 : call wrtout(std_out,msg,'COLL')
1407 : : !
1408 : : !=== Analyze symmetry operations ===
1409 [ - + ]: 6 : if (Dtset%symmorphi==0) then ! Old (Obsolete) implementation: Suppress inversion from symmetries list:
1410 : 0 : nullify(symrel2,tnons2)
1411 : 0 : call remove_inversion(nsym,Dtset%symrel,Dtset%tnons,nsym2,symrel2,tnons2,pinv)
1412 [ # # # # : 0 : if (ANY(ABS(tnons2(:,1:nsym2))>tol8)) then
# # # # ]
1413 : : write(msg,'(3a)')&
1414 : 0 : & ' Non-symmorphic operations still remain in the symmetries list ',ch10,&
1415 : 0 : & ' Program does not stop but _KSS file will not be created...'
1416 : 0 : ABI_WARNING(msg)
1417 : 0 : ierr=ierr+1 ; RETURN
1418 : : end if
1419 [ + - ]: 6 : else if (Dtset%symmorphi==1) then
1420 : : ! If in the input file symmorphi==1 all the symmetry operations are retained:
1421 : : ! both identity and inversion (if any) as well as non-symmorphic operations.
1422 : 6 : nsym2=nsym ; pinv=1
1423 [ - + + - : 18 : ABI_MALLOC(symrel2,(3,3,nsym))
+ - + - -
+ - + ]
1424 [ - + + - : 18 : ABI_MALLOC(tnons2,(3,nsym))
+ - + - -
+ - + ]
1425 [ + + + + : 3750 : symrel2(:,:,:)=Dtset%symrel(:,:,1:nsym)
+ + ]
1426 [ + + + + ]: 1158 : tnons2(:,:) =Dtset%tnons(:,1:nsym)
1427 : : else
1428 : : write(msg,'(a,i4,3a)')&
1429 : 0 : & ' symmorphi = ',Dtset%symmorphi,' while it must be 0 or 1',ch10,&
1430 : 0 : & ' Program does not stop but KSS file will not be created...'
1431 : 0 : ABI_WARNING(msg)
1432 : 0 : ierr=ierr+1 ; RETURN
1433 : : end if
1434 : : !
1435 : : !===================================================================
1436 : : !==== Merge the set of k-centered G-spheres into a big set gbig ====
1437 : : !===================================================================
1438 : : !* Vectors in gbig are ordered by shells
1439 : : !
1440 : 6 : nullify(gbig,shlim)
1441 [ - + ]: 6 : call merge_and_sort_kg(nkpt,Dtset%kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim)
1442 : :
1443 : 6 : nbase = SIZE(shlim) ! Number of independent G in the big sphere.
1444 : 6 : maxpw = shlim(nbase) ! Total number of G"s in the big sphere.
1445 : : !
1446 : : !* Determine optimal number of bands and G"s to be written.
1447 : 6 : npwkss=Dtset%npwkss
1448 [ + - ]: 6 : if ((npwkss==0).or.(npwkss>=maxpw)) then
1449 : 6 : npwkss=maxpw
1450 : : write(msg,'(5a)')&
1451 : 6 : & ' Since the number of g''s to be written on file',ch10,&
1452 : 6 : & ' was 0 or too large, it has been set to the max. value.,',ch10,&
1453 : 12 : & ' computed from the union of the sets of G vectors for the different k-points.'
1454 : 6 : call wrtout(std_out,msg,'COLL')
1455 : : end if
1456 : :
1457 : 6 : ishm=0
1458 [ + + ]: 212 : do ii=1,nbase
1459 [ + - ]: 212 : if (shlim(ii)<=npwkss) then
1460 : 206 : ishm=ii
1461 : : else
1462 : : EXIT
1463 : : end if
1464 : : end do
1465 : :
1466 [ - + ]: 6 : if (shlim(ishm)/=npwkss) then
1467 : 0 : nrst1=shlim(ishm)
1468 : 0 : nrst2=MIN0(shlim(MIN0(ishm+1,nbase)),maxpw)
1469 [ # # ]: 0 : if (IABS(npwkss-nrst2)<IABS(npwkss-nrst1)) nrst1=nrst2
1470 : 0 : npwkss=nrst1
1471 [ # # ]: 0 : if (shlim(ishm)<npwkss) ishm=ishm+1
1472 : : write(msg,'(3a)')&
1473 : 0 : & ' The number of G''s to be written on file is not a whole number of stars ',ch10,&
1474 : 0 : & ' the program set it to the nearest star limit.'
1475 : 0 : call wrtout(std_out,msg,'COLL')
1476 : : end if
1477 : :
1478 : 6 : write(msg,'(a,i5)')' Number of g-vectors written on file is: ',npwkss
1479 : 6 : call wrtout(std_out,msg,'COLL')
1480 : : !
1481 : : !=== Check on the number of stored bands ===
1482 [ + + ]: 6 : if (do_diago) then
1483 : :
1484 [ + - - + ]: 1 : if (Dtset%nbandkss==-1.or.Dtset%nbandkss>=maxpw) then
1485 [ # # ]: 0 : nbandkssk(1:nkpt)=npwarr(1:nkpt)
1486 : 0 : write(msg,'(6a)')ch10,&
1487 : 0 : & ' Since the number of bands to be computed was (-1) or',ch10,&
1488 : 0 : & ' too large, it has been set to the max. value. allowed for each k,',ch10,&
1489 : 0 : & ' thus, the minimum of the number of plane waves for each k point.'
1490 : 0 : call wrtout(std_out,msg,'COLL')
1491 : : else
1492 [ + + ]: 11 : nbandkssk(1:nkpt)=Dtset%nbandkss
1493 : 1 : found=.FALSE.
1494 [ + + ]: 11 : do ikpt=1,nkpt
1495 [ - + ]: 11 : if (Dtset%nbandkss>npwarr(ikpt)) then
1496 : 0 : nbandkssk(ikpt)=npwarr(ikpt)
1497 : 0 : found=.TRUE.
1498 : : end if
1499 : : end do
1500 [ - + ]: 1 : if (found) then
1501 : : write(msg,'(7a)')&
1502 : 0 : & ' The value choosen for the number of bands in file',ch10,&
1503 : 0 : & ' (nbandkss) was greater than at least one number of plane waves ',ch10,&
1504 : 0 : & ' for a given k-point (npw_k).',ch10,' It has been modified consequently.'
1505 : 0 : ABI_WARNING(msg)
1506 : : end if
1507 : : end if
1508 : 1 : found=.FALSE.
1509 [ + + ]: 11 : do ikpt=1,nkpt
1510 [ - + ]: 11 : if (nbandkssk(ikpt)>npwkss) then
1511 : 0 : nbandkssk(ikpt)=npwkss
1512 : 0 : found=.TRUE.
1513 : : end if
1514 : : end do
1515 [ - + ]: 1 : if (found) then
1516 : : write(msg,'(5a)')&
1517 : 0 : & ' The number of bands to be computed (for one k) was',ch10,&
1518 : 0 : & ' greater than the number of g-vectors to be written.',ch10,&
1519 : 0 : & ' It has been modified consequently.'
1520 : 0 : ABI_WARNING(msg)
1521 : : end if
1522 [ + + ]: 11 : nbandksseff=MINVAL(nbandkssk)
1523 : :
1524 : : else ! .not. do_diago
1525 [ + + ]: 59 : do ikpt=1,nkpt
1526 [ + + ]: 113 : do isppol=1,nsppol
1527 : 108 : nbandkssk(ikpt)=Dtset%nband(ikpt+(isppol-1)*nkpt)
1528 : : end do
1529 : : end do
1530 [ + + ]: 59 : nbandksseff=MINVAL(nbandkssk)
1531 [ + - - + ]: 5 : if (Dtset%nbandkss>0 .and. Dtset%nbandkss<nbandksseff) then
1532 : : write(msg,'(a,i5,a,i5,2a)')&
1533 : 0 : & ' Number of bands calculated=',nbandksseff,', greater than nbandkss=',Dtset%nbandkss,ch10,&
1534 : 0 : & ' will write nbandkss bands on the KSS file'
1535 : 0 : ABI_COMMENT(msg)
1536 : 0 : nbandksseff=Dtset%nbandkss
1537 : : end if
1538 : : end if
1539 : :
1540 : 6 : write(msg,'(a,i5)')' Number of bands written on file is: ',nbandksseff
1541 : 6 : call wrtout(std_out,msg,'COLL')
1542 : :
1543 [ + - - + ]: 6 : found= ANY(nbandkssk(1:nkpt)<npwarr(1:nkpt))
1544 : :
1545 [ + + ]: 6 : if (do_diago) then
1546 [ + - ]: 1 : if (found) then
1547 : 1 : write(msg,'(6a)')ch10,&
1548 : 1 : & ' Since the number of bands to be computed',ch10,&
1549 : 1 : & ' is less than the number of G-vectors found,',ch10,&
1550 : 2 : & ' the program will perform partial diagonalizations.'
1551 : : else
1552 : 0 : write(msg,'(6a)')ch10,&
1553 : 0 : & ' Since the number of bands to be computed',ch10,&
1554 : 0 : & ' is equal to the nb of G-vectors found for each k-pt,',ch10,&
1555 : 0 : & ' the program will perform complete diagonalizations.'
1556 : : end if
1557 : 1 : call wrtout(std_out,msg,'COLL')
1558 : : end if
1559 : : !
1560 : : !==========================================================================
1561 : : !=== Open KSS file for output, write header with dimensions and kb sign ===
1562 : : !==========================================================================
1563 : : !
1564 : : !* Output required disk space.
1565 : 6 : call dsksta(ishm,Psps%usepaw,nbandksseff,mpsang,natom,ntypat,npwkss,nkpt,dtset%nspinor,nsppol,nsym2,dimlmn)
1566 : :
1567 [ + + ]: 6 : if (my_rank==master) then
1568 : : call write_kss_header(dtfil%fnameabo_kss,npwkss,ishm,nbandksseff,mband,nsym2,symrel2,tnons2,occ,gbig,shlim,&
1569 [ - + - + : 3 : & crystal,Dtset,Hdr,Psps,iomode,untkss)
- + - + ]
1570 : : end if
1571 : :
1572 [ - + ]: 6 : ABI_FREE(shlim)
1573 : :
1574 [ + + ]: 6 : if ( do_diago) msg = ' Diagonalized eigenvalues'
1575 : 5 : if (.not.do_diago) msg = ' Conjugate gradient eigenvalues'
1576 : 6 : call wrtout(ab_out,msg,'COLL')
1577 : :
1578 [ - + ]: 6 : if (Dtset%enunit==1) then
1579 : 0 : msg=' k eigenvalues [eV]'
1580 : : else
1581 : 6 : msg=' k eigenvalues [Hartree]'
1582 : : end if
1583 : 6 : call wrtout(ab_out,msg,'COLL')
1584 : : !
1585 : : !=== In case of PAW distributed atomic sites, need to retrieve the full paw_ij%dij ===
1586 [ + + - + : 6 : if (do_diago.and.Psps%usepaw==1.and.MPI_enreg%nproc_atom>1) then
- - ]
1587 [ # # # # : 0 : ABI_MALLOC(Paw_ij_all,(dtset%natom))
# # # # #
# # # #
# ]
1588 : 0 : call paw_ij_gather(Paw_ij,Paw_ij_all,-1,MPI_enreg%comm_atom)
1589 : : else
1590 : 6 : paw_ij_all => paw_ij
1591 : : end if
1592 : :
1593 : :
1594 : 6 : call timab(934,2,tsec) ! outkss(Gsort+hd)
1595 : : !
1596 : :
1597 : 6 : k_index=0; bdtot_index=0; ibg=0
1598 : :
1599 [ + + ]: 12 : do isppol=1,nsppol ! Loop over spins
1600 : : !
1601 [ + + ]: 76 : do ikpt=1,nkpt ! Loop over k-points.
1602 : 64 : call timab(935,1,tsec) ! outkss(k-Loop)
1603 : :
1604 : 64 : nband_k =Dtset%nband(ikpt+(isppol-1)*nkpt)
1605 : 64 : npw_k =npwarr(ikpt)
1606 : 64 : istwf_k =Dtset%istwfk(ikpt)
1607 [ + + ]: 256 : kpoint =Dtset%kptns(:,ikpt)
1608 : 64 : nbandkss_k=nbandkssk(ikpt)
1609 : 64 : mtag =5*(ikpt+(isppol-1)*nkpt)
1610 : :
1611 : :
1612 : : ! Get G-vectors, for this k-point.
1613 [ - + ]: 64 : call get_kg(kpoint,istwf_k,ecut_eff,gmet,onpw_k,kg_k)
1614 [ - + ]: 64 : ABI_CHECK(onpw_k==npw_k,"Mismatch in npw_k")
1615 : : !
1616 : : ! ============================================
1617 : : ! ==== Parallelism over k-points and spin ====
1618 : : ! ============================================
1619 [ + + ]: 64 : if (MPI_enreg%proc_distrb(ikpt,1,isppol)==my_rank) then
1620 : :
1621 : 31 : write(msg,'(2a,i3,3x,a)')ch10,' k-point ',ikpt,stag(isppol)
1622 : 31 : call wrtout(std_out, msg)
1623 : :
1624 [ + + ]: 31 : if (do_diago) then
1625 : : ! Direct diagonalization of the KS Hamiltonian.
1626 [ - + ]: 10 : ABI_SFREE_PTR(eig_ene)
1627 [ - + ]: 10 : ABI_SFREE_PTR(eig_vec)
1628 : 10 : comm_self = xmpi_comm_self
1629 : :
1630 : 10 : call timab(936,1,tsec)
1631 : :
1632 : : call init_ddiago_ctl(Diago_ctl,"Vectors",isppol,dtset%nspinor,ecut_eff,Dtset%kptns(:,ikpt),Dtset%nloalg,gmet,&
1633 : 10 : nband_k=nbandkssk(ikpt),effmass_free=Dtset%effmass_free,istwf_k=Dtset%istwfk(ikpt),prtvol=Dtset%prtvol)
1634 : :
1635 : : call ksdiago(Diago_ctl,nbandkssk(ikpt),Dtset%nfft,mgfft,Dtset%ngfft,natom,&
1636 : : Dtset%typat,nfft,dtset%nspinor,nspden,nsppol,Pawtab,Pawfgr,Paw_ij_all,&
1637 [ - + ]: 10 : Psps,rprimd,vtrial,xred,onband_diago,eig_ene,eig_vec,Cprj_diago_k,comm_self,ierr)
1638 : :
1639 : 10 : call timab(936,2,tsec)
1640 : : end if
1641 : :
1642 : : end if ! END of kpt+spin parallelism.
1643 : : !
1644 : : ! ===========================================================
1645 : : ! ==== Transfer data between master and the working proc ====
1646 : : ! ===========================================================
1647 : 64 : call timab(937,1,tsec) !outkss(MPI_exch)
1648 [ - + - + ]: 64 : ABI_MALLOC(Cprjnk_k,(0,0))
1649 [ + + ]: 64 : if (nprocs==1) then
1650 : :
1651 [ - + ]: 20 : if (Psps%usepaw==1) then ! Copy projectors for this k-point
1652 : 0 : n2dim=min(nbandksseff*dtset%nspinor,onband_diago)
1653 [ # # ]: 0 : if (kssform==3) n2dim=nband_k*dtset%nspinor
1654 [ # # ]: 0 : ABI_FREE(Cprjnk_k)
1655 [ # # # # : 0 : ABI_MALLOC(Cprjnk_k,(natom,n2dim))
# # # # #
# # # # #
# # # # #
# ]
1656 : 0 : call pawcprj_alloc(Cprjnk_k,0,dimlmn)
1657 [ # # ]: 0 : if (kssform==3) then
1658 : 0 : call pawcprj_copy(Cprj(:,ibg+1:ibg+dtset%nspinor*nband_k),Cprjnk_k)
1659 : : else
1660 : : !ABI_WARNING("Here I have to use onband_diago") !FIXME
1661 : 0 : call pawcprj_copy(Cprj_diago_k(:,1:n2dim),Cprjnk_k)
1662 : : end if
1663 : : end if
1664 : :
1665 : : else
1666 : : !parallel case
1667 : :
1668 : 44 : receiver=master; sender=MPI_enreg%proc_distrb(ikpt,1,isppol)
1669 : :
1670 : 44 : bufsz=nbandksseff/bufnb; bufrt=nbandksseff-bufnb*bufsz
1671 : :
1672 [ + + ]: 44 : if (my_rank==receiver.or.my_rank==sender) then
1673 : :
1674 [ + - - - ]: 19 : if (do_diago.and.(my_rank==receiver.and.my_rank/=sender)) then ! Alloc arrays if not done yet.
1675 [ # # # # : 0 : ABI_MALLOC(eig_ene,(npw_k*dtset%nspinor))
# # # # #
# # # ]
1676 [ # # # # : 0 : ABI_MALLOC(eig_vec,(2,npw_k*dtset%nspinor,nbandkssk(ikpt)))
# # # # #
# # # # #
# # ]
1677 : : end if
1678 : :
1679 : : if (.not.do_diago) then
1680 : :
1681 [ - + + - : 76 : ABI_MALLOC(eig_vec,(2,npw_k*dtset%nspinor,nbandkssk(ikpt)))
- + + - +
- + - - +
- + ]
1682 : :
1683 [ + + ]: 19 : if (my_rank==sender) then
1684 [ + + ]: 231 : do ib=1,nbandksseff
1685 : 220 : shift = k_index + (ib-1)*npw_k*dtset%nspinor
1686 [ + + ]: 102791 : do ig=1,npw_k*dtset%nspinor
1687 [ + + ]: 307900 : eig_vec(:,ig,ib)=cg(:,ig+shift)
1688 : : end do
1689 : : end do
1690 : : end if
1691 : : !
1692 : : ! In case of PAW and kssform==3, retrieve matrix elements of the PAW projectors for this k-point
1693 [ + - ]: 19 : if (Psps%usepaw==1) then
1694 : 19 : n2dim=min(nbandksseff*dtset%nspinor,onband_diago)
1695 [ + - ]: 19 : if (kssform==3) n2dim=nband_k*dtset%nspinor
1696 [ - + ]: 19 : ABI_FREE(Cprjnk_k)
1697 [ - + + - : 836 : ABI_MALLOC(Cprjnk_k,(natom,n2dim))
- + + - +
- + - - +
- + + + +
+ ]
1698 : 19 : call pawcprj_alloc(Cprjnk_k,0,dimlmn)
1699 [ + + ]: 19 : if (my_rank==sender) then
1700 [ + - ]: 11 : if (kssform==3) then
1701 : 11 : call pawcprj_copy(Cprj(:,ibg+1:ibg+dtset%nspinor*nband_k),Cprjnk_k)
1702 : : else
1703 : : !ABI_WARNING("Here I have to use onband_diago") !FIXME
1704 : 0 : call pawcprj_copy(Cprj_diago_k(:,1:n2dim),Cprjnk_k)
1705 : : end if
1706 : : end if
1707 [ + + ]: 19 : if (sender/=receiver) then
1708 : 16 : call pawcprj_mpi_exch(natom,n2dim,dimlmn,0,Cprjnk_k,Cprjnk_k,sender,receiver,spaceComm,mtag+4,ierr)
1709 : : end if
1710 : : end if ! usepaw
1711 : :
1712 : : else ! do_diago
1713 : 0 : call xmpi_exch(eig_ene,nbandksseff,sender,eig_ene,receiver,spaceComm,mtag+1,ierr)
1714 : : end if
1715 : :
1716 : : ! Exchange eigenvectors.
1717 [ + - ]: 19 : if (bufsz>0) then
1718 [ + + ]: 399 : do i=0,bufnb-1
1719 : : call xmpi_exch(eig_vec(:,:,i*bufsz+1:(i+1)*bufsz),2*npw_k*dtset%nspinor*bufsz,&
1720 : 399 : & sender,eig_vec(:,:,i*bufsz+1:(i+1)*bufsz),receiver,spaceComm,mtag+2,ierr)
1721 : : end do
1722 : : end if
1723 [ - + ]: 19 : if (bufrt>0) then
1724 : : call xmpi_exch(eig_vec(:,:,bufnb*bufsz+1:bufnb*bufsz+bufrt),2*npw_k*dtset%nspinor*bufrt,&
1725 : 0 : & sender,eig_vec(:,:,bufnb*bufsz+1:bufnb*bufsz+bufrt),receiver,spaceComm,mtag+3,ierr)
1726 : : end if
1727 : :
1728 : : end if
1729 : : end if !nprocs > 1
1730 : 64 : call timab(937,2,tsec) !outkss(MPI_exch)
1731 : :
1732 : 64 : call timab(938,1,tsec) !outkss(write)
1733 : :
1734 [ + + ]: 64 : if (my_rank==master) then ! Prepare data for writing on disk.
1735 [ - + + - : 93 : ABI_MALLOC(ene,(nbandksseff))
+ - + - -
+ - + ]
1736 [ - + + - : 124 : ABI_MALLOC(wfg,(2,npwkss*dtset%nspinor,nbandksseff))
- + + - +
- + - - +
- + ]
1737 [ + + + + : 2430991 : ene=zero; wfg=zero
+ + + + ]
1738 : :
1739 [ + + ]: 31 : if (.not.do_diago) then
1740 [ + + ]: 341 : ene(1:nbandksseff)=eigen(1+bdtot_index:nbandksseff+bdtot_index)
1741 : :
1742 [ + + ]: 21 : if (nprocs>1) then
1743 : : call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,&
1744 [ - + - + ]: 11 : & MPI_enreg,gbig,wfg,eig_vec=eig_vec)
1745 : : else
1746 : : call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,&
1747 [ - + ]: 10 : & MPI_enreg,gbig,wfg,icg=k_index,cg=cg)
1748 : : end if
1749 : :
1750 : : else ! Direct diagonalization.
1751 [ + + ]: 1010 : ene(1:nbandksseff)=eig_ene(1:nbandksseff)
1752 : :
1753 : : ! FIXME: receiver does not know Diago_ctl%npw_k
1754 : : call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,&
1755 [ - + - + ]: 10 : & MPI_enreg,gbig,wfg,eig_vec=eig_vec)
1756 : :
1757 : : ! * Check diagonalized eigenvalues with respect to conjugate gradient ones
1758 : 10 : ntemp=MIN(nbandksseff,nband_k)
1759 [ + + + - : 100 : if (ANY(ABS(ene(1:ntemp)-eigen(1+bdtot_index:ntemp+bdtot_index))>tol3)) then
- + ]
1760 : : write(msg,'(3a)')&
1761 : 0 : & ' The diagonalized eigenvalues differ by more than 10^-3 Hartree',ch10,&
1762 : 0 : & ' with respect to the conjugated gradient values.'
1763 : 0 : ABI_WARNING(msg)
1764 : : end if
1765 : : end if
1766 : : !
1767 : : ! * Write out energies
1768 [ - + ]: 31 : if (Dtset%enunit==1) then
1769 : 0 : cfact=Ha_eV ; frmt1='(i4,4x,9(1x,f7.2))' ; frmt2='(8x,9(1x,f7.2))'
1770 : 0 : write(msg,'(a,i3,3x,a)')' Eigenvalues in eV for ikpt= ',ikpt,stag(isppol)
1771 : : else
1772 : 31 : cfact=one ; frmt1='(i4,4x,9(1x,f7.4))' ; frmt2='(8x,9(1x,f7.4))'
1773 : 31 : write(msg,'(a,i3,3x,a)')' Eigenvalues in Hartree for ikpt= ',ikpt,stag(isppol)
1774 : : end if
1775 : 31 : call wrtout(std_out,msg,'COLL')
1776 : :
1777 [ + - + + ]: 310 : write(msg,frmt1)ikpt,(ene(ib)*cfact,ib=1,MIN(9,nbandksseff))
1778 : 31 : call wrtout(std_out,msg,'COLL')
1779 : 31 : call wrtout(ab_out,msg,'COLL')
1780 : :
1781 [ + - ]: 31 : if (nbandksseff>9) then
1782 [ - + + + ]: 173 : do j=10,nbandksseff,9
1783 [ + - + + ]: 1183 : write(msg,frmt2) (ene(ib)*cfact,ib=j,MIN(j+8,nbandksseff))
1784 : 142 : call wrtout(std_out,msg,'COLL')
1785 : 173 : call wrtout(ab_out,msg,'COLL')
1786 : : end do
1787 : : end if
1788 : :
1789 : : if (skip_test_ortho) then ! Set this if to FALSE to skip test below
1790 : : einf=one; esup=one; cinf=zero; csup=zero
1791 : : else
1792 : : !
1793 : : ! Test on the normalization of wavefunctions.
1794 : : ibsp=0
1795 [ + + ]: 1351 : do ib=1,nbandksseff
1796 : : norm=zero
1797 [ + + ]: 2640 : do ispinor=1,dtset%nspinor
1798 : 1320 : ibsp=ibsp+1
1799 : 1320 : spinor_shift1=(ispinor-1)*npwkss
1800 : 1320 : ug1 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ib)
1801 : :
1802 : : !ovlp(1) =ddot(npwkss,ug1(1,:),1,ug1(1,:),1) + ddot(npwkss,ug1(2,:),1,ug1(2,:),1)
1803 : 1320 : ovlp(1) = cg_dznrm2(npwkss,ug1)
1804 : 1320 : ovlp(1) = ovlp(1)**2
1805 : 1320 : ovlp(2) = zero
1806 [ + + ]: 1320 : if (Psps%usepaw==1) ovlp = ovlp &
1807 : : & + paw_overlap(Cprjnk_k(:,ibsp:ibsp),Cprjnk_k(:,ibsp:ibsp),Dtset%typat,Pawtab,&
1808 [ + + ]: 660 : & spinor_comm=MPI_enreg%comm_spinor)
1809 : 2640 : norm = norm + DABS(ovlp(1))
1810 : : end do
1811 [ + + ]: 1320 : if (norm<einf) einf=norm
1812 [ + + ]: 1351 : if (norm>esup) esup=norm
1813 : : end do
1814 : : !
1815 : : ! Test on the orthogonalization of wavefunctions.
1816 [ + + ]: 1351 : do ib=1,nbandksseff
1817 : 1320 : pad1=(ib-1)*dtset%nspinor
1818 [ + + ]: 53391 : do ibp=ib+1,nbandksseff
1819 : 52040 : pad2=(ibp-1)*dtset%nspinor
1820 : 52040 : ovlp(:)=zero
1821 [ + + ]: 104080 : do ispinor=1,dtset%nspinor
1822 : 52040 : ibsp1=pad1+ispinor
1823 : 52040 : ibsp2=pad2+ispinor
1824 : 52040 : spinor_shift1=(ispinor-1)*npwkss
1825 : 52040 : ug1 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ib )
1826 : 52040 : ug2 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ibp)
1827 : :
1828 : : !ovlp(1)=ddot(npwkss,ug1(1,:),1,ug2(1,:),1) + ddot(npwkss,ug1(2,:),1,ug2(2,:),1)
1829 : : !ovlp(2)=ddot(npwkss,ug1(1,:),1,ug2(2,:),1) - ddot(npwkss,ug1(2,:),1,ug2(1,:),1)
1830 : 52040 : ovlp = cg_zdotc(npwkss,ug1,ug2)
1831 : :
1832 [ + + ]: 52040 : if (Psps%usepaw==1) ovlp= ovlp &
1833 : : & + paw_overlap(Cprjnk_k(:,ibsp1:ibsp1),Cprjnk_k(:,ibsp2:ibsp2),Dtset%typat,Pawtab,&
1834 [ + + ]: 58310 : & spinor_comm=MPI_enreg%comm_spinor)
1835 : : end do
1836 : 52040 : norm = DSQRT(ovlp(1)**2+ovlp(2)**2)
1837 [ + + ]: 52040 : if (norm<cinf) cinf=norm
1838 [ + + ]: 53360 : if (norm>csup) csup=norm
1839 : : end do
1840 : : end do
1841 : : end if
1842 : :
1843 : 31 : write(msg,'(a,i3,3x,a)')' Writing out eigenvalues/vectors for ikpt=',ikpt,stag(isppol)
1844 : 31 : call wrtout(std_out,msg,'COLL')
1845 : : !
1846 : : ! * Write occupation numbers on std_out.
1847 [ - + + - : 93 : ABI_MALLOC(occ_k,(MAX(nband_k,nbandksseff)))
+ - + - -
+ - + ]
1848 [ + + ]: 441 : occ_k(1:nband_k)=occ(1+bdtot_index:nband_k+bdtot_index)
1849 [ + + + + ]: 941 : if (nband_k < nbandksseff) occ_k(nband_k+1:nbandksseff)=zero
1850 : :
1851 : 31 : write(msg,'(a,i3,3x,a)')' Occupation numbers for ikpt=',ikpt,stag(isppol)
1852 : 31 : call wrtout(std_out,msg,'COLL')
1853 : 31 : write(msg,'(i4,4x,9(1x,f7.4))')ikpt,(occ_k(ib),ib=1,MIN(9,nbandksseff))
1854 : 31 : call wrtout(std_out,msg,'COLL')
1855 [ + - ]: 31 : if (nbandksseff>9) then
1856 [ - + + + ]: 173 : do j=10,nbandksseff,9
1857 : 142 : write(msg,'(8x,9(1x,f7.4))') (occ_k(ib),ib=j,MIN(j+8,nbandksseff))
1858 : 173 : call wrtout(std_out,msg,'COLL')
1859 : : end do
1860 : : end if
1861 : : !
1862 : : ! =================================================================
1863 : : ! ==== Write wavefunctions, KB and PAW matrix elements on disk ====
1864 : : ! =================================================================
1865 : : call write_kss_wfgk(untkss,ikpt,isppol,kpoint,dtset%nspinor,npwkss,&
1866 [ - + ]: 31 : nbandksseff,natom,Psps,ene,occ_k,rprimd,gbig,wfg,Cprjnk_k,iomode)
1867 : :
1868 [ - + ]: 31 : ABI_FREE(occ_k)
1869 [ - + ]: 31 : ABI_FREE(ene)
1870 [ - + ]: 31 : ABI_FREE(wfg)
1871 : :
1872 : : end if ! my_rank==master
1873 : 64 : call timab(938,2,tsec) !outkss(write)
1874 : :
1875 [ + + + + ]: 64 : if (my_rank==master.or.my_rank==MPI_enreg%proc_distrb(ikpt,1,isppol)) then
1876 [ + + ]: 39 : ABI_SFREE_PTR(eig_ene)
1877 [ + + ]: 39 : ABI_SFREE_PTR(eig_vec)
1878 [ + + ]: 39 : if (Psps%usepaw==1) call pawcprj_free(Cprjnk_k)
1879 : : end if
1880 [ + - + + : 444 : ABI_FREE(Cprjnk_k)
- + - + -
+ ]
1881 [ + - ]: 64 : ABI_SFREE(kg_k)
1882 : :
1883 : : ! if (MPI_enreg%paral_compil_kpt==1) then !cannot be used in seq run!
1884 [ + + ]: 64 : if (.not.(proc_distrb_cycle(MPI_enreg%proc_distrb,ikpt,1,nband_k,isppol,my_rank))) then
1885 : 31 : k_index=k_index+npw_k*nband_k*dtset%nspinor
1886 : 31 : ibg=ibg+dtset%nspinor*nband_k
1887 : : end if
1888 : 64 : bdtot_index=bdtot_index+nband_k
1889 : :
1890 : 64 : call xmpi_barrier(spaceComm) ! FIXME this barrier is detrimental in the case of direct diago!
1891 : :
1892 : 70 : call timab(935,2,tsec) !outkss(k-loop)
1893 : : end do ! ! End loop over k-points.
1894 : : end do ! spin
1895 : :
1896 : : write(msg,'(3a,f9.6,2a,f9.6,4a,f9.6,2a,f9.6,a)')&
1897 : 6 : & ' Test on the normalization of the wavefunctions',ch10,&
1898 : 6 : & ' min sum_G |a(n,k,G)| = ',einf,ch10,&
1899 : 6 : & ' max sum_G |a(n,k,G)| = ',esup,ch10,&
1900 : 6 : & ' Test on the orthogonalization of the wavefunctions',ch10,&
1901 : 6 : & ' min sum_G a(n,k,G)a(n'',k,G) = ',cinf,ch10,&
1902 : 12 : & ' max sum_G a(n,k,G)a(n'',k,G) = ',csup,ch10
1903 : 6 : call wrtout(std_out,msg,'COLL')
1904 : 6 : call wrtout(ab_out,msg,'COLL')
1905 : :
1906 [ - + ]: 6 : ABI_FREE(gbig)
1907 [ - + ]: 6 : ABI_FREE(symrel2)
1908 [ - + ]: 6 : ABI_FREE(tnons2)
1909 [ + + ]: 6 : if (Psps%usepaw==1) then
1910 [ - + ]: 4 : ABI_FREE(dimlmn)
1911 [ - + - - ]: 4 : if (do_diago.and.MPI_enreg%nproc_atom>1) then
1912 [ # # ]: 0 : ABI_FREE(Paw_ij_all)
1913 : : end if
1914 : : end if
1915 : : !
1916 : : !* Close file
1917 [ + + ]: 6 : if (my_rank==master) then
1918 [ - + ]: 3 : if (iomode==IO_MODE_FORTRAN) close(unit=untkss)
1919 [ + - ]: 3 : if (iomode==IO_MODE_ETSF) then
1920 [ - + ]: 3 : NCF_CHECK(nf90_close(untkss))
1921 : : end if
1922 : : end if
1923 : :
1924 [ - + ]: 6 : if (associated(Cprj_diago_k)) then
1925 : 0 : call pawcprj_free(Cprj_diago_k)
1926 [ # # ]: 0 : ABI_FREE(Cprj_diago_k)
1927 : : end if
1928 : :
1929 [ - + ]: 6 : if (lhack) then
1930 [ # # ]: 0 : ABI_FREE(MPI_enreg%proc_distrb)
1931 : : end if
1932 : :
1933 : 6 : call wrtout(std_out, "outkss done", "COLL")
1934 : 6 : call xmpi_barrier(spaceComm)
1935 : :
1936 : : DBG_EXIT("COLL")
1937 [ - + + + : 12 : call timab(933,2,tsec) ! outkss
- + - - -
- - - ]
1938 : :
1939 : : contains
1940 : : !!***
1941 : :
1942 : : !!****f* ABINIT/memkss
1943 : : !! NAME
1944 : : !! memkss
1945 : : !!
1946 : : !! FUNCTION
1947 : : !! This routine evaluates the additional amount of memory required
1948 : : !! by routine 'outkss'.
1949 : : !!
1950 : : !! INPUTS
1951 : : !! mband=maximum number of bands
1952 : : !! mgfft=maximum size of 1D FFTs
1953 : : !! mproj=maximum dimension for number of projection operators for each
1954 : : !! angular momentum for nonlocal pseudopotential
1955 : : !! mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
1956 : : !! mpw=maximum dimensioned size of npw.
1957 : : !! natom=number of atoms in cell.
1958 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1959 : : !! nkpt=number of k points.
1960 : : !! nspinor=number of spinorial components of the wavefunctions
1961 : : !! nsym=number of symmetries in space group
1962 : : !! ntypat=number of types of atoms in unit cell.
1963 : : !!
1964 : : !! NOTES
1965 : : !! This routine is not available for paw calculations
1966 : : !!
1967 : : !! SOURCE
1968 : :
1969 : 1 : subroutine memkss(mband,mgfft,mproj,mpsang,mpw,natom,ngfft,nkpt,nspinor,nsym,ntypat)
1970 : :
1971 : : !Arguments ------------------------------------
1972 : : !scalars
1973 : : integer,intent(in) :: mband,mgfft,mproj,mpsang,mpw,natom,nkpt,nspinor
1974 : : integer,intent(in) :: nsym,ntypat
1975 : : !arrays
1976 : : integer,intent(in) :: ngfft(18)
1977 : :
1978 : : !Local variables-------------------------------
1979 : : !scalars
1980 : : integer(i8b) :: isize,memsize
1981 : : character(len=500) :: msg
1982 : :
1983 : : ! *********************************************************************
1984 : : !
1985 : 1 : isize=580+fnlen+4*(81+nkpt+9*nsym)+8*15 !non allocatable var.
1986 : 1 : if(xmpi_paral==1)then
1987 : : isize=isize+4*4 !kpt_distrb
1988 : : end if
1989 : 1 : memsize=isize
1990 : 1 : isize=isize+4*nkpt+12*mpw+20*nkpt*mpw !nbasek,gbasek,cnormk,gcurr
1991 : 1 : memsize=max(memsize,isize)
1992 : : if(xmpi_paral==1)then
1993 : 1 : isize=isize+12*mpw*nkpt !ibuf1,ibuf2,rbuf1
1994 : 1 : memsize=max(memsize,isize)
1995 : 1 : isize=isize-12*mpw*nkpt !ibuf1,ibuf2,rbuf1
1996 : : end if
1997 : 1 : isize=isize+40*mpw !gbase,cnorm
1998 : 1 : memsize=max(memsize,isize)
1999 : 1 : isize=isize-4*nkpt-20*mpw*nkpt !nbasek,gbasek,cnormk
2000 : 1 : isize=isize+4*mpw !insort
2001 : 1 : memsize=max(memsize,isize)
2002 : 1 : isize=isize-16*mpw !cnorm
2003 : 1 : isize=isize+28*mpw+24*nsym !gbig,nshell,gshell
2004 : 1 : memsize=max(memsize,isize)
2005 : 1 : isize=isize+4*mpw !shlim
2006 : 1 : memsize=max(memsize,isize)
2007 : 1 : isize=isize-44*mpw-24*nsym !gcurr,gbase,gshell,insort,nshell
2008 : 1 : isize=isize-4*mpw !shlim
2009 : : isize=isize+8*mpw*nspinor&
2010 : 1 : & +16*mpw*nspinor*(mpw*nspinor+1) !eigval,eigvec
2011 : 1 : memsize=max(memsize,isize)
2012 : : isize=isize+8*mpw+8*ngfft(4)&
2013 : 1 : & *ngfft(5)*ngfft(6) !ts,vlocal
2014 : 1 : memsize=max(memsize,isize)
2015 : 1 : isize=isize+8*mgfft+4+28*mpw !gbound,indpw_k,kg_k
2016 : 1 : memsize=max(memsize,isize)
2017 : : isize=isize+8*natom&
2018 : 1 : & +24*mpw*ntypat*mpsang*mproj !phkxred,ffnl,kinpw
2019 : 1 : memsize=max(memsize,isize)
2020 : 1 : isize=isize+16*mpw*natom !ph3d
2021 : 1 : memsize=max(memsize,isize)
2022 : : isize=isize+48*mpw*nspinor&
2023 : 1 : & +8*mpw*nspinor*(mpw*nspinor+1) !pwave,subghg,gvnlg
2024 [ - + ]: 1 : if (nspinor==2)&
2025 : 0 : & isize=isize+40*mpw*nspinor !pwave_so,subghg_so
2026 : 1 : memsize=max(memsize,isize)
2027 : 1 : isize=isize+8*mpw*nspinor*(mpw*nspinor+1) !ghg
2028 : 1 : memsize=max(memsize,isize)
2029 : 1 : isize=isize+8*ngfft(4)*ngfft(5)*ngfft(6) !work
2030 : 1 : memsize=max(memsize,isize)
2031 : 1 : isize=isize-8*ngfft(4)*ngfft(5)*ngfft(6) !work
2032 : : isize=isize-8*mgfft+4+28*mpw& !gbound,indpw_k,kg_k
2033 : : &-8*natom-24*mpw*ntypat*mpsang*mproj& !phkxred,ffnl,kinpw
2034 : 1 : &-16*mpw*natom !ph3d
2035 : : isize=isize-48*mpw*nspinor&
2036 : 1 : & -8*mpw*nspinor*(mpw*nspinor+1) !pwave,subghg,gvnlg
2037 [ - + ]: 1 : if (nspinor==2)&
2038 : 0 : & isize=isize-40*mpw*nspinor !pwave_so,subghg_so
2039 : :
2040 : 1 : isize=isize+56*mpw*nspinor !cwork,rwork
2041 : 1 : memsize=max(memsize,isize)
2042 : 1 : isize=isize-56*mpw*nspinor !cwork,rwork
2043 : 1 : isize=isize+112*mpw*nspinor !cwork,rwork,iwork,ifail
2044 : 1 : memsize=max(memsize,isize)
2045 : 1 : isize=isize-112*mpw*nspinor !cwork,rwork,iwork,ifail
2046 : 1 : isize=isize-8*mpw*nspinor*(mpw*nspinor+1) !ghg
2047 : 1 : isize=isize+8*mband !occ_k
2048 : 1 : memsize=max(memsize,isize)
2049 : : isize=isize-8*mband !occ_k
2050 : : isize=isize-8*mpw*nspinor&
2051 : 1 : & -16*mpw*nspinor*(mpw*nspinor+1) !eigval,eigvec
2052 : : isize=isize-32*mpw-8*ngfft(4)&
2053 : 1 : & *ngfft(5)*ngfft(6) !gbig,ts,vlocal
2054 : : if(xmpi_paral==1)then
2055 : : isize=isize-4*4 !kpt_distrb
2056 : : end if
2057 : 1 : isize=isize-580-fnlen-4*(81+nkpt+9*nsym)-8*15 !non allocatable var.
2058 : : !
2059 : 1 : write(msg,'(2a,f8.2,a)')ch10,&
2060 : 2 : & ' Additional amount of memory required by "outkss" routine=',memsize*b2Mb,' Mbytes.'
2061 : 1 : call wrtout(std_out,msg,'COLL')
2062 : : !
2063 : 1 : end subroutine memkss
2064 : : !!***
2065 : :
2066 : : !!****f* ABINIT/dsksta
2067 : : !! NAME
2068 : : !! dsksta
2069 : : !!
2070 : : !! FUNCTION
2071 : : !! This routine evaluates the amount of disk space required by the _KSS file.
2072 : : !!
2073 : : !! INPUTS
2074 : : !! dimlmn(natom*usepaw)=Number of nlm partial waves for each atom.
2075 : : !! ishm=Number of G-shells to be saved in _KSS file.
2076 : : !! mpsang=Max angular momentum +1 for pseudos.
2077 : : !! natom=Number of atoms in the unit cell.
2078 : : !! nbandkss=Number of desired bands to be saved in _KSS file
2079 : : !! nkpt=Number of k points.
2080 : : !! npwkss=Number of desired G-vectors to be saved in _KSS file.
2081 : : !! nspinor=Number of spinorial components.
2082 : : !! nsppol=Number of independent spin polarizations.
2083 : : !! ntypat=Number of type of atoms.
2084 : : !! nsym2=Number of symmetries in space group, without INV
2085 : : !! usepaw=1 if PAW.
2086 : : !!
2087 : : !! OUTPUT
2088 : : !! Writes on standard output
2089 : : !!
2090 : : !! SOURCE
2091 : :
2092 : 6 : subroutine dsksta(ishm,usepaw,nbandkss,mpsang,natom,ntypat,npwkss,nkpt,nspinor,nsppol,nsym2,dimlmn)
2093 : :
2094 : : !Arguments ------------------------------------
2095 : : !scalars
2096 : : integer,intent(in) :: usepaw,ishm,nbandkss,mpsang,natom,ntypat,nkpt
2097 : : integer,intent(in) :: npwkss,nspinor,nsppol,nsym2
2098 : : !arrays
2099 : : integer,intent(in) :: dimlmn(natom*usepaw)
2100 : :
2101 : : !Local variables-------------------------------
2102 : : !scalars
2103 : : integer :: bsize_tot,bsize_hdr,bsize_kb,bsize_wf,bsize_cprj
2104 : : character(len=500) :: msg
2105 : :
2106 : : ! *********************************************************************
2107 : :
2108 : : !The Abinit header is not considered.
2109 : : bsize_hdr= 80*2 + & !title
2110 : : &5*4 + & !nsym2,nbandksseff,npwkss,ishm,mpsang
2111 : : &nsym2*9*4 + & !symrel2
2112 : : &nsym2*3*8 + & !tnons
2113 : : &npwkss*3*4 + & !gbig
2114 : 6 : &ishm*4 !shlim
2115 : :
2116 : : !NOTE: vkb does not depend on nsppol, however the elements are written for each spin.
2117 : 6 : bsize_kb=0
2118 [ + + ]: 6 : if (usepaw==0) then
2119 : : bsize_kb= nsppol* &
2120 : : & ( mpsang*ntypat *8 + & !vkbsign
2121 : : & 2*(nkpt*mpsang*ntypat*npwkss*8) & !vkbd,vkbd
2122 : 2 : & )
2123 : : end if
2124 : :
2125 : : bsize_wf= nsppol* &
2126 : : & ( nkpt*nbandkss *8 + & !energies
2127 : : &nkpt*nbandkss*nspinor*npwkss*2*8 & !wfg
2128 : 6 : &)
2129 : :
2130 : : !For PAW add space required by projectors.
2131 : 6 : bsize_cprj=0
2132 [ + + ]: 6 : if (usepaw==1) then
2133 [ + + ]: 8 : bsize_cprj=SUM(dimlmn(:))*(nsppol*nkpt*nspinor*nbandkss*2*8)
2134 : : end if
2135 : :
2136 : 6 : bsize_tot = bsize_hdr + bsize_kb + bsize_wf + bsize_cprj
2137 : 6 : write(msg,'(2a,f8.2,4a,4(a,f8.2,2a))')ch10,&
2138 : 6 : & ' Total amount of disk space required by _KSS file = ',bsize_tot*b2Mb,' Mb.',ch10,&
2139 : 6 : & ' Subdivided into : ',ch10,&
2140 : 6 : & ' Header = ',bsize_hdr *b2Mb,' Mb.',ch10,&
2141 : 6 : & ' KB elements = ',bsize_kb *b2Mb,' Mb.',ch10,&
2142 : 6 : & ' Wavefunctions (PW) = ',bsize_wf *b2Mb,' Mb.',ch10,&
2143 : 12 : & ' PAW projectors = ',bsize_cprj*b2Mb,' Mb.',ch10
2144 : 6 : call wrtout(std_out,msg,'COLL')
2145 : :
2146 : 6 : end subroutine dsksta
2147 : : !!***
2148 : :
2149 : : end subroutine outkss
2150 : : !!***
2151 : :
2152 : : END MODULE m_io_kss
2153 : : !!***
|
