TABLE OF CONTENTS
ABINIT/init_occ_ent [ Functions ]
NAME
init_occ_ent
FUNCTION
COPYRIGHT
Copyright (C) 2009-2017 ABINIT group (the_author) 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
argin(sizein)=description
OUTPUT
argout(sizeout)=description
PARENTS
getnel
CHILDREN
spline
SOURCE
28 #if defined HAVE_CONFIG_H 29 #include "config.h" 30 #endif 31 32 #include "abi_common.h" 33 34 subroutine init_occ_ent(entfun,limit,nptsdiv2,occfun,occopt,option,smdfun,tphysel,tsmear,tsmearinv,xgrid) 35 36 use defs_basis 37 use m_splines 38 use m_profiling_abi 39 use m_errors 40 41 !This section has been created automatically by the script Abilint (TD). 42 !Do not modify the following lines by hand. 43 #undef ABI_FUNC 44 #define ABI_FUNC 'init_occ_ent' 45 !End of the abilint section 46 47 implicit none 48 49 !Arguments ------------------------------------ 50 !scalars 51 integer,intent(in) :: occopt,option 52 real(dp),intent(in) :: tphysel,tsmear 53 integer,intent(inout) :: nptsdiv2 54 real(dp),intent(out) :: limit,tsmearinv 55 real(dp),intent(inout) :: entfun(-nptsdiv2:nptsdiv2,2),occfun(-nptsdiv2:nptsdiv2,2) 56 real(dp),intent(inout) :: smdfun(-nptsdiv2:nptsdiv2,2),xgrid(-nptsdiv2:nptsdiv2) 57 58 59 !Local variables------------------------------- 60 !scalars 61 integer :: algo,ii,jj,nconvd2 62 integer :: nmaxFD,nminFD 63 integer,parameter :: nptsdiv2_def=6000 64 integer,save :: dblsmr,occopt_prev=-9999 65 real(dp),parameter :: maxFDarg=500.0_dp 66 real(dp),save :: convlim,incconv,limit_occ,tphysel_prev=-9999,tsmear_prev=-9999 67 real(dp) :: aa,dsqrpi,encorr,factor 68 real(dp) :: expinc,expx22,expxo2,gauss,increm 69 real(dp) :: resFD1,resFD2,resFD3,resFD4,resmom,resmom1,resmom2 70 real(dp) :: resmom3,resmom4,secmom,smom1,smom2,thdmom,tmom1,tmom2,tmpexpsum 71 real(dp) :: tmpsmdfun,tratio,tt,xx,yp1,ypn 72 character(len=500) :: message 73 !arrays 74 real(dp),save :: entfun_prev(-nptsdiv2_def:nptsdiv2_def,2),occfun_prev(-nptsdiv2_def:nptsdiv2_def,2) 75 real(dp),save :: smdfun_prev(-nptsdiv2_def:nptsdiv2_def,2),xgrid_prev(-nptsdiv2_def:nptsdiv2_def) 76 real(dp),allocatable :: entder(:),occder(:),smd1(:),smd2(:) 77 real(dp),allocatable :: smdder(:),tgrid(:),work(:),workfun(:) 78 79 ! ************************************************************************* 80 81 !Initialize the occupation function and generalized entropy function, 82 !at the beginning, or if occopt changed 83 84 if(option==-1)then 85 nptsdiv2 = nptsdiv2_def 86 return 87 end if 88 89 90 if(occopt_prev/=occopt .or. & 91 & abs(tsmear_prev-tsmear) >tol12 .or. & 92 & abs(tphysel_prev-tphysel)>tol12 ) then 93 ! write(std_out,*) 'INIT_OCC_ENT CHANGE ..........' 94 occopt_prev=occopt 95 tsmear_prev=tsmear 96 tphysel_prev=tphysel 97 98 ! Check whether input values of tphysel tsmear and occopt are consistent 99 dblsmr = 0 100 if (abs(tphysel)>tol12) then 101 ! Use re-smearing scheme 102 if (abs(tsmear)>tol12) then 103 dblsmr = 1 104 ! Use FD occupations (one smearing) only with "physical" temperature tphysel 105 else if (occopt /= 3) then 106 write(message, '(a,i6,a)' )' tphysel /= 0, tsmear == 0, but occopt is not = 3, but ',occopt,'.' 107 MSG_ERROR(message) 108 end if 109 end if 110 ! write(std_out,*) 'getnel : input read.' 111 ! write(std_out,*) ' dblsmr = ', dblsmr 112 ! write(std_out,*) ' tphysel, tsmear = ', tphysel, tsmear 113 114 ABI_ALLOCATE(entder,(-nptsdiv2_def:nptsdiv2_def)) 115 ABI_ALLOCATE(occder,(-nptsdiv2_def:nptsdiv2_def)) 116 ABI_ALLOCATE(smdder,(-nptsdiv2_def:nptsdiv2_def)) 117 ABI_ALLOCATE(workfun,(-nptsdiv2_def:nptsdiv2_def)) 118 ABI_ALLOCATE(work,(-nptsdiv2_def:nptsdiv2_def)) 119 120 ! Prepare the points on the grid 121 ! limit is the value of the argument that will give 0.0 or 1.0 , with 122 ! less than about 1.0d-15 error for 4<=occopt<=8, and less than about 1.0d-12 123 ! error for occopt==3. It is not worth to compute the function beyond 124 ! that point. Even with a less severe requirement, it is significantly 125 ! larger for occopt==3, with an exponential 126 ! tail, than for the other occupation functions, with a Gaussian tail. 127 ! Note that these values are useful in newocc.f also. 128 limit_occ=6.0_dp 129 if(occopt==3)limit_occ=30.0_dp 130 if(dblsmr /= 0) then 131 tratio = tsmear / tphysel 132 limit_occ=30.0_dp + 6.0_dp*tratio 133 end if 134 135 ! With nptsdiv2_def=6000 (thus increm=0.001 for 4<=occopt<=8, 136 ! and increm=0.005 for occopt==3, the O(1/N4) algorithm gives 1.0d-12 137 ! accuracy on the stored values occfun and entfun. These, together 138 ! with smdfun and xgrid_prev, need permanently about 0.67 MB, which is affordable. 139 increm=limit_occ/nptsdiv2_def 140 do ii=-nptsdiv2_def,nptsdiv2_def 141 xgrid_prev(ii)=ii*increm 142 end do 143 144 ! --------------------------------------------------------- 145 ! Ordinary (unique) smearing function 146 ! --------------------------------------------------------- 147 if (dblsmr == 0) then 148 149 ! Compute the unnormalized smeared delta function between -limit_occ and +limit_occ 150 ! (well, they are actually normalized ...) 151 if(occopt==3)then 152 153 ! Fermi-Dirac 154 do ii=0,nptsdiv2_def 155 xx=xgrid_prev(ii) 156 smdfun_prev( ii,1)=0.25_dp/(cosh(xx/2.0_dp)**2) 157 smdfun_prev(-ii,1)=smdfun_prev(ii,1) 158 end do 159 160 else if(occopt==4 .or. occopt==5)then 161 162 ! Cold smearing of Marzari, two values of the "a" parameter being possible 163 ! first value gives minimization of the bump 164 if(occopt==4)aa=-.5634 165 ! second value gives monotonic occupation function 166 if(occopt==5)aa=-.8165 167 168 dsqrpi=1.0_dp/sqrt(pi) 169 do ii=0,nptsdiv2_def 170 xx=xgrid_prev(ii) 171 gauss=dsqrpi*exp(-xx**2) 172 smdfun_prev( ii,1)=gauss*(1.5_dp+xx*(-aa*1.5_dp+xx*(-1.0_dp+aa*xx))) 173 smdfun_prev(-ii,1)=gauss*(1.5_dp+xx*( aa*1.5_dp+xx*(-1.0_dp-aa*xx))) 174 end do 175 176 else if(occopt==6)then 177 178 ! First order Hermite-Gaussian of Paxton and Methfessel 179 dsqrpi=1.0_dp/sqrt(pi) 180 do ii=0,nptsdiv2_def 181 xx=xgrid_prev(ii) 182 smdfun_prev( ii,1)=dsqrpi*(1.5_dp-xx**2)*exp(-xx**2) 183 smdfun_prev(-ii,1)=smdfun_prev(ii,1) 184 end do 185 186 else if(occopt==7)then 187 188 ! Gaussian smearing 189 dsqrpi=1.0_dp/sqrt(pi) 190 do ii=0,nptsdiv2_def 191 xx=xgrid_prev(ii) 192 smdfun_prev( ii,1)=dsqrpi*exp(-xx**2) 193 smdfun_prev(-ii,1)=smdfun_prev(ii,1) 194 end do 195 196 else if(occopt==8)then 197 198 ! Constant value of the delta function over the smearing interval, for testing purposes only. 199 do ii=0,nptsdiv2_def 200 xx=xgrid_prev(ii) 201 if(xx>half+tol8)then 202 smdfun_prev( ii,1)=zero 203 else if(xx<half-tol8)then 204 smdfun_prev( ii,1)=one 205 else 206 smdfun_prev( ii,1)=half 207 end if 208 smdfun_prev(-ii,1)=smdfun_prev(ii,1) 209 end do 210 211 else 212 write(message, '(a,i0,a)' )' Occopt=',occopt,' is not allowed in getnel. ' 213 MSG_BUG(message) 214 end if 215 216 ! --------------------------------------------------------- 217 ! smear FD delta with occopt delta calculated in smdfun_prev 218 ! --------------------------------------------------------- 219 else if (dblsmr /= 0) then 220 221 nconvd2 = 6000 222 convlim = 10.0_dp 223 incconv = convlim / nconvd2 224 225 ! store smearing functions in smd1 and smd2 226 ABI_ALLOCATE(smd1,(-nconvd2:nconvd2)) 227 ABI_ALLOCATE(smd2,(-nconvd2:nconvd2)) 228 ABI_ALLOCATE(tgrid,(-nconvd2:nconvd2)) 229 230 ! FD function in smd1( ii) and second smearing delta in smd2( ii) 231 ! 232 ! smd1(:) contains delta_FD ( x ) 233 do ii=0,nconvd2 234 tgrid(ii)=ii*incconv 235 tgrid(-ii)=-tgrid(ii) 236 tt=tgrid(ii) 237 smd1( ii)=0.25_dp/(cosh(tt/2.0_dp)**2) 238 smd1(-ii)=smd1(ii) 239 end do 240 241 ! check input values of occopt and fill smd2(:) with appropriate data: 242 ! smd2(:) contains delta_resmear ( x ) 243 if(occopt == 3) then 244 write(message, '(a,a)' )& 245 & 'Occopt=3 is not allowed as a re-smearing.', & 246 & 'Use a single FD, or re-smear with a different delta type (faster cutoff). ' 247 MSG_ERROR(message) 248 else if(occopt==4 .or. occopt==5)then 249 ! Cold smearing of Marzari, two values of the "a" parameter being possible 250 ! first value gives minimization of the bump 251 if(occopt==4)aa=-.5634 252 ! second value gives monotonic occupation function 253 if(occopt==5)aa=-.8165 254 255 dsqrpi=1.0_dp/sqrt(pi) 256 do ii=0,nconvd2 257 tt=tgrid(ii) 258 gauss=dsqrpi*exp(-tt**2) 259 smd2( ii)=gauss*(1.5_dp+tt*(-aa*1.5_dp+tt*(-1.0_dp+aa*tt))) 260 smd2(-ii)=gauss*(1.5_dp+tt*( aa*1.5_dp+tt*(-1.0_dp-aa*tt))) 261 end do 262 else if(occopt==6)then 263 dsqrpi=1.0_dp/sqrt(pi) 264 do ii=0,nconvd2 265 tt=tgrid(ii) 266 smd2( ii)=dsqrpi*(1.5_dp-tt**2)*exp(-tt**2) 267 smd2(-ii)=smd2(ii) 268 end do 269 else if(occopt==7)then 270 dsqrpi=1.0_dp/sqrt(pi) 271 do ii=0,nconvd2 272 tt=tgrid(ii) 273 smd2( ii)=dsqrpi*exp(-tt**2) 274 smd2(-ii)=smd2(ii) 275 end do 276 else if(occopt==8)then 277 do ii=0,nconvd2 278 tt=tgrid(ii) 279 if(tt>half+tol8)then 280 smd2( ii)=zero 281 else if(tt<half-tol8)then 282 smd2( ii)=one 283 else 284 smd2( ii)=half 285 end if 286 smd2(-ii)=smd2(ii) 287 end do 288 else 289 write(message, '(a,i0,a)' )' Occopt= ',occopt,' is not allowed in getnel. ' 290 MSG_BUG(message) 291 end if 292 293 294 ! Use O(1/N4) algorithm from Num Rec (see below) 295 ! 296 ! The grid for the convoluted delta is taken (conservatively) 297 ! to be that for the FD delta ie 6000 pts in [-limit_occ;limit_occ] 298 ! Smearing functions are given on [-dbllim;dbllim] and the grid must 299 ! superpose the normal grid on [-limit_occ:limit_occ] 300 ! The maximal interval for integration of the convolution is 301 ! [-dbllim+limit_occ+lim(delta2);dbllim-limit_occ-lim(delta2)] = 302 ! [-dbllim+36;dbllim-36] 303 304 ! test the smdFD function for extreme values: 305 ! do jj=-nptsdiv2_def,-nptsdiv2_def 306 ! do ii=-nconvd2+4,nconvd2 307 ! call smdFD(xgrid_prev(jj) - tgrid(ii)*tratio, resFD) 308 ! write(std_out,*) 'ii jj = ', ii,jj, ' smdFD (', & 309 ! & xgrid_prev(jj) - tgrid(ii)*tratio, ') ', resFD 310 ! end do 311 ! end do 312 313 expinc = exp(half*incconv*tratio) 314 315 ! jj = position of point at which we are calculating smdfun_prev 316 do jj=-nptsdiv2_def,nptsdiv2_def 317 ! Do not care about the 8 boundary points, 318 ! where the values should be extremely small anyway 319 smdfun_prev(jj,1)=0.0_dp 320 ! only add contribution with delta_FD > 1.0d-100 321 nmaxFD = floor (( maxFDarg+xgrid_prev(jj)) / tratio / incconv ) 322 nmaxFD = min (nmaxFD, nconvd2) 323 nminFD = ceiling((-maxFDarg+xgrid_prev(jj)) / tratio / incconv ) 324 nminFD = max (nminFD, -nconvd2) 325 326 ! Calculate the Fermi-Dirac distrib at point xgrid_prev(jj)-tgrid(ii)*tratio 327 expxo2 = exp (-half*(xgrid_prev(jj) - (nminFD)*incconv*tratio)) 328 expx22 = expxo2*expxo2 329 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 330 resFD4 = tmpexpsum * tmpexpsum 331 expxo2 = expxo2*expinc 332 expx22 = expxo2*expxo2 333 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 334 resFD3 = tmpexpsum * tmpexpsum 335 expxo2 = expxo2*expinc 336 expx22 = expxo2*expxo2 337 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 338 resFD2 = tmpexpsum * tmpexpsum 339 expxo2 = expxo2*expinc 340 expx22 = expxo2*expxo2 341 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 342 resFD1 = tmpexpsum * tmpexpsum 343 344 ! core contribution to the integral with constant weight (48) 345 tmpsmdfun = 0.0_dp 346 do ii=nminFD+4,nmaxFD-4 347 expxo2 = expxo2*expinc 348 ! tmpexpsum = 1.0_dp / (expxo2 + 1.0_dp / expxo2 ) 349 expx22 = expxo2*expxo2 350 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 351 tmpsmdfun = tmpsmdfun + smd2(ii) * tmpexpsum * tmpexpsum 352 end do 353 354 ! Add on end contributions for show (both functions smd and smdFD are very small 355 smdfun_prev(jj,1)=smdfun_prev(jj,1) +48.0_dp*tmpsmdfun & 356 & + 31.0_dp*smd2(nminFD+3)*resFD1 -11.0_dp*smd2(nminFD+2)*resFD2 & 357 & + 5.0_dp*smd2(nminFD+1)*resFD3 - smd2(nminFD)*resFD4 358 359 expxo2 = expxo2*expinc 360 expx22 = expxo2*expxo2 361 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 362 resFD1 = tmpexpsum * tmpexpsum 363 expxo2 = expxo2*expinc 364 expx22 = expxo2*expxo2 365 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 366 resFD2 = tmpexpsum * tmpexpsum 367 expxo2 = expxo2*expinc 368 expx22 = expxo2*expxo2 369 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 370 resFD3 = tmpexpsum * tmpexpsum 371 expxo2 = expxo2*expinc 372 expx22 = expxo2*expxo2 373 tmpexpsum = expxo2 / (expx22 + 1.0_dp) 374 resFD4 = tmpexpsum * tmpexpsum 375 376 ! Contribution above 377 smdfun_prev(jj,1)=smdfun_prev(jj,1) & 378 & + 31.0_dp*smd2(nmaxFD-3)*resFD1 -11.0_dp*smd2(nmaxFD-2)*resFD2 & 379 & + 5.0_dp*smd2(nmaxFD-1)*resFD3 - smd2(nmaxFD)*resFD4 380 smdfun_prev(jj,1)=incconv*smdfun_prev(jj,1)/48.0_dp 381 end do 382 383 secmom = 0.0_dp 384 thdmom = 0.0_dp 385 resmom4 = xgrid_prev(-nptsdiv2_def )*xgrid_prev(-nptsdiv2_def )*smdfun_prev(-nptsdiv2_def , 1) 386 resmom3 = xgrid_prev(-nptsdiv2_def+1)*xgrid_prev(-nptsdiv2_def+1)*smdfun_prev(-nptsdiv2_def+1, 1) 387 resmom2 = xgrid_prev(-nptsdiv2_def+2)*xgrid_prev(-nptsdiv2_def+2)*smdfun_prev(-nptsdiv2_def+2, 1) 388 resmom1 = xgrid_prev(-nptsdiv2_def+3)*xgrid_prev(-nptsdiv2_def+3)*smdfun_prev(-nptsdiv2_def+3, 1) 389 resmom = xgrid_prev(-nptsdiv2_def+4)*xgrid_prev(-nptsdiv2_def+4)*smdfun_prev(-nptsdiv2_def+4, 1) 390 do ii=-nptsdiv2_def+4,nptsdiv2_def-1 391 secmom = secmom + & 392 & ( 17.0_dp*xgrid_prev(ii) *xgrid_prev(ii) *smdfun_prev(ii, 1) & 393 & +42.0_dp*xgrid_prev(ii-1)*xgrid_prev(ii-1)*smdfun_prev(ii-1,1) & 394 & -16.0_dp*xgrid_prev(ii-2)*xgrid_prev(ii-2)*smdfun_prev(ii-2,1) & 395 & + 6.0_dp*xgrid_prev(ii-3)*xgrid_prev(ii-3)*smdfun_prev(ii-3,1) & 396 & - xgrid_prev(ii-4)*xgrid_prev(ii-4)*smdfun_prev(ii-4,1) ) 397 resmom4 = resmom3 398 resmom3 = resmom2 399 resmom2 = resmom1 400 resmom1 = resmom 401 resmom = xgrid_prev(ii+1) *xgrid_prev(ii+1) *smdfun_prev(ii+1, 1) 402 end do 403 secmom=increm * secmom / 48.0_dp 404 ! thdmom=increm * thdmom / 48.0_dp 405 ! 406 ! smom1 = second moment of delta in smd1(:) 407 ! smom2 = second moment of delta in smd2(:) 408 ! 409 smom1 = 0.0_dp 410 smom2 = 0.0_dp 411 tmom1 = 0.0_dp 412 tmom2 = 0.0_dp 413 do ii=-nconvd2+4,nconvd2 414 smom1 = smom1+ & 415 & ( 17.0_dp*tgrid(ii) *tgrid(ii) *smd1(ii) & 416 & +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd1(ii-1) & 417 & -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd1(ii-2) & 418 & + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd1(ii-3) & 419 & - tgrid(ii-4)*tgrid(ii-4)*smd1(ii-4) ) 420 smom2 = smom2+ & 421 & ( 17.0_dp*tgrid(ii) *tgrid(ii) *smd2(ii ) & 422 & +42.0_dp*tgrid(ii-1)*tgrid(ii-1)*smd2(ii-1) & 423 & -16.0_dp*tgrid(ii-2)*tgrid(ii-2)*smd2(ii-2) & 424 & + 6.0_dp*tgrid(ii-3)*tgrid(ii-3)*smd2(ii-3) & 425 & - tgrid(ii-4)*tgrid(ii-4)*smd2(ii-4) ) 426 end do 427 smom1 =incconv * smom1 / 48.0_dp 428 smom2 =incconv * smom2 / 48.0_dp 429 ! tmom1 =incconv * tmom1 / 48.0_dp 430 ! tmom2 =incconv * tmom2 / 48.0_dp 431 432 encorr = smom2*tratio*tratio/secmom 433 434 ! DEBUG 435 ! write(std_out,*) ' getnel : debug, secmoms = ', secmom, smom1, smom2 436 ! write(std_out,*) ' getnel : debug, thdmoms = ', thdmom, tmom1, tmom2 437 ! write(std_out,*) ' getnel : encorr = ', encorr 438 ! ENDDEBUG 439 440 ABI_DEALLOCATE(tgrid) 441 ABI_DEALLOCATE(smd1) 442 ABI_DEALLOCATE(smd2) 443 444 end if 445 446 ! -------------------------------------------------------- 447 ! end of smearing function initialisation, dblsmr case 448 ! -------------------------------------------------------- 449 450 451 ! Now that the smeared delta function has been initialized, compute the 452 ! occupation function 453 occfun_prev(-nptsdiv2_def,1)=zero 454 entfun_prev(-nptsdiv2_def,1)=zero 455 456 ! Different algorithms are possible, corresponding to the formulas 457 ! (4.1.11), (4.1.12) and (4.1.14) in Numerical recipes (pp 107 and 108), 458 ! with respective O(1/N2), O(1/N3), O(1/N4) convergence, where N is the 459 ! number of points in the interval. 460 algo=4 461 462 if(algo==2)then 463 464 ! Extended trapezoidal rule (4.1.11), taken in a cumulative way 465 do ii=-nptsdiv2_def+1,nptsdiv2_def 466 occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*(smdfun_prev(ii,1)+smdfun_prev(ii-1,1))/2.0_dp 467 entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*& 468 & ( -xgrid_prev(ii)*smdfun_prev(ii,1) -xgrid_prev(ii-1)*smdfun_prev(ii-1,1) )/2.0_dp 469 end do 470 471 else if(algo==3)then 472 473 ! Derived from (4.1.12). Converges as O(1/N3). 474 ! Do not care about the following points, 475 ! where the values are extremely small anyway 476 occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+1,1)=0.0_dp 477 do ii=-nptsdiv2_def+2,nptsdiv2_def 478 occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*& 479 & ( 5.0_dp*smdfun_prev(ii,1) + 8.0_dp*smdfun_prev(ii-1,1) - smdfun_prev(ii-2,1) )/12.0_dp 480 entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*& 481 & ( 5.0_dp*(-xgrid_prev(ii) )*smdfun_prev(ii,1) & 482 & +8.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)& 483 & - (-xgrid_prev(ii-2))*smdfun_prev(ii-2,1) )/12.0_dp 484 end do 485 486 else if(algo==4)then 487 488 ! Derived from (4.1.14)- alternative extended Simpsons rule. Converges as O(1/N4). 489 ! Do not care about the following points, 490 ! where the values are extremely small anyway 491 occfun_prev(-nptsdiv2_def+1,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+1,1)=0.0_dp 492 occfun_prev(-nptsdiv2_def+2,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+2,1)=0.0_dp 493 occfun_prev(-nptsdiv2_def+3,1)=0.0_dp ; entfun_prev(-nptsdiv2_def+3,1)=0.0_dp 494 do ii=-nptsdiv2_def+4,nptsdiv2_def 495 occfun_prev(ii,1)=occfun_prev(ii-1,1)+increm*& 496 & ( 17.0_dp*smdfun_prev(ii,1) & 497 & +42.0_dp*smdfun_prev(ii-1,1)& 498 & -16.0_dp*smdfun_prev(ii-2,1)& 499 & + 6.0_dp*smdfun_prev(ii-3,1)& 500 & - smdfun_prev(ii-4,1) )/48.0_dp 501 entfun_prev(ii,1)=entfun_prev(ii-1,1)+increm*& 502 & ( 17.0_dp*(-xgrid_prev(ii) )*smdfun_prev(ii,1) & 503 & +42.0_dp*(-xgrid_prev(ii-1))*smdfun_prev(ii-1,1)& 504 & -16.0_dp*(-xgrid_prev(ii-2))*smdfun_prev(ii-2,1)& 505 & + 6.0_dp*(-xgrid_prev(ii-3))*smdfun_prev(ii-3,1)& 506 & - (-xgrid_prev(ii-4))*smdfun_prev(ii-4,1) )/48.0_dp 507 end do 508 509 end if ! End of choice between different algorithms for integration 510 511 ! Normalize the functions (actually not needed for occopt=3..7) 512 factor=1.0_dp/occfun_prev(nptsdiv2_def,1) 513 smdfun_prev(:,1)=smdfun_prev(:,1)*factor 514 occfun_prev(:,1)=occfun_prev(:,1)*factor 515 entfun_prev(:,1)=entfun_prev(:,1)*factor 516 517 ! Compute the cubic spline fitting of the smeared delta function 518 yp1=0.0_dp ; ypn=0.0_dp 519 workfun(:)=smdfun_prev(:,1) 520 call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, smdder) 521 smdfun_prev(:,2)=smdder(:) 522 523 ! Compute the cubic spline fitting of the occupation function 524 yp1=0.0_dp ; ypn=0.0_dp 525 workfun(:)=occfun_prev(:,1) 526 call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, occder) 527 occfun_prev(:,2)=occder(:) 528 529 ! Compute the cubic spline fitting of the entropy function 530 yp1=0.0_dp ; ypn=0.0_dp 531 workfun(:)=entfun_prev(:,1) 532 call spline(xgrid_prev, workfun, (2*nptsdiv2_def+1), yp1, ypn, entder) 533 entfun_prev(:,2)=entder(:) 534 535 ABI_DEALLOCATE(entder) 536 ABI_DEALLOCATE(occder) 537 ABI_DEALLOCATE(smdder) 538 ABI_DEALLOCATE(work) 539 ABI_DEALLOCATE(workfun) 540 541 end if 542 543 if (abs(tphysel)<tol12) then 544 tsmearinv=one/tsmear 545 else 546 tsmearinv=one/tphysel 547 end if 548 549 entfun(:,:) = entfun_prev(:,:) 550 occfun(:,:) = occfun_prev(:,:) 551 smdfun(:,:) = smdfun_prev(:,:) 552 xgrid(:) = xgrid_prev(:) 553 limit = limit_occ 554 nptsdiv2 = nptsdiv2_def 555 556 end subroutine init_occ_ent