LCOV - code coverage report
Current view: top level - src/motion - cell_opt_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 168 182 92.3 %
Date: 2024-04-20 06:29:22 Functions: 5 5 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief contains a functional that calculates the energy and its derivatives
      10             : !>      for the geometry optimizer
      11             : !> \par History
      12             : !>      03.2008 - Teodoro Laino [tlaino] - University of Zurich - Cell Optimization
      13             : ! **************************************************************************************************
      14             : MODULE cell_opt_utils
      15             :    USE cell_types,                      ONLY: &
      16             :         cell_sym_cubic, cell_sym_hexagonal_gamma_120, cell_sym_hexagonal_gamma_60, &
      17             :         cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, cell_sym_orthorhombic, &
      18             :         cell_sym_rhombohedral, cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, &
      19             :         cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type
      20             :    USE cp_files,                        ONLY: close_file,&
      21             :                                               open_file
      22             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23             :                                               cp_logger_create,&
      24             :                                               cp_logger_get_default_unit_nr,&
      25             :                                               cp_logger_release,&
      26             :                                               cp_logger_set,&
      27             :                                               cp_logger_type,&
      28             :                                               cp_to_string
      29             :    USE input_constants,                 ONLY: fix_none,&
      30             :                                               fix_x,&
      31             :                                               fix_xy,&
      32             :                                               fix_xz,&
      33             :                                               fix_y,&
      34             :                                               fix_yz,&
      35             :                                               fix_z
      36             :    USE input_cp2k_global,               ONLY: create_global_section
      37             :    USE input_enumeration_types,         ONLY: enum_i2c,&
      38             :                                               enumeration_type
      39             :    USE input_keyword_types,             ONLY: keyword_get,&
      40             :                                               keyword_type
      41             :    USE input_section_types,             ONLY: section_get_keyword,&
      42             :                                               section_release,&
      43             :                                               section_type,&
      44             :                                               section_vals_type,&
      45             :                                               section_vals_val_get,&
      46             :                                               section_vals_val_set
      47             :    USE kinds,                           ONLY: default_path_length,&
      48             :                                               default_string_length,&
      49             :                                               dp
      50             :    USE mathconstants,                   ONLY: sqrt3
      51             :    USE mathlib,                         ONLY: angle
      52             :    USE message_passing,                 ONLY: mp_para_env_type
      53             : #include "../base/base_uses.f90"
      54             : 
      55             :    IMPLICIT NONE
      56             :    PRIVATE
      57             : 
      58             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      59             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_opt_utils'
      60             : 
      61             :    PUBLIC :: get_dg_dh, gopt_new_logger_create, &
      62             :              gopt_new_logger_release, read_external_press_tensor, &
      63             :              apply_cell_constraints
      64             : 
      65             : CONTAINS
      66             : 
      67             : ! **************************************************************************************************
      68             : !> \brief creates a new logger used for cell optimization algorithm
      69             : !> \param new_logger ...
      70             : !> \param root_section ...
      71             : !> \param para_env ...
      72             : !> \param project_name ...
      73             : !> \param id_run ...
      74             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
      75             : ! **************************************************************************************************
      76         908 :    SUBROUTINE gopt_new_logger_create(new_logger, root_section, para_env, project_name, &
      77             :                                      id_run)
      78             :       TYPE(cp_logger_type), POINTER                      :: new_logger
      79             :       TYPE(section_vals_type), POINTER                   :: root_section
      80             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      81             :       CHARACTER(len=default_string_length), INTENT(OUT)  :: project_name
      82             :       INTEGER, INTENT(IN)                                :: id_run
      83             : 
      84             :       CHARACTER(len=default_path_length)                 :: c_val, input_file_path, output_file_path
      85             :       CHARACTER(len=default_string_length)               :: label
      86             :       INTEGER                                            :: i, lp, unit_nr
      87             :       TYPE(cp_logger_type), POINTER                      :: logger
      88             :       TYPE(enumeration_type), POINTER                    :: enum
      89             :       TYPE(keyword_type), POINTER                        :: keyword
      90             :       TYPE(section_type), POINTER                        :: section
      91             : 
      92         908 :       NULLIFY (new_logger, logger, enum, keyword, section)
      93         908 :       logger => cp_get_default_logger()
      94             : 
      95         908 :       CALL create_global_section(section)
      96         908 :       keyword => section_get_keyword(section, "RUN_TYPE")
      97         908 :       CALL keyword_get(keyword, enum=enum)
      98         908 :       label = TRIM(enum_i2c(enum, id_run))
      99         908 :       CALL section_release(section)
     100             : 
     101             :       ! Redirecting output of the sub_calculation to a different file
     102         908 :       CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
     103         908 :       input_file_path = TRIM(project_name)
     104         908 :       lp = LEN_TRIM(input_file_path)
     105         908 :       i = logger%iter_info%iteration(logger%iter_info%n_rlevel)
     106         908 :       input_file_path(lp + 1:LEN(input_file_path)) = "-"//TRIM(label)//"-"//ADJUSTL(cp_to_string(i))
     107         908 :       lp = LEN_TRIM(input_file_path)
     108         908 :       CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=input_file_path(1:lp))
     109         908 :       CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
     110             : 
     111             :       ! Redirecting output into a new file
     112         908 :       output_file_path = input_file_path(1:lp)//".out"
     113         908 :       IF (para_env%is_source()) THEN
     114             :          CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
     115         454 :                         file_action="WRITE", file_position="APPEND", unit_number=unit_nr)
     116             :       ELSE
     117         454 :          unit_nr = -1
     118             :       END IF
     119             :       CALL cp_logger_create(new_logger, para_env=para_env, default_global_unit_nr=unit_nr, &
     120         908 :                             close_global_unit_on_dealloc=.FALSE.)
     121         908 :       CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
     122         908 :       IF (c_val /= "") THEN
     123         908 :          CALL cp_logger_set(new_logger, local_filename=TRIM(c_val)//"_localLog")
     124             :       END IF
     125         908 :       new_logger%iter_info%project_name = c_val
     126             :       CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
     127         908 :                                 i_val=new_logger%iter_info%print_level)
     128             : 
     129         908 :    END SUBROUTINE gopt_new_logger_create
     130             : 
     131             : ! **************************************************************************************************
     132             : !> \brief releases a new logger used for cell optimization algorithm
     133             : !> \param new_logger ...
     134             : !> \param root_section ...
     135             : !> \param para_env ...
     136             : !> \param project_name ...
     137             : !> \param id_run ...
     138             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     139             : ! **************************************************************************************************
     140         908 :    SUBROUTINE gopt_new_logger_release(new_logger, root_section, para_env, project_name, id_run)
     141             :       TYPE(cp_logger_type), POINTER                      :: new_logger
     142             :       TYPE(section_vals_type), POINTER                   :: root_section
     143             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     144             :       CHARACTER(len=default_string_length), INTENT(IN)   :: project_name
     145             :       INTEGER, INTENT(IN)                                :: id_run
     146             : 
     147             :       INTEGER                                            :: unit_nr
     148             : 
     149         908 :       IF (para_env%is_source()) THEN
     150         454 :          unit_nr = cp_logger_get_default_unit_nr(new_logger)
     151         454 :          CALL close_file(unit_number=unit_nr)
     152             :       END IF
     153         908 :       CALL cp_logger_release(new_logger)
     154         908 :       CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
     155         908 :       CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
     156             : 
     157         908 :    END SUBROUTINE gopt_new_logger_release
     158             : 
     159             : ! **************************************************************************************************
     160             : !> \brief Reads the external pressure tensor
     161             : !> \param geo_section ...
     162             : !> \param cell ...
     163             : !> \param pres_ext ...
     164             : !> \param mtrx ...
     165             : !> \param rot ...
     166             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     167             : ! **************************************************************************************************
     168         210 :    SUBROUTINE read_external_press_tensor(geo_section, cell, pres_ext, mtrx, rot)
     169             :       TYPE(section_vals_type), POINTER                   :: geo_section
     170             :       TYPE(cell_type), POINTER                           :: cell
     171             :       REAL(KIND=dp), INTENT(OUT)                         :: pres_ext
     172             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: mtrx
     173             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: rot
     174             : 
     175             :       INTEGER                                            :: i, ind, j
     176             :       LOGICAL                                            :: check
     177             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pres_ext_tens
     178         210 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pvals
     179             : 
     180         210 :       NULLIFY (pvals)
     181         210 :       mtrx = 0.0_dp
     182         210 :       pres_ext_tens = 0.0_dp
     183         210 :       pres_ext = 0.0_dp
     184         210 :       CALL section_vals_val_get(geo_section, "EXTERNAL_PRESSURE", r_vals=pvals)
     185         210 :       check = (SIZE(pvals) == 1) .OR. (SIZE(pvals) == 9)
     186         210 :       IF (.NOT. check) &
     187           0 :          CPABORT("EXTERNAL_PRESSURE can have 1 or 9 components only!")
     188             : 
     189         210 :       IF (SIZE(pvals) == 9) THEN
     190             :          ind = 0
     191         392 :          DO i = 1, 3
     192        1274 :             DO j = 1, 3
     193         882 :                ind = ind + 1
     194        1176 :                pres_ext_tens(j, i) = pvals(ind)
     195             :             END DO
     196             :          END DO
     197             :          ! Also the pressure tensor must be oriented in the same canonical directions
     198             :          ! of the simulation cell
     199        3920 :          pres_ext_tens = MATMUL(TRANSPOSE(rot), pres_ext_tens)
     200         392 :          DO i = 1, 3
     201         392 :             pres_ext = pres_ext + pres_ext_tens(i, i)
     202             :          END DO
     203          98 :          pres_ext = pres_ext/3.0_dp
     204         392 :          DO i = 1, 3
     205         392 :             pres_ext_tens(i, i) = pres_ext_tens(i, i) - pres_ext
     206             :          END DO
     207             :       ELSE
     208         112 :          pres_ext = pvals(1)
     209             :       END IF
     210             : 
     211        2730 :       IF (ANY(pres_ext_tens > 1.0E-5_dp)) THEN
     212           0 :          mtrx = cell%deth*MATMUL(cell%h_inv, MATMUL(pres_ext_tens, TRANSPOSE(cell%h_inv)))
     213             :       END IF
     214             : 
     215         210 :    END SUBROUTINE read_external_press_tensor
     216             : 
     217             : ! **************************************************************************************************
     218             : !> \brief Computes the derivatives for the cell
     219             : !> \param gradient ...
     220             : !> \param av_ptens ...
     221             : !> \param pres_ext ...
     222             : !> \param cell ...
     223             : !> \param mtrx ...
     224             : !> \param keep_angles ...
     225             : !> \param keep_symmetry ...
     226             : !> \param pres_int ...
     227             : !> \param pres_constr ...
     228             : !> \param constraint_id ...
     229             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     230             : ! **************************************************************************************************
     231        3901 :    SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &
     232             :                         keep_symmetry, pres_int, pres_constr, constraint_id)
     233             : 
     234             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gradient
     235             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: av_ptens
     236             :       REAL(KIND=dp), INTENT(IN)                          :: pres_ext
     237             :       TYPE(cell_type), POINTER                           :: cell
     238             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: mtrx
     239             :       LOGICAL, INTENT(IN), OPTIONAL                      :: keep_angles, keep_symmetry
     240             :       REAL(KIND=dp), INTENT(OUT)                         :: pres_int, pres_constr
     241             :       INTEGER, INTENT(IN), OPTIONAL                      :: constraint_id
     242             : 
     243             :       INTEGER                                            :: i, my_constraint_id
     244             :       LOGICAL                                            :: my_keep_angles, my_keep_symmetry
     245             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: correction, pten_hinv_old, ptens
     246             : 
     247        3901 :       my_keep_angles = .FALSE.
     248        3901 :       IF (PRESENT(keep_angles)) my_keep_angles = keep_angles
     249        3901 :       my_keep_symmetry = .FALSE.
     250        3901 :       IF (PRESENT(keep_symmetry)) my_keep_symmetry = keep_symmetry
     251       27307 :       gradient = 0.0_dp
     252        3901 :       IF (PRESENT(constraint_id)) THEN
     253        3901 :          my_constraint_id = constraint_id
     254             :       ELSE
     255           0 :          my_constraint_id = fix_none
     256             :       END IF
     257             : 
     258       27307 :       gradient = 0.0_dp
     259             : 
     260        3901 :       ptens = av_ptens
     261             : 
     262             :       ! Evaluating the internal pressure
     263        3901 :       pres_int = 0.0_dp
     264       15604 :       DO i = 1, 3
     265       15604 :          pres_int = pres_int + ptens(i, i)
     266             :       END DO
     267        3901 :       pres_int = pres_int/3.0_dp
     268             : 
     269           0 :       SELECT CASE (my_constraint_id)
     270             :       CASE (fix_x)
     271           0 :          pres_constr = ptens(2, 2) + ptens(3, 3)
     272             :       CASE (fix_y)
     273           0 :          pres_constr = ptens(1, 1) + ptens(3, 3)
     274             :       CASE (fix_z)
     275           5 :          pres_constr = ptens(1, 1) + ptens(2, 2)
     276             :       CASE (fix_xy)
     277          10 :          pres_constr = ptens(3, 3)
     278             :       CASE (fix_xz)
     279           0 :          pres_constr = ptens(2, 2)
     280             :       CASE (fix_yz)
     281           0 :          pres_constr = ptens(1, 1)
     282             :       CASE (fix_none)
     283        3901 :          pres_constr = ptens(1, 1) + ptens(2, 2) + ptens(3, 3)
     284             :       END SELECT
     285        3901 :       pres_constr = pres_constr/3.0_dp
     286             : 
     287        3901 :       ptens(1, 1) = av_ptens(1, 1) - pres_ext
     288        3901 :       ptens(2, 2) = av_ptens(2, 2) - pres_ext
     289        3901 :       ptens(3, 3) = av_ptens(3, 3) - pres_ext
     290             : 
     291      202852 :       pten_hinv_old = cell%deth*MATMUL(cell%h_inv, ptens)
     292      156040 :       correction = MATMUL(mtrx, cell%hmat)
     293             : 
     294        3901 :       gradient(1) = pten_hinv_old(1, 1) - correction(1, 1)
     295        3901 :       gradient(2) = pten_hinv_old(2, 1) - correction(2, 1)
     296        3901 :       gradient(3) = pten_hinv_old(2, 2) - correction(2, 2)
     297        3901 :       gradient(4) = pten_hinv_old(3, 1) - correction(3, 1)
     298        3901 :       gradient(5) = pten_hinv_old(3, 2) - correction(3, 2)
     299        3901 :       gradient(6) = pten_hinv_old(3, 3) - correction(3, 3)
     300             : 
     301        3901 :       CALL apply_cell_constraints(gradient, cell, my_keep_angles, my_keep_symmetry, my_constraint_id)
     302             : 
     303       27307 :       gradient = -gradient
     304             : 
     305        3901 :    END SUBROUTINE get_dg_dh
     306             : 
     307             : ! **************************************************************************************************
     308             : !> \brief Apply cell constraints
     309             : !> \param gradient ...
     310             : !> \param cell ...
     311             : !> \param keep_angles ...
     312             : !> \param keep_symmetry ...
     313             : !> \param constraint_id ...
     314             : !> \author Matthias Krack (October 26, 2017, MK)
     315             : ! **************************************************************************************************
     316        3901 :    SUBROUTINE apply_cell_constraints(gradient, cell, keep_angles, keep_symmetry, constraint_id)
     317             : 
     318             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gradient
     319             :       TYPE(cell_type), POINTER                           :: cell
     320             :       LOGICAL, INTENT(IN)                                :: keep_angles, keep_symmetry
     321             :       INTEGER, INTENT(IN)                                :: constraint_id
     322             : 
     323             :       REAL(KIND=dp)                                      :: a, a_length, ab_length, b_length, cosa, &
     324             :                                                             cosah, cosg, deriv_gamma, g, gamma, &
     325             :                                                             norm, norm_b, norm_c, sina, sinah, sing
     326             : 
     327        3901 :       IF (keep_angles) THEN
     328             :          ! If we want to keep the angles constant we have to project out the
     329             :          ! components of the cell angles
     330         352 :          norm_b = DOT_PRODUCT(cell%hmat(:, 2), cell%hmat(:, 2))
     331         264 :          norm = DOT_PRODUCT(cell%hmat(1:2, 2), gradient(2:3))
     332         440 :          gradient(2:3) = cell%hmat(1:2, 2)/norm_b*norm
     333         352 :          norm_c = DOT_PRODUCT(cell%hmat(:, 3), cell%hmat(:, 3))
     334         352 :          norm = DOT_PRODUCT(cell%hmat(1:3, 3), gradient(4:6))
     335         616 :          gradient(4:6) = cell%hmat(1:3, 3)/norm_c*norm
     336             :          ! Retain an exact orthorhombic cell
     337             :          ! (off-diagonal elements must remain zero identically to keep QS fast)
     338          88 :          IF (cell%orthorhombic) THEN
     339          70 :             gradient(2) = 0.0_dp
     340          70 :             gradient(4) = 0.0_dp
     341          70 :             gradient(5) = 0.0_dp
     342             :          END IF
     343             :       END IF
     344             : 
     345        3901 :       IF (keep_symmetry) THEN
     346         568 :          SELECT CASE (cell%symmetry_id)
     347             :          CASE (cell_sym_cubic, &
     348             :                cell_sym_tetragonal_ab, &
     349             :                cell_sym_tetragonal_ac, &
     350             :                cell_sym_tetragonal_bc, &
     351             :                cell_sym_orthorhombic)
     352         100 :             SELECT CASE (cell%symmetry_id)
     353             :             CASE (cell_sym_cubic)
     354         100 :                g = (gradient(1) + gradient(3) + gradient(6))/3.0_dp
     355         100 :                gradient(1) = g
     356         100 :                gradient(3) = g
     357         465 :                gradient(6) = g
     358             :             CASE (cell_sym_tetragonal_ab, &
     359             :                   cell_sym_tetragonal_ac, &
     360             :                   cell_sym_tetragonal_bc)
     361         186 :                SELECT CASE (cell%symmetry_id)
     362             :                CASE (cell_sym_tetragonal_ab)
     363         186 :                   g = 0.5_dp*(gradient(1) + gradient(3))
     364         186 :                   gradient(1) = g
     365         186 :                   gradient(3) = g
     366             :                CASE (cell_sym_tetragonal_ac)
     367          93 :                   g = 0.5_dp*(gradient(1) + gradient(6))
     368          93 :                   gradient(1) = g
     369          93 :                   gradient(6) = g
     370             :                CASE (cell_sym_tetragonal_bc)
     371          86 :                   g = 0.5_dp*(gradient(3) + gradient(6))
     372          86 :                   gradient(3) = g
     373          86 :                   gradient(6) = g
     374             :                END SELECT
     375             :             CASE (cell_sym_orthorhombic)
     376             :                ! Nothing else to do
     377             :             END SELECT
     378         568 :             gradient(2) = 0.0_dp
     379         568 :             gradient(4) = 0.0_dp
     380         568 :             gradient(5) = 0.0_dp
     381             :          CASE (cell_sym_hexagonal_gamma_60)
     382         236 :             g = 0.5_dp*(gradient(1) + 0.5_dp*(gradient(2) + sqrt3*gradient(3)))
     383         236 :             gradient(1) = g
     384         236 :             gradient(2) = 0.5_dp*g
     385         236 :             gradient(3) = sqrt3*gradient(2)
     386         236 :             gradient(4) = 0.0_dp
     387         236 :             gradient(5) = 0.0_dp
     388             :          CASE (cell_sym_hexagonal_gamma_120)
     389         118 :             g = 0.5_dp*(gradient(1) - 0.5_dp*(gradient(2) - sqrt3*gradient(3)))
     390         118 :             gradient(1) = g
     391         118 :             gradient(2) = -0.5_dp*g
     392         118 :             gradient(3) = -sqrt3*gradient(2)
     393         118 :             gradient(4) = 0.0_dp
     394         118 :             gradient(5) = 0.0_dp
     395             :          CASE (cell_sym_rhombohedral)
     396             :             a = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
     397             :                  angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
     398          93 :                  angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
     399          93 :             cosa = COS(a)
     400          93 :             sina = SIN(a)
     401          93 :             cosah = COS(0.5_dp*a)
     402          93 :             sinah = SIN(0.5_dp*a)
     403          93 :             norm = cosa/cosah
     404          93 :             norm_c = SQRT(1.0_dp - norm*norm)
     405             :             g = (gradient(1) + gradient(2)*cosa + gradient(3)*sina + &
     406          93 :                  gradient(4)*cosah*norm + gradient(5)*sinah*norm + gradient(6)*norm_c)/3.0_dp
     407          93 :             gradient(1) = g
     408          93 :             gradient(2) = g*cosa
     409          93 :             gradient(3) = g*sina
     410          93 :             gradient(4) = g*cosah*norm
     411          93 :             gradient(5) = g*sinah*norm
     412          93 :             gradient(6) = g*norm_c
     413             :          CASE (cell_sym_monoclinic)
     414         126 :             gradient(2) = 0.0_dp
     415         126 :             gradient(5) = 0.0_dp
     416             :          CASE (cell_sym_monoclinic_gamma_ab)
     417             :             ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
     418             :             a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
     419             :                             cell%hmat(2, 1)*cell%hmat(2, 1) + &
     420         106 :                             cell%hmat(3, 1)*cell%hmat(3, 1))
     421             :             b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
     422             :                             cell%hmat(2, 2)*cell%hmat(2, 2) + &
     423         106 :                             cell%hmat(3, 2)*cell%hmat(3, 2))
     424         106 :             ab_length = 0.5_dp*(a_length + b_length)
     425         106 :             gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
     426         106 :             cosg = COS(gamma)
     427         106 :             sing = SIN(gamma)
     428             :             ! Here, g is the average derivative of the cell vector length ab_length, and deriv_gamma is the derivative of the angle gamma
     429         106 :             g = 0.5_dp*(gradient(1) + cosg*gradient(2) + sing*gradient(3))
     430         106 :             deriv_gamma = (gradient(3)*cosg - gradient(2)*sing)/b_length
     431         106 :             gradient(1) = g
     432         106 :             gradient(2) = g*cosg - ab_length*sing*deriv_gamma
     433         106 :             gradient(3) = g*sing + ab_length*cosg*deriv_gamma
     434         106 :             gradient(4) = 0.0_dp
     435        1471 :             gradient(5) = 0.0_dp
     436             :          CASE (cell_sym_triclinic)
     437             :             ! Nothing to do
     438             :          END SELECT
     439             :       END IF
     440             : 
     441        3901 :       SELECT CASE (constraint_id)
     442             :       CASE (fix_x)
     443           0 :          gradient(1:2) = 0.0_dp
     444           0 :          gradient(4) = 0.0_dp
     445             :       CASE (fix_y)
     446           0 :          gradient(2:3) = 0.0_dp
     447           0 :          gradient(5) = 0.0_dp
     448             :       CASE (fix_z)
     449          20 :          gradient(4:6) = 0.0_dp
     450             :       CASE (fix_xy)
     451          60 :          gradient(1:5) = 0.0_dp
     452             :       CASE (fix_xz)
     453           0 :          gradient(1:2) = 0.0_dp
     454           0 :          gradient(4:6) = 0.0_dp
     455             :       CASE (fix_yz)
     456        3901 :          gradient(2:6) = 0.0_dp
     457             :       CASE (fix_none)
     458             :          ! Nothing to do
     459             :       END SELECT
     460             : 
     461        3901 :    END SUBROUTINE apply_cell_constraints
     462             : 
     463             : END MODULE cell_opt_utils

Generated by: LCOV version 1.15