LCOV - code coverage report
Current view: top level - src - qs_cdft_scf_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:71c3ab0) Lines: 87.2 % 187 163
Test Date: 2026-07-25 06:35:44 Functions: 100.0 % 6 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Auxiliary routines for performing a constrained DFT SCF run with Quickstep.
      10              : !> \par History
      11              : !>      - Separated some routines from qs_scf (03.2018) [Nico Holmberg]
      12              : !> \author Nico Holmberg (03.2018)
      13              : ! **************************************************************************************************
      14              : MODULE qs_cdft_scf_utils
      15              :    USE cp_control_types,                ONLY: dft_control_type
      16              :    USE cp_files,                        ONLY: close_file,&
      17              :                                               get_unit_number,&
      18              :                                               open_file
      19              :    USE cp_log_handling,                 ONLY: cp_logger_create,&
      20              :                                               cp_logger_type,&
      21              :                                               cp_to_string
      22              :    USE input_constants,                 ONLY: &
      23              :         broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
      24              :         broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
      25              :         jacobian_fd1, jacobian_fd1_backward, jacobian_fd1_central, jacobian_fd2, &
      26              :         jacobian_fd2_backward, outer_scf_optimizer_broyden, outer_scf_optimizer_newton, &
      27              :         outer_scf_optimizer_newton_ls
      28              :    USE kinds,                           ONLY: default_path_length,&
      29              :                                               dp
      30              :    USE mathlib,                         ONLY: invert_matrix
      31              :    USE message_passing,                 ONLY: mp_para_env_type
      32              :    USE qs_environment_types,            ONLY: get_qs_env,&
      33              :                                               qs_environment_type
      34              :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      35              :    USE scf_control_types,               ONLY: scf_control_type
      36              : #include "./base/base_uses.f90"
      37              : 
      38              :    IMPLICIT NONE
      39              : 
      40              :    PRIVATE
      41              : 
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_scf_utils'
      43              : 
      44              :    PUBLIC :: prepare_jacobian_stencil, build_diagonal_jacobian, &
      45              :              restart_inverse_jacobian, print_inverse_jacobian, &
      46              :              create_tmp_logger, initialize_inverse_jacobian
      47              : 
      48              : CONTAINS
      49              : 
      50              : ! **************************************************************************************************
      51              : !> \brief Prepares the finite difference stencil for computing the Jacobian. The constraints
      52              : !>        are re-evaluated by perturbing each constraint.
      53              : !> \param qs_env the qs_env where to build the Jacobian
      54              : !> \param output_unit the output unit number
      55              : !> \param nwork the number of perturbations to take in the negative direction
      56              : !> \param pwork the number of perturbations to take in the positive direction
      57              : !> \param coeff list of coefficients that determine how to sum up the various perturbations
      58              : !> \param step_multiplier list of values that determine how large steps to take for each perturbatio
      59              : !> \param dh total length of the interval to use for computing the finite difference derivatives
      60              : !> \par History
      61              : !>      03.2018 created [Nico Holmberg]
      62              : ! **************************************************************************************************
      63           62 :    SUBROUTINE prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, &
      64              :                                        coeff, step_multiplier, dh)
      65              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      66              :       INTEGER                                            :: output_unit, nwork, pwork
      67              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: coeff, step_multiplier, dh
      68              : 
      69              :       CHARACTER(len=15)                                  :: fmt_code
      70              :       INTEGER                                            :: ivar
      71              :       TYPE(dft_control_type), POINTER                    :: dft_control
      72              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
      73              :       TYPE(scf_control_type), POINTER                    :: scf_control
      74              : 
      75           62 :       NULLIFY (scf_env, scf_control, dft_control)
      76              : 
      77           62 :       CPASSERT(ASSOCIATED(qs_env))
      78              :       CALL get_qs_env(qs_env, scf_env=scf_env, &
      79              :                       scf_control=scf_control, &
      80           62 :                       dft_control=dft_control)
      81              : 
      82            4 :       IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) /= 1 .AND. &
      83           62 :           SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) /= SIZE(scf_env%outer_scf%variables, 1)) THEN
      84              :          CALL cp_abort(__LOCATION__, &
      85              :                        cp_to_string(SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step))// &
      86              :                        " values passed to keyword JACOBIAN_STEP, expected 1 or "// &
      87            0 :                        cp_to_string(SIZE(scf_env%outer_scf%variables, 1)))
      88              :       END IF
      89              : 
      90          186 :       ALLOCATE (dh(SIZE(scf_env%outer_scf%variables, 1)))
      91           62 :       IF (SIZE(dh) /= SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step)) THEN
      92           60 :          DO ivar = 1, SIZE(dh)
      93           60 :             dh(ivar) = scf_control%outer_scf%cdft_opt_control%jacobian_step(1)
      94              :          END DO
      95              :       ELSE
      96           88 :          dh(:) = scf_control%outer_scf%cdft_opt_control%jacobian_step
      97              :       END IF
      98              : 
      99           62 :       SELECT CASE (scf_control%outer_scf%cdft_opt_control%jacobian_type)
     100              :       CASE DEFAULT
     101              :          CALL cp_abort(__LOCATION__, &
     102              :                        "Unknown Jacobian type: "// &
     103            0 :                        cp_to_string(scf_control%outer_scf%cdft_opt_control%jacobian_type))
     104              :       CASE (jacobian_fd1)
     105              :          ! f'(x0) = [ f(x0+h) - f(x0) ] / h
     106           48 :          nwork = 0
     107           48 :          pwork = 1
     108           48 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     109           48 :          coeff(nwork) = -1.0_dp
     110           48 :          coeff(pwork) = 1.0_dp
     111          144 :          step_multiplier = 1.0_dp
     112              :       CASE (jacobian_fd1_backward)
     113              :          ! f'(x0) = [ f(x0) - f(x0-h) ] / h
     114            2 :          nwork = -1
     115            2 :          pwork = 0
     116            2 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     117            2 :          coeff(nwork) = -1.0_dp
     118            2 :          coeff(pwork) = 1.0_dp
     119            6 :          step_multiplier = -1.0_dp
     120              :       CASE (jacobian_fd2)
     121              :          ! f'(x0) = [ -f(x0+2h) + 4f(x0+h) - 3f(x0) ] / 2h
     122            2 :          nwork = 0
     123            2 :          pwork = 2
     124            2 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     125            2 :          coeff(0) = -3.0_dp
     126            2 :          coeff(1) = 4.0_dp
     127            2 :          coeff(2) = -1.0_dp
     128            2 :          step_multiplier(0) = 0.0_dp
     129            2 :          step_multiplier(1) = 1.0_dp
     130            2 :          step_multiplier(2) = 2.0_dp
     131            4 :          dh(:) = 2.0_dp*dh(:)
     132              :       CASE (jacobian_fd2_backward)
     133              :          ! f'(x0) = [ 3f(x0) - 4f(x0-h) + f(x0-2h) ] / 2h
     134            8 :          nwork = -2
     135            8 :          pwork = 0
     136            8 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     137            8 :          coeff(0) = 3.0_dp
     138            8 :          coeff(-1) = -4.0_dp
     139            8 :          coeff(-2) = 1.0_dp
     140            8 :          step_multiplier(0) = 0.0_dp
     141            8 :          step_multiplier(-1) = -1.0_dp
     142            8 :          step_multiplier(-2) = -2.0_dp
     143           16 :          dh(:) = 2.0_dp*dh(:)
     144              :       CASE (jacobian_fd1_central)
     145              :          ! f'(x0) = [ f(x0+h) - f(x0-h) ] / 2h
     146            2 :          nwork = -1
     147            2 :          pwork = 1
     148            2 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     149            2 :          coeff(0) = 0.0_dp
     150            2 :          coeff(nwork) = -1.0_dp
     151            2 :          coeff(pwork) = 1.0_dp
     152            2 :          step_multiplier(0) = 0.0_dp
     153            2 :          step_multiplier(nwork) = -1.0_dp
     154            2 :          step_multiplier(pwork) = 1.0_dp
     155           66 :          dh(:) = 2.0_dp*dh(:)
     156              :       END SELECT
     157              :       ! Print some info
     158           62 :       IF (output_unit > 0) THEN
     159              :          WRITE (output_unit, FMT="(/,A)") &
     160           31 :             " ================================== JACOBIAN CALCULATION ================================="
     161              :          WRITE (output_unit, FMT="(A)") &
     162           31 :             " Evaluating inverse Jacobian using finite differences"
     163              :          WRITE (output_unit, '(A,I10,A,I10)') &
     164           31 :             " Energy evaluation: ", dft_control%qs_control%cdft_control%ienergy, &
     165           62 :             ", CDFT SCF iteration: ", scf_env%outer_scf%iter_count
     166           55 :          SELECT CASE (scf_control%outer_scf%cdft_opt_control%jacobian_type)
     167              :          CASE (jacobian_fd1)
     168           24 :             WRITE (output_unit, '(A)') " Type               : First order forward difference"
     169              :          CASE (jacobian_fd1_backward)
     170            1 :             WRITE (output_unit, '(A)') " Type               : First order backward difference"
     171              :          CASE (jacobian_fd2)
     172            1 :             WRITE (output_unit, '(A)') " Type               : Second order forward difference"
     173              :          CASE (jacobian_fd2_backward)
     174            4 :             WRITE (output_unit, '(A)') " Type               : Second order backward difference"
     175              :          CASE (jacobian_fd1_central)
     176            1 :             WRITE (output_unit, '(A)') " Type               : First order central difference"
     177              :          CASE DEFAULT
     178              :             CALL cp_abort(__LOCATION__, "Unknown Jacobian type: "// &
     179           31 :                           cp_to_string(scf_control%outer_scf%cdft_opt_control%jacobian_type))
     180              :          END SELECT
     181           31 :          IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN
     182           58 :             WRITE (output_unit, '(A,ES12.4)') " Step size          : ", scf_control%outer_scf%cdft_opt_control%jacobian_step
     183              :          ELSE
     184              :             WRITE (output_unit, '(A,ES12.4,A)') &
     185            2 :                " Step sizes         : ", scf_control%outer_scf%cdft_opt_control%jacobian_step(1), ' (constraint 1)'
     186            2 :             IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) < 10) THEN
     187            2 :                fmt_code = '(ES34.4,A,I2,A)'
     188              :             ELSE
     189            0 :                fmt_code = '(ES34.4,A,I3,A)'
     190              :             END IF
     191            4 :             DO ivar = 2, SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step)
     192            4 :                WRITE (output_unit, fmt_code) scf_control%outer_scf%cdft_opt_control%jacobian_step(ivar), " (constraint", ivar, ")"
     193              :             END DO
     194              :          END IF
     195              :       END IF
     196              : 
     197           62 :    END SUBROUTINE prepare_jacobian_stencil
     198              : ! **************************************************************************************************
     199              : !> \brief Builds a strictly diagonal inverse Jacobian from MD/SCF history.
     200              : !> \param qs_env the qs_environment_type where to compute the Jacobian
     201              : !> \param used_history flag that determines if history was actually used to prepare the Jacobian
     202              : !> \par History
     203              : !>      03.2018 created [Nico Holmberg]
     204              : ! **************************************************************************************************
     205           16 :    SUBROUTINE build_diagonal_jacobian(qs_env, used_history)
     206              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     207              :       LOGICAL                                            :: used_history
     208              : 
     209              :       INTEGER                                            :: i, ihistory, nvar, outer_scf_ihistory
     210              :       LOGICAL                                            :: use_md_history
     211              :       REAL(KIND=dp)                                      :: inv_error
     212           16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: jacobian
     213           16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: gradient_history, inv_jacobian, &
     214           16 :                                                             variable_history
     215              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     216              :       TYPE(scf_control_type), POINTER                    :: scf_control
     217              : 
     218           16 :       NULLIFY (scf_control, scf_env)
     219              : 
     220              :       CALL get_qs_env(qs_env, scf_env=scf_env, &
     221              :                       scf_control=scf_control, &
     222           16 :                       outer_scf_ihistory=outer_scf_ihistory)
     223           16 :       ihistory = scf_env%outer_scf%iter_count
     224              : 
     225           16 :       IF (outer_scf_ihistory >= 3 .AND. .NOT. used_history) THEN
     226              :          ! First, lets try using the history from previous energy evaluations
     227              :          CALL get_qs_env(qs_env, gradient_history=gradient_history, &
     228            0 :                          variable_history=variable_history)
     229            0 :          nvar = SIZE(scf_env%outer_scf%variables, 1)
     230            0 :          use_md_history = .TRUE.
     231              :          ! Check that none of the history values are identical in which case we should try something different
     232            0 :          DO i = 1, nvar
     233            0 :             IF (ABS(variable_history(i, 2) - variable_history(i, 1)) < 1.0E-12_dp) THEN
     234            0 :                use_md_history = .FALSE.
     235              :             END IF
     236              :          END DO
     237            0 :          IF (use_md_history) THEN
     238            0 :             ALLOCATE (jacobian(nvar, nvar))
     239            0 :             DO i = 1, nvar
     240              :                jacobian(i, i) = (gradient_history(i, 2) - gradient_history(i, 1))/ &
     241            0 :                                 (variable_history(i, 2) - variable_history(i, 1))
     242              :             END DO
     243            0 :             IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     244            0 :                ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
     245              :             END IF
     246            0 :             inv_jacobian => scf_env%outer_scf%inv_jacobian
     247            0 :             CALL invert_matrix(jacobian, inv_jacobian, inv_error)
     248            0 :             DEALLOCATE (jacobian)
     249              :             ! Mark that an inverse Jacobian was just built and the next outer_loop_optimize should not perform
     250              :             ! a Broyden update of it
     251            0 :             scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
     252              :             ! Mark that the MD history has been used and should not be reused anymore on this energy evaluation
     253            0 :             used_history = .TRUE.
     254              :          END IF
     255              :       END IF
     256           16 :       IF (ihistory >= 2 .AND. .NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     257              :          ! Next, try history from current SCF procedure
     258            8 :          nvar = SIZE(scf_env%outer_scf%variables, 1)
     259            8 :          IF (SIZE(scf_env%outer_scf%gradient, 2) < 3) THEN
     260              :             CALL cp_abort(__LOCATION__, &
     261              :                           "Keyword EXTRAPOLATION_ORDER in section OUTER_SCF must be greater than or equal "// &
     262            0 :                           "to 3 for optimizers that build the Jacobian from SCF history.")
     263              :          END IF
     264           32 :          ALLOCATE (jacobian(nvar, nvar))
     265           16 :          DO i = 1, nvar
     266              :             jacobian(i, i) = (scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1))/ &
     267           16 :                              (scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1))
     268              :          END DO
     269            8 :          IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     270           24 :             ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
     271              :          END IF
     272            8 :          inv_jacobian => scf_env%outer_scf%inv_jacobian
     273            8 :          CALL invert_matrix(jacobian, inv_jacobian, inv_error)
     274            8 :          DEALLOCATE (jacobian)
     275           16 :          scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
     276              :       ELSE
     277              :          ! No history => will fall back to SD optimizer in outer_loop_optimize
     278              :       END IF
     279              : 
     280           32 :    END SUBROUTINE build_diagonal_jacobian
     281              : 
     282              : ! **************************************************************************************************
     283              : !> \brief Restarts the finite difference inverse Jacobian.
     284              : !> \param qs_env the qs_environment_type where to compute the Jacobian
     285              : !> \par History
     286              : !>      03.2018 created [Nico Holmberg]
     287              : ! **************************************************************************************************
     288            6 :    SUBROUTINE restart_inverse_jacobian(qs_env)
     289              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     290              : 
     291              :       INTEGER                                            :: i, iwork, j, nvar
     292            6 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: inv_jacobian
     293              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     294              :       TYPE(scf_control_type), POINTER                    :: scf_control
     295              : 
     296            6 :       NULLIFY (scf_env, scf_control)
     297            0 :       CPASSERT(ASSOCIATED(qs_env))
     298            6 :       CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control)
     299              : 
     300            6 :       CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control%jacobian_vector))
     301            6 :       nvar = SIZE(scf_env%outer_scf%variables, 1)
     302            6 :       IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_vector) /= nvar**2) THEN
     303              :          CALL cp_abort(__LOCATION__, &
     304            0 :                        "Too many or too few values defined for restarting inverse Jacobian.")
     305              :       END IF
     306            6 :       IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     307           24 :          ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
     308              :       END IF
     309            6 :       inv_jacobian => scf_env%outer_scf%inv_jacobian
     310            6 :       iwork = 1
     311           16 :       DO i = 1, nvar
     312           34 :          DO j = 1, nvar
     313           18 :             inv_jacobian(i, j) = scf_control%outer_scf%cdft_opt_control%jacobian_vector(iwork)
     314           28 :             iwork = iwork + 1
     315              :          END DO
     316              :       END DO
     317            6 :       DEALLOCATE (scf_control%outer_scf%cdft_opt_control%jacobian_vector)
     318            6 :       scf_control%outer_scf%cdft_opt_control%jacobian_restart = .FALSE.
     319            6 :       scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
     320            6 :       scf_env%outer_scf%deallocate_jacobian = .FALSE.
     321              : 
     322            6 :    END SUBROUTINE restart_inverse_jacobian
     323              : 
     324              : ! **************************************************************************************************
     325              : !> \brief Prints the finite difference inverse Jacobian to file
     326              : !> \param logger the default IO logger
     327              : !> \param inv_jacobian the inverse Jacobian matrix
     328              : !> \param iter_count the iteration number
     329              : !> \par History
     330              : !>      03.2018 created [Nico Holmberg]
     331              : ! **************************************************************************************************
     332           55 :    SUBROUTINE print_inverse_jacobian(logger, inv_jacobian, iter_count)
     333              :       TYPE(cp_logger_type), POINTER                      :: logger
     334              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     335              :          POINTER                                         :: inv_jacobian
     336              :       INTEGER                                            :: iter_count
     337              : 
     338              :       CHARACTER(len=default_path_length)                 :: project_name
     339              :       INTEGER                                            :: i, j, lp, nvar, output_unit
     340              : 
     341           55 :       nvar = SIZE(inv_jacobian, 1)
     342           55 :       output_unit = get_unit_number()
     343           55 :       project_name = logger%iter_info%project_name
     344           55 :       lp = LEN_TRIM(project_name)
     345           55 :       project_name(lp + 1:LEN(project_name)) = ".inverseJacobian"
     346              :       CALL open_file(file_name=project_name, file_status="UNKNOWN", &
     347              :                      file_action="WRITE", file_position="APPEND", &
     348           55 :                      unit_number=output_unit)
     349           55 :       WRITE (output_unit, FMT="(/,A)") "Inverse Jacobian matrix in row major order"
     350           55 :       WRITE (output_unit, FMT="(A,I10)") "Iteration: ", iter_count
     351          138 :       DO i = 1, nvar
     352          277 :          DO j = 1, nvar
     353          222 :             WRITE (output_unit, *) inv_jacobian(i, j)
     354              :          END DO
     355              :       END DO
     356           55 :       CALL close_file(unit_number=output_unit)
     357              : 
     358           55 :    END SUBROUTINE print_inverse_jacobian
     359              : 
     360              : ! **************************************************************************************************
     361              : !> \brief Creates a temporary logger for redirecting output to a new file
     362              : !> \param para_env the para_env
     363              : !> \param project_name the project basename
     364              : !> \param suffix the suffix
     365              : !> \param output_unit the default unit number for the newly created temporary logger
     366              : !> \param tmp_logger pointer to the newly created temporary logger
     367              : !> \par History
     368              : !>      03.2018 created [Nico Holmberg]
     369              : ! **************************************************************************************************
     370           70 :    SUBROUTINE create_tmp_logger(para_env, project_name, suffix, output_unit, tmp_logger)
     371              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     372              :       CHARACTER(len=*)                                   :: project_name, suffix
     373              :       INTEGER, INTENT(OUT)                               :: output_unit
     374              :       TYPE(cp_logger_type), INTENT(OUT), POINTER         :: tmp_logger
     375              : 
     376              :       INTEGER                                            :: lp
     377              : 
     378           70 :       IF (para_env%is_source()) THEN
     379           35 :          lp = LEN_TRIM(project_name)
     380           35 :          project_name(lp + 1:LEN(project_name)) = suffix
     381              :          CALL open_file(file_name=project_name, file_status="UNKNOWN", &
     382              :                         file_action="WRITE", file_position="APPEND", &
     383           35 :                         unit_number=output_unit)
     384              :       ELSE
     385           35 :          output_unit = -1
     386              :       END IF
     387              :       CALL cp_logger_create(tmp_logger, &
     388              :                             para_env=para_env, &
     389              :                             default_global_unit_nr=output_unit, &
     390           70 :                             close_global_unit_on_dealloc=.FALSE.)
     391              : 
     392           70 :    END SUBROUTINE create_tmp_logger
     393              : 
     394              : ! **************************************************************************************************
     395              : !> \brief Checks if the inverse Jacobian should be calculated and initializes the calculation
     396              : !> \param scf_control         the scf_control that holds the Jacobian settings
     397              : !> \param scf_env             the scf_env that holds the CDFT iteration information
     398              : !> \param explicit_jacobian   flag that determines if the finite difference Jacobian is needed
     399              : !> \param should_build        flag that determines if the Jacobian should be built
     400              : !> \param used_history        flag that determines if SCF history has been used to build a Jacobian
     401              : !> \par History
     402              : !>      03.2018 created [Nico Holmberg]
     403              : ! **************************************************************************************************
     404          118 :    SUBROUTINE initialize_inverse_jacobian(scf_control, scf_env, explicit_jacobian, &
     405              :                                           should_build, used_history)
     406              :       TYPE(scf_control_type), POINTER                    :: scf_control
     407              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     408              :       LOGICAL                                            :: explicit_jacobian, should_build, &
     409              :                                                             used_history
     410              : 
     411          118 :       CPASSERT(ASSOCIATED(scf_control))
     412          118 :       CPASSERT(ASSOCIATED(scf_env))
     413              : 
     414          118 :       SELECT CASE (scf_control%outer_scf%optimizer)
     415              :       CASE DEFAULT
     416            0 :          CPABORT("Noncompatible optimizer requested.")
     417              :       CASE (outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls)
     418           94 :          CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
     419           94 :          scf_control%outer_scf%cdft_opt_control%build_jacobian = .TRUE.
     420           94 :          explicit_jacobian = .TRUE.
     421              :       CASE (outer_scf_optimizer_broyden)
     422           24 :          CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
     423          118 :          SELECT CASE (scf_control%outer_scf%cdft_opt_control%broyden_type)
     424              :          CASE (broyden_type_1, broyden_type_2, broyden_type_1_ls, broyden_type_2_ls)
     425           20 :             scf_control%outer_scf%cdft_opt_control%build_jacobian = .TRUE.
     426           20 :             explicit_jacobian = .FALSE.
     427              :          CASE (broyden_type_1_explicit, broyden_type_2_explicit, broyden_type_1_explicit_ls, broyden_type_2_explicit_ls)
     428            4 :             scf_control%outer_scf%cdft_opt_control%build_jacobian = .TRUE.
     429           24 :             explicit_jacobian = .TRUE.
     430              :          END SELECT
     431              :       END SELECT
     432          118 :       IF (scf_control%outer_scf%cdft_opt_control%build_jacobian) THEN
     433              :          ! Reset counter
     434          118 :          IF (scf_env%outer_scf%iter_count == 1) scf_control%outer_scf%cdft_opt_control%ijacobian(1) = 0
     435              :          ! Check if an old Jacobian can be reused avoiding a rebuild
     436          118 :          IF (ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     437              :             ! Rebuild if number of previous energy evaluations exceeds limit
     438              :             IF (scf_control%outer_scf%cdft_opt_control%ijacobian(2) >= &
     439           62 :                 scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) .AND. .NOT. used_history .AND. &
     440              :                 scf_control%outer_scf%cdft_opt_control%jacobian_freq(2) > 0) THEN
     441            4 :                should_build = .TRUE.
     442              :                ! Zero the corresponding counters
     443           12 :                scf_control%outer_scf%cdft_opt_control%ijacobian(:) = 0
     444              :                ! Rebuild if number of previous SCF iterations exceeds limit (on this energy eval)
     445              :             ELSE IF (scf_control%outer_scf%cdft_opt_control%ijacobian(1) >= &
     446           58 :                      scf_control%outer_scf%cdft_opt_control%jacobian_freq(1) .AND. &
     447              :                      scf_control%outer_scf%cdft_opt_control%jacobian_freq(1) > 0) THEN
     448           24 :                should_build = .TRUE.
     449              :                ! Zero the corresponding counter
     450           24 :                scf_control%outer_scf%cdft_opt_control%ijacobian(1) = 0
     451              :             END IF
     452           62 :             IF (should_build) DEALLOCATE (scf_env%outer_scf%inv_jacobian)
     453              :          ELSE
     454           56 :             should_build = .TRUE.
     455              :             ! Zero the counter
     456          168 :             scf_control%outer_scf%cdft_opt_control%ijacobian(:) = 0
     457              :          END IF
     458              :       END IF
     459              : 
     460          118 :    END SUBROUTINE initialize_inverse_jacobian
     461              : 
     462              : END MODULE qs_cdft_scf_utils
        

Generated by: LCOV version 2.0-1