File: findlambdadu.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (129 lines) | stat: -rw-r--r-- 3,744 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

! Copyright (C) 2009 F. Bultmark, F. Cricchio and L. Nordstrom.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.

!BOP
! !ROUTINE: findlambdadu
! !INTERFACE:
subroutine findlambdadu(is,l,ufix,lambda0,lambda)
! !INPUT/OUTPUT PARAMETERS:
!   is      : species type (in,integer)
!   l       : angular momentum (in,integer)
!   ufix    : fixed U (in,integer)
!   lambda0 : starting value for screening length  (inout,real)
!   lambda  : screening length corresponding to fixed U (out,real)
! !DESCRIPTION:
!   Find the screening length corresponding to a fixed value of $U$ by using the
!   half-interval method in the first few steps and then the more efficient
!   secant method. For $U=0$ the code automatically sets the screening length to
!   ${\rm lambdamax}=50$. This value is enough to get $F^{(k)}\sim 10^{-3}$
!   corresponding to $U\sim 0$ (that perfectly mimics a bare DFT calculation).
!   However the code is stable up to ${\rm lambdamax}\sim 200$ thanks to
!   improvement of {\tt spline}.
!
! !REVISION HISTORY:
!   Created July 2009 (Francesco Cricchio)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: is
integer, intent(in) :: l
real(8), intent(in) :: ufix
real(8), intent(inout) :: lambda0
real(8), intent(out) :: lambda
! local variables
! max iterations in secant algorithm
integer, parameter :: maxit=100
integer it,nit
! if ufix<umin, assume lambda=lambdamax
real(8), parameter :: umin=1.d-4
! if lambda<lambdamin, perform unscreened calculation
real(8), parameter :: lambdamin=1.d-2
! max value of lambda
! lambdamax=50 is enough to get F^(k)~1.d-3 corresponding to U~0
real(8), parameter :: lambdamax=50.d0
! initial step for lambda
real(8), parameter :: dl0=0.5d0
real(8) f,fold,lambdaold,dl,tol
! external functions
real(8) fyukawa,fyukawa0
external fyukawa,fyukawa0
! small U limit
if (ufix.lt.umin) then
  lambda=lambdamax
  write(*,'("Info(findlambdadu): lambda set to lambdamax (",G18.10,")")') &
   lambdamax
  return
end if
! first perform a search of lambda with half-interval method and low accuracy
! initialize values and search upward from lambda0
lambda=lambda0
dl=dl0
fold=1.d0
tol=1.d-1
nit=0
do it=1,maxit
  if (lambda.lt.lambdamin) then
! unscreened Slater parameters
    f=fyukawa0(is,l,0)-ufix
  else
! screened Slater parameters
    f=fyukawa(is,l,0,lambda)-ufix
  end if
  if ((f*fold).lt.0) dl=-0.5d0*dl
  lambdaold=lambda
  lambda=lambda+dl
  fold=f
  nit=nit+1
  if (abs(f).lt.tol) goto 10
end do
10 continue
! use the found value of lambda to continue the search with secant algorithm and
! higher accuracy
tol=1.d-6
! calculate F^(0)-ufix at lambdaold value
if (lambdaold.lt.lambdamin) then
! unscreened Slater parameters
  fold=fyukawa0(is,l,0)-ufix
else
! screened Slater parameters
  fold=fyukawa(is,l,0,lambdaold)-ufix
end if
! start secant algorithm
do it=1,maxit
! calculate F^(0)-ufix
  if (lambda.lt.lambdamin) then
! unscreened Slater parameters
    f=fyukawa0(is,l,0)-ufix
  else
! screened Slater parameters
    f=fyukawa(is,l,0,lambda)-ufix
  end if
! if requested tolerance has been reached exit the loop
  if (abs(f).lt.tol) goto 20
! update lambda with secant algorithm and roll values
  dl=-f*((lambda-lambdaold)/(f-fold))
  lambdaold=lambda
  lambda=lambda+dl
  fold=f
  nit=nit+1
end do
20 continue
if (nit.ge.maxit) then
  write(*,*)
  write(*,'("Error(findlambdadu): max number of iterations to obtain lambda &
   &reached")')
  write(*,*)
  stop
else
! update initial value for lambda for the next iteration in the SC loop
! 0.5*dl0 is enough
  lambda0=lambda-0.5d0*dl0
  write(*,'("Info(findlambdadu): lambda obtained in ",I4," iterations")') nit
end if
return
end subroutine
!EOC