1!--------------------------------------------------------------------------------------------------!
2!  DFTB+: general package for performing fast atomistic simulations                                !
3!  Copyright (C) 2006 - 2019  DFTB+ developers group                                               !
4!                                                                                                  !
5!  See the LICENSE file for terms of usage and distribution.                                       !
6!--------------------------------------------------------------------------------------------------!
7
8#:include 'common.fypp'
9
10!> Extended Lagrangian dynamics
11module dftbp_extlagrangian_module
12  use dftbp_assert
13  use dftbp_accuracy, only : dp
14  use dftbp_message
15  implicit none
16  private
17
18  public ::ExtLagrangian, ExtLagrangian_init
19  public :: ExtLagrangianInp
20
21
22  !> Input for an extended Lagrangian integrator.
23  !!
24  type :: ExtLagrangianInp
25
26
27    !> Nr. of timesteps to consider for the integration.
28    integer :: nTimeSteps
29
30
31    !> Nr. of elements in the vectors which should be propagated.
32    integer :: nElems
33
34  end type ExtLagrangianInp
35
36
37  !> Represents an extended Lagrangian integrator.
38  type :: ExtLagrangian
39    private
40    integer :: phase
41    integer :: nElems
42    integer :: nTimeSteps
43    integer :: iStep, nSteps, nTransientSteps
44    integer :: ind
45    real(dp) :: scale
46    real(dp), allocatable :: precondMtx(:,:)
47    real(dp), allocatable :: auxVectors(:,:)
48    real(dp), allocatable :: auxCoeffs(:)
49    real(dp) :: alpha, kappa
50  contains
51
52    !> turn on the integrator
53    procedure :: turnOn
54
55    !> input for next time steo
56    procedure :: getNextInput
57
58    !> are converged values required for next step
59    procedure :: needsConvergedValues
60
61    !> Preconditioner for integrator
62    procedure :: setPreconditioner
63
64    !> Internal helper function
65    procedure, private :: updatePhaseAndSteps
66  end type ExtLagrangian
67
68
69  !> Internal type for enumerating different states of the integrator.
70  type :: ExtLagrangianPhases
71    integer :: off
72    integer :: fillingUp
73    integer :: interpolating
74    integer :: on
75  end type ExtLagrangianPhases
76
77
78  !> Internal states of the integrator
79  type(ExtLagrangianPhases), parameter :: phases =&
80      & ExtLagrangianPhases(1, 2, 3, 4)
81
82
83  !> Various integrator parameters for integration with 5 time steps
84  real(dp), parameter :: auxCoeffs5(0:5) = &
85      & [-6.0_dp, 14.0_dp, -8.0_dp, -3.0_dp, 4.0_dp, -1.0_dp]
86  real(dp), parameter :: alpha5 = 18e-3_dp
87  real(dp), parameter :: kappa5 = 1.82_dp
88
89
90  !> Various integrator parameters for integration with 6 time steps
91  real(dp), parameter :: auxCoeffs6(0:6) = &
92      & [-14.0_dp, 36.0_dp, -27.0_dp, -2.0_dp, 12.0_dp, -6.0_dp, 1.0_dp]
93  real(dp), parameter :: alpha6 = 5.5e-3_dp
94  real(dp), parameter :: kappa6 = 1.84_dp
95
96
97  !> Various integrator parameters for integration with 7 time steps
98  real(dp), parameter :: auxCoeffs7(0:7) = &
99      & [-36.0_dp, 99.0_dp, -88.0_dp, 11.0_dp, 32.0_dp, -25.0_dp, 8.0_dp,&
100      & -1.0_dp]
101  real(dp), parameter :: alpha7 = 1.6e-3
102  real(dp), parameter :: kappa7 = 1.86_dp
103
104contains
105
106
107  !> Initializes an extended Lagrangian integrator.
108  subroutine ExtLagrangian_init(this, input)
109
110
111    !> Initialized instance at exit.
112    class(ExtLagrangian), intent(inout) :: this
113
114
115    !> Input container.
116    class(ExtLagrangianInp), intent(in) :: input
117
118    this%nElems = input%nElems
119    this%scale = 1.0_dp
120    this%phase = phases%off
121    this%iStep = 1
122    this%nSteps = 0
123    this%nTransientSteps = 0
124
125    this%nTimeSteps = input%nTimeSteps
126    allocate(this%auxCoeffs(0:this%nTimeSteps))
127    allocate(this%auxVectors(this%nElems, this%nTimeSteps + 1))
128    select case (this%nTimeSteps)
129    case (5)
130      this%auxCoeffs(:) = auxCoeffs5
131      this%alpha = alpha5
132      this%kappa = kappa5
133    case (6)
134      this%auxCoeffs(:) = auxCoeffs6
135      this%alpha = alpha6
136      this%kappa = kappa6
137    case(7)
138      this%auxCoeffs(:) = auxCoeffs7
139      this%alpha = alpha7
140      this%kappa = kappa7
141    case default
142      call error("Invalid number of generations for ExtLagrangian")
143    end select
144
145  end subroutine ExtLagrangian_init
146
147
148  !> Turns on the integrator.
149  !!
150  !! The integrator will start of filling up its database with the subsequent
151  !! vectors passed to it. Whether the database filling is complete, can be
152  !! queried by the `needsConvergedValues()` method. When it returns `.false.`
153  !! the integrator is ready.
154  !!
155  subroutine turnOn(this, nTransientSteps)
156
157
158    !> Instance variable.
159    class(ExtLagrangian), intent(inout) :: this
160
161
162    !> Nr. of transient steps to do *additional* to the ones needed to fill up
163    !! the integrator. During those additional steps, the integrator still needs
164    !! converged quantities and will use them to correct its internal database
165    !! with the predicted input charges to enable a smoother change between the
166    !! fully converged calculations and XL predicted ones.
167    integer, intent(in), optional :: nTransientSteps
168
169    integer :: nTransientSteps0
170
171    if (present(nTransientSteps)) then
172      nTransientSteps0 = nTransientSteps
173    else
174      nTransientSteps0 = 0
175    end if
176    this%phase = phases%fillingUp
177    this%iStep = 1
178    this%nSteps = this%nTimeSteps + 1
179    this%nTransientSteps = nTransientSteps0
180
181  end subroutine turnOn
182
183
184  !> Reads the last output and provides the input for the next timestep.
185  !!
186  subroutine getNextInput(this, outLast, inNext)
187
188
189    !> Instance.
190    class(ExtLagrangian), intent(inout) :: this
191
192
193    !> Output quantity of the last iteration.
194    real(dp), intent(in) :: outLast(:)
195
196
197    !> Input quantity for the next iteration
198    real(dp), intent(out) :: inNext(:)
199
200    real(dp), allocatable :: diff(:)
201    integer :: ind, ind0, ind1
202    integer :: ii
203    real(dp) :: cc
204
205    select case (this%phase)
206
207    case (phases%off)
208
209      inNext(:) = outLast
210
211    case (phases%fillingUp)
212
213      this%ind = modIndex(this%ind + 1, this%nTimeSteps)
214      this%auxVectors(:, this%ind) = outLast
215      inNext(:) = outLast
216
217    case (phases%interpolating, phases%on)
218
219      ind0 = this%ind
220
221      ! Override last prediction in our database with an interpolation between
222      ! it and the exact result, which just have been passed.
223      if (this%phase == phases%interpolating) then
224        cc = real(this%iStep, dp) / real(this%nSteps + 1, dp)
225        this%auxVectors(:,ind0) = cc * this%auxVectors(:,ind0)&
226            & + (1.0_dp - cc) * outLast
227      end if
228
229      ! Build new auxiliary vector by integration
230      allocate(diff(size(inNext)))
231      ind1 = modIndex(ind0 - 1, this%nTimeSteps)
232      inNext(:) = 2.0_dp * this%auxVectors(:,ind0) - this%auxVectors(:,ind1)
233      diff(:) = outLast - this%auxVectors(:,ind0)
234      if (allocated(this%precondMtx)) then
235        diff(:) = matmul(this%precondMtx, diff)
236      end if
237      inNext(:) = inNext + this%kappa * this%scale * diff
238
239      ! Add dissipation
240      do ii = 0, this%nTimeSteps
241        ind = modIndex(ind0 - ii, this%nTimeSteps)
242        inNext(:) = inNext&
243            & + this%alpha * this%auxCoeffs(ii) * this%auxVectors(:,ind)
244      end do
245
246      ! Store predicted vector in the database
247      this%ind = modIndex(this%ind + 1, this%nTimeSteps)
248      this%auxVectors(:,this%ind) = inNext
249    end select
250
251    call this%updatePhaseAndSteps()
252
253  end subroutine getNextInput
254
255  !> helper function
256  pure function modIndex(ind, nTimeSteps)
257    integer, intent(in) :: ind
258    integer, intent(in) :: nTimeSteps
259    integer :: modIndex
260
261    modIndex = modulo(ind - 1, nTimeSteps + 1) + 1
262
263  end function modIndex
264
265  !> Whether next output quantity passed to the integrator should still contain
266  !! fully converged values.
267  !!
268  function needsConvergedValues(this) result(needsConverged)
269
270
271    !> Instance.
272    class(ExtLagrangian), intent(in) :: this
273
274
275    !> Whether converged values are needed.
276    logical :: needsConverged
277
278    needsConverged = any(this%phase ==&
279        & [phases%off, phases%fillingUp, phases%interpolating])
280
281  end function needsConvergedValues
282
283
284  !> Sets a preconditioner for the integrator.
285  !!
286  subroutine setPreconditioner(this, scale, precondMtx)
287
288
289    !> Instance variable.
290    class(ExtLagrangian), intent(inout) :: this
291
292
293    !> Scaling factor for the difference vector (e.g. scaling factor for
294    !! SCF-free XLBOMD). Default: 1.0.
295    real(dp), intent(in), optional :: scale
296
297
298    !> Preconditioning matrix for the difference vector (e.g. inverse Jacobian)
299    !! Default: identity matrix.
300    real(dp), intent(in), optional :: precondMtx(:,:)
301
302  #:call ASSERT_CODE
303    if (present(precondMtx)) then
304      @:ASSERT(all(shape(precondMtx) == [this%nElems, this%nElems]))
305    end if
306  #:endcall ASSERT_CODE
307
308    if (present(scale)) then
309      this%scale = scale
310    end if
311    if (present(precondMtx)) then
312      if (allocated(this%precondMtx)) then
313        deallocate(this%precondMtx)
314        allocate(this%precondMtx(this%nElems, this%nElems))
315      end if
316      this%precondMtx(:,:) = precondMtx
317    end if
318
319  end subroutine setPreconditioner
320
321  ! Private methods
322
323
324  !> helper function
325  subroutine updatePhaseAndSteps(this)
326    class(ExtLagrangian), intent(inout) :: this
327
328    if (any(this%phase == [phases%off, phases%on])) then
329      return
330    end if
331
332    this%iStep = this%iStep + 1
333    if (this%phase == phases%fillingUp .and. this%iStep > this%nSteps) then
334      this%phase = phases%interpolating
335      this%iStep = 1
336      this%nSteps = this%nTransientSteps
337    end if
338    if (this%phase == phases%interpolating .and. this%iStep > this%nSteps) then
339      this%phase = phases%on
340    end if
341
342  end subroutine updatePhaseAndSteps
343
344end module dftbp_extlagrangian_module
345