# Programming auto-correlation function calculation

I have a set of data (heat flow) values in x, y and z directions for which I have to calculate the auto-correlation function to get the thermal conductivity. I have written a small Fortran90 code but something is wrong because I get irrelevant results. I am not sure whether this a problem of units or a problem of coding the algorithm.

10000      ! number of data
3.0000     ! time step (dt)
1281.4610  ! cell volume (V)
2000.0     ! temperature (T)
64         ! number of atoms
9.746633E-03 -1.185578E-02 -2.107013E-02
3.909507E-02 -5.904692E-05 -3.044154E-03
7.439500E-02 2.270370E-02 3.536638E-03
1.027890E-01 3.875051E-02 1.377188E-02
1.218158E-01 4.606921E-02 2.570092E-02
1.306354E-01 4.482340E-02 3.692641E-02
...    !   10000 such lines


The function to calculate is:

$$\frac{1}{3Vk_BT^{3/2}}\int_0^{+\infty}{\langle q(0)q(t)\rangle dt}$$

The Fortran code I wrote is the following:

program kappa

implicit none

integer, parameter :: dp=selected_real_kind(15,307)

real(dp), parameter :: kB=1.380649d-23
real(dp), parameter :: eV2J = 1.602d-19
real(dp), parameter :: ang2m = 1.0d-10
real(dp), parameter :: fs2s = 1.0d-15

real(dp), dimension(:,:), allocatable :: q ! heat flow
real(dp), dimension(3) :: q_0, q_t, int_q ! reference heat flow, integrated heat flow at t, and integrated heat flow
real(dp) :: dt ! time step
real(dp) :: V, T ! Volume and temperature
integer :: count_q ! number of steps (corresponds to number of line of HEAT)
integer :: natoms
integer :: ios, i, j

int_q(:) = 0.0_dp

allocate(q(count_q,3))

open(20, file='thermal_cond.dat', status='unknown', action='write')

! store all the heat flow values in q(:,:)
lp1: do i = 1, count_q
read( 10, *, iostat=ios ) q(i,:)
if (ios < 0) exit lp1 ! unexpected end of file: this should not happen
end do lp1

! for each q_0(:) up to t/2 (t=total time of the simulation)
! calculate the autocorrelation function <q(0)q(t)>
do i = 1, int(count_q / 2)
q_0(:) = q(i,:) ! shift the reference heat flow by one step
q_t(:) = 0.0_dp ! initialize the value of <q(0)q(t)>

do j = i + 1, int(count_q / 2) + (i - 1) ! shift the right end of the interval to have it constant (length of count_q/2)
q_t(:) = q_t(:) + q(j,:) * q_0(:)
end do

q_t(:) = q_t(:) * eV2J**2 / ( int(count_q / 2) * natoms * 3.0 * V * ang2m**3 * kB * T**2 ) * dt * fs2s
int_q(:) = int_q(:) + q_t(:)

write(20, '(4F15.9)') (i - 1) * dt, q_t(:)

end do

write(6,*) 'Kappa(x) = ', int_q(1)
write(6,*) 'Kappa(y) = ', int_q(2)
write(6,*) 'Kappa(z) = ', int_q(3)
write(6,*) '< Kappa > = ', sqrt(int_q(1)**2 + int_q(2)**2 + int_q(3)**2)

close(20)

end program kappa_ml


Could anyone tell me if something is wrong in the program? Thank you

• Can you include the autocorrelation function you are getting and a comparison with a rough image of what you expect it to look like?
– Tyberius
Commented Jun 1, 2022 at 16:20
• Somehow providing the input data would really help - we can then compile and run it. Can you provide a small case which illustrates the problem? Commented Jun 1, 2022 at 16:51
• Well I don't know if it is your problem but you could help yourself by working in more natural units for the problem - there are a number of quantities like eV2J**2 and ang2m**3 * kB which will be tiny and may make the floating point you are using struggle to represent the huge range of numbers you are trying to use. Commented Jun 1, 2022 at 17:21
• I try to enclose the full data file. To answer Ian, the energy is given in eV so I assumed I should multiply by eV2J2, but that's still a problem I have not fully soved. As to ang2m3, this is to transform the volume from anstrom to meter. I am sure for the volume. Commented Jun 1, 2022 at 19:17
• Here is a link where the files can be retrieved: Renater Commented Jun 1, 2022 at 19:25