LCOV - code coverage report
Current view: top level - src - qs_outer_scf.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 89.3 % 290 259
Test Date: 2025-07-25 12:55:17 Functions: 85.7 % 7 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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         4978 :    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         4978 :       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         4978 :          res = 0
      91              :       END SELECT
      92              : 
      93         4978 :    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         5635 :    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         5635 :       CALL timeset(routineN, handle)
     118              : 
     119              :       CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, &
     120         5635 :                       dft_control=dft_control, energy=energy)
     121         5635 :       CPASSERT(scf_control%outer_scf%have_scf)
     122              : 
     123         5635 :       ihistory = scf_env%outer_scf%iter_count
     124         5635 :       CPASSERT(ihistory <= SIZE(scf_env%outer_scf%energy, 1))
     125              : 
     126         5635 :       scf_env%outer_scf%energy(ihistory) = energy%total
     127              : 
     128        10568 :       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         4933 :          scf_env%outer_scf%variables(1, ihistory) = scf_env%iter_delta
     132         4933 :          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         5635 :          CPABORT("")
     169              : 
     170              :       END SELECT
     171              : 
     172         5635 :       CALL timestop(handle)
     173              : 
     174         5635 :    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          976 :    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          976 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ev
     196          976 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: a, b, f, x
     197          976 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: inv_jacobian
     198              : 
     199          976 :       CALL timeset(routineN, handle)
     200              : 
     201          976 :       ihistory = scf_env%outer_scf%iter_count
     202          976 :       optimizer_type = scf_control%outer_scf%optimizer
     203          976 :       NULLIFY (inv_jacobian)
     204              : 
     205          976 :       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         1480 :                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          256 :                   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           24 :                   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         1016 :                CPABORT("")
     369              :             END SELECT
     370              :             EXIT
     371              :          END DO
     372              :       END IF
     373              : 
     374          976 :       CALL timestop(handle)
     375              : 
     376         1952 :    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         1090 :    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         1090 :       CALL timeset(routineN, handle)
     401              : 
     402         1090 :       CALL get_qs_env(qs_env=qs_env, scf_control=scf_control, dft_control=dft_control)
     403         1090 :       ihistory = scf_env%outer_scf%iter_count
     404              : 
     405         1140 :       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         1090 :          CPABORT("")
     426              :       END SELECT
     427              : 
     428         1090 :       CALL timestop(handle)
     429              : 
     430         1090 :    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         3925 :    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         3925 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: extrapolation
     452         3925 :       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         3925 :       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         3925 :                       scf_control=scf_control, dft_control=dft_control)
     462              : 
     463         3925 :       nvariables = SIZE(outer_scf_history, 1)
     464         3925 :       nhistory = SIZE(outer_scf_history, 2)
     465        11775 :       ALLOCATE (extrapolation(nvariables))
     466         3925 :       CPASSERT(nhistory > 0)
     467              : 
     468              :       ! add the current version of qs_env to the history
     469         3925 :       outer_scf_ihistory = outer_scf_ihistory + 1
     470         3925 :       ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory)
     471         7498 :       SELECT CASE (scf_control%outer_scf%type)
     472              :       CASE (outer_scf_none)
     473         3573 :          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         3925 :          CPABORT("")
     493              :       END SELECT
     494         3925 :       CALL set_qs_env(qs_env, outer_scf_ihistory=outer_scf_ihistory)
     495              :       ! multilinear extrapolation
     496         3925 :       nvec = MIN(nhistory, outer_scf_ihistory)
     497         3925 :       alpha = nvec
     498         3925 :       ivec = 1 + MODULO(outer_scf_ihistory - 1, nhistory)
     499         7880 :       extrapolation(:) = alpha*outer_scf_history(:, ivec)
     500         8685 :       DO ihis = 2, nvec
     501         4760 :          alpha = -1.0_dp*alpha*REAL(nvec - ihis + 1, dp)/REAL(ihis, dp)
     502         4760 :          ivec = 1 + MODULO(outer_scf_ihistory - ihis, nhistory)
     503        13453 :          extrapolation(:) = extrapolation + alpha*outer_scf_history(:, ivec)
     504              :       END DO
     505              : 
     506              :       ! update qs_env to use this extrapolation
     507         3951 :       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         3925 :          CPABORT("")
     520              :       END SELECT
     521              : 
     522         3925 :       DEALLOCATE (extrapolation)
     523              : 
     524         3925 :       CALL timestop(handle)
     525              : 
     526         3925 :    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 2.0-1