(git:0de0cc2)
splines_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
15 
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  PURE SUBROUTINE init_splinexy(spl, nn)
48 
49  TYPE(spline_data_type), INTENT(INOUT) :: spl
50  INTEGER, INTENT(IN) :: nn
51 
52  spl%n = nn
53 
54  IF (ASSOCIATED(spl%y)) THEN
55  DEALLOCATE (spl%y)
56  NULLIFY (spl%y)
57  END IF
58 
59  IF (ASSOCIATED(spl%y2)) THEN
60  DEALLOCATE (spl%y2)
61  NULLIFY (spl%y2)
62  END IF
63 
64  ALLOCATE (spl%y(1:nn))
65  ALLOCATE (spl%y2(1:nn))
66 
67  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  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  REAL(kind=dp), POINTER :: ww(:)
108 
109  n = spl%n
110  spl%xn = spl%x1 + (n - 1)*dx
111  spl%h = dx
112  spl%invh = 1.0_dp/dx
113  spl%h26 = dx**2/6.0_dp
114  ALLOCATE (ww(1:n))
115  IF (PRESENT(y1a)) THEN
116  spl%y2(1) = -0.5_dp
117  ww(1) = 3.0_dp*((spl%y(2) - spl%y(1))/dx - y1a)/dx
118  ELSE
119  spl%y2(1) = 0.0_dp
120  ww(1) = 0.0_dp
121  END IF
122  DO i = 2, n - 1
123  s = 0.5_dp
124  p = 0.5_dp*spl%y2(i - 1) + 2.0_dp
125  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  - 0.5_dp*ww(i - 1))/p
128  END DO
129  IF (PRESENT(y1b)) THEN
130  spl%y2(n) = (3.0_dp*(y1b - (spl%y(n) - spl%y(n - 1))/dx)/dx - &
131  0.5_dp*ww(n - 1))/(0.5_dp*spl%y2(n - 1) + 1.0_dp)
132  ELSE
133  spl%y2(n) = 0.0_dp
134  END IF
135  DO i = n - 1, 1, -1
136  spl%y2(i) = spl%y2(i)*spl%y2(i + 1) + ww(i)
137  END DO
138  DEALLOCATE (ww)
139 
140  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  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  xx0 = 1.0_dp/xxi
176  xx = spl_f%rscale(1)*xx0
177  x4 = xx*xx
178  h26 = spl_p(1)%spline_data%h26
179  invh = spl_p(1)%spline_data%invh
180  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  output_unit = cp_logger_get_default_unit_nr(logger)
185  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  '" not in the spline range. Using =",F11.6,T80,"*",/,80("*"))') sqrt(1.0_dp/xx), sqrt(1.0_dp/yy)
188  xx = yy
189  END IF
190  i = int((xx - spl_p(1)%spline_data%x1)*invh + 1)
191  a = (spl_p(1)%spline_data%x1 - xx)*invh + real(i, kind=dp)
192  b = 1.0_dp - a
193 
194  ylo = spl_p(1)%spline_data%y(i)
195  yhi = spl_p(1)%spline_data%y(i + 1)
196  y2lo = spl_p(1)%spline_data%y2(i)
197  y2hi = spl_p(1)%spline_data%y2(i + 1)
198  potential_s = (a*ylo + b*yhi - ((a + 1.0_dp)*y2lo + (b + 1.0_dp)*y2hi)*a*b*h26)*spl_f%fscale(1)
199  y1 = invh*((yhi - ylo) + ((f13 - a*a)*y2lo - (f13 - b*b)*y2hi)*3.0_dp*h26)
200  y1 = 2.0_dp*y1*x4*spl_f%dscale(1)
201 
202  potential_s = potential_s + spl_f%cutoff
203  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  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  h26 = spl%h26
230  invh = spl%invh
231  i = int((xx - spl%x1)*invh + 1)
232 
233  a = (spl%x1 - xx)*invh + real(i, kind=dp)
234  b = 1.0_dp - a
235  ylo = spl%y(i)
236  yhi = spl%y(i + 1)
237  y2lo = spl%y2(i)
238  y2hi = spl%y2(i + 1)
239  spline_value = a*ylo + b*yhi - ((a + 1.0_dp)*y2lo + (b + 1.0_dp)*y2hi)*a*b*h26
240  IF (PRESENT(y1)) y1 = invh*((yhi - ylo) + &
241  ((f13 - a*a)*y2lo - (f13 - b*b)*y2hi)*3.0_dp*h26)
242  END FUNCTION spline_value
243 
244 END MODULE splines_methods
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
routines for handling splines
pure subroutine, public init_splinexy(spl, nn)
allocates storage for function table to be interpolated both x and y are allocated
pure subroutine, public init_spline(spl, dx, y1a, y1b)
allocates storage for y2 table calculates y2 table and other spline parameters
real(kind=dp) function, public potential_s(spl_p, xxi, y1, spl_f, logger)
calculates the potential interpolated with splines value at a given point and the first derivative....
real(kind=dp) function, public spline_value(spl, xx, y1)
calculates the spline value at a given point (and possibly the first derivative) WITHOUT checks and w...
routines for handling splines_types
Definition: splines_types.F:14