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 types used for a Dimer Method calculations
      10             : !> \par History
      11             : !>     -Luca Bellucci 11.2017 - CNR-NANO, Pisa
      12             : !>      add K-DIMER from doi:10.1063/1.4898664
      13             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      14             : ! **************************************************************************************************
      15             : MODULE dimer_methods
      16             :    USE bibliography,                    ONLY: Henkelman2014,&
      17             :                                               cite_reference
      18             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      19             :                                               cp_logger_type
      20             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      21             :                                               cp_print_key_unit_nr
      22             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      23             :                                               cp_subsys_type
      24             :    USE dimer_types,                     ONLY: dimer_env_type,&
      25             :                                               dimer_fixed_atom_control
      26             :    USE dimer_utils,                     ONLY: get_theta
      27             :    USE force_env_methods,               ONLY: force_env_calc_energy_force
      28             :    USE force_env_types,                 ONLY: force_env_get
      29             :    USE geo_opt,                         ONLY: cp_rot_opt
      30             :    USE gopt_f_types,                    ONLY: gopt_f_type
      31             :    USE input_constants,                 ONLY: do_first_rotation_step,&
      32             :                                               do_second_rotation_step,&
      33             :                                               do_third_rotation_step
      34             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      35             :                                               section_vals_type
      36             :    USE kinds,                           ONLY: dp
      37             :    USE motion_utils,                    ONLY: rot_ana,&
      38             :                                               thrs_motion
      39             :    USE particle_list_types,             ONLY: particle_list_type
      40             : #include "../base/base_uses.f90"
      41             : 
      42             :    IMPLICIT NONE
      43             :    PRIVATE
      44             : 
      45             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      46             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dimer_methods'
      47             : 
      48             :    PUBLIC :: cp_eval_at_ts
      49             : 
      50             : CONTAINS
      51             : 
      52             : ! **************************************************************************************************
      53             : !> \brief Computes the dimer energy/gradients (including the rotation of the dimer)
      54             : !> \param gopt_env ...
      55             : !> \param x ...
      56             : !> \param f ...
      57             : !> \param gradient ...
      58             : !> \param calc_force ...
      59             : !> \date  01.2008
      60             : !> \par History
      61             : !>      none
      62             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
      63             : ! **************************************************************************************************
      64        1816 :    RECURSIVE SUBROUTINE cp_eval_at_ts(gopt_env, x, f, gradient, calc_force)
      65             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
      66             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x
      67             :       REAL(KIND=dp), INTENT(OUT)                         :: f
      68             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gradient
      69             :       LOGICAL, INTENT(IN)                                :: calc_force
      70             : 
      71             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_eval_at_ts'
      72             : 
      73             :       INTEGER                                            :: handle, iw
      74             :       LOGICAL                                            :: eval_analytical
      75             :       REAL(KIND=dp)                                      :: angle1, angle2, f1, gm1, gm2, norm, swf
      76             :       TYPE(cp_logger_type), POINTER                      :: logger
      77             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
      78             :       TYPE(section_vals_type), POINTER                   :: print_section
      79             : 
      80        1816 :       NULLIFY (dimer_env, logger, print_section)
      81        1816 :       CALL timeset(routineN, handle)
      82        1816 :       logger => cp_get_default_logger()
      83             : 
      84        1816 :       CPASSERT(ASSOCIATED(gopt_env))
      85        1816 :       dimer_env => gopt_env%dimer_env
      86        1816 :       CPASSERT(ASSOCIATED(dimer_env))
      87        1816 :       iw = cp_print_key_unit_nr(logger, gopt_env%geo_section, "PRINT%PROGRAM_RUN_INFO", extension=".log")
      88             :       ! Possibly rotate Dimer or just compute Gradients of point 0 for Translation
      89        1816 :       IF (gopt_env%dimer_rotation) THEN
      90             :          IF (debug_this_module .AND. (iw > 0)) THEN
      91             :             WRITE (iw, '(A)') "NVEC:"
      92             :             WRITE (iw, '(3F15.9)') dimer_env%nvec
      93             :          END IF
      94        2154 :          SELECT CASE (dimer_env%rot%rotation_step)
      95             :          CASE (do_first_rotation_step, do_third_rotation_step)
      96         726 :             eval_analytical = .TRUE.
      97         726 :             IF ((dimer_env%rot%rotation_step == do_third_rotation_step) .AND. (dimer_env%rot%interpolate_gradient)) THEN
      98             :                eval_analytical = .FALSE.
      99             :             END IF
     100             :             IF (eval_analytical) THEN
     101             :                ! Compute energy, gradients and rotation vector for R1
     102         132 :                CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%rot%g1)
     103             :             ELSE
     104         594 :                angle1 = dimer_env%rot%angle1
     105         594 :                angle2 = dimer_env%rot%angle2
     106             :                dimer_env%rot%g1 = SIN(angle1 - angle2)/SIN(angle1)*dimer_env%rot%g1 + &
     107             :                                   SIN(angle2)/SIN(angle1)*dimer_env%rot%g1p + &
     108       10932 :                                   (1.0_dp - COS(angle2) - SIN(angle2)*TAN(angle1/2.0_dp))*dimer_env%rot%g0
     109             :             END IF
     110             : 
     111             :             ! Determine the theta vector (i.e. the search direction for line minimization)
     112       24954 :             gradient = -2.0_dp*(dimer_env%rot%g1 - dimer_env%rot%g0)
     113             :             IF (debug_this_module .AND. (iw > 0)) THEN
     114             :                WRITE (iw, '(A)') "G1 vector:"
     115             :                WRITE (iw, '(3F15.9)') dimer_env%rot%g1
     116             :                WRITE (iw, '(A)') "-2*(G1-G0) vector:"
     117             :                WRITE (iw, '(3F15.9)') gradient
     118             :             END IF
     119         726 :             CALL get_theta(gradient, dimer_env, norm)
     120         726 :             f = norm
     121         726 :             dimer_env%cg_rot%norm_theta_old = dimer_env%cg_rot%norm_theta
     122         726 :             dimer_env%cg_rot%norm_theta = norm
     123             : 
     124             :             IF (debug_this_module .AND. (iw > 0)) THEN
     125             :                WRITE (iw, '(A,F20.10)') "Rotational Force step (1,3): module:", norm
     126             :                WRITE (iw, '(3F15.9)') gradient
     127             :             END IF
     128             : 
     129             :             ! Compute curvature and derivative of the curvature w.r.t. the rotational angle
     130       12840 :             dimer_env%rot%curvature = DOT_PRODUCT(dimer_env%rot%g1 - dimer_env%rot%g0, dimer_env%nvec)/dimer_env%dr
     131       12840 :             dimer_env%rot%dCdp = 2.0_dp*DOT_PRODUCT(dimer_env%rot%g1 - dimer_env%rot%g0, gradient)/dimer_env%dr
     132             : 
     133         726 :             dimer_env%rot%rotation_step = do_second_rotation_step
     134       12840 :             gradient = -gradient
     135             :          CASE (do_second_rotation_step)
     136             :             ! Compute energy, gradients and rotation vector for R1
     137         702 :             CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%rot%g1p)
     138       12708 :             dimer_env%rot%curvature = DOT_PRODUCT(dimer_env%rot%g1p - dimer_env%rot%g0, dimer_env%nvec)/dimer_env%dr
     139         702 :             dimer_env%rot%rotation_step = do_third_rotation_step
     140             : 
     141             :             ! Determine the theta vector (i.e. the search direction for line minimization)
     142             :             ! This is never used for getting a new theta but is consistent in order to
     143             :             ! give back the right value of f
     144       24714 :             gradient = -2.0_dp*(dimer_env%rot%g1p - dimer_env%rot%g0)
     145         702 :             CALL get_theta(gradient, dimer_env, norm)
     146         702 :             f = norm
     147             : 
     148        2856 :             IF (debug_this_module .AND. (iw > 0)) THEN
     149             :                WRITE (iw, '(A)') "Rotational Force step (1,3):"
     150             :                WRITE (iw, '(3F15.9)') gradient
     151             :             END IF
     152             :          END SELECT
     153             :       ELSE
     154             :          ! Compute energy, gradients and rotation vector for R0
     155         388 :          CALL cp_eval_at_ts_low(gopt_env, x, 0, dimer_env, calc_force, f, dimer_env%rot%g0)
     156             : 
     157             :          ! The dimer is rotated only when we are out of the translation line search
     158         388 :          IF (.NOT. gopt_env%do_line_search) THEN
     159             :             IF (debug_this_module .AND. (iw > 0)) THEN
     160             :                WRITE (iw, '(A)') "Entering the rotation module"
     161             :                WRITE (iw, '(A)') "G0 vector:"
     162             :                WRITE (iw, '(3F15.9)') dimer_env%rot%g0
     163             :             END IF
     164             :             CALL cp_rot_opt(gopt_env%gopt_dimer_env, x, gopt_env%gopt_dimer_param, &
     165         132 :                             gopt_env%gopt_dimer_env%geo_section)
     166         132 :             dimer_env%rot%rotation_step = do_first_rotation_step
     167             :          END IF
     168             : 
     169         388 :          print_section => section_vals_get_subs_vals(gopt_env%gopt_dimer_env%geo_section, "PRINT")
     170             : 
     171             :          ! Correcting gradients for Translation K-DIMER or STANDARD
     172         388 :          IF (dimer_env%kdimer) THEN
     173           6 :             CALL cite_reference(Henkelman2014)
     174             :             ! K-DIMER
     175           6 :             IF (iw > 0) THEN
     176           3 :                WRITE (iw, '(T2,A)') "DIMER| Correcting gradients for Translation with K-DIMER method"
     177             :             END IF
     178           6 :             swf = 1.0_dp + EXP(dimer_env%beta*dimer_env%rot%curvature)
     179           6 :             gm2 = 1.0_dp - (1.0_dp/swf)
     180           6 :             gm1 = (2.0_dp/swf) - 1.0_dp
     181             :             gradient = gm2*(dimer_env%rot%g0 - 2.0_dp*DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec) &
     182         168 :                        - gm1*(DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec)
     183           6 :             CALL remove_rot_transl_component(gopt_env, gradient, print_section)
     184             :             IF (debug_this_module .AND. (iw > 0)) WRITE (iw, *) "K-DIMER", dimer_env%beta, dimer_env%rot%curvature, &
     185             :                dimer_env%rot%dCdp, gm1, gm2, swf
     186             :          ELSE
     187         382 :             IF (iw > 0) THEN
     188         191 :                WRITE (iw, '(T2,A)') "DIMER| Correcting gradients for Translation with standard method"
     189             :             END IF
     190         382 :             IF (dimer_env%rot%curvature > 0) THEN
     191        4518 :                gradient = -DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec
     192          90 :                CALL remove_rot_transl_component(gopt_env, gradient, print_section)
     193             :             ELSE
     194       11578 :                gradient = dimer_env%rot%g0 - 2.0_dp*DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec
     195         292 :                CALL remove_rot_transl_component(gopt_env, gradient, print_section)
     196             :             END IF
     197             :          END IF
     198             :          IF (debug_this_module .AND. (iw > 0)) THEN
     199             :             WRITE (iw, *) "final gradient:", gradient
     200             :             WRITE (iw, '(A,F20.10)') "norm gradient:", SQRT(DOT_PRODUCT(gradient, gradient))
     201             :          END IF
     202         388 :          IF (.NOT. gopt_env%do_line_search) THEN
     203        1908 :             f = SQRT(DOT_PRODUCT(gradient, gradient))
     204             :          ELSE
     205        3772 :             f = -DOT_PRODUCT(gradient, dimer_env%tsl%tls_vec)
     206             :          END IF
     207             :       END IF
     208        1816 :       CALL cp_print_key_finished_output(iw, logger, gopt_env%geo_section, "PRINT%PROGRAM_RUN_INFO")
     209        1816 :       CALL timestop(handle)
     210        1816 :    END SUBROUTINE cp_eval_at_ts
     211             : 
     212             : ! **************************************************************************************************
     213             : !> \brief This function removes translational forces after project of the gradient
     214             : !> \param gopt_env ...
     215             : !> \param gradient ...
     216             : !> \param print_section ...
     217             : !> \par History
     218             : !>      2016/03/02 [LTong] added fixed atom constraint for gradient
     219             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     220             : ! **************************************************************************************************
     221         388 :    SUBROUTINE remove_rot_transl_component(gopt_env, gradient, print_section)
     222             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     223             :       REAL(KIND=dp), DIMENSION(:)                        :: gradient
     224             :       TYPE(section_vals_type), POINTER                   :: print_section
     225             : 
     226             :       CHARACTER(len=*), PARAMETER :: routineN = 'remove_rot_transl_component'
     227             : 
     228             :       INTEGER                                            :: dof, handle, i, j, natoms
     229             :       REAL(KIND=dp)                                      :: norm, norm_gradient_old
     230         388 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: D, mat
     231             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     232             :       TYPE(particle_list_type), POINTER                  :: particles
     233             : 
     234         388 :       CALL timeset(routineN, handle)
     235         388 :       NULLIFY (mat)
     236         388 :       CALL force_env_get(gopt_env%force_env, subsys=subsys)
     237         388 :       CALL cp_subsys_get(subsys, particles=particles)
     238             : 
     239         388 :       natoms = particles%n_els
     240        5680 :       norm_gradient_old = SQRT(DOT_PRODUCT(gradient, gradient))
     241         388 :       IF (norm_gradient_old > 0.0_dp) THEN
     242         388 :          IF (natoms > 1) THEN
     243             :             CALL rot_ana(particles%els, mat, dof, print_section, keep_rotations=.FALSE., &
     244         350 :                          mass_weighted=.FALSE., natoms=natoms)
     245             : 
     246             :             ! Orthogonalize gradient with respect to the full set of Roto-Trasl vectors
     247        1400 :             ALLOCATE (D(3*natoms, dof))
     248             :             ! Check First orthogonality in the first element of the basis set
     249        2406 :             DO i = 1, dof
     250       63664 :                D(:, i) = mat(:, i)
     251        7436 :                DO j = i + 1, dof
     252       81380 :                   norm = DOT_PRODUCT(mat(:, i), mat(:, j))
     253        7086 :                   CPASSERT(ABS(norm) < thrs_motion)
     254             :                END DO
     255             :             END DO
     256        2406 :             DO i = 1, dof
     257       32860 :                norm = DOT_PRODUCT(gradient, D(:, i))
     258       33210 :                gradient = gradient - norm*D(:, i)
     259             :             END DO
     260         350 :             DEALLOCATE (D)
     261         350 :             DEALLOCATE (mat)
     262             :          END IF
     263             :       END IF
     264             :       ! apply constraint
     265         388 :       CALL dimer_fixed_atom_control(gradient, subsys)
     266         388 :       CALL timestop(handle)
     267         388 :    END SUBROUTINE remove_rot_transl_component
     268             : 
     269             : ! **************************************************************************************************
     270             : !> \brief ...
     271             : !> \param gopt_env ...
     272             : !> \param x ...
     273             : !> \param dimer_index ...
     274             : !> \param dimer_env ...
     275             : !> \param calc_force ...
     276             : !> \param f ...
     277             : !> \param gradient ...
     278             : !> \par History
     279             : !>      none
     280             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     281             : ! **************************************************************************************************
     282        2444 :    SUBROUTINE cp_eval_at_ts_low(gopt_env, x, dimer_index, dimer_env, calc_force, &
     283        1222 :                                 f, gradient)
     284             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     285             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x
     286             :       INTEGER, INTENT(IN)                                :: dimer_index
     287             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     288             :       LOGICAL, INTENT(IN)                                :: calc_force
     289             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: f
     290             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: gradient
     291             : 
     292             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'cp_eval_at_ts_low'
     293             : 
     294             :       INTEGER                                            :: handle, idg, idir, ip
     295             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     296             :       TYPE(particle_list_type), POINTER                  :: particles
     297             : 
     298        1222 :       CALL timeset(routineN, handle)
     299        1222 :       idg = 0
     300        1222 :       CALL force_env_get(gopt_env%force_env, subsys=subsys)
     301        1222 :       CALL cp_subsys_get(subsys, particles=particles)
     302        7580 :       DO ip = 1, particles%n_els
     303       26654 :          DO idir = 1, 3
     304       19074 :             idg = idg + 1
     305       25432 :             particles%els(ip)%r(idir) = x(idg) + REAL(dimer_index, KIND=dp)*dimer_env%nvec(idg)*dimer_env%dr
     306             :          END DO
     307             :       END DO
     308             : 
     309             :       ! Compute energy and forces
     310        1222 :       CALL force_env_calc_energy_force(gopt_env%force_env, calc_force=calc_force)
     311             : 
     312             :       ! Possibly take the potential energy
     313        1222 :       IF (PRESENT(f)) THEN
     314        1222 :          CALL force_env_get(gopt_env%force_env, potential_energy=f)
     315             :       END IF
     316             : 
     317             :       ! Possibly take the gradients
     318        1222 :       IF (PRESENT(gradient)) THEN
     319        1222 :          idg = 0
     320        1222 :          CALL cp_subsys_get(subsys, particles=particles)
     321        7580 :          DO ip = 1, particles%n_els
     322       26654 :             DO idir = 1, 3
     323       19074 :                idg = idg + 1
     324       19074 :                CPASSERT(SIZE(gradient) >= idg)
     325       25432 :                gradient(idg) = -particles%els(ip)%f(idir)
     326             :             END DO
     327             :          END DO
     328             :       END IF
     329        1222 :       CALL timestop(handle)
     330        1222 :    END SUBROUTINE cp_eval_at_ts_low
     331             : 
     332             : END MODULE dimer_methods