LCOV - code coverage report
Current view: top level - src/motion - dimer_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:dc34ec9) Lines: 99 110 90.0 %
Date: 2023-03-24 20:09:49 Functions: 4 8 50.0 %

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

Generated by: LCOV version 1.15