LCOV - code coverage report
Current view: top level - src - qs_cdft_scf_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:9ccafda) Lines: 87.1 % 186 162
Test Date: 2026-10-15 00:55: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              :                                               open_file
      18              :    USE cp_log_handling,                 ONLY: cp_logger_create,&
      19              :                                               cp_logger_type,&
      20              :                                               cp_to_string
      21              :    USE input_constants,                 ONLY: &
      22              :         broyden_type_1, broyden_type_1_explicit, broyden_type_1_explicit_ls, broyden_type_1_ls, &
      23              :         broyden_type_2, broyden_type_2_explicit, broyden_type_2_explicit_ls, broyden_type_2_ls, &
      24              :         jacobian_fd1, jacobian_fd1_backward, jacobian_fd1_central, jacobian_fd2, &
      25              :         jacobian_fd2_backward, outer_scf_optimizer_broyden, outer_scf_optimizer_newton, &
      26              :         outer_scf_optimizer_newton_ls
      27              :    USE kinds,                           ONLY: default_path_length,&
      28              :                                               dp
      29              :    USE mathlib,                         ONLY: invert_matrix
      30              :    USE message_passing,                 ONLY: mp_para_env_type
      31              :    USE qs_environment_types,            ONLY: get_qs_env,&
      32              :                                               qs_environment_type
      33              :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      34              :    USE scf_control_types,               ONLY: scf_control_type
      35              : #include "./base/base_uses.f90"
      36              : 
      37              :    IMPLICIT NONE
      38              : 
      39              :    PRIVATE
      40              : 
      41              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_scf_utils'
      42              : 
      43              :    PUBLIC :: prepare_jacobian_stencil, build_diagonal_jacobian, &
      44              :              restart_inverse_jacobian, print_inverse_jacobian, &
      45              :              create_tmp_logger, initialize_inverse_jacobian
      46              : 
      47              : CONTAINS
      48              : 
      49              : ! **************************************************************************************************
      50              : !> \brief Prepares the finite difference stencil for computing the Jacobian. The constraints
      51              : !>        are re-evaluated by perturbing each constraint.
      52              : !> \param qs_env the qs_env where to build the Jacobian
      53              : !> \param output_unit the output unit number
      54              : !> \param nwork the number of perturbations to take in the negative direction
      55              : !> \param pwork the number of perturbations to take in the positive direction
      56              : !> \param coeff list of coefficients that determine how to sum up the various perturbations
      57              : !> \param step_multiplier list of values that determine how large steps to take for each perturbatio
      58              : !> \param dh total length of the interval to use for computing the finite difference derivatives
      59              : !> \par History
      60              : !>      03.2018 created [Nico Holmberg]
      61              : ! **************************************************************************************************
      62           62 :    SUBROUTINE prepare_jacobian_stencil(qs_env, output_unit, nwork, pwork, &
      63              :                                        coeff, step_multiplier, dh)
      64              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      65              :       INTEGER                                            :: output_unit, nwork, pwork
      66              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: coeff, step_multiplier, dh
      67              : 
      68              :       CHARACTER(len=15)                                  :: fmt_code
      69              :       INTEGER                                            :: ivar
      70              :       TYPE(dft_control_type), POINTER                    :: dft_control
      71              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
      72              :       TYPE(scf_control_type), POINTER                    :: scf_control
      73              : 
      74           62 :       NULLIFY (scf_env, scf_control, dft_control)
      75              : 
      76           62 :       CPASSERT(ASSOCIATED(qs_env))
      77              :       CALL get_qs_env(qs_env, scf_env=scf_env, &
      78              :                       scf_control=scf_control, &
      79           62 :                       dft_control=dft_control)
      80              : 
      81            4 :       IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) /= 1 .AND. &
      82           62 :           SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) /= SIZE(scf_env%outer_scf%variables, 1)) THEN
      83              :          CALL cp_abort(__LOCATION__, &
      84              :                        cp_to_string(SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step))// &
      85              :                        " values passed to keyword JACOBIAN_STEP, expected 1 or "// &
      86            0 :                        cp_to_string(SIZE(scf_env%outer_scf%variables, 1)))
      87              :       END IF
      88              : 
      89          186 :       ALLOCATE (dh(SIZE(scf_env%outer_scf%variables, 1)))
      90           62 :       IF (SIZE(dh) /= SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step)) THEN
      91           60 :          DO ivar = 1, SIZE(dh)
      92           60 :             dh(ivar) = scf_control%outer_scf%cdft_opt_control%jacobian_step(1)
      93              :          END DO
      94              :       ELSE
      95           88 :          dh(:) = scf_control%outer_scf%cdft_opt_control%jacobian_step
      96              :       END IF
      97              : 
      98           62 :       SELECT CASE (scf_control%outer_scf%cdft_opt_control%jacobian_type)
      99              :       CASE DEFAULT
     100              :          CALL cp_abort(__LOCATION__, &
     101              :                        "Unknown Jacobian type: "// &
     102            0 :                        cp_to_string(scf_control%outer_scf%cdft_opt_control%jacobian_type))
     103              :       CASE (jacobian_fd1)
     104              :          ! f'(x0) = [ f(x0+h) - f(x0) ] / h
     105           48 :          nwork = 0
     106           48 :          pwork = 1
     107           48 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     108           48 :          coeff(nwork) = -1.0_dp
     109           48 :          coeff(pwork) = 1.0_dp
     110          144 :          step_multiplier = 1.0_dp
     111              :       CASE (jacobian_fd1_backward)
     112              :          ! f'(x0) = [ f(x0) - f(x0-h) ] / h
     113            2 :          nwork = -1
     114            2 :          pwork = 0
     115            2 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     116            2 :          coeff(nwork) = -1.0_dp
     117            2 :          coeff(pwork) = 1.0_dp
     118            6 :          step_multiplier = -1.0_dp
     119              :       CASE (jacobian_fd2)
     120              :          ! f'(x0) = [ -f(x0+2h) + 4f(x0+h) - 3f(x0) ] / 2h
     121            2 :          nwork = 0
     122            2 :          pwork = 2
     123            2 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     124            2 :          coeff(0) = -3.0_dp
     125            2 :          coeff(1) = 4.0_dp
     126            2 :          coeff(2) = -1.0_dp
     127            2 :          step_multiplier(0) = 0.0_dp
     128            2 :          step_multiplier(1) = 1.0_dp
     129            2 :          step_multiplier(2) = 2.0_dp
     130            4 :          dh(:) = 2.0_dp*dh(:)
     131              :       CASE (jacobian_fd2_backward)
     132              :          ! f'(x0) = [ 3f(x0) - 4f(x0-h) + f(x0-2h) ] / 2h
     133            8 :          nwork = -2
     134            8 :          pwork = 0
     135            8 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     136            8 :          coeff(0) = 3.0_dp
     137            8 :          coeff(-1) = -4.0_dp
     138            8 :          coeff(-2) = 1.0_dp
     139            8 :          step_multiplier(0) = 0.0_dp
     140            8 :          step_multiplier(-1) = -1.0_dp
     141            8 :          step_multiplier(-2) = -2.0_dp
     142           16 :          dh(:) = 2.0_dp*dh(:)
     143              :       CASE (jacobian_fd1_central)
     144              :          ! f'(x0) = [ f(x0+h) - f(x0-h) ] / 2h
     145            2 :          nwork = -1
     146            2 :          pwork = 1
     147            2 :          ALLOCATE (coeff(nwork:pwork), step_multiplier(nwork:pwork))
     148            2 :          coeff(0) = 0.0_dp
     149            2 :          coeff(nwork) = -1.0_dp
     150            2 :          coeff(pwork) = 1.0_dp
     151            2 :          step_multiplier(0) = 0.0_dp
     152            2 :          step_multiplier(nwork) = -1.0_dp
     153            2 :          step_multiplier(pwork) = 1.0_dp
     154           66 :          dh(:) = 2.0_dp*dh(:)
     155              :       END SELECT
     156              :       ! Print some info
     157           62 :       IF (output_unit > 0) THEN
     158              :          WRITE (output_unit, FMT="(/,A)") &
     159           31 :             " ================================== JACOBIAN CALCULATION ================================="
     160              :          WRITE (output_unit, FMT="(A)") &
     161           31 :             " Evaluating inverse Jacobian using finite differences"
     162              :          WRITE (output_unit, '(A,I10,A,I10)') &
     163           31 :             " Energy evaluation: ", dft_control%qs_control%cdft_control%ienergy, &
     164           62 :             ", CDFT SCF iteration: ", scf_env%outer_scf%iter_count
     165           55 :          SELECT CASE (scf_control%outer_scf%cdft_opt_control%jacobian_type)
     166              :          CASE (jacobian_fd1)
     167           24 :             WRITE (output_unit, '(A)') " Type               : First order forward difference"
     168              :          CASE (jacobian_fd1_backward)
     169            1 :             WRITE (output_unit, '(A)') " Type               : First order backward difference"
     170              :          CASE (jacobian_fd2)
     171            1 :             WRITE (output_unit, '(A)') " Type               : Second order forward difference"
     172              :          CASE (jacobian_fd2_backward)
     173            4 :             WRITE (output_unit, '(A)') " Type               : Second order backward difference"
     174              :          CASE (jacobian_fd1_central)
     175            1 :             WRITE (output_unit, '(A)') " Type               : First order central difference"
     176              :          CASE DEFAULT
     177              :             CALL cp_abort(__LOCATION__, "Unknown Jacobian type: "// &
     178           31 :                           cp_to_string(scf_control%outer_scf%cdft_opt_control%jacobian_type))
     179              :          END SELECT
     180           31 :          IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) == 1) THEN
     181           58 :             WRITE (output_unit, '(A,ES12.4)') " Step size          : ", scf_control%outer_scf%cdft_opt_control%jacobian_step
     182              :          ELSE
     183              :             WRITE (output_unit, '(A,ES12.4,A)') &
     184            2 :                " Step sizes         : ", scf_control%outer_scf%cdft_opt_control%jacobian_step(1), ' (constraint 1)'
     185            2 :             IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step) < 10) THEN
     186            2 :                fmt_code = '(ES34.4,A,I2,A)'
     187              :             ELSE
     188            0 :                fmt_code = '(ES34.4,A,I3,A)'
     189              :             END IF
     190            4 :             DO ivar = 2, SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_step)
     191            4 :                WRITE (output_unit, fmt_code) scf_control%outer_scf%cdft_opt_control%jacobian_step(ivar), " (constraint", ivar, ")"
     192              :             END DO
     193              :          END IF
     194              :       END IF
     195              : 
     196           62 :    END SUBROUTINE prepare_jacobian_stencil
     197              : ! **************************************************************************************************
     198              : !> \brief Builds a strictly diagonal inverse Jacobian from MD/SCF history.
     199              : !> \param qs_env the qs_environment_type where to compute the Jacobian
     200              : !> \param used_history flag that determines if history was actually used to prepare the Jacobian
     201              : !> \par History
     202              : !>      03.2018 created [Nico Holmberg]
     203              : ! **************************************************************************************************
     204           14 :    SUBROUTINE build_diagonal_jacobian(qs_env, used_history)
     205              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     206              :       LOGICAL                                            :: used_history
     207              : 
     208              :       INTEGER                                            :: i, ihistory, nvar, outer_scf_ihistory
     209              :       LOGICAL                                            :: use_md_history
     210              :       REAL(KIND=dp)                                      :: inv_error
     211           14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: jacobian
     212           14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: gradient_history, inv_jacobian, &
     213           14 :                                                             variable_history
     214              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     215              :       TYPE(scf_control_type), POINTER                    :: scf_control
     216              : 
     217           14 :       NULLIFY (scf_control, scf_env)
     218              : 
     219              :       CALL get_qs_env(qs_env, scf_env=scf_env, &
     220              :                       scf_control=scf_control, &
     221           14 :                       outer_scf_ihistory=outer_scf_ihistory)
     222           14 :       ihistory = scf_env%outer_scf%iter_count
     223              : 
     224           14 :       IF (outer_scf_ihistory >= 3 .AND. .NOT. used_history) THEN
     225              :          ! First, lets try using the history from previous energy evaluations
     226              :          CALL get_qs_env(qs_env, gradient_history=gradient_history, &
     227            0 :                          variable_history=variable_history)
     228            0 :          nvar = SIZE(scf_env%outer_scf%variables, 1)
     229            0 :          use_md_history = .TRUE.
     230              :          ! Check that none of the history values are identical in which case we should try something different
     231            0 :          DO i = 1, nvar
     232            0 :             IF (ABS(variable_history(i, 2) - variable_history(i, 1)) < 1.0E-12_dp) THEN
     233            0 :                use_md_history = .FALSE.
     234              :             END IF
     235              :          END DO
     236            0 :          IF (use_md_history) THEN
     237            0 :             ALLOCATE (jacobian(nvar, nvar))
     238            0 :             DO i = 1, nvar
     239              :                jacobian(i, i) = (gradient_history(i, 2) - gradient_history(i, 1))/ &
     240            0 :                                 (variable_history(i, 2) - variable_history(i, 1))
     241              :             END DO
     242            0 :             IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     243            0 :                ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
     244              :             END IF
     245            0 :             inv_jacobian => scf_env%outer_scf%inv_jacobian
     246            0 :             CALL invert_matrix(jacobian, inv_jacobian, inv_error)
     247            0 :             DEALLOCATE (jacobian)
     248              :             ! Mark that an inverse Jacobian was just built and the next outer_loop_optimize should not perform
     249              :             ! a Broyden update of it
     250            0 :             scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
     251              :             ! Mark that the MD history has been used and should not be reused anymore on this energy evaluation
     252            0 :             used_history = .TRUE.
     253              :          END IF
     254              :       END IF
     255           14 :       IF (ihistory >= 2 .AND. .NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     256              :          ! Next, try history from current SCF procedure
     257            7 :          nvar = SIZE(scf_env%outer_scf%variables, 1)
     258            7 :          IF (SIZE(scf_env%outer_scf%gradient, 2) < 3) THEN
     259              :             CALL cp_abort(__LOCATION__, &
     260              :                           "Keyword EXTRAPOLATION_ORDER in section OUTER_SCF must be greater than or equal "// &
     261            0 :                           "to 3 for optimizers that build the Jacobian from SCF history.")
     262              :          END IF
     263           28 :          ALLOCATE (jacobian(nvar, nvar))
     264           15 :          DO i = 1, nvar
     265              :             jacobian(i, i) = (scf_env%outer_scf%gradient(i, ihistory) - scf_env%outer_scf%gradient(i, ihistory - 1))/ &
     266           15 :                              (scf_env%outer_scf%variables(i, ihistory) - scf_env%outer_scf%variables(i, ihistory - 1))
     267              :          END DO
     268            7 :          IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     269           21 :             ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
     270              :          END IF
     271            7 :          inv_jacobian => scf_env%outer_scf%inv_jacobian
     272            7 :          CALL invert_matrix(jacobian, inv_jacobian, inv_error)
     273            6 :          DEALLOCATE (jacobian)
     274           13 :          scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
     275              :       ELSE
     276              :          ! No history => will fall back to SD optimizer in outer_loop_optimize
     277              :       END IF
     278              : 
     279           27 :    END SUBROUTINE build_diagonal_jacobian
     280              : 
     281              : ! **************************************************************************************************
     282              : !> \brief Restarts the finite difference inverse Jacobian.
     283              : !> \param qs_env the qs_environment_type where to compute the Jacobian
     284              : !> \par History
     285              : !>      03.2018 created [Nico Holmberg]
     286              : ! **************************************************************************************************
     287            6 :    SUBROUTINE restart_inverse_jacobian(qs_env)
     288              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     289              : 
     290              :       INTEGER                                            :: i, iwork, j, nvar
     291            6 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: inv_jacobian
     292              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     293              :       TYPE(scf_control_type), POINTER                    :: scf_control
     294              : 
     295            6 :       NULLIFY (scf_env, scf_control)
     296            0 :       CPASSERT(ASSOCIATED(qs_env))
     297            6 :       CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control)
     298              : 
     299            6 :       CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control%jacobian_vector))
     300            6 :       nvar = SIZE(scf_env%outer_scf%variables, 1)
     301            6 :       IF (SIZE(scf_control%outer_scf%cdft_opt_control%jacobian_vector) /= nvar**2) THEN
     302              :          CALL cp_abort(__LOCATION__, &
     303            0 :                        "Too many or too few values defined for restarting inverse Jacobian.")
     304              :       END IF
     305            6 :       IF (.NOT. ASSOCIATED(scf_env%outer_scf%inv_jacobian)) THEN
     306           24 :          ALLOCATE (scf_env%outer_scf%inv_jacobian(nvar, nvar))
     307              :       END IF
     308            6 :       inv_jacobian => scf_env%outer_scf%inv_jacobian
     309            6 :       iwork = 1
     310           16 :       DO i = 1, nvar
     311           34 :          DO j = 1, nvar
     312           18 :             inv_jacobian(i, j) = scf_control%outer_scf%cdft_opt_control%jacobian_vector(iwork)
     313           28 :             iwork = iwork + 1
     314              :          END DO
     315              :       END DO
     316            6 :       DEALLOCATE (scf_control%outer_scf%cdft_opt_control%jacobian_vector)
     317            6 :       scf_control%outer_scf%cdft_opt_control%jacobian_restart = .FALSE.
     318            6 :       scf_control%outer_scf%cdft_opt_control%broyden_update = .FALSE.
     319            6 :       scf_env%outer_scf%deallocate_jacobian = .FALSE.
     320              : 
     321            6 :    END SUBROUTINE restart_inverse_jacobian
     322              : 
     323              : ! **************************************************************************************************
     324              : !> \brief Prints the finite difference inverse Jacobian to file
     325              : !> \param logger the default IO logger
     326              : !> \param inv_jacobian the inverse Jacobian matrix
     327              : !> \param iter_count the iteration number
     328              : !> \par History
     329              : !>      03.2018 created [Nico Holmberg]
     330              : ! **************************************************************************************************
     331          110 :    SUBROUTINE print_inverse_jacobian(logger, inv_jacobian, iter_count)
     332              :       TYPE(cp_logger_type), POINTER                      :: logger
     333              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     334              :          POINTER                                         :: inv_jacobian
     335              :       INTEGER                                            :: iter_count
     336              : 
     337              :       CHARACTER(len=default_path_length)                 :: project_name
     338              :       INTEGER                                            :: i, j, lp, nvar, output_unit
     339              : 
     340           55 :       nvar = SIZE(inv_jacobian, 1)
     341           55 :       project_name = logger%iter_info%project_name
     342           55 :       lp = LEN_TRIM(project_name)
     343           55 :       project_name(lp + 1:LEN(project_name)) = ".inverseJacobian"
     344              :       CALL open_file(file_name=project_name, &
     345              :                      file_status="UNKNOWN", &
     346              :                      file_action="WRITE", &
     347              :                      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          105 :       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          114 :    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          114 :       CPASSERT(ASSOCIATED(scf_control))
     412          114 :       CPASSERT(ASSOCIATED(scf_env))
     413              : 
     414          114 :       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           20 :          CPASSERT(ASSOCIATED(scf_control%outer_scf%cdft_opt_control))
     423          114 :          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           16 :             scf_control%outer_scf%cdft_opt_control%build_jacobian = .TRUE.
     426           16 :             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           20 :             explicit_jacobian = .TRUE.
     430              :          END SELECT
     431              :       END SELECT
     432          114 :       IF (scf_control%outer_scf%cdft_opt_control%build_jacobian) THEN
     433              :          ! Reset counter
     434          114 :          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          114 :          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           60 :                 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           56 :                      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           60 :             IF (should_build) DEALLOCATE (scf_env%outer_scf%inv_jacobian)
     453              :          ELSE
     454           54 :             should_build = .TRUE.
     455              :             ! Zero the counter
     456          162 :             scf_control%outer_scf%cdft_opt_control%ijacobian(:) = 0
     457              :          END IF
     458              :       END IF
     459              : 
     460          114 :    END SUBROUTINE initialize_inverse_jacobian
     461              : 
     462              : END MODULE qs_cdft_scf_utils
        

Generated by: LCOV version 2.0-1