1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for performing an outer scf loop 8!> \par History 9!> Created [2006.03] 10!> \author Joost VandeVondele 11! ************************************************************************************************** 12MODULE qs_outer_scf 13 USE cp_control_types, ONLY: ddapc_restraint_type,& 14 dft_control_type,& 15 s2_restraint_type 16 USE cp_log_handling, ONLY: cp_to_string 17 USE input_constants, ONLY: & 18 broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, & 19 broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, & 20 cdft2ot, do_ddapc_constraint, do_s2_constraint, ot2cdft, outer_scf_basis_center_opt, & 21 outer_scf_cdft_constraint, outer_scf_ddapc_constraint, outer_scf_none, & 22 outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, outer_scf_optimizer_diis, & 23 outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, outer_scf_optimizer_none, & 24 outer_scf_optimizer_sd, outer_scf_optimizer_secant, outer_scf_s2_constraint 25 USE kinds, ONLY: dp 26 USE mathlib, ONLY: diamat_all 27 USE qs_basis_gradient, ONLY: qs_basis_center_gradient,& 28 qs_update_basis_center_pos,& 29 return_basis_center_gradient_norm 30 USE qs_cdft_opt_types, ONLY: cdft_opt_type_copy,& 31 cdft_opt_type_release 32 USE qs_cdft_types, ONLY: cdft_control_type 33 USE qs_energy_types, ONLY: qs_energy_type 34 USE qs_environment_types, ONLY: get_qs_env,& 35 qs_environment_type,& 36 set_qs_env 37 USE qs_scf_types, ONLY: qs_scf_env_type 38 USE scf_control_types, ONLY: scf_control_type 39#include "./base/base_uses.f90" 40 41 IMPLICIT NONE 42 43 PRIVATE 44 45! *** Global parameters *** 46 47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_outer_scf' 48 49! *** Public subroutines *** 50 51 PUBLIC :: outer_loop_gradient, outer_loop_optimize, outer_loop_update_qs_env, & 52 outer_loop_variables_count, outer_loop_extrapolate, & 53 outer_loop_switch, outer_loop_purge_history 54 55CONTAINS 56 57! ************************************************************************************************** 58!> \brief returns the number of variables that is employed in the outer loop. with a CDFT constraint 59!> this value is returned by the cdft_control type 60!> \param scf_control the outer loop control type 61!> \param cdft_control the cdft loop control type 62!> \return the number of variables 63!> \par History 64!> 03.2006 created [Joost VandeVondele] 65! ************************************************************************************************** 66 FUNCTION outer_loop_variables_count(scf_control, cdft_control) RESULT(res) 67 TYPE(scf_control_type), POINTER :: scf_control 68 TYPE(cdft_control_type), INTENT(IN), OPTIONAL, & 69 POINTER :: cdft_control 70 INTEGER :: res 71 72 SELECT CASE (scf_control%outer_scf%type) 73 CASE (outer_scf_ddapc_constraint) 74 res = 1 75 CASE (outer_scf_s2_constraint) 76 res = 1 77 CASE (outer_scf_cdft_constraint) 78 IF (PRESENT(cdft_control)) THEN 79 res = SIZE(cdft_control%target) 80 ELSE 81 res = 1 82 END IF 83 CASE (outer_scf_basis_center_opt) 84 res = 1 85 CASE (outer_scf_none) ! just needed to communicate the gradient criterion 86 res = 1 87 CASE DEFAULT 88 res = 0 89 END SELECT 90 91 END FUNCTION outer_loop_variables_count 92 93! ************************************************************************************************** 94!> \brief computes the gradient wrt to the outer loop variables 95!> \param qs_env ... 96!> \param scf_env ... 97!> \par History 98!> 03.2006 created [Joost VandeVondele] 99! ************************************************************************************************** 100 SUBROUTINE outer_loop_gradient(qs_env, scf_env) 101 TYPE(qs_environment_type), POINTER :: qs_env 102 TYPE(qs_scf_env_type), POINTER :: scf_env 103 104 CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_gradient', & 105 routineP = moduleN//':'//routineN 106 107 INTEGER :: handle, ihistory, ivar, n 108 LOGICAL :: is_constraint 109 TYPE(cdft_control_type), POINTER :: cdft_control 110 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control 111 TYPE(dft_control_type), POINTER :: dft_control 112 TYPE(qs_energy_type), POINTER :: energy 113 TYPE(s2_restraint_type), POINTER :: s2_restraint_control 114 TYPE(scf_control_type), POINTER :: scf_control 115 116 CALL timeset(routineN, handle) 117 118 CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, & 119 dft_control=dft_control, energy=energy) 120 CPASSERT(scf_control%outer_scf%have_scf) 121 122 ihistory = scf_env%outer_scf%iter_count 123 CPASSERT(ihistory <= SIZE(scf_env%outer_scf%energy, 1)) 124 125 scf_env%outer_scf%energy(ihistory) = energy%total 126 127 SELECT CASE (scf_control%outer_scf%type) 128 CASE (outer_scf_none) 129 ! just pass the inner loop scf criterion to the outer loop one 130 scf_env%outer_scf%variables(1, ihistory) = scf_env%iter_delta 131 scf_env%outer_scf%gradient(1, ihistory) = scf_env%iter_delta 132 CASE (outer_scf_ddapc_constraint) 133 CPASSERT(dft_control%qs_control%ddapc_restraint) 134 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control) 135 NULLIFY (ddapc_restraint_control) 136 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)%ddapc_restraint_control 137 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint) 138 IF (is_constraint) EXIT 139 END DO 140 CPASSERT(is_constraint) 141 142 scf_env%outer_scf%variables(:, ihistory) = ddapc_restraint_control%strength 143 scf_env%outer_scf%gradient(:, ihistory) = ddapc_restraint_control%ddapc_order_p - & 144 ddapc_restraint_control%target 145 CASE (outer_scf_s2_constraint) 146 CPASSERT(dft_control%qs_control%s2_restraint) 147 s2_restraint_control => dft_control%qs_control%s2_restraint_control 148 is_constraint = (s2_restraint_control%functional_form == do_s2_constraint) 149 CPASSERT(is_constraint) 150 151 scf_env%outer_scf%variables(:, ihistory) = s2_restraint_control%strength 152 scf_env%outer_scf%gradient(:, ihistory) = s2_restraint_control%s2_order_p - & 153 s2_restraint_control%target 154 CASE (outer_scf_cdft_constraint) 155 CPASSERT(dft_control%qs_control%cdft) 156 cdft_control => dft_control%qs_control%cdft_control 157 DO ivar = 1, SIZE(scf_env%outer_scf%gradient, 1) 158 scf_env%outer_scf%variables(ivar, ihistory) = cdft_control%strength(ivar) 159 scf_env%outer_scf%gradient(ivar, ihistory) = cdft_control%value(ivar) - & 160 cdft_control%target(ivar) 161 END DO 162 CASE (outer_scf_basis_center_opt) 163 CALL qs_basis_center_gradient(qs_env) 164 scf_env%outer_scf%gradient(:, ihistory) = return_basis_center_gradient_norm(qs_env) 165 166 CASE DEFAULT 167 CPABORT("") 168 169 END SELECT 170 171 CALL timestop(handle) 172 173 END SUBROUTINE outer_loop_gradient 174 175! ************************************************************************************************** 176!> \brief optimizes the parameters of the outer_scf 177!> \param scf_env the scf_env where to optimize the parameters 178!> \param scf_control control parameters for the optimization 179!> \par History 180!> 03.2006 created [Joost VandeVondele] 181!> 01.2017 added Broyden and Newton optimizers [Nico Holmberg] 182!> \note 183!> ought to be general, and independent of the actual kind of variables 184! ************************************************************************************************** 185 SUBROUTINE outer_loop_optimize(scf_env, scf_control) 186 TYPE(qs_scf_env_type), POINTER :: scf_env 187 TYPE(scf_control_type), POINTER :: scf_control 188 189 CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_optimize', & 190 routineP = moduleN//':'//routineN 191 192 INTEGER :: handle, i, ibuf, ihigh, ihistory, ilow, & 193 j, jbuf, nb, nvar, optimizer_type 194 REAL(KIND=dp) :: interval, scale, tmp 195 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev 196 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, f, x 197 REAL(KIND=dp), DIMENSION(:, :), POINTER :: inv_jacobian 198 199 CALL timeset(routineN, handle) 200 201 ihistory = scf_env%outer_scf%iter_count 202 optimizer_type = scf_control%outer_scf%optimizer 203 NULLIFY (inv_jacobian) 204 205 IF (scf_control%outer_scf%type == outer_scf_basis_center_opt) THEN 206 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) 207 ELSE 208 DO WHILE (.TRUE.) ! if we need a different run type we'll restart here 209 210 SELECT CASE (optimizer_type) 211 CASE (outer_scf_optimizer_bisect) ! bisection on the gradient, needs to be 1D 212 CPASSERT(SIZE(scf_env%outer_scf%gradient(:, 1)) == 1) 213 ! find the pair of points that bracket a zero of the gradient, with the smallest interval possible 214 ilow = -1 215 ihigh = -1 216 interval = HUGE(interval) 217 DO i = 1, ihistory 218 DO j = i + 1, ihistory 219 ! distrust often used points 220 IF (scf_env%outer_scf%count(i) .GT. scf_control%outer_scf%bisect_trust_count) CYCLE 221 IF (scf_env%outer_scf%count(j) .GT. scf_control%outer_scf%bisect_trust_count) CYCLE 222 223 ! if they bracket a zero use them 224 IF (scf_env%outer_scf%gradient(1, i)* & 225 scf_env%outer_scf%gradient(1, j) < 0.0_dp) THEN 226 tmp = ABS(scf_env%outer_scf%variables(1, i) - scf_env%outer_scf%variables(1, j)) 227 IF (tmp < interval) THEN 228 ilow = i 229 ihigh = j 230 interval = tmp 231 ENDIF 232 ENDIF 233 ENDDO 234 ENDDO 235 IF (ilow == -1) THEN ! we didn't bracket a minimum yet, try something else 236 optimizer_type = outer_scf_optimizer_diis 237 CYCLE 238 ENDIF 239 scf_env%outer_scf%count(ilow) = scf_env%outer_scf%count(ilow) + 1 240 scf_env%outer_scf%count(ihigh) = scf_env%outer_scf%count(ihigh) + 1 241 scf_env%outer_scf%variables(:, ihistory + 1) = 0.5_dp*(scf_env%outer_scf%variables(:, ilow) + & 242 scf_env%outer_scf%variables(:, ihigh)) 243 CASE (outer_scf_optimizer_none) 244 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) 245 CASE (outer_scf_optimizer_sd) 246 ! Notice that we are just trying to find a stationary point 247 ! e.g. the ddpac_constraint, one maximizes the function, so the stepsize might have 248 ! to be negative 249 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - & 250 scf_control%outer_scf%step_size*scf_env%outer_scf%gradient(:, ihistory) 251 CASE (outer_scf_optimizer_diis) 252 CPASSERT(scf_control%outer_scf%diis_buffer_length > 0) 253 ! set up DIIS matrix 254 nb = MIN(ihistory, scf_control%outer_scf%diis_buffer_length) 255 IF (nb < 2) THEN 256 optimizer_type = outer_scf_optimizer_sd 257 CYCLE 258 ELSE 259 ALLOCATE (b(nb + 1, nb + 1), a(nb + 1, nb + 1), ev(nb + 1)) 260 DO I = 1, nb 261 DO J = I, nb 262 ibuf = ihistory - nb + i 263 jbuf = ihistory - nb + j 264 b(I, J) = DOT_PRODUCT(scf_env%outer_scf%gradient(:, ibuf), & 265 scf_env%outer_scf%gradient(:, jbuf)) 266 b(J, I) = b(I, J) 267 ENDDO 268 ENDDO 269 b(nb + 1, :) = -1.0_dp 270 b(:, nb + 1) = -1.0_dp 271 b(nb + 1, nb + 1) = 0.0_dp 272 273 CALL diamat_all(b, ev) 274 a(:, :) = b 275 DO I = 1, nb + 1 276 IF (ABS(ev(I)) .LT. 1.0E-12_dp) THEN 277 a(:, I) = 0.0_dp 278 ELSE 279 a(:, I) = a(:, I)/ev(I) 280 ENDIF 281 END DO 282 ev(:) = -MATMUL(a, b(nb + 1, :)) 283 284 scf_env%outer_scf%variables(:, ihistory + 1) = 0.0_dp 285 DO i = 1, nb 286 ibuf = ihistory - nb + i 287 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory + 1) + & 288 ev(i)*scf_env%outer_scf%variables(:, ibuf) 289 ENDDO 290 DEALLOCATE (a, b, ev) 291 ENDIF 292 CASE (outer_scf_optimizer_secant) 293 CPASSERT(SIZE(scf_env%outer_scf%gradient, 2) >= 3) 294 CPASSERT(SIZE(scf_env%outer_scf%gradient, 1) == 1) 295 nvar = SIZE(scf_env%outer_scf%gradient, 1) 296 IF (ihistory < 2) THEN 297 ! Need two history values to use secant, switch to sd 298 optimizer_type = outer_scf_optimizer_sd 299 CYCLE 300 END IF 301 ! secant update 302 scf_env%outer_scf%variables(1, ihistory + 1) = scf_env%outer_scf%variables(1, ihistory) - & 303 (scf_env%outer_scf%variables(1, ihistory) - & 304 scf_env%outer_scf%variables(1, ihistory - 1))/ & 305 (scf_env%outer_scf%gradient(1, ihistory) - & 306 scf_env%outer_scf%gradient(1, ihistory - 1))* & 307 scf_env%outer_scf%gradient(1, ihistory) 308 CASE (outer_scf_optimizer_broyden) 309 IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN 310 ! Inverse Jacobian not yet built, switch to sd 311 optimizer_type = outer_scf_optimizer_sd 312 CYCLE 313 END IF 314 inv_jacobian => scf_env%outer_scf%inv_jacobian 315 IF (ihistory < 2) THEN 316 ! Cannot perform a Broyden update without enough SCF history on this energy evaluation 317 scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE. 318 END IF 319 IF (scf_control%outer_scf%cdft_opt_control%broyden_update) THEN 320 ! Perform a Broyden update of the inverse Jacobian J^(-1) 321 IF (SIZE(scf_env%outer_scf%gradient, 2) .LT. 3) & 322 CALL cp_abort(__LOCATION__, & 323 "Keyword EXTRAPOLATION_ORDER in section OUTER_SCF "// & 324 "must be greater than or equal to 3 for Broyden optimizers.") 325 nvar = SIZE(scf_env%outer_scf%gradient, 1) 326 ALLOCATE (f(nvar, 1), x(nvar, 1)) 327 DO i = 1, nvar 328 f(i, 1) = scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1) 329 x(i, 1) = scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1) 330 END DO 331 SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type) 332 CASE (broyden_type_1, broyden_type_1_explicit, broyden_type_1_ls, broyden_type_1_explicit_ls) 333 ! Broyden's 1st method 334 ! Denote: dx_n = \delta x_n; df_n = \delta f_n 335 ! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(dx_n^T * J_n^(-1))/(dx_n^T * J_n^(-1) * df_n) 336 scale = SUM(MATMUL(TRANSPOSE(x), MATMUL(inv_jacobian, f))) 337 scale = 1.0_dp/scale 338 IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp 339 inv_jacobian = inv_jacobian + scale*MATMUL((x - MATMUL(inv_jacobian, f)), & 340 MATMUL(TRANSPOSE(x), inv_jacobian)) 341 CASE (broyden_type_2, broyden_type_2_explicit, broyden_type_2_ls, broyden_type_2_explicit_ls) 342 ! Broyden's 2nd method 343 ! J_(n+1)^(-1) = J_n^(-1) + (dx_n - J_n^(-1)*df_n)*(df_n^T)/(||df_n||^2) 344 scale = SUM(MATMUL(TRANSPOSE(f), f)) 345 scale = 1.0_dp/scale 346 IF (scale < 1.0E-12_dp) scale = 1.0E-12_dp 347 inv_jacobian = inv_jacobian + scale*MATMUL((x - MATMUL(inv_jacobian, f)), TRANSPOSE(inv_jacobian)) 348 CASE DEFAULT 349 CALL cp_abort(__LOCATION__, & 350 "Unknown Broyden type: "// & 351 cp_to_string(scf_control%outer_scf%cdft_opt_control%broyden_type)) 352 END SELECT 353 ! Clean up 354 DEALLOCATE (f, x) 355 END IF 356 ! Update variables x_(n+1) = x_n - J^(-1)*f(x_n) 357 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - & 358 scf_control%outer_scf%cdft_opt_control%newton_step* & 359 MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory)) 360 scf_control%outer_scf%cdft_opt_control%broyden_update = .TRUE. 361 CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls) 362 CPASSERT(ASSOCIATED(scf_env%outer_scf%inv_jacobian)) 363 inv_jacobian => scf_env%outer_scf%inv_jacobian 364 scf_env%outer_scf%variables(:, ihistory + 1) = scf_env%outer_scf%variables(:, ihistory) - & 365 scf_control%outer_scf%cdft_opt_control%newton_step* & 366 MATMUL(inv_jacobian, scf_env%outer_scf%gradient(:, ihistory)) 367 CASE DEFAULT 368 CPABORT("") 369 END SELECT 370 EXIT 371 ENDDO 372 END IF 373 374 CALL timestop(handle) 375 376 END SUBROUTINE outer_loop_optimize 377 378! ************************************************************************************************** 379!> \brief propagates the updated variables to wherever they need to be set in 380!> qs_env 381!> \param qs_env ... 382!> \param scf_env ... 383!> \par History 384!> 03.2006 created [Joost VandeVondele] 385! ************************************************************************************************** 386 SUBROUTINE outer_loop_update_qs_env(qs_env, scf_env) 387 TYPE(qs_environment_type), POINTER :: qs_env 388 TYPE(qs_scf_env_type), POINTER :: scf_env 389 390 CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_update_qs_env', & 391 routineP = moduleN//':'//routineN 392 393 INTEGER :: handle, ihistory, n 394 LOGICAL :: is_constraint 395 TYPE(cdft_control_type), POINTER :: cdft_control 396 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control 397 TYPE(dft_control_type), POINTER :: dft_control 398 TYPE(s2_restraint_type), POINTER :: s2_restraint_control 399 TYPE(scf_control_type), POINTER :: scf_control 400 401 CALL timeset(routineN, handle) 402 403 CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control) 404 ihistory = scf_env%outer_scf%iter_count 405 406 SELECT CASE (scf_control%outer_scf%type) 407 CASE (outer_scf_none) 408 ! do nothing 409 CASE (outer_scf_ddapc_constraint) 410 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control) 411 NULLIFY (ddapc_restraint_control) 412 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)%ddapc_restraint_control 413 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint) 414 IF (is_constraint) EXIT 415 END DO 416 ddapc_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1) 417 CASE (outer_scf_s2_constraint) 418 s2_restraint_control => dft_control%qs_control%s2_restraint_control 419 s2_restraint_control%strength = scf_env%outer_scf%variables(1, ihistory + 1) 420 CASE (outer_scf_cdft_constraint) 421 cdft_control => dft_control%qs_control%cdft_control 422 cdft_control%strength(:) = scf_env%outer_scf%variables(:, ihistory + 1) 423 CASE (outer_scf_basis_center_opt) 424 CALL qs_update_basis_center_pos(qs_env) 425 CASE DEFAULT 426 CPABORT("") 427 END SELECT 428 429 CALL timestop(handle) 430 431 END SUBROUTINE outer_loop_update_qs_env 432 433! ************************************************************************************************** 434!> \brief uses the outer_scf_history to extrapolate new values for the variables 435!> and updates their value in qs_env accordingly 436!> \param qs_env the qs_environment_type where to update the variables 437!> \par History 438!> 03.2006 created [Joost VandeVondele] 439!> \note 440!> it assumes that the current value of qs_env still needs to be added to the history 441!> simple multilinear extrapolation is employed 442! ************************************************************************************************** 443 SUBROUTINE outer_loop_extrapolate(qs_env) 444 TYPE(qs_environment_type), POINTER :: qs_env 445 446 CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_extrapolate', & 447 routineP = moduleN//':'//routineN 448 449 INTEGER :: handle, ihis, ivec, n, nhistory, & 450 nvariables, nvec, outer_scf_ihistory 451 LOGICAL :: is_constraint 452 REAL(kind=dp) :: alpha 453 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: extrapolation 454 REAL(kind=dp), DIMENSION(:, :), POINTER :: outer_scf_history 455 TYPE(ddapc_restraint_type), POINTER :: ddapc_restraint_control 456 TYPE(dft_control_type), POINTER :: dft_control 457 TYPE(scf_control_type), POINTER :: scf_control 458 459 CALL timeset(routineN, handle) 460 461 CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, & 462 outer_scf_ihistory=outer_scf_ihistory, & 463 scf_control=scf_control, dft_control=dft_control) 464 465 nvariables = SIZE(outer_scf_history, 1) 466 nhistory = SIZE(outer_scf_history, 2) 467 ALLOCATE (extrapolation(nvariables)) 468 CPASSERT(nhistory > 0) 469 470 ! add the current version of qs_env to the history 471 outer_scf_ihistory = outer_scf_ihistory + 1 472 ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory) 473 SELECT CASE (scf_control%outer_scf%type) 474 CASE (outer_scf_none) 475 outer_scf_history(1, ivec) = 0.0_dp 476 CASE (outer_scf_ddapc_constraint) 477 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control) 478 NULLIFY (ddapc_restraint_control) 479 ddapc_restraint_control => dft_control%qs_control%ddapc_restraint_control(n)%ddapc_restraint_control 480 is_constraint = (ddapc_restraint_control%functional_form == do_ddapc_constraint) 481 IF (is_constraint) EXIT 482 END DO 483 outer_scf_history(1, ivec) = & 484 ddapc_restraint_control%strength 485 CASE (outer_scf_s2_constraint) 486 outer_scf_history(1, ivec) = & 487 dft_control%qs_control%s2_restraint_control%strength 488 CASE (outer_scf_cdft_constraint) 489 outer_scf_history(:, ivec) = & 490 dft_control%qs_control%cdft_control%strength(:) 491 CASE (outer_scf_basis_center_opt) 492 outer_scf_history(1, ivec) = 0.0_dp 493 CASE DEFAULT 494 CPABORT("") 495 END SELECT 496 CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory) 497 ! multilinear extrapolation 498 nvec = MIN(nhistory, outer_scf_ihistory) 499 alpha = nvec 500 ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory) 501 extrapolation(:) = alpha*outer_scf_history(:, ivec) 502 DO ihis = 2, nvec 503 alpha = -1.0_dp*alpha*REAL(nvec - ihis + 1, dp)/REAL(ihis, dp) 504 ivec = 1 + MODULO(outer_scf_ihistory - ihis, nhistory) 505 extrapolation(:) = extrapolation + alpha*outer_scf_history(:, ivec) 506 ENDDO 507 508 ! update qs_env to use this extrapolation 509 SELECT CASE (scf_control%outer_scf%type) 510 CASE (outer_scf_none) 511 ! nothing 512 CASE (outer_scf_ddapc_constraint) 513 ddapc_restraint_control%strength = extrapolation(1) 514 CASE (outer_scf_s2_constraint) 515 dft_control%qs_control%s2_restraint_control%strength = extrapolation(1) 516 CASE (outer_scf_cdft_constraint) 517 dft_control%qs_control%cdft_control%strength(:) = extrapolation(:) 518 CASE (outer_scf_basis_center_opt) 519 ! nothing to do 520 CASE DEFAULT 521 CPABORT("") 522 END SELECT 523 524 DEALLOCATE (extrapolation) 525 526 CALL timestop(handle) 527 528 END SUBROUTINE outer_loop_extrapolate 529 530! ************************************************************************************************** 531!> \brief switch between two outer_scf envs stored in cdft_control 532!> \param scf_env the scf_env where values need to be updated using cdft_control 533!> \param scf_control the scf_control where values need to be updated using cdft_control 534!> \param cdft_control container for the second outer_scf env 535!> \param dir determines what switching operation to perform 536!> \par History 537!> 12.2015 created [Nico Holmberg] 538! ************************************************************************************************** 539 540 SUBROUTINE outer_loop_switch(scf_env, scf_control, cdft_control, dir) 541 TYPE(qs_scf_env_type), POINTER :: scf_env 542 TYPE(scf_control_type), POINTER :: scf_control 543 TYPE(cdft_control_type), POINTER :: cdft_control 544 INTEGER, INTENT(IN) :: dir 545 546 CHARACTER(len=*), PARAMETER :: routineN = 'outer_loop_switch', & 547 routineP = moduleN//':'//routineN 548 549 INTEGER :: nvariables 550 551 SELECT CASE (dir) 552 CASE (cdft2ot) 553 ! Constraint -> OT 554 ! Switch data in scf_control: first save values that might have changed 555 IF (ASSOCIATED(scf_control%outer_scf%cdft_opt_control)) THEN 556 CPASSERT(ASSOCIATED(cdft_control%constraint_control%cdft_opt_control)) 557 CALL cdft_opt_type_copy(cdft_control%constraint_control%cdft_opt_control, & 558 scf_control%outer_scf%cdft_opt_control) 559 ! OT SCF does not need cdft_opt_control 560 CALL cdft_opt_type_release(scf_control%outer_scf%cdft_opt_control) 561 END IF 562 ! Now switch 563 scf_control%outer_scf%have_scf = cdft_control%ot_control%have_scf 564 scf_control%outer_scf%max_scf = cdft_control%ot_control%max_scf 565 scf_control%outer_scf%eps_scf = cdft_control%ot_control%eps_scf 566 scf_control%outer_scf%step_size = cdft_control%ot_control%step_size 567 scf_control%outer_scf%type = cdft_control%ot_control%type 568 scf_control%outer_scf%optimizer = cdft_control%ot_control%optimizer 569 scf_control%outer_scf%diis_buffer_length = cdft_control%ot_control%diis_buffer_length 570 scf_control%outer_scf%bisect_trust_count = cdft_control%ot_control%bisect_trust_count 571 ! Switch data in scf_env: first save current values for constraint 572 cdft_control%constraint%iter_count = scf_env%outer_scf%iter_count 573 cdft_control%constraint%energy = scf_env%outer_scf%energy 574 cdft_control%constraint%variables = scf_env%outer_scf%variables 575 cdft_control%constraint%gradient = scf_env%outer_scf%gradient 576 cdft_control%constraint%count = scf_env%outer_scf%count 577 cdft_control%constraint%deallocate_jacobian = scf_env%outer_scf%deallocate_jacobian 578 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN 579 nvariables = SIZE(scf_env%outer_scf%inv_jacobian, 1) 580 IF (.NOT. ASSOCIATED(cdft_control%constraint%inv_jacobian)) & 581 ALLOCATE (cdft_control%constraint%inv_jacobian(nvariables, nvariables)) 582 cdft_control%constraint%inv_jacobian = scf_env%outer_scf%inv_jacobian 583 END IF 584 ! Now switch 585 IF (ASSOCIATED(scf_env%outer_scf%energy)) & 586 DEALLOCATE (scf_env%outer_scf%energy) 587 ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1)) 588 scf_env%outer_scf%energy = 0.0_dp 589 IF (ASSOCIATED(scf_env%outer_scf%variables)) & 590 DEALLOCATE (scf_env%outer_scf%variables) 591 ALLOCATE (scf_env%outer_scf%variables(1, scf_control%outer_scf%max_scf + 1)) 592 scf_env%outer_scf%variables = 0.0_dp 593 IF (ASSOCIATED(scf_env%outer_scf%gradient)) & 594 DEALLOCATE (scf_env%outer_scf%gradient) 595 ALLOCATE (scf_env%outer_scf%gradient(1, scf_control%outer_scf%max_scf + 1)) 596 scf_env%outer_scf%gradient = 0.0_dp 597 IF (ASSOCIATED(scf_env%outer_scf%count)) & 598 DEALLOCATE (scf_env%outer_scf%count) 599 ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1)) 600 scf_env%outer_scf%count = 0 601 ! OT SCF does not need Jacobian 602 scf_env%outer_scf%deallocate_jacobian = .TRUE. 603 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) & 604 DEALLOCATE (scf_env%outer_scf%inv_jacobian) 605 CASE (ot2cdft) 606 ! OT -> constraint 607 scf_control%outer_scf%have_scf = cdft_control%constraint_control%have_scf 608 scf_control%outer_scf%max_scf = cdft_control%constraint_control%max_scf 609 scf_control%outer_scf%eps_scf = cdft_control%constraint_control%eps_scf 610 scf_control%outer_scf%step_size = cdft_control%constraint_control%step_size 611 scf_control%outer_scf%type = cdft_control%constraint_control%type 612 scf_control%outer_scf%optimizer = cdft_control%constraint_control%optimizer 613 scf_control%outer_scf%diis_buffer_length = cdft_control%constraint_control%diis_buffer_length 614 scf_control%outer_scf%bisect_trust_count = cdft_control%constraint_control%bisect_trust_count 615 CALL cdft_opt_type_copy(scf_control%outer_scf%cdft_opt_control, & 616 cdft_control%constraint_control%cdft_opt_control) 617 nvariables = SIZE(cdft_control%constraint%variables, 1) 618 IF (ASSOCIATED(scf_env%outer_scf%energy)) & 619 DEALLOCATE (scf_env%outer_scf%energy) 620 ALLOCATE (scf_env%outer_scf%energy(scf_control%outer_scf%max_scf + 1)) 621 scf_env%outer_scf%energy = cdft_control%constraint%energy 622 IF (ASSOCIATED(scf_env%outer_scf%variables)) & 623 DEALLOCATE (scf_env%outer_scf%variables) 624 ALLOCATE (scf_env%outer_scf%variables(nvariables, scf_control%outer_scf%max_scf + 1)) 625 scf_env%outer_scf%variables = cdft_control%constraint%variables 626 IF (ASSOCIATED(scf_env%outer_scf%gradient)) & 627 DEALLOCATE (scf_env%outer_scf%gradient) 628 ALLOCATE (scf_env%outer_scf%gradient(nvariables, scf_control%outer_scf%max_scf + 1)) 629 scf_env%outer_scf%gradient = cdft_control%constraint%gradient 630 IF (ASSOCIATED(scf_env%outer_scf%count)) & 631 DEALLOCATE (scf_env%outer_scf%count) 632 ALLOCATE (scf_env%outer_scf%count(scf_control%outer_scf%max_scf + 1)) 633 scf_env%outer_scf%count = cdft_control%constraint%count 634 scf_env%outer_scf%iter_count = cdft_control%constraint%iter_count 635 scf_env%outer_scf%deallocate_jacobian = cdft_control%constraint%deallocate_jacobian 636 IF (ASSOCIATED(cdft_control%constraint%inv_jacobian)) THEN 637 IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) & 638 DEALLOCATE (scf_env%outer_scf%inv_jacobian) 639 ALLOCATE (scf_env%outer_scf%inv_jacobian(nvariables, nvariables)) 640 scf_env%outer_scf%inv_jacobian = cdft_control%constraint%inv_jacobian 641 END IF 642 CASE DEFAULT 643 CPABORT("") 644 END SELECT 645 646 END SUBROUTINE outer_loop_switch 647 648! ************************************************************************************************** 649!> \brief purges outer_scf_history zeroing everything except 650!> the latest value of the outer_scf variable stored in qs_control 651!> \param qs_env the qs_environment_type where to purge 652!> \par History 653!> 05.2016 created [Nico Holmberg] 654! ************************************************************************************************** 655 SUBROUTINE outer_loop_purge_history(qs_env) 656 TYPE(qs_environment_type), POINTER :: qs_env 657 658 CHARACTER(LEN=*), PARAMETER :: routineN = 'outer_loop_purge_history', & 659 routineP = moduleN//':'//routineN 660 661 INTEGER :: handle, outer_scf_ihistory 662 REAL(kind=dp), DIMENSION(:, :), POINTER :: gradient_history, outer_scf_history, & 663 variable_history 664 665 CALL timeset(routineN, handle) 666 667 CALL get_qs_env(qs_env, outer_scf_history=outer_scf_history, & 668 outer_scf_ihistory=outer_scf_ihistory, & 669 gradient_history=gradient_history, & 670 variable_history=variable_history) 671 CPASSERT(SIZE(outer_scf_history, 2) > 0) 672 outer_scf_ihistory = 0 673 outer_scf_history = 0.0_dp 674 gradient_history = 0.0_dp 675 variable_history = 0.0_dp 676 CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory) 677 678 CALL timestop(handle) 679 680 END SUBROUTINE outer_loop_purge_history 681 682END MODULE qs_outer_scf 683