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