LCOV - code coverage report
Current view: top level - src - qs_outer_scf.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 259 290 89.3 %
Date: 2024-03-29 07:50:05 Functions: 6 7 85.7 %

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

Generated by: LCOV version 1.15