1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8module m_fire 9 10use precision, only: dp 11use parallel, only: Node 12 13implicit none 14 15 real(dp), parameter, public :: FIRE_DEF_dt_inc = 1.1_dp 16 real(dp), parameter, public :: FIRE_DEF_dt_dec = 0.5_dp 17 real(dp), parameter, public :: FIRE_DEF_alphamax = 0.1_dp 18 real(dp), parameter, public :: FIRE_DEF_alpha_factor = 0.99_dp 19 real(dp), parameter, public :: FIRE_DEF_dtmax = 10.0_dp 20 integer, parameter, public :: FIRE_DEF_nmin = 5 21 22type, public :: fire_t 23 integer :: n 24 integer :: step = 0 25 integer :: npos = 0 26 real(dp) :: alpha 27 real(dp) :: dt 28 integer :: nmin = FIRE_DEF_nmin 29 real(dp) :: dt_inc = FIRE_DEF_dt_inc 30 real(dp) :: dt_dec = FIRE_DEF_dt_dec 31 real(dp) :: alphamax = FIRE_DEF_alphamax 32 real(dp) :: alpha_factor = FIRE_DEF_alpha_factor 33 real(dp) :: dtmax = FIRE_DEF_dtmax 34 real(dp), pointer :: v(:) => null() ! velocities 35 logical :: debug = .false. 36end type 37 38private 39 40public :: fire_setup, fire_step 41 42CONTAINS 43 44 subroutine fire_setup(b,n,dt,nmin,dt_inc,dt_dec, & 45 alphamax,alpha_factor,debug) 46 47 type(fire_t), intent(inout) :: b 48 integer, intent(in) :: n ! size of problem 49 real(dp), intent(in) :: dt 50 logical, intent(in), optional :: debug 51 integer, intent(in), optional :: nmin 52 real(dp), intent(in), optional :: dt_inc 53 real(dp), intent(in), optional :: dt_dec 54 real(dp), intent(in), optional :: alphamax 55 real(dp), intent(in), optional :: alpha_factor 56 57 if (present(debug)) then 58 b%debug = debug 59 endif 60 if (present(nmin)) then 61 b%nmin = nmin 62 endif 63 if (present(dt_inc)) then 64 b%dt_inc = dt_inc 65 endif 66 if (present(dt_dec)) then 67 b%dt_dec = dt_dec 68 endif 69 if (present(alphamax)) then 70 b%alphamax = alphamax 71 endif 72 if (present(alpha_factor)) then 73 b%alpha_factor = alpha_factor 74 endif 75 76 b%alpha = b%alphamax 77 b%n = n 78 b%dt = dt 79 b%npos = 0 80 b%step = 0 81 82 allocate(b%v(n)) 83 b%v = 0.0_dp 84 85 end subroutine fire_setup 86 87!---------------------------------------------------------- 88 89 subroutine fire_step(b,f,x,maxstep) 90 type(fire_t), intent(inout) :: b 91 real(dp), intent(in) :: f(:) 92 real(dp), intent(inout) :: x(:) 93 real(dp), intent(in) :: maxstep(:) 94 95 integer :: i 96 logical :: debug 97 98! For a general problem, we should choose dt and the masses 99! so that the first step does not travel too far compared to some 100! "typical length". Say, delta_x = 0.1 bohr 101 102 real(dp) dt, norm_v, norm_f, power 103 104 real(dp), allocatable :: dx(:) 105 integer, allocatable :: count(:) 106 107 allocate(dx(size(x))) 108 allocate(count(size(x))) 109 count = 0 110 111 debug = (b%debug .and. Node .eq. 0) 112 113 dt = b%dt 114 b%step = b%step + 1 115 116 Power = sum(f*b%v) 117 If (Power < 0) then 118 if (debug) print *, "Power < 0" 119 b%v = 0.0_dp 120 b%dt = b%dt * b%dt_dec 121 b%alpha = b%alphamax 122 b%npos = 0 123 else 124 if (debug) print *, "Power > 0" 125 norm_v = sqrt(sum(b%v*b%v)) 126 norm_f = sqrt(sum(f*f)) 127 b%v = (1-b%alpha)*b%v + b%alpha * norm_v * f / norm_f 128 if (b%npos > b%nmin) then 129 b%dt = min(b%dt * b%dt_inc, b%dtmax) 130 b%alpha = b%alpha * b%alpha_factor 131 endif 132 b%npos = b%npos + 1 133 endif ! Power 134 135 ! Euler with midpoint correction 136 dx = b%v*dt/2 ! half of the mid-point rule (v(0)) 137 b%v = b%v + f*dt ! v at dt 138 dx = dx + b%v*dt/2 ! the other half (v(t)) 139 140 where (abs(dx) <= maxstep) 141 x = x + dx 142 elsewhere 143 x = x + sign(maxstep,dx) 144 count = count + 1 145 end where 146 147 if (debug) then 148 write(6,"(a)") 'FIRE: step, dt, npos, alpha, magv, magf, long' 149 write(6,"(a, i6, f10.4, i4, f10.4, 2f12.6, i5)") & 150 "FIREdat:", b%step, b%dt, b%npos, b%alpha, & 151 sqrt(sum(b%v*b%v)), sqrt(sum(f*f)), sum(count) 152 if (b%n <=100) then 153 write(6,"(a)") "FIRE: nvar, x, v, f" 154 do i = 1, b%n 155 write(6,*) i, x(i), b%v(i), f(i) 156 enddo 157 endif 158 endif 159 160 deallocate(dx,count) 161 162 end subroutine fire_step 163 164 end module m_fire 165 166