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 utilities for a Dimer Method calculations
      10             : !> \par History
      11             : !>      none
      12             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      13             : ! **************************************************************************************************
      14             : MODULE dimer_utils
      15             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      16             :    USE dimer_types,                     ONLY: dimer_env_type
      17             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      18             :                                               section_vals_remove_values,&
      19             :                                               section_vals_type,&
      20             :                                               section_vals_val_set
      21             :    USE kinds,                           ONLY: dp
      22             :    USE memory_utilities,                ONLY: reallocate
      23             : #include "../base/base_uses.f90"
      24             : 
      25             :    IMPLICIT NONE
      26             :    PRIVATE
      27             : 
      28             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      29             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dimer_utils'
      30             : 
      31             :    PUBLIC :: rotate_dimer, update_dimer_vec, get_theta
      32             :    REAL(KIND=dp), PARAMETER, PUBLIC     :: dimer_thrs = EPSILON(0.0_dp)*1.0E4_dp
      33             : 
      34             : CONTAINS
      35             : 
      36             : ! **************************************************************************************************
      37             : !> \brief Performs a rotation of the unit dimer vector
      38             : !> \param nvec ...
      39             : !> \param theta ...
      40             : !> \param dt ...
      41             : !> \par History
      42             : !>      none
      43             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      44             : ! **************************************************************************************************
      45        1998 :    SUBROUTINE rotate_dimer(nvec, theta, dt)
      46             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: nvec, theta
      47             :       REAL(KIND=dp)                                      :: dt
      48             : 
      49             :       INTEGER                                            :: output_unit
      50             :       LOGICAL                                            :: check
      51             : 
      52        1998 :       output_unit = cp_logger_get_default_io_unit()
      53             : 
      54             :       ! Orthogonality check for the rotation..
      55       36348 :       check = ABS(DOT_PRODUCT(nvec, theta)) < MAX(1.0E-9_dp, dimer_thrs)
      56        1998 :       IF (.NOT. check .AND. (output_unit > 0)) THEN
      57           0 :          WRITE (output_unit, *) "NVEC and THETA should be orthogonal! Residue: ", &
      58           0 :             ABS(DOT_PRODUCT(nvec, theta))
      59             :       END IF
      60        1998 :       CPASSERT(check)
      61       70698 :       nvec = nvec*COS(dt) + theta*SIN(dt)
      62             : 
      63        1998 :    END SUBROUTINE rotate_dimer
      64             : 
      65             : ! **************************************************************************************************
      66             : !> \brief Updates the orientation of the dimer vector in the input file
      67             : !> \param dimer_env ...
      68             : !> \param motion_section ...
      69             : !> \par History
      70             : !>      none
      71             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      72             : ! **************************************************************************************************
      73         858 :    SUBROUTINE update_dimer_vec(dimer_env, motion_section)
      74             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
      75             :       TYPE(section_vals_type), POINTER                   :: motion_section
      76             : 
      77             :       INTEGER                                            :: i, i_rep_val, isize, j, size_array
      78         858 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: array
      79             :       TYPE(section_vals_type), POINTER                   :: nvec_section
      80             : 
      81             :       nvec_section => section_vals_get_subs_vals(motion_section, &
      82        1716 :                                                  "GEO_OPT%TRANSITION_STATE%DIMER%DIMER_VECTOR")
      83             :       ! Clean the content of the section first..
      84         858 :       CALL section_vals_remove_values(nvec_section)
      85             :       ! Fill in the section with the present values..
      86         858 :       size_array = 6
      87         858 :       isize = 0
      88         858 :       i_rep_val = 0
      89        2354 :       Main_Loop: DO i = 1, SIZE(dimer_env%nvec), size_array
      90        2354 :          ALLOCATE (array(size_array))
      91        2354 :          i_rep_val = i_rep_val + 1
      92       15386 :          DO j = 1, size_array
      93       13890 :             isize = isize + 1
      94       13890 :             array(j) = dimer_env%nvec(isize)
      95       15386 :             IF (isize == SIZE(dimer_env%nvec)) THEN
      96         858 :                CALL reallocate(array, 1, j)
      97             :                CALL section_vals_val_set(nvec_section, "_DEFAULT_KEYWORD_", r_vals_ptr=array, &
      98         858 :                                          i_rep_val=i_rep_val)
      99         858 :                EXIT Main_Loop
     100             :             END IF
     101             :          END DO
     102             :          CALL section_vals_val_set(nvec_section, "_DEFAULT_KEYWORD_", r_vals_ptr=array, &
     103        1496 :                                    i_rep_val=i_rep_val)
     104             :       END DO Main_Loop
     105         858 :       CPASSERT(isize == SIZE(dimer_env%nvec))
     106         858 :    END SUBROUTINE update_dimer_vec
     107             : 
     108             : ! **************************************************************************************************
     109             : !> \brief This function orthonormalize the vector for the rotational search
     110             : !> \param gradient ...
     111             : !> \param dimer_env ...
     112             : !> \param norm ...
     113             : !> \par History
     114             : !>      none
     115             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     116             : ! **************************************************************************************************
     117        1428 :    SUBROUTINE get_theta(gradient, dimer_env, norm)
     118             :       REAL(KIND=dp), DIMENSION(:)                        :: gradient
     119             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     120             :       REAL(KIND=dp), INTENT(OUT)                         :: norm
     121             : 
     122       49668 :       gradient = gradient - DOT_PRODUCT(gradient, dimer_env%nvec)*dimer_env%nvec
     123       25548 :       norm = SQRT(DOT_PRODUCT(gradient, gradient))
     124        1428 :       IF (norm < EPSILON(0.0_dp)) THEN
     125             :          ! This means that NVEC is totally aligned with minimum curvature mode
     126           0 :          gradient = 0.0_dp
     127             :       ELSE
     128       25548 :          gradient = gradient/norm
     129             :       END IF
     130        1428 :    END SUBROUTINE get_theta
     131             : 
     132             : END MODULE dimer_utils