1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module logrid_oct_m
22  use global_oct_m
23  use messages_oct_m
24  use profiling_oct_m
25
26  implicit none
27
28  private
29  public ::         &
30    logrid_t,       &
31    logrid_init,    &
32    logrid_end,     &
33    logrid_copy,    &
34    logrid_index,   &
35    logrid_radius,  &
36    derivative_in_log_grid
37
38  integer, parameter, public :: &
39    LOGRID_PSF  = 1, & !< log grid used in Troullier-Martins code
40    LOGRID_CPI  = 2    !< log grid used in FHI code
41
42  type logrid_t
43    ! Components are public by default
44    integer  :: flavor
45
46    FLOAT    :: a, b
47    integer  :: nrval
48
49    FLOAT, pointer :: rofi(:)  !< r value of the point i
50    FLOAT, pointer :: r2ofi(:) !< r value of the point i
51    FLOAT, pointer :: drdi(:)  !< Jacobian, i.e., the derivative of r in terms of i
52    FLOAT, pointer :: s(:)     !< sqrt of drdi
53  end type logrid_t
54
55contains
56
57  ! ---------------------------------------------------------
58  subroutine logrid_init(grid, flavor, aa, bb, nrval)
59    type(logrid_t), intent(out) :: grid
60    integer,        intent(in)  :: flavor
61    FLOAT,          intent(in)  :: aa, bb
62    integer,        intent(in)  :: nrval
63
64    FLOAT :: rpb, ea
65    integer  :: ir
66
67    PUSH_SUB(logrid_init)
68
69    ASSERT(flavor==LOGRID_PSF.or.flavor==LOGRID_CPI)
70
71    grid%flavor = flavor
72    grid%a = aa
73    grid%b = bb
74    grid%nrval = nrval
75
76    SAFE_ALLOCATE(grid%rofi(1:nrval))
77    SAFE_ALLOCATE(grid%r2ofi(1:nrval))
78    SAFE_ALLOCATE(grid%drdi(1:nrval))
79    SAFE_ALLOCATE(grid%s(1:nrval))
80
81    select case(grid%flavor)
82    case(LOGRID_PSF)
83      rpb = bb
84      ea  = exp(aa)
85      do ir = 1, nrval
86        grid%drdi(ir) = aa*rpb
87        rpb           = rpb*ea
88        grid%rofi(ir) = bb*(exp(aa*(ir-1)) - M_ONE)
89      end do
90
91    case(LOGRID_CPI)
92      grid%rofi(1) = M_ZERO
93      grid%drdi(1) = M_ZERO
94
95      rpb = log(aa)
96      grid%rofi(2) = bb
97      grid%drdi(2) = bb*rpb
98      do ir = 3, grid%nrval
99        grid%rofi(ir) = grid%rofi(ir-1)*aa
100        grid%drdi(ir) = grid%rofi(ir)*rpb
101      end do
102    end select
103
104    ! calculate sqrt(drdi)
105    do ir = 1, grid%nrval
106      grid%s(ir)     = sqrt(grid%drdi(ir))
107      grid%r2ofi(ir) = grid%rofi(ir)**2
108    end do
109
110    POP_SUB(logrid_init)
111  end subroutine logrid_init
112
113
114  ! ---------------------------------------------------------
115  subroutine logrid_end(grid)
116    type(logrid_t), intent(inout) :: grid
117
118    PUSH_SUB(logrid_end)
119
120    SAFE_DEALLOCATE_P(grid%rofi)
121    SAFE_DEALLOCATE_P(grid%r2ofi)
122    SAFE_DEALLOCATE_P(grid%drdi)
123    SAFE_DEALLOCATE_P(grid%s)
124
125    POP_SUB(logrid_end)
126  end subroutine logrid_end
127
128
129  ! ---------------------------------------------------------
130  subroutine logrid_copy(grid_in, grid_out)
131    type(logrid_t), intent(in)  :: grid_in
132    type(logrid_t), intent(out) :: grid_out
133
134    PUSH_SUB(logrid_copy)
135
136    grid_out%flavor = grid_in%flavor
137    grid_out%a      = grid_in%a
138    grid_out%b      = grid_in%b
139    grid_out%nrval  = grid_in%nrval
140
141    SAFE_ALLOCATE(grid_out%rofi (1:grid_out%nrval))
142    SAFE_ALLOCATE(grid_out%r2ofi(1:grid_out%nrval))
143    SAFE_ALLOCATE(grid_out%drdi (1:grid_out%nrval))
144    SAFE_ALLOCATE(grid_out%s    (1:grid_out%nrval))
145
146    grid_out%rofi(:)  = grid_in%rofi(:)
147    grid_out%r2ofi(:) = grid_in%r2ofi(:)
148    grid_out%drdi(:)  = grid_in%drdi(:)
149    grid_out%s(:)     = grid_in%s(:)
150
151    POP_SUB(logrid_copy)
152  end subroutine logrid_copy
153
154
155  ! ---------------------------------------------------------
156  integer function logrid_index(grid, rofi) result(ii)
157    type(logrid_t), intent(in) :: grid
158    FLOAT,          intent(in) :: rofi
159
160    integer :: ir
161
162    PUSH_SUB(logrid_index)
163
164    ii = 0
165    do ir = 1, grid%nrval-1
166
167      if(rofi >= grid%rofi(ir).and.rofi < grid%rofi(ir+1)) then
168        if(abs(rofi-grid%rofi(ir)) < abs(rofi-grid%rofi(ir+1))) then
169          ii = ir
170        else
171          ii = ir + 1
172        end if
173        exit
174      end if
175
176    end do
177
178    POP_SUB(logrid_index)
179  end function logrid_index
180
181
182  ! ---------------------------------------------------------
183  subroutine derivative_in_log_grid(grid, ff, dfdr)
184    type(logrid_t), intent(in)   :: grid
185    FLOAT,          intent(in)   :: ff(:)
186    FLOAT,          intent(out)  :: dfdr(:)
187
188    integer :: ii
189
190    PUSH_SUB(derivative_in_log_grid)
191
192    dfdr(1) = (ff(2) - ff(1))/(grid%rofi(2) - grid%rofi(1))
193    do ii = 2, grid%nrval-1
194      dfdr(ii) = (ff(ii+1) - ff(ii-1))/(grid%rofi(ii+1) - grid%rofi(ii-1))
195    end do
196    dfdr(grid%nrval) = (ff(grid%nrval) - ff(grid%nrval-1))/(grid%rofi(grid%nrval) - grid%rofi(grid%nrval-1))
197
198    POP_SUB(derivative_in_log_grid)
199  end subroutine derivative_in_log_grid
200
201  ! ----------------------------------------------------------
202  FLOAT pure function logrid_radius(grid) result(radius)
203    type(logrid_t), intent(in)   :: grid
204
205    radius = grid%rofi(grid%nrval)
206  end function logrid_radius
207
208end module logrid_oct_m
209
210!! Local Variables:
211!! mode: f90
212!! coding: utf-8
213!! End:
214