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