TABLE OF CONTENTS
ABINIT/newton [ 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