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