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!> Implements the Extended Lagrangian Born-Oppenheimer MD.
11!> Aradi et al. Extended lagrangian density functional tight-binding molecular dynamics for
12!> molecules and solids. J. Chem. Theory Comput. 11:3357-3363, 2015
13module dftbp_xlbomd_module
14  use dftbp_assert
15  use dftbp_accuracy
16  use dftbp_io
17  use dftbp_message
18  use dftbp_extlagrangian_module
19  implicit none
20  private
21
22  public :: XlbomdInp, Xlbomd, Xlbomd_init
23
24
25  !> File for reading the inverse of the Jacobian matrix if needed
26  character(*), parameter :: JacobianKernelFile = "neginvjac.dat"
27
28
29  !> Input for the Xlbomd driver.
30  type :: XlbomdInp
31
32    !> Number of generation to consider during the integration (5, 6, 7)
33    integer :: nKappa
34
35    !> Scaling factor (only for fast Xlbomd, otherwise 1.0)
36    real(dp) :: scale
37
38    !> SCC parameter override: minimal SCC iterations per timestep
39    integer :: minSCCIter
40
41    !> SCC parameter override: maximal SCC iterations per timestep
42    integer :: maxSCCIter
43
44    !> SCC parameter override: scc tolerance
45    real(dp) :: sccTol
46
47    !> Number of time steps to precede the actual start of the XL integrator
48    integer :: nPreSteps
49
50    !> Number of full SCC steps after the XL integrator has been started. Those
51    !> steps are used to fill up the integrator and to average the Jacobian.
52    integer :: nFullSCCSteps
53
54    !> Number of transient steps (during which prediction is corrected)
55    integer :: nTransientSteps
56
57    !> Whether Xlbomd should use inverse Jacobian instead of scaling
58    logical :: useInverseJacobian
59
60    !> Whether Jacobian should be read from disk.
61    logical :: readInverseJacobian
62
63  end type XlbomdInp
64
65
66  !> Contains the data for the Xlbomd driver.
67  type :: Xlbomd
68    private
69    type(ExtLagrangian) :: extLagr
70    integer :: nKappa
71    integer :: minSCCIter, maxSCCIter, minSCCIter0, maxSCCIter0
72    real(dp) :: sccTol, sccTol0
73    integer :: iStep
74    integer :: iStartXl, nPreSteps, nTransientSteps, nFullSCCSteps
75    logical :: useInverseJacobian, readInverseJacobian
76    real(dp), allocatable :: invJacobian(:,:)
77  contains
78    procedure :: setDefaultSCCParameters
79    procedure :: isActive
80    procedure :: getSCCParameters
81    procedure :: getNextCharges
82    procedure :: needsInverseJacobian
83    procedure :: setInverseJacobian
84    procedure, private :: readJacobianKernel
85  end type Xlbomd
86
87contains
88
89
90  !> Initializes the Xlbomd instance.
91  subroutine Xlbomd_init(this, input, nElems)
92
93    !> Instance.
94    type(Xlbomd), intent(out) :: this
95
96    !> Basic input parameters.
97    type(XlbomdInp), intent(in) :: input
98
99    !> Nr. of elements in the charge vector
100    integer, intent(in) :: nElems
101
102    type(ExtLagrangianInp) :: extLagrInp
103
104    @:ASSERT(input%scale >= 0.0_dp)
105    @:ASSERT(input%minSCCIter >= 1)
106    @:ASSERT(input%maxSCCIter >= 1 .and. input%maxSCCIter >= input%minSCCIter)
107    @:ASSERT(input%sccTol > 0.0_dp)
108
109    extLagrInp%nTimeSteps = input%nKappa
110    extLagrInp%nElems = nElems
111    call ExtLagrangian_init(this%extLagr, extLagrInp)
112
113    call this%extLagr%setPreconditioner(scale=input%scale)
114
115    this%nKappa = input%nKappa
116    this%iStep = 1
117    this%nTransientSteps = input%nTransientSteps
118    this%nPreSteps = input%nPreSteps
119    this%nFullSCCSteps = input%nFullSCCSteps
120    this%iStartXl = this%nPreSteps + input%nFullSCCSteps - input%nKappa
121
122    this%minSCCIter = input%minSCCIter
123    this%maxSCCIter = input%maxSCCIter
124    this%sccTol = input%sccTol
125
126    this%useInverseJacobian = input%useInverseJacobian
127    this%readInverseJacobian = input%readInverseJacobian
128    if (this%useInverseJacobian) then
129      allocate(this%invJacobian(nElems, nElems))
130      if (this%readInverseJacobian) then
131        call this%readJacobianKernel()
132      else
133        this%invJacobian(:,:) = 0.0_dp
134      end if
135    end if
136
137  end subroutine Xlbomd_init
138
139
140  !> Delivers charges for the next time step.
141  subroutine getNextCharges(this, qCurrent, qNext)
142
143    !> Instance.
144    class(Xlbomd), intent(inout) :: this
145
146    !> Charges for current time step.
147    real(dp), intent(in) :: qCurrent(:)
148
149    !> Charges for next time step.
150    real(dp), intent(out) :: qNext(:)
151
152    if (this%iStep == this%iStartXl) then
153      call this%extLagr%turnOn(this%nTransientSteps)
154    end if
155    call this%extLagr%getNextInput(qCurrent, qNext)
156    this%iStep = this%iStep + 1
157
158  end subroutine getNextCharges
159
160
161  !> Sets default SCC parameters to return, if no override is done.
162  subroutine setDefaultSCCParameters(this, minSCCIter, maxSCCIter, sccTol)
163
164    !> Instance.
165    class(Xlbomd), intent(inout) :: this
166
167    !> Minimal number of SCC iterations.
168    integer, intent(in) :: minSCCIter
169
170    !> Maximal number of SCC iterations.
171    integer, intent(in) :: maxSCCIter
172
173    !> SCC tolerance.
174    real(dp), intent(in) :: sccTol
175
176    this%minSCCIter0 = minSCCIter
177    this%maxSCCIter0 = maxSCCIter
178    this%sccTol0 = sccTol
179
180  end subroutine setDefaultSCCParameters
181
182
183  !> Signals whether XLBOMD integration is active.
184  function isActive(this)
185
186    !> Instance.
187    class(Xlbomd), intent(in) :: this
188
189    !> True if XLBOMD integration is active (no SCC convergence needed)
190    logical :: isActive
191
192    isActive = this%extLagr%needsConvergedValues()
193
194  end function isActive
195
196
197  !> Returns the SCC parameters to be used when the integrator is active.
198  subroutine getSCCParameters(this, minSCCIter, maxSCCIter, sccTol)
199
200    !> Instance.
201    class(Xlbomd), intent(in) :: this
202
203    !> Minimal number of SCC cycles.
204    integer, intent(out) :: minSCCIter
205
206    !> Maximal number of SCC cycles
207    integer, intent(out) :: maxSCCIter
208
209    !> Tolerance for SCC convergence.
210    real(dp), intent(out) :: sccTol
211
212    if (this%extLagr%needsConvergedValues()) then
213      minSCCIter = this%minSCCIter0
214      maxSCCIter = this%maxSCCIter0
215      sccTol = this%sccTol0
216    else
217      minSCCIter = this%minSCCIter
218      maxSCCIter = this%maxSCCIter
219      sccTol = this%sccTol
220    end if
221
222  end subroutine getSCCParameters
223
224
225  !> Whether integrator needs the inverse Jacobian at the current step.
226  function needsInverseJacobian(this)
227
228    !> Instance.
229    class(Xlbomd), intent(in) :: this
230
231    !> True, if a inverse Jacobian is needed. It should be passed via the setInverseJacobian()
232    !> procedure.
233    logical :: needsInverseJacobian
234
235    integer :: iStart, iEnd
236
237    needsInverseJacobian = .false.
238    if (this%useInverseJacobian .and. .not. this%readInverseJacobian) then
239      iStart = this%nPreSteps + 1
240      iEnd = this%nPreSteps + this%nFullSCCSteps
241      needsInverseJacobian = (this%iStep >= iStart .and. this%iStep <= iEnd)
242    end if
243
244  end function needsInverseJacobian
245
246
247  !> Sets the inverse Jacobian for driving the Xlbomd simulation.
248  subroutine setInverseJacobian(this, invJacobian)
249
250    !> Instance.
251    class(Xlbomd), intent(inout) :: this
252
253    !> Inverse Jacobian.
254    real(dp), intent(in) :: invJacobian(:,:)
255
256    real(dp) :: coeffOld, coeffNew, normFactor
257    integer :: nn
258
259    if (this%needsInverseJacobian()) then
260      nn = this%iStep - this%nPreSteps
261      normFactor = 1.0_dp / sum(invJacobian(:,1))
262      coeffNew = 1.0_dp / real(nn, dp)
263      coeffOld = real(nn - 1, dp) * coeffNew
264      this%invJacobian(:,:) = this%invJacobian * coeffOld &
265          & + normFactor * invJacobian * coeffNew
266      call this%extLagr%setPreconditioner(precondMtx=this%invJacobian)
267    end if
268
269  end subroutine setInverseJacobian
270
271
272  !> Read inverse Jacobian from disc
273  subroutine readJacobianKernel(this)
274
275    !> Instance.
276    class(Xlbomd), intent(inout) :: this
277
278    open(12, file=JacobianKernelFile, status="old", action="read")
279    read(12, *) this%invJacobian
280    close(12)
281    this%invJacobian = transpose(this%invJacobian)
282    write(stdout, "(A,A,A)") "Negative inverse Jacobian read from '", &
283        & JacobianKernelFile, "'"
284
285  end subroutine readJacobianKernel
286
287end module dftbp_xlbomd_module
288