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