TABLE OF CONTENTS


ABINIT/prteigrs [ Functions ]

[ Top ] [ Functions ]

NAME

 prteigrs

FUNCTION

 Print out eigenvalues band by band and k point by k point.
 If option=1, do it in a standard way, for self-consistent calculations.
 If option=2, print out residuals and eigenvalues, in a format
 adapted for nonself-consistent calculations, within the loops.
 If option=3, print out eigenvalues, in a format
 adapted for nonself-consistent calculations, at the end of the job.
 If option=4, print out derivatives of eigenvalues (same format as option==3, except header that is printed)
 If option=5, print out Fan contribution to zero-point motion correction to eigenvalues (averaged)
                  (same format as option==3, except header that is printed)
 If option=6, print out DDW contribution to zero-point motion correction to eigenvalues (averaged)
                  (same format as option==3, except header that is printed)
 If option=7, print out Fan+DDW contribution to zero-point motion correction to eigenvalues (averaged)
                  (same format as option==3, except header that is printed)

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  eigen(mband*nkpt*nsppol)=eigenvalues (hartree)
   or, if option==4, diagonal of derivative of eigenvalues
   or, if option==5...7, zero-point motion correction to eigenvalues (averaged)
  enunit=choice parameter: 0=>output in hartree; 1=>output in eV;
   2=> output in both hartree and eV
  fermie=fermi energy (Hartree)
  fname_eig=filename for printing of the eigenenergies
  iout=unit number for formatted output file
  iscf=option for self-consistency
  kptns(3,nkpt)=k points in reduced coordinates
  kptopt=option for the generation of k points
  mband=maximum number of bands
  nband(nkpt)=number of bands at each k point
  nkpt=number of k points
  nnsclo_now=number of non-self-consistent loops for the current vtrial
    (often 1 for SCF calculation, =nstep for non-SCF calculations)
  nsppol=1 for unpolarized, 2 for spin-polarized
  occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point
  occopt=option for occupancies
  option= (see above)
  prteig=control print eigenenergies
  prtvol=control print volume and debugging
  resid(mband*nkpt*nsppol)=residuals (hartree**2)
  tolwfr=tolerance on band residual of wf, hartrees**2 (needed when option=2)
  vxcavg=average of vxc potential
  wtk(nkpt)=k-point weights

OUTPUT

  (only writing)

PARENTS

      clnup1,dfpt_looppert,respfn,scprqt,vtorho

CHILDREN

      wrtout

SOURCE

 68 #if defined HAVE_CONFIG_H
 69 #include "config.h"
 70 #endif
 71 
 72 #include "abi_common.h"
 73 
 74 
 75 subroutine prteigrs(eigen,enunit,fermie,fname_eig,iout,iscf,kptns,kptopt,mband,nband,&
 76 &  nkpt,nnsclo_now,nsppol,occ,occopt,option,prteig,prtvol,resid,tolwfr,vxcavg,wtk)
 77 
 78  use defs_basis
 79  use m_profiling_abi
 80  use m_errors
 81 
 82  use m_io_tools,  only : open_file
 83 
 84 !This section has been created automatically by the script Abilint (TD).
 85 !Do not modify the following lines by hand.
 86 #undef ABI_FUNC
 87 #define ABI_FUNC 'prteigrs'
 88  use interfaces_14_hidewrite
 89 !End of the abilint section
 90 
 91  implicit none
 92 
 93 !Arguments ------------------------------------
 94 !scalars
 95  integer,intent(in) :: enunit,iout,iscf,kptopt,mband,nkpt,nnsclo_now,nsppol
 96  integer,intent(in) :: occopt,option,prteig,prtvol
 97  real(dp),intent(in) :: fermie,tolwfr,vxcavg
 98  character(len=*),intent(in) :: fname_eig
 99 !arrays
100  integer,intent(in) :: nband(nkpt*nsppol)
101  real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kptns(3,nkpt)
102  real(dp),intent(in) :: occ(mband*nkpt*nsppol),resid(mband*nkpt*nsppol)
103  real(dp),intent(in) :: wtk(nkpt)
104 
105 !Local variables-------------------------------
106 !scalars
107  integer,parameter :: nkpt_max=50
108  integer :: band_index,iband,ienunit,ii,ikpt,isppol,nband_index,nband_k,nkpt_eff,tmagnet,tmetal,temp_unit
109  real(dp) :: convrt,magnet,residk,rhodn,rhoup
110  character(len=2) :: ibnd_fmt,ikpt_fmt
111  character(len=7) :: strunit1,strunit2
112  character(len=39) :: kind_of_output
113  character(len=500) :: message
114 
115 ! *************************************************************************
116 
117  if (nsppol<1.or.nsppol>2) then
118    write(message, '(a,i0)' )' nsppol must be 1 or 2. Argument was ',nsppol
119    MSG_BUG(message)
120  end if
121 
122  if (enunit<0.or.enunit>2) then
123    write(message, '(a,i0)' )' enunit must be 0, 1 or 2. Argument was ',enunit
124    MSG_BUG(message)
125  end if
126 
127  if (prteig > 0) then
128    write(message, '(a,a)' ) ' prteigrs : about to open file ',TRIM(fname_eig)
129    call wrtout(iout,message,'COLL')
130    if (open_file(fname_eig, message, newunit=temp_unit, status='unknown', form='formatted') /= 0) then
131      MSG_ERROR(message)
132    end if
133    rewind(temp_unit) ! always rewind disk file and print latest eigenvalues
134  end if
135 
136  kind_of_output=              ' Eigenvalues                          '
137  if(option==4) kind_of_output=' Expectation of eigenvalue derivatives'
138  if(option==5) kind_of_output=' Fan corrections to eigenvalues at T=0'
139  if(option==6) kind_of_output=' DDW corrections to eigenvalues at T=0'
140  if(option==7) kind_of_output=' Fan+DDW corrs   to eigenvalues at T=0'
141 
142  nkpt_eff=nkpt
143 
144 !DEBUG
145 !write(message,'(a,5i5)')' prtvol,iscf,kptopt,nkpt_eff,nkpt_max ',prtvol,iscf,kptopt,nkpt_eff,nkpt_max
146 !call wrtout(iout,message,'COLL')
147 !ENDDEBUG
148 
149  if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max
150  if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>1 .and. iout==ab_out)nkpt_eff=1
151 
152  if(option==1 .or. (option>=3 .and. option<=7))then
153 
154    do ienunit=0,1
155 
156      if (enunit==1 .and. ienunit==0)cycle
157      if (enunit==0 .and. ienunit==1)cycle
158 !  Print eigenvalues in hartree for enunit=0 or 2
159 !  The definition of two different strings is quite ridiculous. Historical reasons ...
160      
161      if (ienunit==0)then
162        convrt=one
163        strunit1='hartree'
164        strunit2='hartree'
165      end if
166      if (ienunit==1)then
167        convrt=Ha_eV
168        strunit1='   eV  '
169        strunit2='eV     '
170      end if
171 
172      band_index=0
173 
174      if(ienunit==0)then  ! XG20140730 I do not know why this is only done when ienunit==0
175        tmetal=0
176        if(option==1 .and. occopt>=3 .and. occopt<=8)tmetal=1
177        tmagnet=0
178        if(tmetal==1 .and. nsppol==2)then
179          tmagnet=1
180          rhoup = 0._dp
181          rhodn = 0._dp
182          nband_index = 1
183          do isppol=1,nsppol
184            do ikpt=1,nkpt
185              nband_k=nband(ikpt+(isppol-1)*nkpt)
186              do iband=1,nband_k
187                if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index)
188                if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index)
189                nband_index = nband_index + 1
190              end do
191            end do
192          end do
193          magnet = abs(rhoup - rhodn)
194        end if
195      end if
196      
197      if(iscf>=0 .and. (ienunit==0 .or. option==1))then
198        write(message, '(3a,f10.5,3a,f10.5)' ) &
199 &       ' Fermi (or HOMO) energy (',trim(strunit2),') =',convrt*fermie,'   Average Vxc (',trim(strunit2),')=',convrt*vxcavg
200        call wrtout(iout,message,'COLL')
201        if (prteig > 0) then
202          call wrtout(temp_unit,message,'COLL')
203        end if
204      end if
205 
206      
207 !    if( (iscf>=0 .or. iscf==-3) .and. ienunit==0)then     ! This is the most correct
208      if(iscf>=0 .and. ienunit==0)then ! For historical reasons
209        if(tmagnet==1)then
210          write(message, '(a,es16.8,a,a,es16.8,a,es16.8)' )&
211 &         ' Magnetization (Bohr magneton)=',magnet,ch10,&
212 &         ' Total spin up =',rhoup,'   Total spin down =',rhodn
213          call wrtout(iout,message,'COLL')
214          if (prteig > 0) then
215            call wrtout(temp_unit,message,'COLL')
216          end if
217        end if
218      end if
219 
220 !    Loop over spins (suppress spin data if nsppol not 2)
221      do isppol=1,nsppol
222 
223        ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
224        if (nsppol==2.and.isppol==1) then
225          write(message, '(4a,'//ikpt_fmt//',2x,a)' ) &
226 &         trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN UP:'
227        else if (nsppol==2.and.isppol==2) then
228          write(message, '(4a,'//ikpt_fmt//',2x,a)' ) &
229 &         trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN DOWN:'
230        else
231          write(message, '(4a,'//ikpt_fmt//',2x,a)' ) &
232 &         trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points:'
233        end if
234        call wrtout(iout,message,'COLL')
235        if (prteig > 0) then
236          call wrtout(temp_unit,message,'COLL')
237        end if
238 
239        if(ienunit==0)then
240          if(option>=4 .and. option<=7)then
241            message = '  (in case of degenerate eigenvalues, averaged derivative)'
242            call wrtout(iout,message,'COLL')
243            if (prteig > 0) then
244              call wrtout(temp_unit,message,'COLL')
245            end if
246          end if
247        end if
248 
249        do ikpt=1,nkpt
250          nband_k=nband(ikpt+(isppol-1)*nkpt)
251          ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
252          ibnd_fmt="i3" ; if(nband_k>=1000)ibnd_fmt="i6" ; if(nband_k>=1000000)ibnd_fmt="i9"
253          if(ikpt<=nkpt_eff)then
254            write(message, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
255 &           ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',&
256 &           kptns(1:3,ikpt)+tol10,' (reduced coord)'
257            call wrtout(iout,message,'COLL')
258            if (prteig > 0) then
259              call wrtout(temp_unit,message,'COLL')
260            end if
261            do ii=0,(nband_k-1)/8
262 !            write(message, '(8f15.10)' ) (convrt*eigen(iband+band_index),&
263              write(message, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index),&
264 &             iband=1+ii*8,min(nband_k,8+ii*8))
265              call wrtout(iout,message,'COLL')
266              if (prteig > 0) then
267                call wrtout(temp_unit,message,'COLL')
268              end if
269            end do
270            if(ienunit==0 .and. option==1 .and. occopt>=3 .and. occopt<=8)then
271              write(message, '(5x,a,'//ikpt_fmt//')' )  ' occupation numbers for kpt#',ikpt
272              call wrtout(iout,message,'COLL')
273              do ii=0,(nband_k-1)/8
274                write(message, '(8(f10.5,1x))' ) (occ(iband+band_index),&
275 &               iband=1+ii*8,min(nband_k,8+ii*8))
276                call wrtout(iout,message,'COLL')
277              end do
278            end if
279 
280          else
281            if(ikpt==nkpt_eff+1)then
282              write(message, '(a,a)' ) &
283 &             ' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10
284              call wrtout(iout,message,'COLL')
285            end if
286            if (prteig > 0) then
287              write(message, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
288 &             ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',&
289 &             kptns(1:3,ikpt)+tol10,' (reduced coord)'
290              call wrtout(temp_unit,message,'COLL')
291              do ii=0,(nband_k-1)/8
292                write(message, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index),&
293 &               iband=1+ii*8,min(nband_k,8+ii*8))
294                call wrtout(temp_unit,message,'COLL')
295              end do
296            end if
297          end if
298          band_index=band_index+nband_k
299        end do ! do ikpt=1,nkpt
300      end do ! do isppol=1,nsppol
301 
302    end do ! End loop over Hartree or eV
303 
304  else if(option==2)then
305 
306    band_index=0
307    do isppol=1,nsppol
308 
309      if(nsppol==2)then
310        if(isppol==1)write(message, '(2a)' ) ch10,' SPIN UP channel '
311        if(isppol==2)write(message, '(2a)' ) ch10,' SPIN DOWN channel '
312        call wrtout(iout,message,'COLL')
313        if(prteig>0) then
314          call wrtout(temp_unit,message,'COLL')
315        end if
316      end if
317 
318      do ikpt=1,nkpt
319        nband_k=nband(ikpt+(isppol-1)*nkpt)
320        ikpt_fmt="i5" ; if(nkpt>=10000)ikpt_fmt="i7" ; if(nkpt>=1000000)ikpt_fmt="i9"
321 
322        if(ikpt<=nkpt_eff)then
323          write(message, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) &
324 &         'Non-SCF case, kpt',ikpt,' (',(kptns(ii,ikpt),ii=1,3),'), residuals and eigenvalues='
325          call wrtout(iout,message,'COLL')
326          if (prteig > 0) then
327            write(message, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) &
328 &           'Non-SCF case, kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') '
329            call wrtout(temp_unit,message,'COLL')
330          end if
331          do ii=0,(nband_k-1)/8
332            write(message, '(1p,8e10.2)' )(resid(iband+band_index),iband=1+8*ii,min(8+8*ii,nband_k))
333            call wrtout(iout,message,'COLL')
334          end do
335          do ii=0,(nband_k-1)/6
336            write(message, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
337            call wrtout(iout,message,'COLL')
338            if (prteig > 0) then
339              call wrtout(temp_unit,message,'COLL')
340            end if
341          end do
342        else
343          if(ikpt==nkpt_eff+1)then
344            write(message, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10
345            call wrtout(iout,message,'COLL')
346          end if
347          if (prteig > 0) then
348            write(message, '(1x,a,i5,a,f9.5,2f9.5,a)' ) &
349 &           'Non-SCF kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') '
350            call wrtout(temp_unit,message,'COLL')
351            do ii=0,(nband_k-1)/6
352              write(message, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
353              call wrtout(temp_unit,message,'COLL')
354            end do
355          end if
356        end if
357 
358        residk=maxval(resid(band_index+1:band_index+nband_k))
359        if (residk>tolwfr) then
360          write(message, '(1x,a,2i5,a,1p,e13.5)' ) &
361 &         ' prteigrs : nnsclo,ikpt=',nnsclo_now,ikpt,' max resid (incl. the buffer)=',residk
362          call wrtout(iout,message,'COLL')
363        end if
364 
365        band_index=band_index+nband_k
366      end do
367    end do
368    call wrtout(iout," ",'COLL')
369 
370  else
371    write(message, '(a,i0,a)' )' option = ',option,', is not an allowed value.'
372    MSG_BUG(message)
373  end if
374 
375  if (prteig > 0) close (temp_unit)
376 
377 end subroutine prteigrs