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 357 : 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 357 : integer,pointer :: gbig(:,:)
781 : character(len=500) :: msg
782 : !arrays
783 357 : integer,pointer :: symrel2(:,:,:),shlim(:)
784 357 : real(dp),pointer :: tnons2(:,:)
785 : ! *********************************************************************
786 :
787 357 : ierr = 0
788 357 : write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.'
789 357 : 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 357 : 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 342 : 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 342 : nsym2=nsym ; pinv=1
812 1026 : ABI_MALLOC(symrel2,(3,3,nsym))
813 1026 : ABI_MALLOC(tnons2,(3,nsym))
814 185098 : symrel2(:,:,:)=symrel(:,:,1:nsym)
815 57190 : 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 357 : nullify(gbig,shlim)
830 357 : call merge_and_sort_kg(nkpt,kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim)
831 :
832 357 : nbase = SIZE(shlim) ! Number of independent G in the big sphere.
833 357 : 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 357 : if ((npwkss==0).or.(npwkss>=maxpw)) then
838 357 : npwkss=maxpw
839 : write(msg,'(5a)')&
840 357 : & ' Since the number of g''s to be written on file',ch10,&
841 357 : & ' was 0 or too large, it has been set to the max. value.,',ch10,&
842 714 : & ' computed from the union of the sets of G vectors for the different k-points.'
843 357 : call wrtout(std_out,msg,'COLL')
844 : end if
845 :
846 357 : ishm=0
847 20668 : do ii=1,nbase
848 20668 : if (shlim(ii)<=npwkss) then
849 20311 : ishm=ii
850 : else
851 : EXIT
852 : end if
853 : end do
854 : !ishm=bisect(shlim,npwkss)
855 :
856 357 : 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 357 : write(msg,'(a,i5)')' Number of G-vectors is: ',npwkss
869 357 : call wrtout(std_out,msg,'COLL')
870 :
871 1071 : ABI_MALLOC(gvec_kss,(3,npwkss))
872 2069597 : gvec_kss = gbig(:,1:npwkss)
873 :
874 357 : ABI_FREE(gbig)
875 357 : ABI_FREE(symrel2)
876 357 : ABI_FREE(tnons2)
877 357 : ABI_FREE(shlim)
878 :
879 357 : 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 : !!***
|
