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