Line data Source code
1 : module m_intgrf
2 : implicit none
3 : TYPE :: intgrf_out
4 : REAL :: value ! value of the integration
5 : INTEGER :: ierror ! error code
6 : END TYPE intgrf_out
7 :
8 : ! error and warning codes for intgrf function
9 : INTEGER, PARAMETER :: NO_ERROR = 0
10 : INTEGER, PARAMETER :: NEGATIVE_EXPONENT_WARNING = 1
11 : INTEGER, PARAMETER :: NEGATIVE_EXPONENT_ERROR = 2
12 : CONTAINS
13 :
14 : !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
15 : !
16 : ! Integrates function f numerically (Lagrange and Simpson integration)
17 : ! on grid(itype) and is much faster than intgr.
18 : ! (Only normal outward integration.)
19 : ! Before first use of this function it has to be initialized with
20 : ! intgrf_init.
21 :
22 148744 : FUNCTION intgrf(f, atoms, itype, gridf)
23 :
24 : use m_juDFT
25 : use m_types_setup
26 : USE m_constants
27 :
28 : IMPLICIT NONE
29 :
30 : REAL :: intgrf
31 : type(t_atoms) :: atoms
32 : INTEGER, INTENT(IN) :: itype
33 : REAL, INTENT(IN) :: gridf(atoms%jmtd, atoms%ntype)
34 : REAL, INTENT(IN) :: f(:)
35 : ! - local -
36 : TYPE(intgrf_out) :: fct_res
37 :
38 : fct_res = pure_intgrf(f, atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, &
39 148744 : atoms%ntype, itype, gridf)
40 148744 : IF (fct_res%ierror == NEGATIVE_EXPONENT_WARNING) THEN
41 : write (oUnit, *) 'intgrf: Warning!'// &
42 0 : 'Negative exponent x in extrapolation a+c*r**x'
43 148744 : ELSEIF (fct_res%ierror == NEGATIVE_EXPONENT_ERROR) THEN
44 : write (oUnit, *) &
45 0 : 'intgrf: Negative exponent x in extrapolation a+c*r**x'
46 : CALL juDFT_error( &
47 0 : 'intgrf: Negative exponent x in extrapolation a+c*r**x')
48 : END IF
49 148744 : intgrf = fct_res%value
50 :
51 148744 : END FUNCTION intgrf
52 :
53 : ! pure wrapper for intgrf with same functionality
54 : ! can be used within forall loops
55 :
56 148744 : PURE FUNCTION pure_intgrf(f, jri, jmtd, rmsh, dx, ntype, itype, gridf)
57 :
58 : IMPLICIT NONE
59 :
60 : TYPE(intgrf_out) :: pure_intgrf
61 :
62 : INTEGER, INTENT(IN) :: itype, ntype, jmtd
63 : INTEGER, INTENT(IN) :: jri(ntype)
64 : REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
65 : REAL, INTENT(IN) :: gridf(jmtd, ntype)
66 : REAL, INTENT(IN) :: f(:)
67 : ! - local -
68 : INTEGER :: n
69 : REAL :: r1, h, a, x
70 :
71 148744 : n = jri(itype)
72 148744 : r1 = rmsh(1, itype)
73 148744 : h = dx(itype)
74 :
75 148744 : pure_intgrf%ierror = NO_ERROR
76 :
77 : ! integral from 0 to r1 approximated by leading term in power series expansion
78 148744 : IF (f(1)*f(2) > 1e-10 .AND. h > 0) THEN
79 1726 : IF (f(2) == f(1)) THEN
80 0 : pure_intgrf%value = r1*f(1)
81 : ELSE
82 1726 : x = (f(3) - f(2))/(f(2) - f(1))
83 1726 : a = (f(2) - x*f(1))/(1 - x)
84 1726 : x = log(x)/h
85 1726 : IF (x < 0) THEN
86 0 : IF (x > -1) THEN
87 : pure_intgrf%ierror = NEGATIVE_EXPONENT_WARNING
88 : ELSE IF (x <= -1) THEN
89 0 : pure_intgrf%ierror = NEGATIVE_EXPONENT_ERROR
90 0 : RETURN
91 : END IF
92 : END IF
93 :
94 1726 : pure_intgrf%value = r1*(f(1) + x*a)/(x + 1)
95 :
96 : ! x = f(2) / f(1)
97 : ! x = log(x)/h
98 : ! IF(x.lt.0) THEN
99 : ! IF(x.gt.-1) write(oUnit,'(A,ES9.1)') 'intgrf: Warning!&
100 : ! & Negative exponent x in&
101 : ! & extrapolation c*r**x:',x
102 : ! IF(x.le.-1) write(oUnit,'(A,ES9.1)') 'intgrf: Negative exponent&
103 : ! & x in extrapolation&
104 : ! & c*r**x:',x
105 : ! IF(x.le.-1) call juDFT_error('intgrf: Negative exponent&
106 : ! & x in extrapolation &
107 : ! & c*r**x')
108 : ! END IF
109 : ! intgrf = (r1*f(1))/(x+1)
110 :
111 : END IF
112 : ELSE
113 : pure_intgrf%value = 0
114 : END IF
115 :
116 : ! integrate from r(1) to r(n) by multiplying with gridf
117 : pure_intgrf%value = pure_intgrf%value &
118 117574152 : + dot_product(gridf(:n, itype), f(:n))
119 :
120 148744 : END FUNCTION pure_intgrf
121 :
122 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123 :
124 : ! Initializes fast numerical integration intgrf (see below).
125 :
126 36 : SUBROUTINE intgrf_init(ntype, jmtd, jri, dx, rmsh, gridf)
127 :
128 : IMPLICIT NONE
129 :
130 : INTEGER, INTENT(IN) :: ntype, jmtd
131 : INTEGER, INTENT(IN) :: jri(ntype)
132 : REAL, INTENT(IN) :: dx(ntype), rmsh(jmtd, ntype)
133 : REAL, INTENT(OUT), ALLOCATABLE :: gridf(:, :)
134 :
135 : ! - local -
136 : INTEGER :: i, j, itype
137 : INTEGER :: n, nstep, n0 = 6
138 : INTEGER, PARAMETER :: simpson(7) = (/41, 216, 27, 272, 27, 216, 41/)
139 : REAL :: r1, h, dr
140 : REAL :: r(7)
141 : REAL, PARAMETER :: lagrange(7, 6) = reshape( &
142 : (/19087., 65112., -46461., 37504., -20211., 6312., -863., &
143 : -863., 25128., 46989., -16256., 7299., -2088., 271., &
144 : & 271., -2760., 30819., 37504., -6771., 1608., -191., &
145 : -191., 1608., -6771., 37504., 30819., -2760., 271., &
146 : & 271., -2088., 7299., -16256., 46989., 25128., -863., &
147 : -863., 6312., -20211., 37504., -46461., 65112., 19087./), & ! The last row is actually never used.
148 : (/7, 6/))
149 :
150 36 : n = jmtd
151 144 : ALLOCATE (gridf(n, ntype))
152 :
153 49044 : gridf = 0
154 :
155 96 : DO itype = 1, ntype
156 :
157 60 : n = jri(itype)
158 60 : r1 = rmsh(1, itype)
159 60 : h = dx(itype)
160 :
161 60 : nstep = (n - 1)/6
162 60 : n0 = n - 6*nstep
163 60 : dr = exp(h)
164 :
165 : ! Calculate Lagrange-integration coefficients from r(1) to r(n0)
166 60 : r(1) = r1
167 60 : IF (n0 > 1) THEN
168 252 : DO i = 2, 7
169 252 : r(i) = r(i - 1)*dr
170 : END DO
171 288 : DO i = 1, 7
172 792 : gridf(i, itype) = h/60480*r(i)*sum(lagrange(i, 1:n0 - 1))
173 : END DO
174 36 : r(1) = r(n0)
175 : END IF
176 :
177 : ! Calculate Simpson-integration coefficients from r(n0) to r(n)
178 8064 : DO i = 1, nstep
179 55776 : DO j = 2, 7
180 55776 : r(j) = r(j - 1)*dr
181 : END DO
182 63744 : DO j = n0, n0 + 6
183 : gridf(j, itype) = gridf(j, itype) + h/140*r(j - n0 + 1)* &
184 63744 : simpson(j - n0 + 1)
185 : END DO
186 7968 : n0 = n0 + 6
187 8028 : r(1) = r(7)
188 : END DO
189 :
190 : END DO
191 :
192 36 : END SUBROUTINE intgrf_init
193 0 : end module m_intgrf
|
