1#include "f_defs.h"
2
3!-----------------------------------------------------------------------
4subroutine input_kernel(crys,gvec,kg,kgq,kp,syms,xct,flagbz, &
5  intwfnv,intwfnc)
6!-----------------------------------------------------------------------
7!
8!     Read parameters from file WFN_co
9!     Initialize k-points sampling, kg type
10!     Initialize G-space, gvec
11!
12!     input: xct type
13!
14!     output: crys,gvec,syms,kg types
15!             INT_VWFN_* and INT_CWFN_* files
16!
17  use global_m
18  use checkbz_m
19  use fullbz_m
20  use input_utils_m
21  use misc_m
22  use sort_m
23  use wfn_rho_vxc_io_m
24  use read_rho_vxc_m
25  use epsread_hdf5_m
26#ifdef HDF5
27  use hdf5
28#endif
29  implicit none
30
31  type (crystal), intent(out) :: crys
32  type (gspace), intent(out) :: gvec
33  type (grid), intent(out) :: kg,kgq
34  type (kpoints), intent(out) :: kp
35  type (symmetry), intent(out) :: syms
36  type (xctinfo), intent(inout) :: xct
37  integer, intent(in) :: flagbz
38  type (int_wavefunction), intent(out) :: intwfnc,intwfnv
39
40
41  type (wavefunction) :: wfnc,wfnv
42  character :: tmpfn*16
43  integer :: iwritev,iwritec,iwritek
44  integer :: ickmem,irk
45  integer :: ii,jj,is,isp,ig,ikq,ik,umk
46  integer :: irks,ivband,icband
47
48  real(DP) :: diffvol,vcell,kt(3),div,tol,delta,qq_temp(3)
49
50  real(DP), allocatable :: ek_tmp(:)
51  integer, allocatable :: index(:),indxk(:),indxkq(:),k_tmp(:,:)
52  integer, allocatable :: isrti(:)
53  SCALAR, allocatable :: cg(:,:), cgarray(:)
54
55  character(len=3) :: sheader
56  integer :: iflavor
57  type(gspace) :: gvec_kpt
58  logical :: skip_checkbz
59
60
61  PUSH_SUB(input_kernel)
62
63  call logit('input_kernel:  reading WFN_co')
64  if (peinf%inode == 0) call open_file(25,file='WFN_co',form='unformatted',status='old')
65
66  sheader = 'WFN'
67  iflavor = 0
68  call read_binary_header_type(25, sheader, iflavor, kp, gvec, syms, crys, dont_warn_kgrid=xct%patched_sampling_co)
69  call check_trunc_kpts(xct%icutv, kp)
70
71  if(any(kp%shift > TOL_Zero) .and. peinf%inode == 0) then
72    write(0,'(a)') "WARNING: WFN_co has a shift. This is not recommended."
73  endif
74
75  call logit('input_kernel:  reading gvec info')
76  SAFE_ALLOCATE(gvec%components, (3, gvec%ng))
77  call read_binary_gvectors(25, gvec%ng, gvec%ng, gvec%components)
78
79  call get_volume(vcell,crys%bdot)
80  diffvol=abs(crys%celvol-vcell)
81  if (diffvol.gt.TOL_Small) then
82    call die('volume mismatch', only_root_writes = .true.)
83  endif
84
85  ! there is no fine grid to set the Fermi level in kernel
86  call find_efermi(xct%rfermi, xct%efermi, xct%efermi_input, kp, kp%mnband, 1, &
87    "coarse grid", should_search = .true., should_update = .true., write7 = .false.)
88
89  call assess_degeneracies(kp, kp%el(kp%mnband, :, :), kp%mnband - 1, xct%efermi, TOL_Degeneracy)
90
91  xct%nspin=kp%nspin
92  ! consistency check on status of spinor...
93  if (xct%nspinor .eq. 2 .and. kp%nspinor .ne. 2) &
94    call die("Keyword SPINOR used in absorption.inp but WFN is not of spinor type", only_root_writes = .true.)
95
96  if(any(kp%ifmax(:,:) == 0)) &
97    call die("BSE codes cannot handle a system where some k-points have no occupied bands.", only_root_writes = .true.)
98
99  kp%nvband=minval(kp%ifmax(:,:)-kp%ifmin(:,:))+1
100  kp%ncband=kp%mnband-maxval(kp%ifmax(:,:))
101
102!----------------------------------------------------------------
103! (gsm) check whether the requested number of bands
104!       is available in the wavefunction file
105
106  if(xct%nvb_co.gt.kp%nvband) then
107    call die("The requested number of valence bands is not available in WFN_co.", &
108             only_root_writes=.true.)
109  endif
110  if(xct%ncb_co.gt.kp%ncband) then
111    call die("The requested number of conduction bands is not available in WFN_co.", &
112             only_root_writes=.true.)
113  endif
114
115! DAS: degenerate subspace check
116
117  if (peinf%inode.eq.0) then
118    if(xct%ncb_co.eq.kp%ncband) then
119      call die("You must provide one more conduction band in WFN_co in order to assess degeneracy.", &
120               only_root_writes=.true.)
121    endif
122    do jj = 1, kp%nspin
123      do ii = 1, kp%nrk
124        if(kp%ifmax(ii, jj) - xct%nvb_co > 0) then
125          ! no need to compare against band 0 if all valence are included
126          if(abs(kp%el(kp%ifmax(ii, jj) - xct%nvb_co + 1, ii, jj) &
127            - kp%el(kp%ifmax(ii, jj) - xct%nvb_co, ii, jj)) .lt. TOL_Degeneracy) then
128            if(xct%degeneracy_check_override) then
129              write(0,'(a)') &
130                "WARNING: Selected number of valence bands breaks degenerate subspace in WFN_co. " // &
131                "Run degeneracy_check.x for allowable numbers."
132              write(0,*)
133            else
134              write(0,'(a)') &
135                "Run degeneracy_check.x for allowable numbers, or use keyword " // &
136                "degeneracy_check_override to run anyway (at your peril!)."
137              call die("Selected number of valence bands breaks degenerate subspace in WFN_co.", &
138                       only_root_writes=.true.)
139            endif
140          endif
141        endif
142        if(abs(kp%el(kp%ifmax(ii, jj) + xct%ncb_co, ii, jj) &
143          - kp%el(kp%ifmax(ii, jj) + xct%ncb_co + 1, ii, jj)) .lt. TOL_Degeneracy) then
144          if(xct%degeneracy_check_override) then
145            write(0,'(a)') &
146              "WARNING: Selected number of conduction bands breaks degenerate subspace in WFN_co. " // &
147              "Run degeneracy_check.x for allowable numbers."
148            write(0,*)
149          else
150            write(0,'(a)') &
151              "Run degeneracy_check.x for allowable numbers, or use keyword " // &
152              "degeneracy_check_override to run anyway (at your peril!)."
153            call die("Selected number of conduction bands breaks degenerate subspace in WFN_co.",&
154                     only_root_writes=.true.)
155          endif
156        endif
157      enddo
158    enddo
159  endif
160
161!-----------------------------------------------------------------------
162!     Read the k-point sampling from kpoints (if it exists) or from
163!     WFN_co
164
165  if (xct%read_kpoints) then
166    if (peinf%inode.eq.0) then
167      call open_file(9,file='kpoints_co',form='formatted',status='old')
168      read(9,*) kg%nr
169      SAFE_ALLOCATE(kg%r, (3,kg%nr))
170      do ii=1,kg%nr
171        read(9,*) (kg%r(jj,ii),jj=1,3),div
172        kg%r(:,ii) = kg%r(:,ii)/div
173      enddo
174      call close_file(9)
175    endif
176#ifdef MPI
177    call MPI_BCAST(kg%nr,   1,     MPI_INTEGER,0,MPI_COMM_WORLD,mpierr)
178    if(peinf%inode.ne.0) then
179      SAFE_ALLOCATE(kg%r, (3,kg%nr))
180    endif
181    call MPI_BCAST(kg%r,    3*kg%nr,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr)
182#endif
183
184!----------------------------------------------------------------
185!     indxk : stores the correspondence between k-points kg%r and kp%rk
186!     (it is used to select the set of wavefunctions to be stored)
187!     tol : tolerance in the coordinates of k-points
188
189    tol = 1.d-4
190    SAFE_ALLOCATE(indxk, (kg%nr))
191    indxk=0
192    do jj=1,kg%nr
193      do ii=1,kp%nrk
194        kt(:) = kg%r(:,jj) - kp%rk(:,ii)
195        if ((abs(kt(1)).lt.tol).and.(abs(kt(2)).lt.tol) &
196          .and.(abs(kt(3)).lt.tol)) then
197          if (indxk(jj).ne.0) write(0,*) 'WARNING: multiple definition of k-point',jj,indxk(jj),kg%r(:,jj)
198          indxk(jj)=ii
199        endif
200      enddo
201      if (indxk(jj).eq.0) write(0,*) 'WARNING: could not find vector ',kg%r(:,jj),' in WFN_co'
202!
203!     no need to stop here; if indxk.eq.0, the job will stop in genwf
204!
205    enddo
206  else
207    kg%nr=kp%nrk
208    SAFE_ALLOCATE(kg%r, (3,kg%nr))
209    kg%r(1:3,1:kg%nr)=kp%rk(1:3,1:kp%nrk)
210    SAFE_ALLOCATE(indxk, (kg%nr))
211    do ii=1,kg%nr
212      indxk(ii) = ii
213    enddo
214  endif
215
216!-----------------------------------------------------------------------
217!     Order g-vectors with respect to their kinetic energy
218!
219
220  call logit('input_kernel:  reordering gvecs')
221
222  ! FHJ: Figure out ecute and ecutg
223  call get_ecut()
224
225  SAFE_ALLOCATE(index, (gvec%ng))
226  SAFE_ALLOCATE(gvec%ekin, (gvec%ng))
227  call kinetic_energies(gvec, crys%bdot, gvec%ekin)
228  call sortrx(gvec%ng, gvec%ekin, index, gvec = gvec%components)
229
230  xct%ng   = gcutoff(gvec%ng, gvec%ekin, index, xct%ecutg)
231  if (xct%theory.eq.0) &
232  xct%neps = gcutoff(gvec%ng, gvec%ekin, index, xct%ecute) ! cuteness energy??
233
234  if(xct%ecute > xct%ecutg .and. xct%theory.eq.0) then
235    write(0,*) 'ecute = ', xct%ecute, ' ecutg = ', xct%ecutg
236    call die("The screened_coulomb_cutoff cannot be greater than the bare_coulomb_cutoff.", only_root_writes = .true.)
237  endif
238
239  SAFE_ALLOCATE(ek_tmp, (gvec%ng))
240  ek_tmp = gvec%ekin
241  SAFE_ALLOCATE(k_tmp, (3,gvec%ng))
242  k_tmp = gvec%components
243  do ii=1,gvec%ng
244    gvec%ekin(ii) = ek_tmp(index(ii))
245    gvec%components(:,ii) = k_tmp(:,index(ii))
246  enddo
247
248  call gvec_index(gvec)
249
250  ! If we are not doing just purely TDHF then
251  ! read the charge density/fxc
252  if ((xct%theory == 1) .and. ((1.0d0 - xct%coulomb_mod%long_range_frac_fock > TOL_SMALL) .or. &
253     (1.0d0 - xct%coulomb_mod%short_range_frac_fock > TOL_SMALL))) then
254    xct%coul_mod_flag=.true.
255    SAFE_ALLOCATE(isrti, (gvec%ng))
256    do ii=1,gvec%ng
257      isrti(index(ii)) = ii
258    enddo
259    call read_rho(xct%wpg, gvec, kp, syms, crys, isrti, index, 'WFN_co')
260    SAFE_DEALLOCATE(isrti)
261  endif
262
263  SAFE_DEALLOCATE(index)
264  SAFE_DEALLOCATE(ek_tmp)
265  SAFE_DEALLOCATE(k_tmp)
266
267
268!-----------------------------------------------------------------------
269!     Generate full brillouin zone from irreducible wedge, rk -> fk
270!
271!     If flagbz.eq.1, only Identity will be used as
272!     symmetry operation. In this case, kg%r (irreducible BZ) and kg%f
273!     (full BZ) will be identical.
274!
275  if (flagbz.eq.0.and.peinf%inode.eq.0) write(6,801)
276  if (flagbz.eq.1.and.peinf%inode.eq.0) write(6,802)
277801 format(1x,'Using symmetries to expand the coarse grid sampling')
278802 format(1x,'No symmetries used in the coarse grid sampling')
279!
280  call timacc(7,1)
281  if (flagbz.eq.1) then
282    call fullbz(crys,syms,kg,1,skip_checkbz,wigner_seitz=.true.,paranoid=.true.)
283  else
284    call fullbz(crys,syms,kg,syms%ntran,skip_checkbz,wigner_seitz=.true.,paranoid=.true.)
285  endif
286  tmpfn='WFN_co'
287  if (.not. skip_checkbz.and..not.xct%patched_sampling_co) then
288    call checkbz(kg%nf,kg%f,kp%kgrid,kp%shift,crys%bdot, &
289      tmpfn,'k',.true.,xct%freplacebz,xct%fwritebz)
290  endif
291  call timacc(7,2)
292  xct%nkpt_co=kg%nf
293
294  if (peinf%verb_high .and. peinf%inode==0) then
295    write(6,'(/1x,a6,14x,a7,12x,2(1x,a6),3x,a3)') 'i', 'k-point', 'indr', 'itran', 'kg0'
296    write(6,'(1x,6("-"),1x,32("-"),2(1x,6("-")),1x,8("-"))')
297    do ii=1,kg%nf
298      write(6,'(1x,i6,3(1x,f10.6),2(1x,i6),3(1x,i2))') &
299        ii, kg%f(:,ii), kg%indr(ii), kg%itran(ii), kg%kg0(:,ii)
300    enddo
301  endif
302
303!------------------------------------------------------------------------
304! If there is a finite center-of-mass momentum, Q, find mapping between k
305! and k+Q
306
307  SAFE_ALLOCATE(xct%indexq,(kg%nf))
308  if (xct%qflag.eq.1) then
309    do ik=1,kg%nf
310      xct%indexq(ik) = ik
311    enddo
312  endif
313
314  SAFE_ALLOCATE(kgq%f,(3,kg%nf))
315  SAFE_ALLOCATE(kgq%kg0,(3,kg%nf))
316  SAFE_ALLOCATE(kgq%indr,(kg%nf))
317
318!------------------------------------------------
319! Distribute vcks-quadruplets among the PEs
320
321  call logit('input_kernel:  calling distrib_kernel')
322  if (xct%qflag.eq.1) then
323    call distrib_kernel(xct,kp%ngkmax,kg,kg,gvec)
324  endif
325  if (peinf%verb_debug .and. peinf%inode==0) then
326    write(6,'(/1x,a,3(1x,i0)/)') "Allocating isort", peinf%iownwfv(peinf%inode+1), peinf%iownwfc(peinf%inode+1), gvec%ng
327  endif
328  if (xct%qflag.eq.1) then
329    SAFE_ALLOCATE(intwfnv%cg, (kp%ngkmax,peinf%iownwfv(peinf%inode+1),kp%nspin*kp%nspinor))
330  endif
331  SAFE_ALLOCATE(intwfnc%cg, (kp%ngkmax,peinf%iownwfc(peinf%inode+1),kp%nspin*kp%nspinor))
332  if (xct%qflag.eq.1) then
333    SAFE_ALLOCATE(intwfnv%isort, (gvec%ng,peinf%iownwfk(peinf%inode+1)))
334    SAFE_ALLOCATE(intwfnv%ng, (peinf%iownwfk(peinf%inode+1)))
335  endif
336  SAFE_ALLOCATE(intwfnc%isort, (gvec%ng,peinf%iownwfk(peinf%inode+1)))
337  SAFE_ALLOCATE(intwfnc%ng, (peinf%iownwfk(peinf%inode+1)))
338  intwfnv%nspin=kp%nspin
339  intwfnv%nspinor=kp%nspinor
340  intwfnc%nspin=kp%nspin
341  intwfnc%nspinor=kp%nspinor
342
343#ifdef DEBUG
344  if (peinf%inode.eq.0) then
345    write(6,*) 'kp%ngkmax',kp%ngkmax
346  endif
347#endif
348
349!------------------------------------------------
350! Begin loop that distributes wave functions
351
352  SAFE_ALLOCATE(wfnv%isort, (gvec%ng))
353  SAFE_ALLOCATE(wfnc%isort, (gvec%ng))
354  wfnv%nspin=kp%nspin
355  wfnv%nspinor=kp%nspinor
356  wfnc%nspin=kp%nspin
357  wfnc%nspinor=kp%nspinor
358
359  do irk=1,kp%nrk
360    irks = 0
361    do ii=1,kg%nr
362      if (irk == indxk(ii)) then
363        irks=ii
364        exit
365      endif
366    enddo
367
368    SAFE_ALLOCATE(gvec_kpt%components, (3, kp%ngk(irk)))
369    call read_binary_gvectors(25, kp%ngk(irk), kp%ngk(irk), gvec_kpt%components)
370
371    SAFE_ALLOCATE(cg, (kp%ngk(irk),kp%nspin*kp%nspinor))
372    if(irks > 0) then
373      do ii = 1, kp%ngk(irk)
374        call findvector(wfnv%isort(ii), gvec_kpt%components(:, ii), gvec)
375        if(wfnv%isort(ii) == 0) call die('could not find gvec', only_root_writes=.true.)
376      enddo
377
378      wfnv%ng=kp%ngk(irk)
379      wfnc%ng=kp%ngk(irk)
380      SAFE_ALLOCATE(cgarray, (kp%ngk(irk)))
381      wfnc%isort(1:gvec%ng)=wfnv%isort(1:gvec%ng)
382      ickmem=0
383    endif
384
385!        write(6,*) peinf%inode, 'loop wfnup', irk
386
387! Loop Over Bands
388
389    do ii=1,kp%mnband
390
391      call read_binary_data(25, kp%ngk(irk), kp%ngk(irk), kp%nspin*kp%nspinor, cg,bcast=.false.)
392
393! If we do not need this band, skip it...
394      if(irks == 0) cycle
395!
396      do is=1, kp%nspin
397        if (ii .ge. (kp%ifmax(irk,is)-xct%nvb_co+1) .and. ii .le. (kp%ifmax(irk,is)+xct%ncb_co)) then
398          do isp=1, kp%nspinor
399            if (peinf%inode.eq.0) then
400              ! Check normalization of this band
401              call checknorm('WFN_co',ii,irks,kp%ngk(irk),kp%nspin,kp%nspinor,cg(:,:))
402              do ig=1, kp%ngk(irk)
403                cgarray(ig)=cg(ig, is*isp)
404              end do
405            end if
406!           write(6,*) peinf%inode, 'Read', irk, ii
407
408!-----------------------
409! If ii is one of the selected valence band...
410          if((ii.le.kp%ifmax(irk,is)).and. &
411            (ii.ge.(kp%ifmax(irk,is)-xct%nvb_co+1))) then
412
413            ! Skip current valence band if valence bands will be read from WFNq_co
414            if (xct%qflag.eq.0) cycle
415
416!            write(6,*) peinf%inode, 'in val',irk,ii
417
418#ifdef MPI
419              call MPI_BCAST(cgarray,kp%ngk(irk),MPI_SCALAR,0, &
420                MPI_COMM_WORLD,mpierr)
421#endif
422
423            ivband=kp%ifmax(irk,is)-ii+1
424            iwritev=peinf%ipev(peinf%inode+1,ivband,irk)
425            if (xct%qflag.eq.1) then
426              iwritek=peinf%ipek(peinf%inode+1,irk)
427            else if (xct%qflag.eq.2) then
428              iwritek=peinf%ipekq(peinf%inode+1,irk)
429            endif
430
431!             write(6,*) peinf%inode, 'bcast val',irk,ii,ivband,
432!     >       iwritev
433              if(iwritev.ne.0) then
434                intwfnv%cg(1:wfnv%ng,iwritev,is*isp)=cgarray(1:kp%ngk(irk))
435              endif
436              if(iwritek.ne.0) then
437                intwfnv%isort(1:gvec%ng,iwritek)=wfnv%isort(1:gvec%ng)
438                intwfnv%ng(iwritek)=kp%ngk(irk)
439              endif
440
441!             write(*,'(a, 4i4)') 'valence', peinf%inode+1, ivband, irk, iwritev
442!             write(*, '(a, 6i5)') "valence", ivband, ii, is, peinf%inode, irk, iwritev
443!             write(6,*) peinf%inode, 'wrote val',irk,ii
444
445            endif !ii is one of the selected valence band
446
447!-----------------------
448! If ii is one of the selected conduction band...
449
450            if((ii.ge.(kp%ifmax(irk,is)+1)) &
451              .and.(ii.le.(kp%ifmax(irk,is)+xct%ncb_co))) then
452
453#ifdef MPI
454              call MPI_BCAST(cgarray,kp%ngk(irk),MPI_SCALAR,0, &
455                MPI_COMM_WORLD,mpierr)
456#endif
457
458              icband=ii-kp%ifmax(irk,is)
459              iwritec=peinf%ipec(peinf%inode+1,icband,irk)
460              iwritek=peinf%ipek(peinf%inode+1,irk)
461!             write(6,*) 'icband',ii,icband,iwritec
462              if(iwritec.ne.0) then
463                intwfnc%cg(1:wfnc%ng,iwritec,is*isp)=cgarray(1:kp%ngk(irk))
464              endif
465              if(iwritek.ne.0) then
466                intwfnc%isort(1:gvec%ng,iwritek)=wfnc%isort(1:gvec%ng)
467                intwfnc%ng(iwritek)=kp%ngk(irk)
468              endif
469
470!             write(*,'(a, 4i4)') 'conduction', peinf%inode+1, icband, irk, iwritec
471!             write(*, '(a, 6i5)') "conduction", icband, ii, is, peinf%inode, irk, iwritec
472
473            endif ! ii is one of the selected conduction bands
474          end do ! isp
475        end if
476      end do ! is
477      if (ii>maxval(kp%ifmax)+xct%ncb_co .and. irk==kp%nrk) then
478        exit
479      endif
480    enddo ! ii (loop on bands)
481
482    SAFE_DEALLOCATE(cg)
483    SAFE_DEALLOCATE(cgarray)
484
485
486  enddo !end loop over k-points
487
488! Write out info about xtal
489
490  if (peinf%inode.eq.0) then
491    write(6,'(/1x,a)') 'Crystal wavefunctions read from file WFN_co:'
492    write(6,'(1x,a,i0)') '- Number of k-points in WFN_co: ', kp%nrk
493    write(6,'(1x,a,i0)') '- Number of k-points in the full BZ of WFN_co: ', kg%nf
494    if (peinf%verb_high) then
495      write(6,'(1x,a)') '- K-points:'
496      write(6,'(1(2x,3(1x,f10.6)))') kg%r(1:3,1:kg%nr)
497    endif
498    call close_file(25)
499  endif ! node 0
500
501
502  SAFE_DEALLOCATE_P(wfnv%isort)
503  SAFE_DEALLOCATE_P(wfnc%isort)
504  SAFE_DEALLOCATE_P(gvec_kpt%components)
505  SAFE_DEALLOCATE(indxk)
506  SAFE_DEALLOCATE(indxkq)
507
508  POP_SUB(input_kernel)
509
510  return
511
512contains
513
514  ! FHJ: Figure out xct%ecute and xct%ecutg, if the user didn`t specify them.
515  subroutine get_ecut()
516#ifdef HDF5
517    integer :: ng, nq, nfreq, nfreq_imag, nmtx_max, qgrid(3), freq_dep
518#endif
519    real(DP) :: ecuts
520    !real(DP) :: raux
521
522    PUSH_SUB(input_kernel.get_ecut)
523
524    if (xct%theory/=0) then
525      xct%ecute = 0
526      if (xct%ecutg < TOL_ZERO) xct%ecutg = kp%ecutwfc
527      return
528    endif
529    if (peinf%inode==0) then
530#ifdef HDF5
531      if(xct%use_hdf5) then
532        call read_eps_grid_sizes_hdf5(ng, nq, ecuts, nfreq, &
533          nfreq_imag, nmtx_max, qgrid, freq_dep, 'eps0mat.h5')
534      else
535#endif
536        call open_file(10, file='eps0mat', form='unformatted', status='old')
537        read(10)
538        read(10)
539        read(10)
540        read(10)
541        read(10)
542        read(10)
543        read(10) ecuts
544        call close_file(10)
545      endif
546#ifdef HDF5
547    endif
548#endif
549#ifdef MPI
550    call MPI_Bcast(ecuts, 1, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr)
551#endif
552    if (xct%ecute < TOL_ZERO) xct%ecute = ecuts
553    !raux = maxval(crys%bdot)
554    !if (xct%ecutg < TOL_ZERO) xct%ecutg = xct%ecute + raux
555    if (xct%ecutg < TOL_ZERO) xct%ecutg = kp%ecutwfc
556
557    POP_SUB(input_kernel.get_ecut)
558
559  end subroutine get_ecut
560
561end subroutine input_kernel
562