Programming auto-correlation function calculation

Question

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.

The data look like this (I am adding comments):

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
    
    
    read( 10, * ) count_q
    allocate(q(count_q,3))
    
    read( 10, * ) dt
    read( 10, * ) V
    read( 10, * ) T
    read( 10, * ) natoms
    
    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

Answer 1

I cannot comment on the Fortran code, but in general, calculating the integrated autocorrelation correctly is a subtle exercise. There is increasing uncertainty in estimating correlation at increasing time lags, as you will see if you plot your integrand as a function of t, which means that the more of your autocorrelation series you try to integrate over the noisier your final estimate will be.

One thing you need to do is to censor your integral -- instead of taking it to infinite time, stop it at some T which is a suitable multiple of the correlation time. A typical value is between 5 to 10 correlation times. For more details on why this works, and an implementation in Python, see: Autocorrelation time estimation. (The autocorrelation integral you need is simply the autocorrelation time multiplied by the variance.)