TABLE OF CONTENTS


ABINIT/m_occ [ Modules ]

[ Top ] [ Modules ]

NAME

 m_occ

FUNCTION

  Low-level functions for occupation factors.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (XG, AF)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_occ
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_splines
33  use m_xmpi
34 
35  use m_time,         only : timab
36  use m_fstrings,     only : sjoin, itoa
37  use defs_abitypes,  only : MPI_type
38  use m_mpinfo,       only : proc_distrb_cycle
39 
40  implicit none
41 
42  private

m_abinit/getnel [ Functions ]

[ Top ] [ Functions ]

NAME

 getnel

FUNCTION

 Option=1 :
 Get the total number of electrons nelect, given a trial fermienergy fermie.
 For this, compute new occupation numbers at each k point,
 from eigenenergies eigen, according to the
 smearing scheme defined by occopt (and smearing width tsmear or tphysel).

 Option=2 :
 Compute and output the smeared density of states, and the integrated density
 of states, then write these data

 Warning : this routine assumes checks have been done in the calling
 routine, and that the values of the arguments are sensible

 NOTE
 in order to speed the calculation, it would be easy to
 compute the entropy only when the fermi energy is well converged

INPUTS

 dosdeltae= DOS delta of Energy (needed if Option=2)
 eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
 fermie= fermi energy (Hartree)
 maxocc=asymptotic maximum occupation number per band
 mband=maximum number of bands
 nband(nkpt*nsppol)=number of bands at each k point
 nkpt=number of k points
 nsppol=1 for unpolarized, 2 for spin-polarized
 occopt=option for occupancies, or re-smearing scheme if dblsmr /= 0
 option=see above
 tphysel="physical" electronic temperature with FD occupations
 tsmear=smearing width (or temperature)
 unitdos=unit number of output of the DOS. Not needed if option==1
 wtk(nkpt)=k point weights

OUTPUT

 doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy for each band and k point.
 entropy= entropy associated with the smearing (adimensional)
 nelect=number of electrons per unit cell
 occ(mband*nkpt*nsppol)=occupancies for each band and k point.

NOTES

 Modified beginning 23/11/2000 by MV
 Add an additional smearing on top of a FD type, in order to improve k-point
 convergence: tsmear = 0 and tphysel ~= 2.e-3 corresponds to a small (300K)
 temperature on the electrons insufficient for convergence purposes.
 Feed re-smeared "Dirac delta" to the rest of ABINIT with only one parameter,
 tphysel, which is the physical temperature.
 encorr = correction to energy for terms of order tsmear^2:
       $  E_{phys} = E_{free} - encorr*(E_{int}-E_{free}) + O(tsmear^3)  $

PARENTS

      cchi0q0_intraband,clnup1,conducti_nc,dfpt_looppert,m_ebands,newocc

CHILDREN

      dos_hdr_write,init_occ_ent,splfit,wrtout

SOURCE

116 subroutine getnel(doccde,dosdeltae,eigen,entropy,fermie,maxocc,mband,nband,&
117 &  nelect,nkpt,nsppol,occ,occopt,option,tphysel,tsmear,unitdos,wtk)
118 
119 
120 !This section has been created automatically by the script Abilint (TD).
121 !Do not modify the following lines by hand.
122 #undef ABI_FUNC
123 #define ABI_FUNC 'getnel'
124 !End of the abilint section
125 
126  implicit none
127 
128 !Arguments ------------------------------------
129 !scalars
130  integer,intent(in) :: mband,nkpt,nsppol,occopt,option,unitdos
131  real(dp),intent(in) :: dosdeltae,fermie,maxocc,tphysel,tsmear
132  real(dp),intent(out) :: entropy,nelect
133 !arrays
134  integer,intent(in) :: nband(nkpt*nsppol)
135  real(dp),intent(in) :: eigen(mband*nkpt*nsppol),wtk(nkpt)
136  real(dp),intent(out) :: doccde(mband*nkpt*nsppol) !vz_i
137  real(dp),intent(inout) :: occ(mband*nkpt*nsppol) !vz_i
138 
139 !Local variables-------------------------------
140 ! nptsdiv2 is the number of integration points, divided by 2.
141 ! tratio  = ratio tsmear/tphysel for convoluted smearing function
142 ! save values so we can impose recalculation of smdfun when
143 ! the smearing or electronic temperature change between
144 ! datasets
145 ! corresponds roughly to delta_FD (maxFDarg) = 1.0d-100
146 !
147 ! return fermi-dirac smearing function analytically
148 !
149 ! real(dp) :: smdFD
150 ! smdFD (tt) = 1.0_dp / (exp(-tt/2.0_dp) + exp(tt/2.0_dp))**2
151 !scalars
152 ! TODO: This parameter is defined in init_occ_ent but we cannot call the
153 ! routine to get this value since the same variable is used to dimension the
154 ! arrays! This Constants should be stored somewhere in a module.
155  integer,parameter :: nptsdiv2_def=6000
156  integer,parameter :: prtdos1=1
157  integer :: bantot,iband,iene,ikpt,index,index_start,isppol
158  integer :: nene,nptsdiv2
159  real(dp) :: buffer,deltaene,dosdbletot,doshalftot,dostot
160  real(dp) :: enemax,enemin,enex,intdostot,limit,tsmearinv
161  character(len=500) :: message
162 !arrays
163  real(dp),allocatable :: entfun(:,:),occfun(:,:)
164  real(dp),allocatable :: smdfun(:,:),xgrid(:)
165  real(dp),allocatable :: arg(:),derfun(:),dos(:),dosdble(:),doshalf(:),ent(:)
166  real(dp),allocatable :: intdos(:)
167 
168 ! *************************************************************************
169 
170  DBG_ENTER("COLL")
171 
172  if(option/=1 .and. option/=2)then
173    MSG_BUG(sjoin('Option must be either 1 or 2. It is:', itoa(option)))
174  end if
175 
176 !Initialize the occupation function and generalized entropy function,
177 !at the beginning, or if occopt changed
178 
179 !Just get the number nptsdiv2 and allocate entfun, occfun,
180 !smdfun and xgrid accordingly
181  nptsdiv2 = nptsdiv2_def
182 
183 ! call init_occ_ent(entfun, limit, &
184 !& nptsdiv2, occfun, occopt, -1, smdfun, tphysel, &
185 !& tsmear, tsmearinv, xgrid)
186 
187  ABI_ALLOCATE(entfun,(-nptsdiv2:nptsdiv2,2))
188  ABI_ALLOCATE(occfun,(-nptsdiv2:nptsdiv2,2))
189  ABI_ALLOCATE(smdfun,(-nptsdiv2:nptsdiv2,2))
190  ABI_ALLOCATE(xgrid,(-nptsdiv2:nptsdiv2))
191 
192 !Call to init_occ_ent
193  call init_occ_ent(entfun, limit, &
194 & nptsdiv2, occfun, occopt, option, smdfun, tphysel, &
195 & tsmear, tsmearinv, xgrid)
196 
197 !The initialisation of occfun and entfun is done
198 
199 !---------------------------------------------------------------------
200 
201 !write(std_out,*)' getnel : debug  tphysel, tsmear = ', tphysel, tsmear
202  bantot=sum(nband(:))
203 
204  ABI_ALLOCATE(arg,(bantot))
205  ABI_ALLOCATE(derfun,(bantot))
206  ABI_ALLOCATE(ent,(bantot))
207 
208  if(option==1)then
209    !normal evaluation of occupations and entropy
210 
211 !  Compute the arguments of the occupation and entropy functions
212    arg(:)=(fermie-eigen(1:bantot))*tsmearinv
213 
214 !  Compute the values of the occupation function, and the entropy function
215 !  Note : splfit also takes care of the points outside of the interval,
216 !  and assign to them the value of the closest extremal point,
217 !  which is what is needed here.
218 
219    call splfit(xgrid,doccde,occfun,1,arg,occ,(2*nptsdiv2+1),bantot)
220    call splfit(xgrid,derfun,entfun,0,arg,ent,(2*nptsdiv2+1),bantot)
221 
222 !  Normalize occ and ent, and sum number of electrons and entropy
223    nelect=zero; entropy=zero
224    index=0
225    do isppol=1,nsppol
226      do ikpt=1,nkpt
227        do iband=1,nband(ikpt+nkpt*(isppol-1))
228          index=index+1
229          ent(index)=ent(index)*maxocc
230          occ(index)=occ(index)*maxocc
231          doccde(index)=-doccde(index)*maxocc*tsmearinv
232          entropy=entropy+wtk(ikpt)*ent(index)
233          nelect=nelect+wtk(ikpt)*occ(index)
234        end do
235      end do
236    end do
237 
238 !  write(std_out,*) ' getnel : debug   wtk, occ, eigen = ', wtk, occ, eigen
239 !  write(std_out,*)xgrid(-nptsdiv2),xgrid(nptsdiv2)
240 !  write(std_out,*)'fermie',fermie
241 !  do ii=1,bantot
242 !  write(std_out,*)ii,arg(ii),doccde(ii)
243 !  end do
244 !  write(std_out,*)'eigen',eigen(:)
245 !  write(std_out,*)'arg',arg(:)
246 !  write(std_out,*)'occ',occ(:)
247 !  write(std_out,*)'nelect',nelect
248 
249  else if(option==2)then
250   ! evaluate DOS for smearing, half smearing, and double.
251 
252    buffer=limit/tsmearinv*.5_dp
253 
254   ! A Similar section is present is dos_calcnwrite. Should move all DOS stuff to m_ebands
255   ! Choose the lower and upper energies
256    enemax=maxval(eigen(1:bantot))+buffer
257    enemin=minval(eigen(1:bantot))-buffer
258 
259   ! Extend the range to a nicer value
260    enemax=0.1_dp*ceiling(enemax*10._dp)
261    enemin=0.1_dp*floor(enemin*10._dp)
262 
263   ! Choose the energy increment
264    if(abs(dosdeltae)<tol10)then
265      deltaene=0.001_dp
266      if(prtdos1>=2)deltaene=0.0005_dp ! Higher resolution possible (and wanted) for tetrahedron
267    else
268      deltaene=dosdeltae
269    end if
270    nene=nint((enemax-enemin)/deltaene)+1
271 
272 !  Write the header of the DOS file, and also decides the energy range and increment
273    call dos_hdr_write(deltaene,eigen,enemax,enemin,fermie,mband,nband,nene,&
274 &   nkpt,nsppol,occopt,prtdos1,tphysel,tsmear,unitdos)
275 
276    ABI_ALLOCATE(dos,(bantot))
277    ABI_ALLOCATE(dosdble,(bantot))
278    ABI_ALLOCATE(doshalf,(bantot))
279    ABI_ALLOCATE(intdos,(bantot))
280 
281    do isppol=1,nsppol
282 
283      if (nsppol==2) then
284        if(isppol==1) write(message,'(a,16x,a)')  '#','Spin-up DOS'
285        if(isppol==2) write(message,'(2a,16x,a)')  ch10,'#','Spin-dn DOS '
286        call wrtout(unitdos,message,'COLL')
287      end if
288      index_start=0
289      if(isppol==2)then
290        do ikpt=1,nkpt
291          index_start=index_start+nband(ikpt)
292        end do
293      end if
294 
295      enex=enemin
296      do iene=1,nene
297 
298 !      Compute the arguments of the dos and occupation function
299        arg(:)=(enex-eigen(1:bantot))*tsmearinv
300 
301        call splfit(xgrid,derfun,smdfun,0,arg,dos,(2*nptsdiv2+1),bantot)
302        call splfit(xgrid,derfun,occfun,0,arg,intdos,(2*nptsdiv2+1),bantot)
303 !      Also compute the dos with tsmear halved and doubled
304        arg(:)=arg(:)*2.0_dp
305        call splfit(xgrid,derfun,smdfun,0,arg,doshalf,(2*nptsdiv2+1),bantot)
306 !      Since arg was already doubled, must divide by four
307        arg(:)=arg(:)*0.25_dp
308        call splfit(xgrid,derfun,smdfun,0,arg,dosdble,(2*nptsdiv2+1),bantot)
309 
310 !      Now, accumulate the contribution from each eigenenergy
311        dostot=zero
312        intdostot=zero
313        doshalftot=zero
314        dosdbletot=zero
315        index=index_start
316 
317 !      write(std_out,*)' eigen, arg, dos, intdos, doshalf, dosdble'
318        do ikpt=1,nkpt
319          do iband=1,nband(ikpt+nkpt*(isppol-1))
320            index=index+1
321            dostot=dostot+wtk(ikpt)*maxocc*dos(index)*tsmearinv
322            intdostot=intdostot+wtk(ikpt)*maxocc*intdos(index)
323            doshalftot=doshalftot+wtk(ikpt)*maxocc*doshalf(index)*tsmearinv*2.0_dp
324            dosdbletot=dosdbletot+wtk(ikpt)*maxocc*dosdble(index)*tsmearinv*0.5_dp
325          end do
326        end do
327 
328 !      Print the data for this energy
329        write(unitdos, '(f8.3,2f14.6,2f14.3)' )enex,dostot,intdostot,doshalftot,dosdbletot
330 
331        enex=enex+deltaene
332      end do ! iene
333    end do ! isppol
334 
335    ABI_DEALLOCATE(dos)
336    ABI_DEALLOCATE(dosdble)
337    ABI_DEALLOCATE(doshalf)
338    ABI_DEALLOCATE(intdos)
339 
340 !  MG: It does not make sense to close the unit here since the routines
341 !  did not open the file here!
342 !  Close the DOS file
343    close(unitdos)
344  end if
345 
346  ABI_DEALLOCATE(arg)
347  ABI_DEALLOCATE(derfun)
348  ABI_DEALLOCATE(ent)
349  ABI_DEALLOCATE(entfun)
350  ABI_DEALLOCATE(occfun)
351  ABI_DEALLOCATE(smdfun)
352  ABI_DEALLOCATE(xgrid)
353 
354  DBG_EXIT("COLL")
355 
356 end subroutine getnel

m_occ/dos_hdr_write [ Functions ]

[ Top ] [ m_occ ] [ Functions ]

NAME

 dos_hdr_write

FUNCTION

 Write the header of the DOS files, for both smearing and tetrahedron methods.

INPUTS

 deltaene=increment of DOS energy arguments
 enemax=maximal value of the DOS energy argument
 enemin=minimal value of the DOS energy argument
 nene=number of DOS energy argument
 eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
 fermie=fermi energy useful for band alignment...
 mband=maximum number of bands
 nband(nkpt*nsppol)=number of bands at each k point
 nkpt=number of k points
 nsppol=1 for unpolarized, 2 for spin-polarized
 occopt=option for occupancies, or re-smearing scheme if dblsmr /= 0
 prtdos=1 for smearing technique, 2 or 3 for tetrahedron technique
 tphysel="physical" electronic temperature with FD occupations
 tsmear=smearing width (or temperature)
 unitdos=unit number of output of the DOS.

OUTPUT

   Only writing.

PARENTS

      getnel,m_epjdos

CHILDREN

      wrtout

SOURCE

1466 subroutine dos_hdr_write(deltaene,&
1467 &  eigen,enemax,enemin,fermie,mband,nband,nene,&
1468 &  nkpt,nsppol,occopt,prtdos,tphysel,tsmear,unitdos)
1469 
1470 
1471 !This section has been created automatically by the script Abilint (TD).
1472 !Do not modify the following lines by hand.
1473 #undef ABI_FUNC
1474 #define ABI_FUNC 'dos_hdr_write'
1475 !End of the abilint section
1476 
1477  implicit none
1478 
1479 !Arguments ------------------------------------
1480 !scalars
1481  integer,intent(in) :: mband,nkpt,nsppol,occopt,prtdos,unitdos,nene
1482  real(dp),intent(in) :: fermie,tphysel,tsmear
1483  real(dp),intent(in) :: deltaene,enemax,enemin
1484 !arrays
1485  integer,intent(in) :: nband(nkpt*nsppol)
1486  real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
1487 
1488 !Local variables-------------------------------
1489 !scalars
1490  character(len=500) :: msg
1491 
1492 ! *************************************************************************
1493 
1494 !Write the DOS file
1495  write(msg, '(7a,i2,a,i5,a,i4)' ) "#",ch10, &
1496 & '# ABINIT package : DOS file  ',ch10,"#",ch10,&
1497 & '# nsppol =',nsppol,', nkpt =',nkpt,', nband(1)=',nband(1)
1498  call wrtout(unitdos,msg,'COLL')
1499 
1500  if (any(prtdos== [1,4])) then
1501    write(msg, '(a,i2,a,f6.3,a,f6.3,a)' )  &
1502 &   '# Smearing technique, occopt =',occopt,', tsmear=',tsmear,' Hartree, tphysel=',tphysel,' Hartree'
1503  else
1504    write(msg, '(a)' ) '# Tetrahedron method '
1505  end if
1506  call wrtout(unitdos,msg,'COLL')
1507 
1508  if (mband*nkpt*nsppol>=3) then
1509    write(msg, '(a,3f8.3,2a)' )'# For identification : eigen(1:3)=',eigen(1:3),ch10,"#"
1510  else
1511    write(msg, '(a,3f8.3)' ) '# For identification : eigen=',eigen
1512    write(msg, '(3a)')trim(msg),ch10,"#"
1513  end if
1514  call wrtout(unitdos,msg,'COLL')
1515 
1516  write(msg, '(a,f16.8)' ) '# Fermi energy : ', fermie
1517  call wrtout(unitdos,msg,'COLL')
1518 
1519  if (prtdos==1) then
1520    write(msg, '(5a)' ) "#",ch10,&
1521 &   '# The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell),',&
1522 &   ch10,'# as well as the DOS with tsmear halved and doubled, are computed,'
1523 
1524  else if (prtdos==2)then
1525    write(msg, '(3a)' ) "#",ch10,&
1526 &   '# The DOS (in electrons/Hartree/cell) and integrated DOS (in electrons/cell) are computed,'
1527 
1528  else if (any(prtdos == [3, 4])) then
1529    write(msg, '(5a)' ) "#",ch10,&
1530 &   '# The local DOS (in electrons/Hartree for one atomic sphere)',ch10,&
1531 &   '# and integrated local DOS (in electrons for one atomic sphere) are computed.'
1532 
1533  else if (prtdos==5)then
1534    write(msg, '(9a)' ) "#",ch10,&
1535 &   '# The spin component DOS (in electrons/Hartree/cell)',ch10,&
1536 &   '# and integrated spin component DOS (in electrons/cell) are computed.',ch10,&
1537 &   '# Remember that the wf are eigenstates of S_z and S^2, not S_x and S_y',ch10,&
1538 &   '#   so the latter will not always sum to 0 for paired electronic states.'
1539  end if
1540  call wrtout(unitdos,msg,'COLL')
1541 
1542  write(msg, '(a,i5,a,a,a,f9.4,a,f9.4,a,f8.5,a,a,a)' )&
1543 & '# at ',nene,' energies (in Hartree) covering the interval ',ch10,&
1544 & '# between ',enemin,' and ',enemax,' Hartree by steps of ',deltaene,' Hartree.',ch10,"#"
1545  call wrtout(unitdos,msg,'COLL')
1546 
1547  if (prtdos==1) then
1548    write(msg, '(a,a)' )&
1549 &   '#       energy        DOS       Integr. DOS   ','     DOS           DOS    '
1550    call wrtout(unitdos,msg,'COLL')
1551 
1552    write(msg, '(a)' )&
1553 &   '#                                              (tsmear/2)    (tsmear*2) '
1554    call wrtout(unitdos,msg,'COLL')
1555  else
1556    write(msg, '(a)' ) '#       energy        DOS '
1557  end if
1558 
1559 end subroutine dos_hdr_write

m_occ/init_occ_ent [ Functions ]

[ Top ] [ m_occ ] [ Functions ]

NAME

 init_occ_ent

FUNCTION

INPUTS

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

PARENTS

      getnel

CHILDREN

      spline

SOURCE

 779 subroutine init_occ_ent(entfun,limit,nptsdiv2,occfun,occopt,option,smdfun,tphysel,tsmear,tsmearinv,xgrid)
 780 
 781 
 782 !This section has been created automatically by the script Abilint (TD).
 783 !Do not modify the following lines by hand.
 784 #undef ABI_FUNC
 785 #define ABI_FUNC 'init_occ_ent'
 786 !End of the abilint section
 787 
 788  implicit none
 789 
 790 !Arguments ------------------------------------
 791 !scalars
 792  integer,intent(in) :: occopt,option
 793  real(dp),intent(in) :: tphysel,tsmear
 794  integer,intent(inout) :: nptsdiv2
 795  real(dp),intent(out) :: limit,tsmearinv
 796  real(dp),intent(inout) :: entfun(-nptsdiv2:nptsdiv2,2),occfun(-nptsdiv2:nptsdiv2,2)
 797  real(dp),intent(inout) :: smdfun(-nptsdiv2:nptsdiv2,2),xgrid(-nptsdiv2:nptsdiv2)
 798 
 799 
 800 !Local variables-------------------------------
 801 !scalars
 802  integer :: algo,ii,jj,nconvd2
 803  integer :: nmaxFD,nminFD
 804  integer,parameter :: nptsdiv2_def=6000
 805  integer,save :: dblsmr,occopt_prev=-9999
 806  real(dp),parameter :: maxFDarg=500.0_dp
 807  real(dp),save :: convlim,incconv,limit_occ,tphysel_prev=-9999,tsmear_prev=-9999
 808  real(dp) :: aa,dsqrpi,encorr,factor
 809  real(dp) :: expinc,expx22,expxo2,gauss,increm
 810  real(dp) :: resFD1,resFD2,resFD3,resFD4,resmom,resmom1,resmom2
 811  real(dp) :: resmom3,resmom4,secmom,smom1,smom2,thdmom,tmom1,tmom2,tmpexpsum
 812  real(dp) :: tmpsmdfun,tratio,tt,xx,yp1,ypn
 813  character(len=500) :: message
 814 !arrays
 815  real(dp),save :: entfun_prev(-nptsdiv2_def:nptsdiv2_def,2),occfun_prev(-nptsdiv2_def:nptsdiv2_def,2)
 816  real(dp),save :: smdfun_prev(-nptsdiv2_def:nptsdiv2_def,2),xgrid_prev(-nptsdiv2_def:nptsdiv2_def)
 817  real(dp),allocatable :: entder(:),occder(:),smd1(:),smd2(:)
 818  real(dp),allocatable :: smdder(:),tgrid(:),work(:),workfun(:)
 819 
 820 ! *************************************************************************
 821 
 822 !Initialize the occupation function and generalized entropy function,
 823 !at the beginning, or if occopt changed
 824 
 825  if(option==-1)then
 826    nptsdiv2 = nptsdiv2_def
 827    return
 828  end if
 829 
 830 
 831  if(occopt_prev/=occopt           .or. &
 832 & abs(tsmear_prev-tsmear)  >tol12 .or. &
 833 & abs(tphysel_prev-tphysel)>tol12       ) then
 834 !  write(std_out,*) 'INIT_OCC_ENT CHANGE ..........'
 835    occopt_prev=occopt
 836    tsmear_prev=tsmear
 837    tphysel_prev=tphysel
 838 
 839 !  Check whether input values of tphysel tsmear and occopt are consistent
 840    dblsmr = 0
 841    if (abs(tphysel)>tol12) then
 842 !    Use re-smearing scheme
 843      if (abs(tsmear)>tol12) then
 844        dblsmr = 1
 845 !      Use FD occupations (one smearing) only with "physical" temperature tphysel
 846      else if (occopt /= 3) then
 847        write(message, '(a,i6,a)' )' tphysel /= 0, tsmear == 0, but occopt is not = 3, but ',occopt,'.'
 848        MSG_ERROR(message)
 849      end if
 850    end if
 851 !  write(std_out,*) 'getnel : input read.'
 852 !  write(std_out,*) '  dblsmr = ', dblsmr
 853 !  write(std_out,*) '  tphysel, tsmear = ', tphysel, tsmear
 854 
 855    ABI_ALLOCATE(entder,(-nptsdiv2_def:nptsdiv2_def))
 856    ABI_ALLOCATE(occder,(-nptsdiv2_def:nptsdiv2_def))
 857    ABI_ALLOCATE(smdder,(-nptsdiv2_def:nptsdiv2_def))
 858    ABI_ALLOCATE(workfun,(-nptsdiv2_def:nptsdiv2_def))
 859    ABI_ALLOCATE(work,(-nptsdiv2_def:nptsdiv2_def))
 860 
 861 !  Prepare the points on the grid
 862 !  limit is the value of the argument that will give 0.0 or 1.0 , with
 863 !  less than about 1.0d-15 error for 4<=occopt<=8, and less than about 1.0d-12
 864 !  error for occopt==3. It is not worth to compute the function beyond
 865 !  that point. Even with a less severe requirement, it is significantly
 866 !  larger for occopt==3, with an exponential
 867 !  tail, than for the other occupation functions, with a Gaussian tail.
 868 !  Note that these values are useful in newocc.f also.
 869    limit_occ=6.0_dp
 870    if(occopt==3)limit_occ=30.0_dp
 871    if(dblsmr /= 0) then
 872      tratio = tsmear / tphysel
 873      limit_occ=30.0_dp + 6.0_dp*tratio
 874    end if
 875 
 876 !  With nptsdiv2_def=6000 (thus increm=0.001 for 4<=occopt<=8,
 877 !  and increm=0.005 for occopt==3, the O(1/N4) algorithm gives 1.0d-12
 878 !  accuracy on the stored values occfun and entfun. These, together
 879 !  with smdfun and xgrid_prev, need permanently about 0.67 MB, which is affordable.
 880    increm=limit_occ/nptsdiv2_def
 881    do ii=-nptsdiv2_def,nptsdiv2_def
 882      xgrid_prev(ii)=ii*increm
 883    end do
 884 
 885 !  ---------------------------------------------------------
 886 !  Ordinary (unique) smearing function
 887 !  ---------------------------------------------------------
 888    if (dblsmr == 0) then
 889 
 890 !    Compute the unnormalized smeared delta function between -limit_occ and +limit_occ
 891 !    (well, they are actually normalized ...)
 892      if(occopt==3)then
 893 
 894 !      Fermi-Dirac
 895        do ii=0,nptsdiv2_def
 896          xx=xgrid_prev(ii)
 897          smdfun_prev( ii,1)=0.25_dp/(cosh(xx/2.0_dp)**2)
 898          smdfun_prev(-ii,1)=smdfun_prev(ii,1)
 899        end do
 900 
 901      else if(occopt==4 .or. occopt==5)then
 902 
 903 !      Cold smearing of Marzari, two values of the "a" parameter being possible
 904 !      first value gives minimization of the bump
 905        if(occopt==4)aa=-.5634
 906 !      second value gives monotonic occupation function
 907        if(occopt==5)aa=-.8165
 908 
 909        dsqrpi=1.0_dp/sqrt(pi)
 910        do ii=0,nptsdiv2_def
 911          xx=xgrid_prev(ii)
 912          gauss=dsqrpi*exp(-xx**2)
 913          smdfun_prev( ii,1)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx)))
 914          smdfun_prev(-ii,1)=gauss*(1.5_dp+xx*( aa*1.5_dp+xx*(-1.0_dp-aa*xx)))
 915        end do
 916 
 917      else if(occopt==6)then
 918 
 919 !      First order Hermite-Gaussian of Paxton and Methfessel
 920        dsqrpi=1.0_dp/sqrt(pi)
 921        do ii=0,nptsdiv2_def
 922          xx=xgrid_prev(ii)
 923          smdfun_prev( ii,1)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2)
 924          smdfun_prev(-ii,1)=smdfun_prev(ii,1)
 925        end do
 926 
 927      else if(occopt==7)then
 928 
 929 !      Gaussian smearing
 930        dsqrpi=1.0_dp/sqrt(pi)
 931        do ii=0,nptsdiv2_def
 932          xx=xgrid_prev(ii)
 933          smdfun_prev( ii,1)=dsqrpi*exp(-xx**2)
 934          smdfun_prev(-ii,1)=smdfun_prev(ii,1)
 935        end do
 936 
 937      else if(occopt==8)then
 938 
 939 !      Constant value of the delta function over the smearing interval, for testing purposes only.
 940        do ii=0,nptsdiv2_def
 941          xx=xgrid_prev(ii)
 942          if(xx>half+tol8)then
 943            smdfun_prev( ii,1)=zero
 944          else if(xx<half-tol8)then
 945            smdfun_prev( ii,1)=one
 946          else
 947            smdfun_prev( ii,1)=half
 948          end if
 949          smdfun_prev(-ii,1)=smdfun_prev(ii,1)
 950        end do
 951 
 952      else
 953        write(message, '(a,i0,a)' )' Occopt=',occopt,' is not allowed in getnel. '
 954        MSG_BUG(message)
 955      end if
 956 
 957 !    ---------------------------------------------------------
 958 !    smear FD delta with occopt delta calculated in smdfun_prev
 959 !    ---------------------------------------------------------
 960    else if (dblsmr /= 0) then
 961 
 962      nconvd2 = 6000
 963      convlim = 10.0_dp
 964      incconv = convlim / nconvd2
 965 
 966 !    store smearing functions in smd1 and smd2
 967      ABI_ALLOCATE(smd1,(-nconvd2:nconvd2))
 968      ABI_ALLOCATE(smd2,(-nconvd2:nconvd2))
 969      ABI_ALLOCATE(tgrid,(-nconvd2:nconvd2))
 970 
 971 !    FD function in smd1( ii) and second smearing delta in smd2( ii)
 972 !
 973 !    smd1(:) contains delta_FD ( x )
 974      do ii=0,nconvd2
 975        tgrid(ii)=ii*incconv
 976        tgrid(-ii)=-tgrid(ii)
 977        tt=tgrid(ii)
 978        smd1( ii)=0.25_dp/(cosh(tt/2.0_dp)**2)
 979        smd1(-ii)=smd1(ii)
 980      end do
 981 
 982 !    check input values of occopt and fill smd2(:) with appropriate data:
 983 !    smd2(:) contains delta_resmear ( x )
 984      if(occopt == 3) then
 985        write(message, '(a,a)' )&
 986 &       'Occopt=3 is not allowed as a re-smearing.', &
 987 &       'Use a single FD, or re-smear with a different delta type (faster cutoff). '
 988        MSG_ERROR(message)
 989      else if(occopt==4 .or. occopt==5)then
 990 !      Cold smearing of Marzari, two values of the "a" parameter being possible
 991 !      first value gives minimization of the bump
 992        if(occopt==4)aa=-.5634
 993 !      second value gives monotonic occupation function
 994        if(occopt==5)aa=-.8165
 995 
 996        dsqrpi=1.0_dp/sqrt(pi)
 997        do ii=0,nconvd2
 998          tt=tgrid(ii)
 999          gauss=dsqrpi*exp(-tt**2)
1000          smd2( ii)=gauss*(1.5_dp+tt*(-aa*1.5_dp+tt*(-1.0_dp+aa*tt)))
1001          smd2(-ii)=gauss*(1.5_dp+tt*( aa*1.5_dp+tt*(-1.0_dp-aa*tt)))
1002        end do
1003      else if(occopt==6)then
1004        dsqrpi=1.0_dp/sqrt(pi)
1005        do ii=0,nconvd2
1006          tt=tgrid(ii)
1007          smd2( ii)=dsqrpi*(1.5_dp-tt**2)*exp(-tt**2)
1008          smd2(-ii)=smd2(ii)
1009        end do
1010      else if(occopt==7)then
1011        dsqrpi=1.0_dp/sqrt(pi)
1012        do ii=0,nconvd2
1013          tt=tgrid(ii)
1014          smd2( ii)=dsqrpi*exp(-tt**2)
1015          smd2(-ii)=smd2(ii)
1016        end do
1017      else if(occopt==8)then
1018        do ii=0,nconvd2
1019          tt=tgrid(ii)
1020          if(tt>half+tol8)then
1021            smd2( ii)=zero
1022          else if(tt<half-tol8)then
1023            smd2( ii)=one
1024          else
1025            smd2( ii)=half
1026          end if
1027          smd2(-ii)=smd2(ii)
1028        end do
1029      else
1030        write(message, '(a,i0,a)' )' Occopt= ',occopt,' is not allowed in getnel. '
1031        MSG_BUG(message)
1032      end if
1033 
1034 
1035 !    Use O(1/N4) algorithm from Num Rec (see below)
1036 !
1037 !    The grid for the convoluted delta is taken (conservatively)
1038 !    to be that for the FD delta ie 6000 pts in [-limit_occ;limit_occ]
1039 !    Smearing functions are given on [-dbllim;dbllim] and the grid must
1040 !    superpose the normal grid on [-limit_occ:limit_occ]
1041 !    The maximal interval for integration of the convolution is
1042 !    [-dbllim+limit_occ+lim(delta2);dbllim-limit_occ-lim(delta2)] =
1043 !    [-dbllim+36;dbllim-36]
1044 
1045 !    test the smdFD function for extreme values:
1046 !    do jj=-nptsdiv2_def,-nptsdiv2_def
1047 !    do ii=-nconvd2+4,nconvd2
1048 !    call smdFD(xgrid_prev(jj) - tgrid(ii)*tratio, resFD)
1049 !    write(std_out,*) 'ii jj = ', ii,jj, ' smdFD (', &
1050 !    &    xgrid_prev(jj) - tgrid(ii)*tratio, ') ', resFD
1051 !    end do
1052 !    end do
1053 
1054      expinc = exp(half*incconv*tratio)
1055 
1056 !    jj = position of point at which we are calculating smdfun_prev
1057      do jj=-nptsdiv2_def,nptsdiv2_def
1058 !      Do not care about the 8 boundary points,
1059 !      where the values should be extremely small anyway
1060        smdfun_prev(jj,1)=0.0_dp
1061 !      only add contribution with delta_FD > 1.0d-100
1062        nmaxFD = floor  (( maxFDarg+xgrid_prev(jj)) / tratio / incconv )
1063        nmaxFD = min (nmaxFD, nconvd2)
1064        nminFD = ceiling((-maxFDarg+xgrid_prev(jj)) / tratio / incconv )
1065        nminFD = max (nminFD, -nconvd2)
1066 
1067 !      Calculate the Fermi-Dirac distrib at point xgrid_prev(jj)-tgrid(ii)*tratio
1068        expxo2 = exp (-half*(xgrid_prev(jj) - (nminFD)*incconv*tratio))
1069        expx22 = expxo2*expxo2
1070        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1071        resFD4 = tmpexpsum * tmpexpsum
1072        expxo2 = expxo2*expinc
1073        expx22 = expxo2*expxo2
1074        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1075        resFD3 = tmpexpsum * tmpexpsum
1076        expxo2 = expxo2*expinc
1077        expx22 = expxo2*expxo2
1078        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1079        resFD2 = tmpexpsum * tmpexpsum
1080        expxo2 = expxo2*expinc
1081        expx22 = expxo2*expxo2
1082        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1083        resFD1 = tmpexpsum * tmpexpsum
1084 
1085 !      core contribution to the integral with constant weight (48)
1086        tmpsmdfun = 0.0_dp
1087        do ii=nminFD+4,nmaxFD-4
1088          expxo2 = expxo2*expinc
1089 !        tmpexpsum = 1.0_dp / (expxo2 + 1.0_dp / expxo2 )
1090          expx22 = expxo2*expxo2
1091          tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1092          tmpsmdfun = tmpsmdfun + smd2(ii) * tmpexpsum * tmpexpsum
1093        end do
1094 
1095 !      Add on end contributions for show (both functions smd and smdFD are very small
1096        smdfun_prev(jj,1)=smdfun_prev(jj,1)       +48.0_dp*tmpsmdfun             &
1097 &       + 31.0_dp*smd2(nminFD+3)*resFD1 -11.0_dp*smd2(nminFD+2)*resFD2 &
1098 &       +  5.0_dp*smd2(nminFD+1)*resFD3 -       smd2(nminFD)*resFD4
1099 
1100        expxo2 = expxo2*expinc
1101        expx22 = expxo2*expxo2
1102        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1103        resFD1 = tmpexpsum * tmpexpsum
1104        expxo2 = expxo2*expinc
1105        expx22 = expxo2*expxo2
1106        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1107        resFD2 = tmpexpsum * tmpexpsum
1108        expxo2 = expxo2*expinc
1109        expx22 = expxo2*expxo2
1110        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1111        resFD3 = tmpexpsum * tmpexpsum
1112        expxo2 = expxo2*expinc
1113        expx22 = expxo2*expxo2
1114        tmpexpsum = expxo2 / (expx22 + 1.0_dp)
1115        resFD4 = tmpexpsum * tmpexpsum
1116 
1117 !      Contribution above
1118        smdfun_prev(jj,1)=smdfun_prev(jj,1)                                      &
1119 &       + 31.0_dp*smd2(nmaxFD-3)*resFD1  -11.0_dp*smd2(nmaxFD-2)*resFD2 &
1120 &       +  5.0_dp*smd2(nmaxFD-1)*resFD3  -       smd2(nmaxFD)*resFD4
1121        smdfun_prev(jj,1)=incconv*smdfun_prev(jj,1)/48.0_dp
1122      end do
1123 
1124      secmom = 0.0_dp
1125      thdmom = 0.0_dp
1126      resmom4 = xgrid_prev(-nptsdiv2_def  )*xgrid_prev(-nptsdiv2_def  )*smdfun_prev(-nptsdiv2_def  ,  1)
1127      resmom3 = xgrid_prev(-nptsdiv2_def+1)*xgrid_prev(-nptsdiv2_def+1)*smdfun_prev(-nptsdiv2_def+1,  1)
1128      resmom2 = xgrid_prev(-nptsdiv2_def+2)*xgrid_prev(-nptsdiv2_def+2)*smdfun_prev(-nptsdiv2_def+2,  1)
1129      resmom1 = xgrid_prev(-nptsdiv2_def+3)*xgrid_prev(-nptsdiv2_def+3)*smdfun_prev(-nptsdiv2_def+3,  1)
1130      resmom  = xgrid_prev(-nptsdiv2_def+4)*xgrid_prev(-nptsdiv2_def+4)*smdfun_prev(-nptsdiv2_def+4,  1)
1131      do ii=-nptsdiv2_def+4,nptsdiv2_def-1
1132        secmom = secmom +                                   &
1133 &       ( 17.0_dp*xgrid_prev(ii)  *xgrid_prev(ii)  *smdfun_prev(ii,  1)   &
1134 &       +42.0_dp*xgrid_prev(ii-1)*xgrid_prev(ii-1)*smdfun_prev(ii-1,1)   &
1135 &       -16.0_dp*xgrid_prev(ii-2)*xgrid_prev(ii-2)*smdfun_prev(ii-2,1)   &
1136 &       + 6.0_dp*xgrid_prev(ii-3)*xgrid_prev(ii-3)*smdfun_prev(ii-3,1)   &
1137 &       -       xgrid_prev(ii-4)*xgrid_prev(ii-4)*smdfun_prev(ii-4,1)  )
1138        resmom4 = resmom3
1139        resmom3 = resmom2
1140        resmom2 = resmom1
1141        resmom1 = resmom
1142        resmom  = xgrid_prev(ii+1)  *xgrid_prev(ii+1)  *smdfun_prev(ii+1,  1)
1143      end do
1144      secmom=increm * secmom / 48.0_dp
1145 !    thdmom=increm * thdmom / 48.0_dp
1146 !
1147 !    smom1  = second moment of delta in smd1(:)
1148 !    smom2  = second moment of delta in smd2(:)
1149 !
1150      smom1  = 0.0_dp
1151      smom2  = 0.0_dp
1152      tmom1  = 0.0_dp
1153      tmom2  = 0.0_dp
1154      do ii=-nconvd2+4,nconvd2
1155        smom1 = smom1+                                       &
1156 &       ( 17.0_dp*tgrid(ii)  *tgrid(ii)  *smd1(ii)         &
1157 &       +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd1(ii-1)       &
1158 &       -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd1(ii-2)       &
1159 &       + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd1(ii-3)       &
1160 &       -       tgrid(ii-4)*tgrid(ii-4)*smd1(ii-4)  )
1161        smom2 = smom2+                                       &
1162 &       ( 17.0_dp*tgrid(ii)  *tgrid(ii)  *smd2(ii  )     &
1163 &       +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd2(ii-1)     &
1164 &       -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd2(ii-2)     &
1165 &       + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd2(ii-3)     &
1166 &       -       tgrid(ii-4)*tgrid(ii-4)*smd2(ii-4)  )
1167      end do
1168      smom1 =incconv * smom1  / 48.0_dp
1169      smom2 =incconv * smom2  / 48.0_dp
1170 !    tmom1 =incconv * tmom1  / 48.0_dp
1171 !    tmom2 =incconv * tmom2  / 48.0_dp
1172 
1173      encorr =  smom2*tratio*tratio/secmom
1174 
1175 !    DEBUG
1176 !    write(std_out,*) ' getnel : debug, secmoms = ', secmom, smom1, smom2
1177 !    write(std_out,*) ' getnel : debug, thdmoms = ', thdmom, tmom1, tmom2
1178 !    write(std_out,*) ' getnel : encorr = ', encorr
1179 !    ENDDEBUG
1180 
1181      ABI_DEALLOCATE(tgrid)
1182      ABI_DEALLOCATE(smd1)
1183      ABI_DEALLOCATE(smd2)
1184 
1185    end if
1186 
1187 !  --------------------------------------------------------
1188 !  end of smearing function initialisation, dblsmr case
1189 !  --------------------------------------------------------
1190 
1191 
1192 !  Now that the smeared delta function has been initialized, compute the
1193 !  occupation function
1194    occfun_prev(-nptsdiv2_def,1)=zero
1195    entfun_prev(-nptsdiv2_def,1)=zero
1196 
1197 !  Different algorithms are possible, corresponding to the formulas
1198 !  (4.1.11), (4.1.12) and (4.1.14) in Numerical recipes (pp 107 and 108),
1199 !  with respective O(1/N2), O(1/N3), O(1/N4) convergence, where N is the
1200 !  number of points in the interval.
1201    algo=4
1202 
1203    if(algo==2)then
1204 
1205 !    Extended trapezoidal rule (4.1.11), taken in a cumulative way
1206      do ii=-nptsdiv2_def+1,nptsdiv2_def
1207        occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*(smdfun_prev(ii,1)+smdfun_prev(ii-1,1))/2.0_dp
1208        entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1209 &       ( -xgrid_prev(ii)*smdfun_prev(ii,1) -xgrid_prev(ii-1)*smdfun_prev(ii-1,1) )/2.0_dp
1210      end do
1211 
1212    else if(algo==3)then
1213 
1214 !    Derived from (4.1.12). Converges as O(1/N3).
1215 !    Do not care about the following points,
1216 !    where the values are extremely small anyway
1217      occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ;   entfun_prev(-nptsdiv2_def+1,1)=0.0_dp
1218      do ii=-nptsdiv2_def+2,nptsdiv2_def
1219        occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*&
1220 &       ( 5.0_dp*smdfun_prev(ii,1) + 8.0_dp*smdfun_prev(ii-1,1) - smdfun_prev(ii-2,1) )/12.0_dp
1221        entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1222 &       ( 5.0_dp*(-xgrid_prev(ii)  )*smdfun_prev(ii,1)  &
1223 &       +8.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)&
1224 &       -      (-xgrid_prev(ii-2))*smdfun_prev(ii-2,1) )/12.0_dp
1225      end do
1226 
1227    else if(algo==4)then
1228 
1229 !    Derived from (4.1.14)- alternative extended Simpsons rule. Converges as O(1/N4).
1230 !    Do not care about the following points,
1231 !    where the values are extremely small anyway
1232      occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ;   entfun_prev(-nptsdiv2_def+1,1)=0.0_dp
1233      occfun_prev(-nptsdiv2_def+2,1)=0.0_dp ;   entfun_prev(-nptsdiv2_def+2,1)=0.0_dp
1234      occfun_prev(-nptsdiv2_def+3,1)=0.0_dp ;   entfun_prev(-nptsdiv2_def+3,1)=0.0_dp
1235      do ii=-nptsdiv2_def+4,nptsdiv2_def
1236        occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*&
1237 &       ( 17.0_dp*smdfun_prev(ii,1)  &
1238 &       +42.0_dp*smdfun_prev(ii-1,1)&
1239 &       -16.0_dp*smdfun_prev(ii-2,1)&
1240 &       + 6.0_dp*smdfun_prev(ii-3,1)&
1241 &       -       smdfun_prev(ii-4,1) )/48.0_dp
1242        entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*&
1243 &       ( 17.0_dp*(-xgrid_prev(ii)  )*smdfun_prev(ii,1)  &
1244 &       +42.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)&
1245 &       -16.0_dp*(-xgrid_prev(ii-2))*smdfun_prev(ii-2,1)&
1246 &       + 6.0_dp*(-xgrid_prev(ii-3))*smdfun_prev(ii-3,1)&
1247 &       -       (-xgrid_prev(ii-4))*smdfun_prev(ii-4,1) )/48.0_dp
1248      end do
1249 
1250    end if ! End of choice between different algorithms for integration
1251 
1252 !  Normalize the functions (actually not needed for occopt=3..7)
1253    factor=1.0_dp/occfun_prev(nptsdiv2_def,1)
1254    smdfun_prev(:,1)=smdfun_prev(:,1)*factor
1255    occfun_prev(:,1)=occfun_prev(:,1)*factor
1256    entfun_prev(:,1)=entfun_prev(:,1)*factor
1257 
1258 !  Compute the cubic spline fitting of the smeared delta function
1259    yp1=0.0_dp ; ypn=0.0_dp
1260    workfun(:)=smdfun_prev(:,1)
1261    call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, smdder)
1262    smdfun_prev(:,2)=smdder(:)
1263 
1264 !  Compute the cubic spline fitting of the occupation function
1265    yp1=0.0_dp ; ypn=0.0_dp
1266    workfun(:)=occfun_prev(:,1)
1267    call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, occder)
1268    occfun_prev(:,2)=occder(:)
1269 
1270 !  Compute the cubic spline fitting of the entropy function
1271    yp1=0.0_dp ; ypn=0.0_dp
1272    workfun(:)=entfun_prev(:,1)
1273    call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, entder)
1274    entfun_prev(:,2)=entder(:)
1275 
1276    ABI_DEALLOCATE(entder)
1277    ABI_DEALLOCATE(occder)
1278    ABI_DEALLOCATE(smdder)
1279    ABI_DEALLOCATE(work)
1280    ABI_DEALLOCATE(workfun)
1281 
1282  end if
1283 
1284  if (abs(tphysel)<tol12) then
1285    tsmearinv=one/tsmear
1286  else
1287    tsmearinv=one/tphysel
1288  end if
1289 
1290  entfun(:,:) = entfun_prev(:,:)
1291  occfun(:,:) = occfun_prev(:,:)
1292  smdfun(:,:) = smdfun_prev(:,:)
1293  xgrid(:) = xgrid_prev(:)
1294  limit = limit_occ
1295  nptsdiv2 = nptsdiv2_def
1296 
1297 end subroutine init_occ_ent

m_occ/newocc [ Functions ]

[ Top ] [ m_occ ] [ Functions ]

NAME

 newocc

FUNCTION

 Compute new occupation numbers at each k point,
 from eigenenergies eigen, according to the
 smearing scheme defined by occopt (smearing width tsmear and
 physical temperature tphysel),
 with the constraint of number of valence electrons per unit cell nelect.

INPUTS

  eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), hartree
  spinmagntarget=if differ from -99.99_dp, fix the magnetic moment (in Bohr magneton)
  mband=maximum number of bands
  nband(nkpt)=number of bands at each k point
  nelect=number of electrons per unit cell
  nkpt=number of k points
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  occopt=option for occupancies
  prtvol=control print volume and debugging output
  stmbias=if non-zero, compute occupation numbers for STM (non-zero around the Fermi energy)
   NOTE : in this case, only fermie and occ are meaningful outputs.
  tphysel="physical" electronic temperature with FD occupations
  tsmear=smearing width (or temperature)
  wtk(nkpt)=k point weights

OUTPUT

  doccde(maxval(nband(:))*nkpt*nsppol)=derivative of occupancies wrt
           the energy for each band and k point
  entropy= entropy associated with the smearing (adimensional)
  fermie= fermi energy (Hartree)
  occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point

PARENTS

      gstate,m_ebands,respfn,vtorho

CHILDREN

      getnel,timab,wrtout

SOURCE

402 subroutine newocc(doccde,eigen,entropy,fermie,spinmagntarget,mband,nband,&
403 &  nelect,nkpt,nspinor,nsppol,occ,occopt,prtvol,stmbias,tphysel,tsmear,wtk)
404 
405 
406 !This section has been created automatically by the script Abilint (TD).
407 !Do not modify the following lines by hand.
408 #undef ABI_FUNC
409 #define ABI_FUNC 'newocc'
410 !End of the abilint section
411 
412  implicit none
413 
414 !Arguments ------------------------------------
415 !scalars
416  integer,intent(in) :: mband,nkpt,nspinor,nsppol,occopt,prtvol
417  real(dp),intent(in) :: spinmagntarget,nelect,stmbias,tphysel,tsmear
418  real(dp),intent(out) :: entropy,fermie
419 !arrays
420  integer,intent(in) :: nband(nkpt*nsppol)
421  real(dp),intent(in) :: eigen(mband*nkpt*nsppol),wtk(nkpt)
422  real(dp),intent(out) :: doccde(mband*nkpt*nsppol)
423  real(dp),intent(inout) :: occ(mband*nkpt*nsppol) !vz_i
424 
425 !Local variables-------------------------------
426  integer,parameter :: niter_max=120,nkpt_max=50,fake_unit=-666,option1=1
427  integer :: cnt,cnt2,cnt3,ib,ii,ik,ikpt,is,isppol,nkpt_eff
428  integer :: sign
429  integer,allocatable :: nbandt(:)
430  real(dp) :: dosdeltae,entropy_tmp,fermihi,fermilo,fermimid,fermimid_tmp
431  real(dp) :: fermi_biased,maxocc
432  real(dp) :: nelect_tmp,nelecthi,nelectlo,nelectmid,nelect_biased
433  real(dp) :: entropyt(2),fermihit(2),fermilot(2),fermimidt(2),nelecthit(2)
434  real(dp) :: nelectlot(2),nelectt(2),tsec(2)
435  real(dp),allocatable :: doccdet(:),eigent(:),occt(:)
436  character(len=500) :: message
437 
438 ! *************************************************************************
439 
440  DBG_ENTER("COLL")
441 
442  call timab(74,1,tsec)
443 
444 !Here treat the case where occopt does not correspond to a metallic occupation scheme
445  if (occopt<3 .or. occopt>8) then
446    write(message,'(a,i0,a)')' occopt= ',occopt,', a value not allowed in newocc.'
447    MSG_BUG(message)
448  end if ! test of metallic occopt
449 
450 !Check whether nband is a constant for all k point and spin-pol
451  do isppol=1,nsppol
452    do ikpt=1,nkpt
453      if(nband(ikpt+(isppol-1)*nkpt)/=nband(1))then
454        write(message,'(3a,i0,a,i0,a,i0,a)')&
455 &       'The number of bands must be the same for all k-points ',ch10,&
456 &       'but nband(1)= ',nband(1),' is different of nband(',&
457 &       ikpt+(isppol-1)*nkpt,') = ',nband(ikpt+(isppol-1)*nkpt),'.'
458        MSG_BUG(message)
459      end if
460    end do
461  end do
462 
463 !Check whether nelect is strictly positive
464  if(nelect<=zero)then
465    write(message,'(3a,es16.8,a)')&
466 &   'nelect must be a positive number, while ',ch10,&
467 &   'the calling routine ask nelect=',nelect,'.'
468    MSG_BUG(message)
469  end if
470 
471  maxocc=two/(nsppol*nspinor)
472 !Check whether nelect is coherent with nband (nband(1) is enough,
473 !since it was checked that nband is independent of k-point and spin-pol
474  if( nelect > nband(1)*nsppol*maxocc )then
475    write(message,'(3a,es16.8,a,i0,a,es16.8,a)' )&
476 &   'nelect must be smaller than nband*maxocc, while ',ch10,&
477 &   'the calling routine gives nelect= ',nelect,', nband= ',nband(1),' and maxocc= ',maxocc,'.'
478    MSG_BUG(message)
479  end if
480 
481 !Use bissection algorithm to find fermi energy
482 !This choice is due to the fact that it will always give sensible
483 !result (because the answer is bounded, even if the smearing function
484 !is non-monotonic (which is the case for occopt=4 or 6)
485 !Might speed up it, if needed !
486 
487 !Lowest and largest trial fermi energies, and corresponding number of electrons
488 !They are obtained from the smallest or largest eigenenergy, plus a range of
489 !energy that allows for complete occupation of all bands, or, on the opposite,
490 !for zero occupation of all bands (see getnel.f)
491  dosdeltae=zero  ! the DOS is not computed, with option=1
492  fermilo=minval(eigen(1:nband(1)*nkpt*nsppol))-6.001_dp*tsmear
493  if(occopt==3)fermilo=fermilo-24.0_dp*tsmear
494 
495  call getnel(doccde,dosdeltae,eigen,entropy,fermilo,maxocc,mband,nband,&
496 & nelectlo,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk)
497 
498  fermihi=maxval(eigen(1:nband(1)*nkpt*nsppol))+6.001_dp*tsmear
499 !safety value
500  fermihi = min(fermihi, 1.e6_dp)
501  if(occopt==3)fermihi=fermihi+24.0_dp*tsmear
502 
503  call getnel(doccde,dosdeltae,eigen,entropy,fermihi,maxocc,mband,nband,&
504 & nelecthi,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk)
505 
506 !Prepare fixed moment calculation
507  if(abs(spinmagntarget+99.99_dp)>1.0d-10)then
508    sign = 1
509    do is = 1, nsppol
510      fermihit(is) = fermihi
511      fermilot(is) = fermilo
512      nelectt(is) = half*(nelect+sign*spinmagntarget)
513      sign = -sign
514      nelecthit(is) = nelecthi
515      nelectlot(is) = nelectlo
516    end do
517  end if
518 
519 !If the target nelect is not between nelectlo and nelecthi, exit
520  if(nelect<nelectlo .or. nelect>nelecthi)then
521    write(message, '(a,a,a,a,d16.8,a,a,d16.8,a,d16.8,a,a,d16.8,a,d16.8)') ch10,&
522 &   ' newocc : ',ch10,&
523 &   '  The calling routine gives nelect=',nelect,ch10,&
524 &   '  The lowest bound is ',fermilo,', with nelect=',nelectlo,ch10,&
525 &   '  The highest bound is ',fermihi,', with nelect=',nelecthi
526    call wrtout(std_out,message,'COLL')
527 
528    write(message, '(11a)' )&
529 &   'In order to get the right number of electrons,',ch10,&
530 &   'it seems that the Fermi energy must be outside the range',ch10,&
531 &   'of eigenenergies, plus 6 or 30 times the smearing, which is strange.',ch10,&
532 &   'It might be that your number of bands (nband) corresponds to the strictly',ch10,&
533 &   'minimum number of bands to accomodate your electrons (so, OK for an insulator),',ch10,&
534 &   'while you are trying to describe a metal. In this case, increase nband, otherwise ...'
535    MSG_BUG(message)
536  end if
537 
538  if( abs(spinmagntarget+99.99_dp) < tol10 ) then
539 
540 !  Usual bissection loop
541    do ii=1,niter_max
542      fermimid=(fermihi+fermilo)*half
543 !    Produce nelectmid from fermimid
544      call getnel(doccde,dosdeltae,eigen,entropy,fermimid,maxocc,mband,nband,&
545 &     nelectmid,nkpt,nsppol,occ,occopt,option1,tphysel,tsmear,fake_unit,wtk)
546 !    write(std_out,'(a,es24.16,a,es24.16)' )' newocc : from fermi=',fermimid,', getnel gives nelect=',nelectmid
547      if(nelectmid>nelect*(one-tol14))then
548        fermihi=fermimid
549        nelecthi=nelectmid
550      end if
551      if(nelectmid<nelect*(one+tol14))then
552        fermilo=fermimid
553        nelectlo=nelectmid
554      end if
555      if( abs(nelecthi-nelectlo) <= nelect*two*tol14 .or. &
556 &     abs(fermihi-fermilo) <= tol14*abs(fermihi+fermilo) ) exit
557      if(ii==niter_max)then
558        write(message,'(a,i0,3a,es22.14,a,es22.14,a)')&
559 &       'It was not possible to find Fermi energy in ',niter_max,' bissections.',ch10,&
560 &       'nelecthi = ',nelecthi,', and nelectlo = ',nelectlo,'.'
561        MSG_BUG(message)
562      end if
563    end do ! End of bissection loop
564 
565    fermie=fermimid
566    write(message, '(a,f14.6,a,f14.6,a,a,i4)' ) &
567 &   ' newocc: new Fermi energy is ',fermie,' , with nelect=',nelectmid,ch10,&
568 &   '  Number of bissection calls =',ii
569    call wrtout(std_out,message,'COLL')
570 
571 !  Compute occupation numbers for prtstm/=0, close to the Fermi energy
572    if(abs(stmbias)>tol10)then
573      fermi_biased=fermie-stmbias
574      ABI_ALLOCATE(occt,(mband*nkpt*nsppol))
575      call getnel(doccde,dosdeltae,eigen,entropy,fermi_biased,maxocc,mband,nband,&
576 &     nelect_biased,nkpt,nsppol,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk)
577      occ(:)=occ(:)-occt(:)
578      nelect_biased=abs(nelectmid-nelect_biased)
579 !    Here, arrange to have globally positive occupation numbers,
580 !    irrespective of the stmbias sign
581      if(-stmbias>tol10)occ(:)=-occ(:)
582      ABI_DEALLOCATE(occt)
583 
584      write(message,'(a,f14.6)')' newocc : the number of electrons in the STM range is nelect_biased=',nelect_biased
585      call wrtout(std_out,message,'COLL')
586 
587    end if
588 
589  else ! Calculations with a specified moment
590 
591 !  Bissection loop
592    cnt2=0
593    cnt3=0
594    entropy=zero
595    maxocc=one
596    ABI_ALLOCATE(doccdet,(nkpt*mband))
597    ABI_ALLOCATE(eigent,(nkpt*mband))
598    ABI_ALLOCATE(occt,(nkpt*mband))
599    ABI_ALLOCATE(nbandt,(nkpt))
600 
601    do is = 1, nsppol
602      nelect_tmp = nelectt(is)
603      fermihi = fermihit(is)
604      fermilo = fermilot(is)
605      nelecthi = nelecthit(is)
606      nelectlo = nelectlot(is)
607 !    DEBUG
608 !    write(std_out,'(a,i1,3(f8.4,1x))') "Spin, N(spin):", is, nelect, fermihi, fermilo
609 !    write(std_out,'(a,2(f8.4,1x))') "Hi, lo:", nelecthi, nelectlo
610 !    ENDDEBUG
611 
612      do ii=1,niter_max
613        fermimid_tmp=(fermihi+fermilo)/2.0_dp
614 !      temporary arrays
615        cnt = 0
616        do ik = 1, nkpt
617          nbandt(ik) = mband
618          do ib = 1, mband
619            cnt = cnt + 1
620            eigent(cnt) = eigen(cnt+cnt2)
621            occt(cnt) = occ(cnt+cnt2)
622            doccdet(cnt) = doccde(cnt+cnt2)
623          end do
624        end do
625 
626 !      Produce nelectmid from fermimid
627        call getnel(doccdet,dosdeltae,eigent,entropy_tmp,fermimid_tmp,maxocc,mband,nbandt,&
628 &       nelectmid,nkpt,1,occt,occopt,option1,tphysel,tsmear,fake_unit,wtk)
629        entropyt(is) = entropy_tmp
630        fermimidt(is) = fermimid_tmp
631        fermimid = fermimidt(is)
632 !      temporary arrays
633        cnt = 0
634        do ik = 1, nkpt
635          do ib = 1, mband
636            cnt = cnt + 1
637            occ(cnt+cnt2) = occt(cnt)
638            doccde(cnt+cnt2) = doccdet(cnt)
639          end do
640        end do
641 
642 !      DEBUG
643 !      write(std_out,'(a,es24.16,a,es24.16)' )&
644 !      &    ' newocc : from fermi=',fermimid,', getnel gives nelect=',nelectmid
645 !      ENDDEBUG
646 
647        if(nelectmid>=nelect_tmp)then
648          fermihi=fermimid_tmp
649          nelecthi=nelectmid
650        else
651          fermilo=fermimid_tmp
652          nelectlo=nelectmid
653        end if
654        if( abs(nelecthi-nelectlo) <= 1.0d-13 .or. abs(fermihi-fermilo) <= 0.5d-14*abs(fermihi+fermilo) ) exit
655 
656        if(ii==niter_max)then
657          write(message,'(a,i3,3a,es22.14,a,es22.14,a)')&
658 &         '  It was not possible to find Fermi energy in ',niter_max,' bissections.',ch10,&
659 &         '  nelecthi=',nelecthi,', and nelectlo=',nelectlo,'.'
660          MSG_BUG(message)
661        end if
662      end do ! End of bissection loop
663 
664      cnt2 = cnt2 + nkpt*mband
665      entropy = entropy + entropyt(is)
666      fermie=fermimid
667      write(message, '(a,i2,a,f14.6,a,f14.6,a,a,i4)' ) &
668 &     ' newocc : new Fermi energy for spin ', is, ' is ',fermie,' , with nelect=',nelectmid,ch10,&
669 &     '  Number of bissection calls =',ii
670      call wrtout(std_out,message,'COLL')
671 
672    end do ! spin
673 
674    ABI_DEALLOCATE(doccdet)
675    ABI_DEALLOCATE(eigent)
676    ABI_DEALLOCATE(nbandt)
677    ABI_DEALLOCATE(occt)
678 
679  end if !  End of logical on fixed moment calculations
680 
681 !write(std_out,*) "kT*Entropy:", entropy*tsmear
682 
683  nkpt_eff=nkpt
684  if(prtvol==0)nkpt_eff=min(nkpt_max,nkpt)
685 
686  if(nsppol==1)then
687    write(message, '(a,i0,a)' ) &
688 &   ' newocc : computed new occ. numbers for occopt= ',occopt,' , spin-unpolarized case. '
689    call wrtout(std_out,message,'COLL')
690    do ikpt=1,nkpt_eff
691      write(message,'(a,i4,a)' ) ' k-point number ',ikpt,' :'
692      do ii=0,(nband(1)-1)/12
693        write(message,'(12f6.3)') &
694 &       occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
695        call wrtout(std_out,message,'COLL')
696      end do
697    end do
698    if(nkpt/=nkpt_eff)then
699      call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information','COLL')
700    end if
701 
702 !  DEBUG
703 !  write(message, '(a)' ) &
704 !  &   ' newocc : corresponding derivatives are '
705 !  call wrtout(std_out,message,'COLL')
706 !  do ikpt=1,nkpt_eff
707 !  write(message,'(a,i4,a)' ) ' k-point number ',ikpt,' :'
708 !  do ii=0,(nband(1)-1)/12
709 !  write(message,'(12f6.1)') &
710 !  &    doccde(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
711 !  call wrtout(std_out,message,'COLL')
712 !  end do
713 !  end do
714 !  if(nkpt/=nkpt_eff)then
715 !  write(message,'(a)') &
716 !  &    ' newocc : prtvol=0, stop printing more k-point informations'
717 !  call wrtout(std_out,message,'COLL')
718 !  end if
719 !  ENDDEBUG
720  else
721    write(message, '(a,i0,a,a)' ) &
722 &   ' newocc : computed new occupation numbers for occopt= ',occopt,&
723 &   ch10,'  (1) spin up   values  '
724    call wrtout(std_out,message,'COLL')
725    do ikpt=1,nkpt_eff
726      write(message,'(a,i0,a)' ) ' k-point number ',ikpt,':'
727      do ii=0,(nband(1)-1)/12
728        write(message,'(12f6.3)') &
729 &       occ(1+ii*12+(ikpt-1)*nband(1):min(12+ii*12,nband(1))+(ikpt-1)*nband(1))
730        call wrtout(std_out,message,'COLL')
731      end do
732    end do
733    if(nkpt/=nkpt_eff)then
734      call wrtout(std_out,'newocc: prtvol=0, stop printing more k-point information','COLL')
735    end if
736 
737    call wrtout(std_out,'  (2) spin down values  ','COLL')
738    do ikpt=1,nkpt_eff
739      do ii=0,(nband(1)-1)/12
740        write(message,'(12f6.3)') &
741 &       occ( 1+ii*12+(ikpt-1+nkpt)*nband(1) : &
742 &       min(12+ii*12,nband(1))+(ikpt-1+nkpt)*nband(1) )
743        call wrtout(std_out,message,'COLL')
744      end do
745    end do
746    if(nkpt/=nkpt_eff)then
747      call wrtout(std_out,' newocc: prtvol=0, stop printing more k-point information','COLL')
748    end if
749 
750  end if !  End choice based on spin
751 
752  call timab(74,2,tsec)
753 
754  DBG_EXIT("COLL")
755 
756 end subroutine newocc

m_occ/occeig [ Functions ]

[ Top ] [ m_occ ] [ Functions ]

NAME

 occeig

FUNCTION

 For each pair of active bands (m,n), generates ratios
 that depend on the difference between occupation numbers
 and eigenvalues.

INPUTS

  doccde_k(nband_k)=derivative of occ_k wrt the energy
  doccde_kq(nband_k)=derivative of occ_kq wrt the energy
  eig0_k(nband_k)=GS eigenvalues at k
  eig0_kq(nband_k)=GS eigenvalues at k+q
  nband_k=number of bands
  occopt=option for occupancies
  occ_k(nband_k)=occupation number for each band at k
  occ_kq(nband_k)=occupation number for each band at k+q

OUTPUT

  rocceig(nband_k,nband_k)$= (occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$,
   if this ratio has been attributed to the band n, 0.0_dp otherwise

NOTES

 Supposing the occupations numbers differ:
 if $abs(occ_{k,q}(m)) < abs(occ_k(n))$
  $rocceig(m,n)=(occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n)) $
 if $abs(occ_{k,q}(m))>abs(occ_k(n))$
  rocceig(m,n)=0.0_dp

 If the occupation numbers are close enough, then
 if the eigenvalues are also close, take the derivative
  $ rocceig(m,n)=\frac{1}{2}*docc/deig0 $
 otherwise,
  $ rocceig(m,n)=\frac{1}{2}*(occ_{k,q}(m)-occ_k(n))/(eig0_{k,q}(m)-eig0_k(n))$

PARENTS

      dfpt_nstpaw,dfpt_rhofermi,dfpt_vtorho

CHILDREN

SOURCE

1343 subroutine occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,occopt,occ_k,occ_kq,rocceig)
1344 
1345 
1346 !This section has been created automatically by the script Abilint (TD).
1347 !Do not modify the following lines by hand.
1348 #undef ABI_FUNC
1349 #define ABI_FUNC 'occeig'
1350 !End of the abilint section
1351 
1352  implicit none
1353 
1354 !Arguments ------------------------------------
1355 !scalars
1356  integer,intent(in) :: nband_k,occopt
1357 !arrays
1358  real(dp),intent(in) :: doccde_k(nband_k),doccde_kq(nband_k),eig0_k(nband_k)
1359  real(dp),intent(in) :: eig0_kq(nband_k),occ_k(nband_k),occ_kq(nband_k)
1360  real(dp),intent(out) :: rocceig(nband_k,nband_k)
1361 
1362 !Local variables-------------------------------
1363 !scalars
1364  integer :: ibandk,ibandkq
1365  real(dp) :: diffabsocc,diffeig,diffocc,ratio,sumabsocc
1366  character(len=500) :: message
1367 
1368 ! *************************************************************************
1369 
1370 !The parameter tol5 defines the treshhold for degeneracy, and the width of the step function
1371 
1372  rocceig(:,:)=0.0_dp
1373 
1374  do ibandk=1,nband_k
1375    do ibandkq=1,nband_k
1376      diffeig=eig0_kq(ibandkq)-eig0_k(ibandk)
1377      diffocc=occ_kq(ibandkq)-occ_k(ibandk)
1378 
1379      if( abs(diffeig) > tol5 ) then
1380        ratio=diffocc/diffeig
1381      else
1382        if(occopt<3)then
1383 !        In a non-metallic case, if the eigenvalues are degenerate,
1384 !        the occupation numbers must also be degenerate, in which
1385 !        case there is no contribution from this pair of bands
1386          if( abs(diffocc) > tol5 ) then
1387            write(message,'(a,a,a,a,a,a,a,2(a,i4,a,es16.6,a,es16.6,a,a),a)' ) &
1388 &           'In a non-metallic case (occopt<3), for a RF calculation,',ch10,&
1389 &           'if the eigenvalues are degenerate,',&
1390 &           ' the occupation numbers must also be degenerate.',ch10,&
1391 &           'However, the following pair of states gave :',ch10,&
1392 &           'k -state, band number',ibandk,', occ=',occ_k(ibandk),&
1393 &           'eigenvalue=',eig0_k(ibandk),',',ch10,&
1394 &           ' kq-state, band number',ibandkq,', occ=',occ_kq(ibandkq),&
1395 &           ', eigenvalue=',eig0_kq(ibandkq),'.',ch10,&
1396 &           'Action : change occopt, consistently, in GS and RF calculations.'
1397            MSG_ERROR(message)
1398          end if
1399          ratio=0.0_dp
1400        else
1401 !        In the metallic case, one can compute a better approximation of the
1402 !        ratio by using derivatives doccde
1403          ratio=0.5_dp*(doccde_kq(ibandkq)+doccde_k(ibandk))
1404 !        DEBUG
1405 !        write(std_out,*)' occeig : ibandkq,doccde_kq(ibandkq)',ibandkq,doccde_kq(ibandkq)
1406 !        write(std_out,*)'          ibandk ,doccde_k (ibandk )',ibandk,doccde_k(ibandk)
1407 !        ENDDEBUG
1408        end if
1409      end if
1410 
1411 !    Here, must pay attention to the smallness of some coefficient
1412      diffabsocc=abs(occ_k(ibandk))-abs(occ_kq(ibandkq))
1413      sumabsocc=abs(occ_k(ibandk))+abs(occ_kq(ibandkq))
1414      if(sumabsocc>tol8)then
1415        if( diffabsocc > sumabsocc*tol5 ) then
1416          rocceig(ibandkq,ibandk)=ratio
1417        else if ( diffabsocc >= -sumabsocc*tol5 ) then
1418          rocceig(ibandkq,ibandk)=0.5_dp*ratio
1419        else
1420          rocceig(ibandkq,ibandk)=0.0_dp
1421        end if
1422      end if
1423 
1424    end do ! ibandkq
1425  end do ! ibandk
1426 
1427 end subroutine occeig

m_occ/pareigocc [ Functions ]

[ Top ] [ m_occ ] [ Functions ]

NAME

 pareigocc

FUNCTION

 This subroutine transmit to all processors, using MPI:
   - the eigenvalues and,
   - if ground-state, the occupation numbers
     (In fact, in the present status of the routine,
     occupation numbers are NOT transmitted)
     transmit_occ = 2 is used in case the occ should be transmitted.
     Yet the code is not already written.

INPUTS

  formeig=format of eigenvalues (0 for GS, 1 for RF)
  localrdwf=(for parallel case) if 1, the eig and occ initial values
            are local to each machine, if 0, they are on proc me=0.
  mband=maximum number of bands of the output wavefunctions
  mpi_enreg=information about MPI parallelization
  nband(nkpt*nsppol)=desired number of bands at each k point
  nkpt=number of k points
  nsppol=1 for unpolarized, 2 for spin-polarized, output wf file processors,
         Warning : defined only when paralbd=1
  transmit_occ/=2 transmit only eigenvalues, =2 for transmission of occ also
         (yet transmit_occ=2 is not safe or finished at all)

OUTPUT

  (see side effects)

SIDE EFFECTS

  eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
  occ(mband*nkpt*nsppol)=occupation (input or init to 0.0)  NOT USED NOW

NOTES

 * The case paralbd=1 with formeig=0 is implemented, but not yet used.

 * The transmission of occ is not activated yet !

 * The routine takes the eigenvalues in the eigen array on one of the
   processors that possess the wavefunctions, and transmit it to all procs.
   If localrdwf==0, me=0 has the full array at start,
   If localrdwf==1, the transfer might be more complex.

 * This routine should not be used for RF wavefunctions, since
   it does not treat the eigenvalues as a matrix.

PARENTS

      newkpt,wfsinp

CHILDREN

      timab,xmpi_bcast,xmpi_sum

SOURCE

1616 subroutine pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,transmit_occ)
1617 
1618 
1619 !This section has been created automatically by the script Abilint (TD).
1620 !Do not modify the following lines by hand.
1621 #undef ABI_FUNC
1622 #define ABI_FUNC 'pareigocc'
1623 !End of the abilint section
1624 
1625  implicit none
1626 
1627 !Arguments ------------------------------------
1628 !scalars
1629  integer,intent(in) :: formeig,localrdwf,mband,nkpt,nsppol,transmit_occ
1630  type(MPI_type),intent(in) :: mpi_enreg
1631 !arrays
1632  integer,intent(in) :: nband(nkpt*nsppol)
1633  real(dp),intent(inout) :: eigen(mband*(2*mband)**formeig*nkpt*nsppol)
1634  real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
1635 
1636 !Local variables-------------------------------
1637 !scalars
1638  integer :: band_index,iband,ierr,ikpt,isppol,me,nbks,spaceComm
1639  !character(len=500) :: message
1640 !arrays
1641  real(dp) :: tsec(2)
1642  real(dp),allocatable :: buffer1(:),buffer2(:)
1643 
1644 ! *************************************************************************
1645 
1646  if(xmpi_paral==1)then
1647 
1648 !  Init mpi_comm
1649    spaceComm=mpi_enreg%comm_cell
1650    if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
1651    if(mpi_enreg%paral_hf==1) spaceComm=mpi_enreg%comm_kpt
1652 !  Init me
1653    me=mpi_enreg%me_kpt
1654 
1655    if(localrdwf==0)then
1656      call xmpi_bcast(eigen,0,spaceComm,ierr)
1657 
1658    else if(localrdwf==1)then
1659 
1660 !    Prepare transmission of eigen (and occ)
1661      ABI_ALLOCATE(buffer1,(2*mband**(formeig+1)*nkpt*nsppol))
1662      ABI_ALLOCATE(buffer2,(2*mband**(formeig+1)*nkpt*nsppol))
1663      buffer1(:)=zero
1664      buffer2(:)=zero
1665 
1666      band_index=0
1667      do isppol=1,nsppol
1668        do ikpt=1,nkpt
1669          nbks=nband(ikpt+(isppol-1)*nkpt)
1670 
1671          if(mpi_enreg%paralbd==0)then
1672 
1673            if(formeig==0)then
1674              buffer1(2*band_index+1:2*band_index+nbks)=&
1675 &             eigen(band_index+1:band_index+nbks)
1676              if(transmit_occ==2) then
1677                buffer1(2*band_index+nbks+1:2*band_index+2*nbks)=&
1678 &               occ(band_index+1:band_index+nbks)
1679              end if
1680              band_index=band_index+nbks
1681            else if(formeig==1)then
1682              buffer1(band_index+1:band_index+2*nbks**2)=&
1683 &             eigen(band_index+1:band_index+2*nbks**2)
1684              band_index=band_index+2*nbks**2
1685            end if
1686 
1687          else if(mpi_enreg%paralbd==1)then
1688 
1689 !          Skip this k-point if not the proper processor
1690            if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbks,isppol,me)) then
1691              if(formeig==0) then
1692                band_index=band_index+nbks
1693              else
1694                band_index=band_index+2*nbks**2
1695              end if
1696              cycle
1697            end if
1698 !          Loop on bands
1699            do iband=1,nbks
1700              if(mpi_enreg%proc_distrb(ikpt, iband,isppol) /= me)cycle
1701              if(formeig==0)then
1702                buffer1(2*band_index+iband)=eigen(band_index+iband)
1703 !              if(transmit_occ==2) then
1704 !              buffer1(2*band_index+iband+nbdks)=occ(band_index+iband)
1705 !              end if
1706              else if (formeig==1)then
1707                buffer1(band_index+(iband-1)*2*nbks+1: &
1708 &               band_index+(iband-1)*2*nbks+2*nbks)=&
1709 &               eigen(band_index+(iband-1)*2*nbks+1: &
1710 &               band_index+(iband-1)*2*nbks+2*nbks)
1711              end if
1712            end do
1713            if(formeig==0)then
1714              band_index=band_index+nbks
1715            else
1716              band_index=band_index+2*nbks**2
1717            end if
1718          end if
1719 
1720        end do
1721      end do
1722 
1723 !    Build sum of everything
1724      call timab(48,1,tsec)
1725 !    call wrtout(std_out,' pareigocc : MPI_ALLREDUCE','COLL')
1726      if(formeig==0)band_index=band_index*2
1727      call xmpi_sum(buffer1,buffer2,band_index,spaceComm,ierr)
1728      call timab(48,2,tsec)
1729 
1730      band_index=0
1731      do isppol=1,nsppol
1732        do ikpt=1,nkpt
1733          nbks=nband(ikpt+(isppol-1)*nkpt)
1734          if(formeig==0)then
1735            eigen(band_index+1:band_index+nbks)=&
1736 &           buffer2(2*band_index+1:2*band_index+nbks)
1737            if(transmit_occ==2) then
1738              occ(band_index+1:band_index+nbks)=&
1739 &             buffer2(2*band_index+nbks+1:2*band_index+2*nbks)
1740            end if
1741            band_index=band_index+nbks
1742          else if(formeig==1)then
1743            eigen(band_index+1:band_index+2*nbks**2)=&
1744 &           buffer1(band_index+1:band_index+2*nbks**2)
1745            band_index=band_index+2*nbks**2
1746          end if
1747        end do
1748      end do
1749 
1750      ABI_DEALLOCATE(buffer1)
1751      ABI_DEALLOCATE(buffer2)
1752 
1753    end if
1754 
1755  end if
1756 
1757 end subroutine pareigocc