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