LCOV - code coverage report
Current view: top level - src/motion - dimer_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 95 106 89.6 %
Date: 2024-04-20 06:29:22 Functions: 4 8 50.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 types used for a Dimer Method calculations
      10             : !> \par History
      11             : !>      Luca Bellucci 11.2017 added kdimer and beta
      12             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      13             : ! **************************************************************************************************
      14             : MODULE dimer_types
      15             : 
      16             :    USE cell_types,                      ONLY: use_perd_x,&
      17             :                                               use_perd_xy,&
      18             :                                               use_perd_xyz,&
      19             :                                               use_perd_xz,&
      20             :                                               use_perd_y,&
      21             :                                               use_perd_yz,&
      22             :                                               use_perd_z
      23             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      24             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      25             :                                               cp_subsys_type
      26             :    USE global_types,                    ONLY: global_environment_type
      27             :    USE input_constants,                 ONLY: do_first_rotation_step
      28             :    USE input_section_types,             ONLY: section_vals_get,&
      29             :                                               section_vals_get_subs_vals,&
      30             :                                               section_vals_type,&
      31             :                                               section_vals_val_get
      32             :    USE kinds,                           ONLY: dp
      33             :    USE molecule_kind_list_types,        ONLY: molecule_kind_list_type
      34             :    USE molecule_kind_types,             ONLY: fixd_constraint_type,&
      35             :                                               get_molecule_kind,&
      36             :                                               molecule_kind_type
      37             : #include "../base/base_uses.f90"
      38             : 
      39             :    IMPLICIT NONE
      40             :    PRIVATE
      41             : 
      42             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      43             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dimer_types'
      44             : 
      45             :    PUBLIC :: dimer_env_type, &
      46             :              dimer_env_create, &
      47             :              dimer_env_retain, &
      48             :              dimer_env_release, &
      49             :              dimer_fixed_atom_control
      50             : 
      51             : ! **************************************************************************************************
      52             : !> \brief Type containing all informations abour the rotation of the Dimer
      53             : !> \par History
      54             : !>      none
      55             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      56             : ! **************************************************************************************************
      57             :    TYPE dimer_rotational_type
      58             :       ! Rotational parameters
      59             :       INTEGER                                    :: rotation_step = 0
      60             :       LOGICAL                                    :: interpolate_gradient = .FALSE.
      61             :       REAL(KIND=dp)                              :: angle_tol = 0.0_dp, angle1 = 0.0_dp, angle2 = 0.0_dp, &
      62             :                                                     dCdp = 0.0_dp, curvature = 0.0_dp
      63             :       REAL(KIND=dp), POINTER, DIMENSION(:)       :: g0 => NULL(), g1 => NULL(), g1p => NULL()
      64             :    END TYPE dimer_rotational_type
      65             : 
      66             : ! **************************************************************************************************
      67             : !> \brief Type containing all informations abour the translation of the Dimer
      68             : !> \par History
      69             : !>      none
      70             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      71             : ! **************************************************************************************************
      72             :    TYPE dimer_translational_type
      73             :       ! Translational parameters
      74             :       REAL(KIND=dp), POINTER, DIMENSION(:)       :: tls_vec => NULL()
      75             :    END TYPE dimer_translational_type
      76             : 
      77             : ! **************************************************************************************************
      78             : !> \brief Conjugate Directions type
      79             : !> \par History
      80             : !>      none
      81             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      82             : ! **************************************************************************************************
      83             :    TYPE dimer_cg_rot_type
      84             :       REAL(KIND=dp)                              :: norm_theta = 0.0_dp, norm_theta_old = 0.0_dp, norm_h = 0.0_dp
      85             :       REAL(KIND=dp), POINTER, DIMENSION(:)       :: nvec_old => NULL()
      86             :    END TYPE dimer_cg_rot_type
      87             : 
      88             : ! **************************************************************************************************
      89             : !> \brief Defines the environment for a Dimer Method calculation
      90             : !> \par History
      91             : !>      none
      92             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
      93             : ! **************************************************************************************************
      94             :    TYPE dimer_env_type
      95             :       INTEGER                                    :: ref_count = 0
      96             :       REAL(KIND=dp)                              :: dr = 0.0_dp
      97             :       REAL(KIND=dp), POINTER, DIMENSION(:)       :: nvec => NULL()
      98             :       REAL(KIND=dp)                              :: beta = 0.0_dp
      99             :       TYPE(dimer_rotational_type)                :: rot = dimer_rotational_type()
     100             :       TYPE(dimer_translational_type)             :: tsl = dimer_translational_type()
     101             :       TYPE(dimer_cg_rot_type)                    :: cg_rot = dimer_cg_rot_type()
     102             :       LOGICAL                                    :: kdimer = .FALSE.
     103             :    END TYPE dimer_env_type
     104             : 
     105             : CONTAINS
     106             : 
     107             : ! **************************************************************************************************
     108             : !> \brief ...
     109             : !> \param dimer_env ...
     110             : !> \param subsys ...
     111             : !> \param globenv ...
     112             : !> \param dimer_section ...
     113             : !> \par History
     114             : !>      Luca Bellucci 11.2017 added K-DIMER and BETA
     115             : !>      2016/03/03 [LTong] changed input natom to subsys
     116             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     117             : ! **************************************************************************************************
     118          22 :    SUBROUTINE dimer_env_create(dimer_env, subsys, globenv, dimer_section)
     119             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     120             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     121             :       TYPE(global_environment_type), POINTER             :: globenv
     122             :       TYPE(section_vals_type), POINTER                   :: dimer_section
     123             : 
     124             :       INTEGER                                            :: i, isize, j, k, n_rep_val, natom, unit_nr
     125             :       LOGICAL                                            :: explicit
     126             :       REAL(KIND=dp)                                      :: norm, xval(3)
     127          22 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: array
     128             :       TYPE(section_vals_type), POINTER                   :: nvec_section
     129             : 
     130          44 :       unit_nr = cp_logger_get_default_io_unit()
     131          22 :       CPASSERT(.NOT. ASSOCIATED(dimer_env))
     132          22 :       ALLOCATE (dimer_env)
     133          22 :       dimer_env%ref_count = 1
     134             :       ! Setup NVEC
     135             :       ! get natom
     136          22 :       CALL cp_subsys_get(subsys=subsys, natom=natom)
     137             :       ! Allocate the working arrays
     138          66 :       ALLOCATE (dimer_env%nvec(natom*3))
     139          66 :       ALLOCATE (dimer_env%rot%g0(natom*3))
     140          66 :       ALLOCATE (dimer_env%rot%g1(natom*3))
     141          66 :       ALLOCATE (dimer_env%rot%g1p(natom*3))
     142             :       ! Check if the dimer vector is available in the input or not..
     143          22 :       nvec_section => section_vals_get_subs_vals(dimer_section, "DIMER_VECTOR")
     144          22 :       CALL section_vals_get(nvec_section, explicit=explicit)
     145          22 :       IF (explicit) THEN
     146           4 :          IF (unit_nr > 0) WRITE (unit_nr, *) "Reading Dimer Vector from file!"
     147           4 :          NULLIFY (array)
     148           4 :          CALL section_vals_val_get(nvec_section, "_DEFAULT_KEYWORD_", n_rep_val=n_rep_val)
     149           4 :          isize = 0
     150          10 :          DO i = 1, n_rep_val
     151           6 :             CALL section_vals_val_get(nvec_section, "_DEFAULT_KEYWORD_", r_vals=array, i_rep_val=i)
     152          34 :             DO j = 1, SIZE(array)
     153          24 :                isize = isize + 1
     154          30 :                dimer_env%nvec(isize) = array(j)
     155             :             END DO
     156             :          END DO
     157           4 :          CPASSERT(isize == SIZE(dimer_env%nvec))
     158             :       ELSE
     159          18 :          CALL globenv%gaussian_rng_stream%fill(dimer_env%nvec)
     160             :       END IF
     161             :       ! Check for translation in the dimer vector and remove them
     162          22 :       IF (natom > 1) THEN
     163          20 :          xval = 0.0_dp
     164          90 :          DO j = 1, natom
     165         300 :             DO k = 1, 3
     166         210 :                i = (j - 1)*3 + k
     167         280 :                xval(k) = xval(k) + dimer_env%nvec(i)
     168             :             END DO
     169             :          END DO
     170             :          ! Subtract net translations
     171          80 :          xval = xval/REAL(natom*3, KIND=dp)
     172          90 :          DO j = 1, natom
     173         300 :             DO k = 1, 3
     174         210 :                i = (j - 1)*3 + k
     175         280 :                dimer_env%nvec(i) = dimer_env%nvec(i) - xval(k)
     176             :             END DO
     177             :          END DO
     178             :       END IF
     179             :       ! set nvec components to zero for the corresponding constraints
     180          22 :       CALL dimer_fixed_atom_control(dimer_env%nvec, subsys)
     181         238 :       norm = SQRT(SUM(dimer_env%nvec**2))
     182          22 :       IF (norm <= EPSILON(0.0_dp)) &
     183           0 :          CPABORT("The norm of the dimer vector is 0! Calculation cannot proceed further.")
     184         238 :       dimer_env%nvec = dimer_env%nvec/norm
     185          22 :       dimer_env%rot%rotation_step = do_first_rotation_step
     186          22 :       CALL section_vals_val_get(dimer_section, "DR", r_val=dimer_env%dr)
     187             :       CALL section_vals_val_get(dimer_section, "INTERPOLATE_GRADIENT", &
     188          22 :                                 l_val=dimer_env%rot%interpolate_gradient)
     189             :       CALL section_vals_val_get(dimer_section, "ANGLE_TOLERANCE", &
     190          22 :                                 r_val=dimer_env%rot%angle_tol)
     191             :       CALL section_vals_val_get(dimer_section, "K-DIMER", &
     192          22 :                                 l_val=dimer_env%kdimer)
     193             :       CALL section_vals_val_get(dimer_section, "BETA", &
     194          22 :                                 r_val=dimer_env%beta)
     195             :       ! initialise values
     196          22 :       dimer_env%cg_rot%norm_h = 1.0_dp
     197         238 :       dimer_env%rot%g0 = 0.0_dp
     198         238 :       dimer_env%rot%g1 = 0.0_dp
     199         238 :       dimer_env%rot%g1p = 0.0_dp
     200          66 :       ALLOCATE (dimer_env%cg_rot%nvec_old(natom*3))
     201          66 :    END SUBROUTINE dimer_env_create
     202             : 
     203             : ! **************************************************************************************************
     204             : !> \brief ...
     205             : !> \param dimer_env ...
     206             : !> \par History
     207             : !>      none
     208             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     209             : ! **************************************************************************************************
     210          22 :    SUBROUTINE dimer_env_retain(dimer_env)
     211             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     212             : 
     213          22 :       CPASSERT(ASSOCIATED(dimer_env))
     214          22 :       CPASSERT(dimer_env%ref_count > 0)
     215          22 :       dimer_env%ref_count = dimer_env%ref_count + 1
     216          22 :    END SUBROUTINE dimer_env_retain
     217             : 
     218             : ! **************************************************************************************************
     219             : !> \brief ...
     220             : !> \param dimer_env ...
     221             : !> \par History
     222             : !>      none
     223             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
     224             : ! **************************************************************************************************
     225        2012 :    SUBROUTINE dimer_env_release(dimer_env)
     226             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     227             : 
     228        2012 :       IF (ASSOCIATED(dimer_env)) THEN
     229          44 :          CPASSERT(dimer_env%ref_count > 0)
     230          44 :          dimer_env%ref_count = dimer_env%ref_count - 1
     231          44 :          IF (dimer_env%ref_count == 0) THEN
     232          22 :             IF (ASSOCIATED(dimer_env%nvec)) THEN
     233          22 :                DEALLOCATE (dimer_env%nvec)
     234             :             END IF
     235          22 :             IF (ASSOCIATED(dimer_env%rot%g0)) THEN
     236          22 :                DEALLOCATE (dimer_env%rot%g0)
     237             :             END IF
     238          22 :             IF (ASSOCIATED(dimer_env%rot%g1)) THEN
     239          22 :                DEALLOCATE (dimer_env%rot%g1)
     240             :             END IF
     241          22 :             IF (ASSOCIATED(dimer_env%rot%g1p)) THEN
     242          22 :                DEALLOCATE (dimer_env%rot%g1p)
     243             :             END IF
     244          22 :             IF (ASSOCIATED(dimer_env%cg_rot%nvec_old)) THEN
     245          22 :                DEALLOCATE (dimer_env%cg_rot%nvec_old)
     246             :             END IF
     247             :             ! No need to deallocate tls_vec (just a pointer to aother local array)
     248          22 :             NULLIFY (dimer_env%tsl%tls_vec)
     249          22 :             DEALLOCATE (dimer_env)
     250             :          END IF
     251             :       END IF
     252        2012 :    END SUBROUTINE dimer_env_release
     253             : 
     254             : ! **************************************************************************************************
     255             : !> \brief Set parts of a given array vec to zero according to fixed atom constraints.
     256             : !>        When atoms are (partially) fixed then the relevant components of
     257             : !>        nvec should be set to zero.  Furthermore, the relevant components
     258             : !>        of the gradient in CG should also be set to zero.
     259             : !> \param vec : vector to be modified
     260             : !> \param subsys : subsys type object used by CP2k
     261             : !> \par History
     262             : !>      2016/03/03 [LTong] created
     263             : !> \author Lianheng Tong [LTong]
     264             : ! **************************************************************************************************
     265         410 :    SUBROUTINE dimer_fixed_atom_control(vec, subsys)
     266             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: vec
     267             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     268             : 
     269             :       INTEGER                                            :: ii, ikind, ind, iparticle, nfixed_atoms, &
     270             :                                                             nkinds
     271         410 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     272             :       TYPE(molecule_kind_list_type), POINTER             :: molecule_kinds
     273         410 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     274             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     275             : 
     276         410 :       NULLIFY (molecule_kinds, molecule_kind, fixd_list)
     277             : 
     278             :       ! need to get constraint information from molecule information
     279             :       CALL cp_subsys_get(subsys=subsys, &
     280         410 :                          molecule_kinds=molecule_kinds)
     281         410 :       molecule_kind_set => molecule_kinds%els
     282             : 
     283             :       ! get total number of fixed atoms
     284             :       ! nkinds is the kinds of molecules, not atoms
     285         410 :       nkinds = molecule_kinds%n_els
     286        1026 :       DO ikind = 1, nkinds
     287         616 :          molecule_kind => molecule_kind_set(ikind)
     288             :          CALL get_molecule_kind(molecule_kind, &
     289             :                                 nfixd=nfixed_atoms, &
     290         616 :                                 fixd_list=fixd_list)
     291        1026 :          IF (ASSOCIATED(fixd_list)) THEN
     292         224 :             DO ii = 1, nfixed_atoms
     293         224 :                IF (.NOT. fixd_list(ii)%restraint%active) THEN
     294          56 :                   iparticle = fixd_list(ii)%fixd
     295          56 :                   ind = (iparticle - 1)*3
     296             :                   ! apply constraint to nvec
     297          56 :                   SELECT CASE (fixd_list(ii)%itype)
     298             :                   CASE (use_perd_x)
     299           0 :                      vec(ind + 1) = 0.0_dp
     300             :                   CASE (use_perd_y)
     301           0 :                      vec(ind + 2) = 0.0_dp
     302             :                   CASE (use_perd_z)
     303           0 :                      vec(ind + 3) = 0.0_dp
     304             :                   CASE (use_perd_xy)
     305           0 :                      vec(ind + 1) = 0.0_dp
     306           0 :                      vec(ind + 2) = 0.0_dp
     307             :                   CASE (use_perd_xz)
     308           0 :                      vec(ind + 1) = 0.0_dp
     309           0 :                      vec(ind + 3) = 0.0_dp
     310             :                   CASE (use_perd_yz)
     311           0 :                      vec(ind + 2) = 0.0_dp
     312           0 :                      vec(ind + 3) = 0.0_dp
     313             :                   CASE (use_perd_xyz)
     314          56 :                      vec(ind + 1) = 0.0_dp
     315          56 :                      vec(ind + 2) = 0.0_dp
     316          56 :                      vec(ind + 3) = 0.0_dp
     317             :                   END SELECT
     318             :                END IF ! .NOT.fixd_list(ii)%restraint%active
     319             :             END DO ! ii
     320             :          END IF ! ASSOCIATED(fixd_list)
     321             :       END DO ! ikind
     322         410 :    END SUBROUTINE dimer_fixed_atom_control
     323             : 
     324           0 : END MODULE dimer_types

Generated by: LCOV version 1.15