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