TABLE OF CONTENTS


ABINIT/newton [ Functions ]

[ Top ] [ Functions ]

NAME

 newton

FUNCTION

  Compute root of a function with newton methods (newton/halley)

COPYRIGHT

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

INPUTS

  x_input      : input of x             
  x_precision  : required precision on x
  max_iter     : maximum number of iterations
  opt_noninter

OUTPUT

  f_precision  : output precision on function F
  ierr_hh      : different from zero if an error occurs

PARENTS

      fermi_green

CHILDREN

      compute_green,integrate_green

SOURCE

 33 #if defined HAVE_CONFIG_H
 34 #include "config.h"
 35 #endif
 36 
 37 
 38 #include "abi_common.h"
 39 
 40 
 41 subroutine newton(cryst_struc,green,paw_dmft,pawang,self,&
 42 & x_input,x_precision,max_iter,f_precision,ierr_hh,opt_noninter,opt_algo)
 43 
 44  use defs_basis
 45  use defs_abitypes
 46  use m_profiling_abi
 47  use m_errors
 48 
 49  use m_pawang,    only : pawang_type
 50  use m_crystal,   only : crystal_t
 51  use m_green,     only : green_type
 52  use m_paw_dmft,  only: paw_dmft_type
 53  use m_self,      only : self_type
 54 
 55 !This section has been created automatically by the script Abilint (TD).
 56 !Do not modify the following lines by hand.
 57 #undef ABI_FUNC
 58 #define ABI_FUNC 'newton'
 59  use interfaces_14_hidewrite
 60 !End of the abilint section
 61 
 62  implicit none
 63 !Arguments ------------------------------------
 64 !scalars
 65  type(crystal_t),intent(in) :: cryst_struc
 66  type(green_type),intent(inout) :: green
 67  !type(MPI_type), intent(in) :: mpi_enreg
 68  type(paw_dmft_type), intent(inout) :: paw_dmft
 69  type(pawang_type),intent(in) :: pawang
 70  type(self_type), intent(inout) :: self
 71  integer,intent(in) :: opt_noninter,max_iter
 72  integer,intent(out) :: ierr_hh
 73  real(dp),intent(inout) :: x_input,x_precision
 74  real(dp),intent(inout) :: f_precision
 75  real(dp),intent(in), optional :: opt_algo
 76 !Local variables-------------------------------
 77  integer iter
 78  real(dp) Fx,Fxprime,Fxdouble,xold,x_before,Fxoptimum
 79  real(dp) :: nb_elec_x
 80  integer option,optalgo
 81  logical l_minus,l_plus
 82  real(dp) :: x_minus,x_plus,x_optimum
 83  character(len=500) :: message
 84 ! *********************************************************************
 85  x_minus=-10_dp
 86  x_plus=-11_dp
 87  xold=-12_dp
 88 
 89  if(present(opt_algo)) then
 90    optalgo=opt_algo
 91  else
 92    optalgo=1 ! newton
 93  end if
 94 
 95  x_input=paw_dmft%fermie
 96  ierr_hh=0
 97  option =2  ! Halley method
 98  option =1  ! Newton method
 99 
100 !write(std_out,*) "ldaprint",opt_noninter
101 
102 !--- Start of iterations
103  write(message,'(a,3a)') "     Fermi level   Charge    Difference"
104  call wrtout(std_out,message,'COLL')
105 !do iter=1, 40
106 !x_input=float(iter)/100_dp
107 !call function_and_deriv(cryst_struc,f_precision,green,iter,mpi_enreg,paw_dmft,pawang,self,&
108 !&  x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
109 !write(std_out,*) x_input,Fx
110 !enddo
111 !call leave_new('COLL')
112 
113  l_minus=.false.
114  l_plus=.false.
115  Fxoptimum=1_dp
116 !========================================
117 !start iteration to find fermi level
118 !========================================
119  do iter=1, max_iter
120 
121 !  ========================================
122 !  If zero is located between two values: apply newton method or dichotomy
123 !  ========================================
124    if(l_minus.and.l_plus) then
125 
126 !    ==============================================
127 !    Compute the function and derivatives for newton
128 !    ==============================================
129      call function_and_deriv(cryst_struc,f_precision,green,iter,paw_dmft,pawang,self,&
130 &     x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
131 
132 !    Apply stop criterion on Fx
133      if(abs(Fx) < f_precision) then
134 !      write(message,'(a,2f12.6)') "Fx,f_precision",Fx,f_precision
135 !      call wrtout(std_out,message,'COLL')
136        x_precision=x_input-xold
137        return
138      end if 
139      if(iter==max_iter) then
140        write(message,'(a,2f12.6)') "   Fermi level could not be found"
141        call wrtout(std_out,message,'COLL')
142        x_input=x_optimum
143        ierr_hh=-123
144        return
145      end if
146 
147 !    Cannot divide by Fxprime if too small
148      if(abs(Fxprime) .le. 1.e-15)then
149        ierr_hh=-314
150        write(message,'(a,f12.7)') "Fxprime=",Fxprime
151        call wrtout(std_out,message,'COLL')
152        return
153      end if
154 
155      x_precision=x_input-xold
156 
157 !    ==============================================
158 !    Newton/Halley's  formula for next iteration
159 !    ==============================================
160      xold=x_input
161      if(option==1) x_input=x_input-Fx/Fxprime
162      if(option==2) x_input=x_input-2*Fx*Fxprime/(2*Fxprime**2-Fx*Fxdouble)
163 
164 !    ==============================================
165 !    If newton does not work well, use dichotomy.
166 !    ==============================================
167      if(x_input<x_minus.or.x_input>x_plus) then
168        call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
169 &       Fx,opt_noninter,nb_elec_x,xold)
170        write(message,'(a,3f12.6)') " ---",x_input,nb_elec_x,Fx
171        call wrtout(std_out,message,'COLL')
172        if(Fx>0) then
173          x_plus=xold
174        else if(Fx<0) then
175          x_minus=xold
176        end if
177        x_input=(x_plus+x_minus)/2.d0
178        
179      end if
180 !    write(std_out,'(a,2f12.6)') " Q(xold) and dQ/dx=",Fx,Fxprime
181 !    write(std_out,'(a,f12.6)') " =>  new Fermi level",x_input
182 !    ========================================
183 !    Locate zero between two values
184 !    ========================================
185    else 
186      call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
187 &     Fx,opt_noninter,nb_elec_x,x_input)
188      write(message,'(a,3f12.6)') "  --",x_input,nb_elec_x,Fx
189 !    Possible improvement for large systems, removed temporarely for
190 !    automatic tests: more study is necessary: might worsen the convergency
191 !    if(iter==1) then
192 !    f_precision=max(abs(Fx/50),f_precision)
193 !    write(message,'(a,4x,a,e12.6)') ch10," Precision required changed to:",f_precision
194 !    call wrtout(std_out,message,'COLL')
195 !    endif
196      call wrtout(std_out,message,'COLL')
197      if(Fx>0) then
198        l_plus=.true.
199        x_plus=x_input
200        x_input=x_input-0.02
201      else if(Fx<0) then
202        l_minus=.true.
203        x_minus=x_input
204        x_input=x_input+0.02
205      end if
206      
207    end if
208 
209    if(abs(Fx)<abs(Fxoptimum)) then
210      Fxoptimum=Fx
211      x_optimum=x_input
212    end if
213 
214 
215 
216 !  if(myid==master) then
217 !  write(std_out,'(a,i4,3f12.6)') "i,xnew,F,Fprime",i,x_input,Fx,Fxprime
218 !  endif
219 
220 
221 !  ! Apply stop criterion on x
222 !  if(abs(x_input-xold) .le. x_input*x_precision) then
223 !  !    write(std_out,'(a,4f12.6)') "x_input-xold, x_precision*x_input   "&
224 !  !&    ,x_input-xold,x_precision*x_input,x_precision
225 !  f_precision=Fx
226 !  return
227 !  endif
228 
229  end do
230 !--- End of iterations
231 
232 
233  ierr_hh=1
234  return
235 
236  CONTAINS  !========================================================================================
237 !-----------------------------------------------------------------------

newton/compute_nb_elec [ Functions ]

[ Top ] [ newton ] [ Functions ]

NAME

 compute_nb_elec

FUNCTION

  Compute nb of electrons as a function of Fermi level

INPUTS

 fermie       : input of energy 
 opt_noninter

 OUTPUTS
 Fx           : Value of F(x).
 nb_elec_x    : Number of electrons for the value of x

PARENTS

      newton

CHILDREN

      compute_green,integrate_green

SOURCE

377 subroutine compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
378 &  Fx,opt_noninter,nb_elec_x,fermie)
379 
380  use m_profiling_abi
381 
382  use defs_basis
383  use defs_abitypes
384  use m_errors
385  use m_crystal, only : crystal_t
386  use m_green, only : green_type,compute_green,integrate_green
387  use m_paw_dmft, only: paw_dmft_type
388  use m_self, only : self_type
389 
390 !This section has been created automatically by the script Abilint (TD).
391 !Do not modify the following lines by hand.
392 #undef ABI_FUNC
393 #define ABI_FUNC 'compute_nb_elec'
394 !End of the abilint section
395 
396  implicit none
397 !Arguments ------------------------------------
398 !scalars
399  type(crystal_t),intent(in) :: cryst_struc
400  type(green_type),intent(inout) :: green
401  !type(MPI_type), intent(in) :: mpi_enreg
402  type(paw_dmft_type), intent(inout) :: paw_dmft
403  type(pawang_type),intent(in) :: pawang
404  type(self_type), intent(inout) :: self
405  integer,intent(in) :: opt_noninter
406  real(dp),intent(in)  :: fermie
407  real(dp),intent(out) :: Fx,nb_elec_x
408 ! *********************************************************************
409    paw_dmft%fermie=fermie
410    call compute_green(cryst_struc,green,paw_dmft,pawang,0,self,opt_self=1,&
411 &   opt_nonxsum=1,opt_nonxsum2=1)
412    call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=0,&
413 &   opt_ksloc=-1) !,opt_nonxsum=1)
414 !  opt_ksloc=-1, compute total charge
415    nb_elec_x=green%charge_ks
416    Fx=nb_elec_x-paw_dmft%nelectval
417 
418    if(opt_noninter==1) then
419    end if
420  end subroutine compute_nb_elec

newton/function_and_deriv [ Functions ]

[ Top ] [ newton ] [ Functions ]

NAME

  function_and_deriv 

FUNCTION

  Compute value of a function and its numerical derivatives

INPUTS

  x_input      : input of x             
  option       : if 1 compute only first derivative
                 if 2 compute the two first derivatives.
  opt_noninter

 OUTPUTS
  Fx           : Value of F(x)
  Fxprime      : Value of F'(x)
  Fxdouble     : Value of F''(x)

PARENTS

      newton

CHILDREN

      compute_green,integrate_green

SOURCE

265 subroutine function_and_deriv(cryst_struc,f_precision,green,iter,paw_dmft,pawang,self&
266 & ,x_input,x_old,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option)
267 
268  use m_profiling_abi
269 
270  use defs_basis
271  use defs_abitypes
272  use m_errors
273  use m_crystal, only : crystal_t
274  use m_green, only : green_type
275  use m_paw_dmft, only: paw_dmft_type
276  use m_self, only : self_type
277 
278 !This section has been created automatically by the script Abilint (TD).
279 !Do not modify the following lines by hand.
280 #undef ABI_FUNC
281 #define ABI_FUNC 'function_and_deriv'
282  use interfaces_14_hidewrite
283 !End of the abilint section
284 
285  implicit none
286 !Arguments ------------------------------------
287 !scalars
288  type(crystal_t),intent(in) :: cryst_struc
289  type(green_type),intent(inout) :: green
290  !type(MPI_type), intent(in) :: mpi_enreg
291  type(paw_dmft_type), intent(inout) :: paw_dmft
292  type(pawang_type),intent(in) :: pawang
293  type(self_type), intent(inout) :: self
294  integer,intent(in) :: iter,opt_noninter,option
295  real(dp),intent(inout)  :: f_precision,x_input,x_precision
296  real(dp),intent(out) :: Fx,Fxprime,Fxdouble
297  real(dp),intent(inout) :: x_old
298 !Local variables-------------------------------
299  real(dp) :: deltax,nb_elec_x,Fxmoins,Fxplus,xmoins,xplus,x0
300  character(len=500) :: message
301 ! *********************************************************************
302 
303 !  choose deltax: for numeric evaluation of derivative
304    if(iter==1) then
305 !    deltax=0.02
306    end if
307 !  deltax=max((x_input-x_old)/10.d0,min(0.00001_dp,x_precision/100_dp))
308    deltax=min(0.00001_dp,x_precision/100_dp)  ! small but efficient
309 !  endif
310 !  write(std_out,*) "iter,x_input,deltax",iter,x_input,deltax
311    x0=x_input
312    xmoins=x0-deltax
313    xplus=x0+deltax
314 
315    call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
316 &   Fx,opt_noninter,nb_elec_x,x0)
317 
318    write(message,'(a,3f12.6)') "  - ",x0,nb_elec_x,Fx
319    call wrtout(std_out,message,'COLL')
320 !  write(std_out,*) "Fx", Fx
321    if(abs(Fx)<f_precision) return
322 
323    call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
324 &   Fxplus,opt_noninter,nb_elec_x,xplus)
325 
326    write(message,'(a,3f12.6)') "  - ",xplus,nb_elec_x,Fxplus
327    call wrtout(std_out,message,'COLL')
328 
329    if(option==2) then
330      call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,&
331 &     Fxmoins,opt_noninter,nb_elec_x,xmoins)
332 
333      write(message,'(a,3f12.6)') "  - ",xmoins,nb_elec_x,Fxmoins
334      call wrtout(std_out,message,'COLL')
335    end if
336 
337    if(option==1) then
338      Fxprime=(Fxplus-Fx)/deltax
339    else if (option==2) then
340      Fxprime=(Fxplus-Fxmoins)/(2*deltax)
341      Fxdouble=(Fxplus+Fxmoins-2*Fx)/(deltax**2)
342    end if
343 !  write(std_out,*) "after computation of Fxprime",myid
344    if(Fxprime<zero) then
345      write(message,'(a,f12.6)') "  Warning: slope of charge versus fermi level is negative !", Fxprime
346      call wrtout(std_out,message,'COLL')
347    end if
348    x_old=x_input
349 
350    return
351  end subroutine function_and_deriv