1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Handling of the Wiener process currently employed in turn of the 8!> Langevin dynamics. 9!> \par History 10!> none 11!> \author Matthias Krack (05.07.2005) 12! ************************************************************************************************** 13MODULE wiener_process 14 15 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 16 USE cp_para_types, ONLY: cp_para_env_type 17 USE cp_subsys_types, ONLY: cp_subsys_get,& 18 cp_subsys_type 19 USE distribution_1d_types, ONLY: distribution_1d_type 20 USE force_env_types, ONLY: force_env_get,& 21 force_env_type 22 USE input_section_types, ONLY: section_vals_get,& 23 section_vals_get_subs_vals,& 24 section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: dp 27 USE md_environment_types, ONLY: get_md_env,& 28 md_environment_type,& 29 need_per_atom_wiener_process 30 USE metadynamics_types, ONLY: meta_env_type 31 USE parallel_rng_types, ONLY: GAUSSIAN,& 32 create_rng_stream,& 33 next_rng_seed,& 34 read_rng_stream,& 35 rng_record_length 36 USE particle_list_types, ONLY: particle_list_type 37 USE simpar_types, ONLY: simpar_type 38 USE string_utilities, ONLY: compress 39#include "../base/base_uses.f90" 40 41 IMPLICIT NONE 42 43 PRIVATE 44 45 ! Global parameters in this module 46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'wiener_process' 47 48 ! Public subroutines 49 PUBLIC :: create_wiener_process, create_wiener_process_cv 50 51CONTAINS 52 53! ************************************************************************************************** 54!> \brief Create a Wiener process for Langevin dynamics and initialize an 55!> independent random number generator for each atom in all force 56!> environment and all the subsystems/fragments therein. 57!> \param md_env ... 58!> \par History 59!> Creation (06.07.2005,MK) 60! ************************************************************************************************** 61 SUBROUTINE create_wiener_process(md_env) 62 63 TYPE(md_environment_type), POINTER :: md_env 64 65 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_wiener_process', & 66 routineP = moduleN//':'//routineN 67 68 CHARACTER(LEN=40) :: name 69 INTEGER :: iparticle, iparticle_kind, & 70 iparticle_local, nparticle, & 71 nparticle_kind, nparticle_local 72 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: seed 73 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 74 TYPE(cp_para_env_type), POINTER :: para_env 75 TYPE(cp_subsys_type), POINTER :: subsys 76 TYPE(distribution_1d_type), POINTER :: local_particles 77 TYPE(force_env_type), POINTER :: force_env 78 TYPE(particle_list_type), POINTER :: particles 79 TYPE(section_vals_type), POINTER :: force_env_section, subsys_section, & 80 work_section 81 TYPE(simpar_type), POINTER :: simpar 82 83 NULLIFY (work_section, force_env) 84 CPASSERT(ASSOCIATED(md_env)) 85 86 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env, & 87 simpar=simpar) 88 89 ![NB] shouldn't the calling process know if it's needed 90 IF (need_per_atom_wiener_process(md_env)) THEN 91 92 CALL force_env_get(force_env, force_env_section=force_env_section, & 93 subsys=subsys) 94 95 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS") 96 97 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, & 98 particles=particles) 99 100 nparticle_kind = atomic_kinds%n_els 101 nparticle = particles%n_els 102 103 ! Allocate the (local) data structures for the Wiener process 104 ALLOCATE (local_particles%local_particle_set(nparticle_kind)) 105 106 DO iparticle_kind = 1, nparticle_kind 107 nparticle_local = local_particles%n_el(iparticle_kind) 108 ALLOCATE (local_particles%local_particle_set(iparticle_kind)%rng(nparticle_local)) 109 DO iparticle_local = 1, nparticle_local 110 NULLIFY (local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream) 111 END DO 112 END DO 113 114 ! Each process generates all seeds. The seed generation should be 115 ! quite fast and in this way a broadcast is avoided. 116 ALLOCATE (seed(3, 2, nparticle)) 117 118 ! Load initial seed (not needed for a restart) 119 seed(:, :, 1) = subsys%seed(:, :) 120 121 DO iparticle = 2, nparticle 122 seed(:, :, iparticle) = next_rng_seed(seed(:, :, iparticle - 1)) 123 END DO 124 125 ! Update initial seed 126 subsys%seed(:, :) = next_rng_seed(seed(:, :, nparticle)) 127 128 ! Create a random number stream (Wiener process) for each particle 129 DO iparticle_kind = 1, nparticle_kind 130 nparticle_local = local_particles%n_el(iparticle_kind) 131 DO iparticle_local = 1, nparticle_local 132 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 133 WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for particle", iparticle 134 CALL compress(name) 135 CALL create_rng_stream(rng_stream=local_particles%local_particle_set(iparticle_kind)% & 136 rng(iparticle_local)%stream, name=name, distribution_type=GAUSSIAN, & 137 extended_precision=.TRUE., seed=seed(:, :, iparticle)) 138 END DO 139 END DO 140 141 DEALLOCATE (seed) 142 143 ! Possibly restart Wiener process 144 NULLIFY (work_section) 145 work_section => section_vals_get_subs_vals(section_vals=subsys_section, & 146 subsection_name="RNG_INIT") 147 CALL init_local_particle_set(distribution_1d=local_particles, & 148 nparticle_kind=nparticle_kind, & 149 work_section=work_section) 150 END IF 151 152 END SUBROUTINE create_wiener_process 153 154! ************************************************************************************************** 155!> \brief Helper routine for create_wiener_process. 156!> \param distribution_1d ... 157!> \param nparticle_kind ... 158!> \param work_section ... 159!> \par History 160!> 01.2014 moved from distribution_1d_types (Ole Schuett) 161! ************************************************************************************************** 162 SUBROUTINE init_local_particle_set(distribution_1d, nparticle_kind, & 163 work_section) 164 165 TYPE(distribution_1d_type), POINTER :: distribution_1d 166 INTEGER, INTENT(in) :: nparticle_kind 167 TYPE(section_vals_type), POINTER :: work_section 168 169 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_local_particle_set', & 170 routineP = moduleN//':'//routineN 171 172 CHARACTER(LEN=rng_record_length) :: rng_record 173 INTEGER :: iparticle, iparticle_kind, & 174 iparticle_local, nparticle_local 175 LOGICAL :: explicit 176 177! ------------------------------------------------------------------------- 178 179 CPASSERT(ASSOCIATED(distribution_1d)) 180 181 IF (ASSOCIATED(work_section)) THEN 182 CALL section_vals_get(work_section, explicit=explicit) 183 IF (explicit) THEN 184 DO iparticle_kind = 1, nparticle_kind 185 nparticle_local = distribution_1d%n_el(iparticle_kind) 186 DO iparticle_local = 1, nparticle_local 187 iparticle = distribution_1d%list(iparticle_kind)%array(iparticle_local) 188 IF (iparticle == distribution_1d%list(iparticle_kind)%array(iparticle_local)) THEN 189 CALL section_vals_val_get(section_vals=work_section, & 190 keyword_name="_DEFAULT_KEYWORD_", & 191 i_rep_val=iparticle, & 192 c_val=rng_record) 193 CALL read_rng_stream(rng_stream=distribution_1d% & 194 local_particle_set(iparticle_kind)% & 195 rng(iparticle_local)%stream, & 196 rng_record=rng_record) 197 END IF 198 END DO 199 END DO 200 END IF 201 END IF 202 203 END SUBROUTINE init_local_particle_set 204 205! ************************************************************************************************** 206!> \brief Create a Wiener process for Langevin dynamics used for 207!> metadynamics and initialize an 208!> independent random number generator for each COLVAR. 209!> \param meta_env ... 210!> \date 01.2009 211!> \author Fabio Sterpone 212!> 213! ************************************************************************************************** 214 SUBROUTINE create_wiener_process_cv(meta_env) 215 216 TYPE(meta_env_type), POINTER :: meta_env 217 218 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_wiener_process_cv', & 219 routineP = moduleN//':'//routineN 220 221 CHARACTER(LEN=40) :: name 222 INTEGER :: i_c 223 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: seed 224 REAL(KIND=dp), DIMENSION(3, 2) :: initial_seed 225 226 IF (.NOT. ASSOCIATED(meta_env)) RETURN 227 228 initial_seed = next_rng_seed() 229 230 DO i_c = 1, meta_env%n_colvar 231 NULLIFY (meta_env%rng(i_c)%stream) 232 END DO 233 234 ! Each process generates all seeds. The seed generation should be 235 ! quite fast and in this way a broadcast is avoided. 236 237 ALLOCATE (seed(3, 2, meta_env%n_colvar)) 238 239 seed(:, :, 1) = initial_seed 240 DO i_c = 2, meta_env%n_colvar 241 seed(:, :, i_c) = next_rng_seed(seed(:, :, i_c - 1)) 242 END DO 243 244 ! Update initial seed 245 initial_seed = next_rng_seed(seed(:, :, meta_env%n_colvar)) 246 247 ! Create a random number stream (Wiener process) for each particle 248 DO i_c = 1, meta_env%n_colvar 249 WRITE (UNIT=name, FMT="(A,I8)") "Wiener process for COLVAR", i_c 250 CALL compress(name) 251 CALL create_rng_stream(rng_stream=meta_env%rng(i_c)%stream, name=name, distribution_type=GAUSSIAN, & 252 extended_precision=.TRUE., seed=seed(:, :, i_c)) 253 END DO 254 DEALLOCATE (seed) 255 256 END SUBROUTINE create_wiener_process_cv 257 258END MODULE wiener_process 259