RMPCDMD
correlator.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2016 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
15 
16 module correlator
17  implicit none
18 
19  private
20 
21  public :: correlator_t
22  public :: axial_correlator_t
23  public :: correlate_block_distsq
24  public :: correlate_block_dot
25  public :: get_n_blocks
26  public :: write_correlator_block
27 
29  integer :: block_length
30  integer :: n_blocks
31  integer :: dim
32  integer :: n
33  double precision, allocatable :: value(:,:,:,:)
34  double precision, allocatable :: correlation(:,:,:,:)
35  integer, allocatable :: count(:)
36  contains
37  procedure, public :: init
38  procedure, public :: add
39  end type correlator_t
40 
42  type(correlator_t) :: msd
43  type(correlator_t) :: oacf
44  type(correlator_t) :: vacf
45  type(correlator_t) :: parallel_vacf
46  type(correlator_t) :: transverse_vacf
47  contains
48  procedure :: init => axial_init
49  procedure :: add => axial_add
50  procedure :: add_fast => axial_add_fast
51  procedure :: write => axial_write
52  end type axial_correlator_t
53 
54 contains
55 
57  subroutine init(this, l, b, dim, n)
58  class(correlator_t), intent(out) :: this
59  integer, intent(in) :: l, b
60  integer, intent(in), optional :: dim
61  integer, intent(in), optional :: n
62 
63  if (present(dim)) then
64  this%dim = dim
65  else
66  this%dim = 1
67  end if
68 
69  if (present(n)) then
70  this%n = n
71  else
72  this%n = 1
73  end if
74 
75  this%block_length = l
76  this%n_blocks = b
77 
78  allocate(this%value(this%dim, this%n, l, b))
79  allocate(this%correlation(this%dim, this%n, l, b))
80  allocate(this%count(b))
81 
82  this%value(:,:,:,:) = 0
83  this%correlation(:,:,:,:) = 0
84  this%count(:) = 0
85 
86  end subroutine init
87 
89  subroutine add(this, i, block_operation, x, x_n, xvec, xvec_n)
90  class(correlator_t), intent(inout) :: this
91  integer, intent(in) :: i
92  double precision, intent(in), optional :: x
93  double precision, intent(in), optional :: x_n(:)
94  double precision, intent(in), optional :: xvec(:)
95  double precision, intent(in), optional :: xvec_n(:,:)
96  interface
97  subroutine block_operation(s, c, idx, l, n, dim)
98  double precision, intent(in) :: s(:,:,:)
99  double precision, intent(inout) :: c(:,:,:)
100  integer, intent(in) :: idx
101  integer, intent(in) :: l, n, dim
102  end subroutine block_operation
103  end interface
104 
105  integer :: l, b
106  integer :: i_block, normed_idx, x_count
107  double precision, allocatable :: x_gen(:,:)
108 
109  x_count = 0
110  allocate(x_gen(this%dim, this%n))
111  if (present(x)) then
112  x_count = x_count+1
113  x_gen = x
114  else if (present(x_n)) then
115  x_count = x_count+1
116  x_gen(1,:) = x_n
117  else if (present(xvec)) then
118  x_count = x_count+1
119  x_gen(:,1) = xvec
120  else if (present(xvec_n)) then
121  x_count = x_count+1
122  x_gen = xvec_n
123  end if
124 
125  if (x_count/=1) stop 'multiple (or no) arguments to correlator_t add'
126 
127  l = this%block_length
128  b = this%n_blocks
129 
130  block_loop: do i_block = 0, b-1
131  if (modulo(i, l**i_block) /= 0) exit block_loop
132  normed_idx = modulo(i / l**i_block, l)
133  this%value(:, :, normed_idx+1, i_block+1) = x_gen
134  if (i > l**(i_block+1)) then
135  ! correlate block b, for point normed_idx
136  call block_operation(this%value(:,:,:,i_block+1), &
137  this%correlation(:,:,:,i_block+1), normed_idx, l, this%n, this%dim)
138  this%count(i_block+1) = this%count(i_block+1) + 1
139  end if
140  end do block_loop
141 
142  deallocate(x_gen)
143 
144  end subroutine add
145 
147  subroutine correlate_block_dot(s, c, idx, l, n, dim)
148  double precision, intent(in) :: s(:,:,:)
149  double precision, intent(inout) :: c(:,:,:)
150  integer, intent(in) :: idx, l, n, dim
151 
152  integer :: i, j, n_idx
153  double precision :: current(dim,n)
154 
155  current = s(:,:,idx+1)
156 
157  do i = 0, l-1
158  j = modulo(idx-i, l) + 1
159  do n_idx = 1, n
160  c(:, n_idx, j) = c(:, n_idx, j) + (current(:,n_idx)*s(:,n_idx,i+1))
161  end do
162  end do
163 
164  end subroutine correlate_block_dot
165 
167  subroutine correlate_block_distsq(s, c, idx, l, n, dim)
168  double precision, intent(in) :: s(:,:,:)
169  double precision, intent(inout) :: c(:,:,:)
170  integer, intent(in) :: idx, l, n, dim
171 
172  integer :: i, j, n_idx
173  double precision :: current(dim, n)
174 
175  current = s(:,:,idx+1)
176 
177  do i = 0, l-1
178  j = modulo(idx-i, l) + 1
179  do n_idx = 1, n
180  c(:, n_idx, j) = c(:, n_idx, j) + (current(:,n_idx)-s(:,n_idx,i+1))**2
181  end do
182  end do
183 
184  end subroutine correlate_block_distsq
185 
187  function get_n_blocks(l, n_blocks_max, n_samples) result(n_blocks)
188  integer, intent(in) :: l, n_blocks_max, n_samples
189  integer :: n_blocks
190 
191  do n_blocks = 1, 7
192  if (l**n_blocks >= n_samples/l) exit
193  end do
194 
195  end function get_n_blocks
196 
198  subroutine axial_init(this, block_length, n_samples, n_samples_fast)
199  class(axial_correlator_t), intent(inout) :: this
200  integer, intent(in) :: block_length, n_samples, n_samples_fast
201 
202  integer :: n_blocks
203 
204  n_blocks = get_n_blocks(block_length, 7, n_samples)
205  call this%msd%init(block_length, n_blocks, dim=3)
206  call this%oacf%init(block_length, n_blocks, dim=3)
207 
208  n_blocks = get_n_blocks(block_length, 8, n_samples_fast)
209  call this%vacf%init(block_length, n_blocks, dim=3)
210  call this%parallel_vacf%init(block_length, n_blocks, dim=3)
211  call this%transverse_vacf%init(block_length, n_blocks, dim=3)
212 
213  end subroutine axial_init
214 
216  subroutine axial_add(this, i, com_pos, unit_r)
217  class(axial_correlator_t), intent(inout) :: this
218  integer, intent(in) :: i
219  double precision, intent(in) :: com_pos(3), unit_r(3)
220 
221  call this%msd%add(i, correlate_block_distsq, xvec=com_pos)
222  call this%oacf%add(i, correlate_block_dot, xvec=unit_r)
223 
224  end subroutine axial_add
225 
227  subroutine axial_add_fast(this, i, v_com, unit_r)
228  class(axial_correlator_t), intent(inout) :: this
229  integer, intent(in) :: i
230  double precision, intent(in) :: v_com(3), unit_r(3)
231 
232  double precision :: parallel_v(3), transverse_v(3)
233 
234  parallel_v = dot_product(v_com, unit_r) * unit_r
235  transverse_v = v_com - parallel_v
236  call this%vacf%add(i, correlate_block_dot, xvec=v_com)
237  call this%parallel_vacf%add(i, correlate_block_dot, xvec=parallel_v)
238  call this%transverse_vacf%add(i, correlate_block_dot, xvec=transverse_v)
239 
240  end subroutine axial_add_fast
241 
243  subroutine axial_write(this, correlator_group, sampling, dt, fast_sampling, fast_dt)
244  use hdf5
245  implicit none
246  class(axial_correlator_t), intent(inout) :: this
247  integer(HID_T) :: correlator_group
248  integer, intent(in) :: sampling, fast_sampling
249  double precision, intent(in) :: dt, fast_dt
250 
251  call write_correlator_block(correlator_group, 'mean_square_displacement', &
252  this%msd, sampling, dt)
253  call write_correlator_block(correlator_group, 'orientation_autocorrelation', &
254  this%oacf, sampling, dt)
255 
256  call write_correlator_block(correlator_group, 'velocity_autocorrelation', &
257  this%vacf, fast_sampling, fast_dt)
258  call write_correlator_block(correlator_group, 'parallel_velocity_autocorrelation', &
259  this%parallel_vacf, fast_sampling, fast_dt)
260  call write_correlator_block(correlator_group, 'transverse_velocity_autocorrelation', &
261  this%transverse_vacf, fast_sampling, fast_dt)
262 
263  end subroutine axial_write
264 
266  subroutine write_correlator_block(loc, name, c, step, time)
267  use hdf5
268  use h5md_module
269  implicit none
270  integer(HID_T), intent(inout) :: loc
271  character(len=*), intent(in) :: name
272  type(correlator_t), intent(in) :: c
273  integer, intent(in) :: step
274  double precision, intent(in) :: time
275 
276  integer(HID_T) :: group
277  integer :: error
278 
279  call h5gcreate_f(loc, name, group, error)
280  call h5md_write_dataset(group, 'value', c%correlation)
281  call h5md_write_dataset(group, 'count', c%count)
282  call h5md_write_dataset(group, 'step', step)
283  call h5md_write_dataset(group, 'time', time)
284  call h5gclose_f(group, error)
285 
286  end subroutine write_correlator_block
287 
288 end module correlator
subroutine, public write_correlator_block(loc, name, c, step, time)
Write the data from a correlator_t variable to a HDF5 group.
Definition: correlator.f90:267
subroutine axial_init(this, block_length, n_samples, n_samples_fast)
Initialize an axial correlator.
Definition: correlator.f90:199
subroutine add(this, i, block_operation, x, x_n, xvec, xvec_n)
Add a data element to the correlator.
Definition: correlator.f90:90
subroutine axial_add_fast(this, i, v_com, unit_r)
Add a data element to an axial correlator (fast time scale)
Definition: correlator.f90:228
subroutine init(this, l, b, dim, n)
Initialize a correlator_t variable.
Definition: correlator.f90:58
subroutine axial_add(this, i, com_pos, unit_r)
Add a data element to an axial correlator (slow time scale)
Definition: correlator.f90:217
integer function, public get_n_blocks(l, n_blocks_max, n_samples)
Return the number of blocks needed for a trajectory of length n_samples.
Definition: correlator.f90:188
subroutine, public correlate_block_dot(s, c, idx, l, n, dim)
Correlation routine for dot product.
Definition: correlator.f90:148
Block correlators.
Definition: correlator.f90:16
subroutine axial_write(this, correlator_group, sampling, dt, fast_sampling, fast_dt)
Write all correlators from an axial correlator variable to a HDF5 group.
Definition: correlator.f90:244
subroutine, public correlate_block_distsq(s, c, idx, l, n, dim)
Correlation routine for squared distance.
Definition: correlator.f90:168