Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief routines for handling splines
10 : !> \par History
11 : !> 2001-09-21-HAF added this doc entry and changed formatting
12 : !> \author various
13 : ! **************************************************************************************************
14 : MODULE splines_methods
15 :
16 : USE cp_log_handling, ONLY: cp_logger_get_default_unit_nr,&
17 : cp_logger_type
18 : USE kinds, ONLY: dp
19 : USE splines_types, ONLY: spline_data_p_type,&
20 : spline_data_type,&
21 : spline_factor_type
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'splines_methods'
28 :
29 : PUBLIC :: init_splinexy ! allocates x and y vectors for splines
30 : PUBLIC :: init_spline ! generate table for spline (allocates y2)
31 : PUBLIC :: potential_s ! return value of spline and 1. derivative
32 : ! without checks (fast routine for pair_potential)
33 : PUBLIC :: spline_value ! return value of spline and 1. derivative
34 : ! without check (without assumption of 1/x^2 grid)
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief allocates storage for function table to be interpolated
40 : !> both x and y are allocated
41 : !> \param spl spline_data structure to be initialized
42 : !> \param nn integer number of datapoints, that the function table will hold
43 : !> \par History
44 : !> 2001-09-21-HAF added this doc entry and changed formatting
45 : !> \author unknown
46 : ! **************************************************************************************************
47 334179 : PURE SUBROUTINE init_splinexy(spl, nn)
48 :
49 : TYPE(spline_data_type), INTENT(INOUT) :: spl
50 : INTEGER, INTENT(IN) :: nn
51 :
52 334179 : spl%n = nn
53 :
54 334179 : IF (ASSOCIATED(spl%y)) THEN
55 302539 : DEALLOCATE (spl%y)
56 : NULLIFY (spl%y)
57 : END IF
58 :
59 334179 : IF (ASSOCIATED(spl%y2)) THEN
60 302539 : DEALLOCATE (spl%y2)
61 : NULLIFY (spl%y2)
62 : END IF
63 :
64 1002537 : ALLOCATE (spl%y(1:nn))
65 1002537 : ALLOCATE (spl%y2(1:nn))
66 :
67 334179 : END SUBROUTINE init_splinexy
68 :
69 : ! **************************************************************************************************
70 : !> \brief allocates storage for y2 table
71 : !> calculates y2 table and other spline parameters
72 : !> \param spl spline_data structure to be initialized
73 : !> spl%y() must hold the function values
74 : !> spl%x() must hold the absissa values in increasing
75 : !> order OR if dx (below) is given, spl%x(1) must hold
76 : !> the starting (left-most) point.
77 : !> \param dx x(i) are assumed to be x(1)+dx*(i-1)
78 : !> (spline evaluations will also be faster)
79 : !> y1a : (OPTIONAL) if present, the 1-deriv of the left endpoint
80 : !> if not present, natural spline condition at this end
81 : !> (2-deriv == 0)
82 : !> y1b : (OPTIONAL) if present, the 1-deriv of the right endpoint
83 : !> if not present, natural spline condition at this end
84 : !> (2-deriv == 0)
85 : !> \param y1a ...
86 : !> \param y1b ...
87 : !> \par Examples
88 : !> CALL init_spline(spline,dx=0.1_dp)
89 : !> CALL init_spline(spline,y1b=0.0_dp)
90 : !> CALL init_spline(spline,0.1_dp,0.0_dp,0.0_dp)
91 : !> \par History
92 : !> 2001-09-21-HAF added this doc entry and changed formatting
93 : !> 2001-09-24-HAF changed interface and re-written
94 : !> \author unknown
95 : !> \note
96 : !> if dx is given, the x array will be used as y2 array instead of
97 : !> allocating a new array. (y2 will become x, and x will be nullified)
98 : ! **************************************************************************************************
99 334429 : PURE SUBROUTINE init_spline(spl, dx, y1a, y1b)
100 :
101 : TYPE(spline_data_type), INTENT(INOUT) :: spl
102 : REAL(KIND=dp), INTENT(IN) :: dx
103 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: y1a, y1b
104 :
105 : INTEGER :: i, n
106 : REAL(KIND=dp) :: p, s
107 334429 : REAL(KIND=dp), POINTER :: ww(:)
108 :
109 334429 : n = spl%n
110 334429 : spl%xn = spl%x1 + (n - 1)*dx
111 334429 : spl%h = dx
112 334429 : spl%invh = 1.0_dp/dx
113 334429 : spl%h26 = dx**2/6.0_dp
114 1003287 : ALLOCATE (ww(1:n))
115 334429 : IF (PRESENT(y1a)) THEN
116 0 : spl%y2(1) = -0.5_dp
117 0 : ww(1) = 3.0_dp*((spl%y(2) - spl%y(1))/dx - y1a)/dx
118 : ELSE
119 334429 : spl%y2(1) = 0.0_dp
120 334429 : ww(1) = 0.0_dp
121 : END IF
122 128972491 : DO i = 2, n - 1
123 128638062 : s = 0.5_dp
124 128638062 : p = 0.5_dp*spl%y2(i - 1) + 2.0_dp
125 128638062 : spl%y2(i) = -0.5_dp/p
126 : ww(i) = (3.0_dp*(spl%y(i + 1) - 2.0_dp*spl%y(i) + spl%y(i - 1))/(dx*dx) &
127 128972491 : - 0.5_dp*ww(i - 1))/p
128 : END DO
129 334429 : IF (PRESENT(y1b)) THEN
130 : spl%y2(n) = (3.0_dp*(y1b - (spl%y(n) - spl%y(n - 1))/dx)/dx - &
131 0 : 0.5_dp*ww(n - 1))/(0.5_dp*spl%y2(n - 1) + 1.0_dp)
132 : ELSE
133 334429 : spl%y2(n) = 0.0_dp
134 : END IF
135 129306920 : DO i = n - 1, 1, -1
136 129306920 : spl%y2(i) = spl%y2(i)*spl%y2(i + 1) + ww(i)
137 : END DO
138 334429 : DEALLOCATE (ww)
139 :
140 334429 : END SUBROUTINE init_spline
141 :
142 : ! **************************************************************************************************
143 : !> \brief calculates the potential interpolated with splines value at a given point
144 : !> and the first derivative. Checks included to avoid just segfaulting!!
145 : !> \param spl_p spline_data structure
146 : !> \param xxi absissa value
147 : !> \param y1 1. derivative at xx
148 : !> \param spl_f ...
149 : !> \param logger ...
150 : !> \return ...
151 : !> \par Output
152 : !> spline interpolated value at xx
153 : !> \par History
154 : !> 2001-09-25-HAF added this doc entry and changed formatting
155 : !> \author unknown
156 : !> \note
157 : !> the spline MUST have uniform x values and xx MUST be
158 : !> in the interpolation interval. No checks are done to ensure
159 : !> either condition.
160 : ! **************************************************************************************************
161 1326074278 : FUNCTION potential_s(spl_p, xxi, y1, spl_f, logger)
162 : TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spl_p
163 : REAL(KIND=dp), INTENT(IN) :: xxi
164 : REAL(KIND=dp), INTENT(OUT) :: y1
165 : TYPE(spline_factor_type), POINTER :: spl_f
166 : TYPE(cp_logger_type), POINTER :: logger
167 : REAL(KIND=dp) :: potential_s
168 :
169 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp
170 :
171 : INTEGER :: i, output_unit
172 : REAL(KIND=dp) :: a, b, h26, invh, x4, xx, xx0, y2hi, &
173 : y2lo, yhi, ylo, yy
174 :
175 1326074278 : xx0 = 1.0_dp/xxi
176 1326074278 : xx = spl_f%rscale(1)*xx0
177 1326074278 : x4 = xx*xx
178 1326074278 : h26 = spl_p(1)%spline_data%h26
179 1326074278 : invh = spl_p(1)%spline_data%invh
180 1326074278 : IF (xx >= spl_p(1)%spline_data%xn) THEN
181 : ! In case the value is not on the spline let's print a warning and give the value
182 : ! for the smaller point available in the spline..
183 : ! This should happen in very few cases though..
184 0 : output_unit = cp_logger_get_default_unit_nr(logger)
185 0 : yy = spl_p(1)%spline_data%xn - spl_p(1)%spline_data%h
186 : WRITE (output_unit, FMT='(/,80("*"),/,"*",1X,"Value of r in Input =",F11.6,'// &
187 0 : '" not in the spline range. Using =",F11.6,T80,"*",/,80("*"))') SQRT(1.0_dp/xx), SQRT(1.0_dp/yy)
188 0 : xx = yy
189 : END IF
190 1326074278 : i = INT((xx - spl_p(1)%spline_data%x1)*invh + 1)
191 1326074278 : a = (spl_p(1)%spline_data%x1 - xx)*invh + REAL(i, kind=dp)
192 1326074278 : b = 1.0_dp - a
193 :
194 1326074278 : ylo = spl_p(1)%spline_data%y(i)
195 1326074278 : yhi = spl_p(1)%spline_data%y(i + 1)
196 1326074278 : y2lo = spl_p(1)%spline_data%y2(i)
197 1326074278 : y2hi = spl_p(1)%spline_data%y2(i + 1)
198 1326074278 : potential_s = (a*ylo + b*yhi - ((a + 1.0_dp)*y2lo + (b + 1.0_dp)*y2hi)*a*b*h26)*spl_f%fscale(1)
199 1326074278 : y1 = invh*((yhi - ylo) + ((f13 - a*a)*y2lo - (f13 - b*b)*y2hi)*3.0_dp*h26)
200 1326074278 : y1 = 2.0_dp*y1*x4*spl_f%dscale(1)
201 :
202 1326074278 : potential_s = potential_s + spl_f%cutoff
203 1326074278 : END FUNCTION potential_s
204 :
205 : ! **************************************************************************************************
206 : !> \brief calculates the spline value at a given point
207 : !> (and possibly the first derivative) WITHOUT checks
208 : !> and without any funny scaling business, or weird
209 : !> 1/x^2 grid assumptions
210 : !>
211 : !> \param spl ...
212 : !> \param xx ...
213 : !> \param y1 ...
214 : !> \return ...
215 : !> \author HAF
216 : ! **************************************************************************************************
217 333560962 : FUNCTION spline_value(spl, xx, y1)
218 : ! Return value
219 : TYPE(spline_data_type), INTENT(IN) :: spl
220 : REAL(KIND=dp), INTENT(IN) :: xx
221 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: y1
222 : REAL(KIND=dp) :: spline_value
223 :
224 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp
225 :
226 : INTEGER :: i
227 : REAL(KIND=dp) :: a, b, h26, invh, y2hi, y2lo, yhi, ylo
228 :
229 333560962 : h26 = spl%h26
230 333560962 : invh = spl%invh
231 333560962 : i = INT((xx - spl%x1)*invh + 1)
232 :
233 333560962 : a = (spl%x1 - xx)*invh + REAL(i, kind=dp)
234 333560962 : b = 1.0_dp - a
235 333560962 : ylo = spl%y(i)
236 333560962 : yhi = spl%y(i + 1)
237 333560962 : y2lo = spl%y2(i)
238 333560962 : y2hi = spl%y2(i + 1)
239 333560962 : spline_value = a*ylo + b*yhi - ((a + 1.0_dp)*y2lo + (b + 1.0_dp)*y2hi)*a*b*h26
240 333560962 : IF (PRESENT(y1)) y1 = invh*((yhi - ylo) + &
241 0 : ((f13 - a*a)*y2lo - (f13 - b*b)*y2hi)*3.0_dp*h26)
242 333560962 : END FUNCTION spline_value
243 :
244 : END MODULE splines_methods
|
