5
$\begingroup$

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

$\endgroup$
11
  • $\begingroup$ Can you include the autocorrelation function you are getting and a comparison with a rough image of what you expect it to look like? $\endgroup$
    – Tyberius
    Commented Jun 1, 2022 at 16:20
  • $\begingroup$ 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? $\endgroup$
    – Ian Bush
    Commented Jun 1, 2022 at 16:51
  • $\begingroup$ 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. $\endgroup$
    – Ian Bush
    Commented Jun 1, 2022 at 17:21
  • $\begingroup$ 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. $\endgroup$
    – Pascal
    Commented Jun 1, 2022 at 19:17
  • $\begingroup$ Here is a link where the files can be retrieved: Renater $\endgroup$
    – Pascal
    Commented Jun 1, 2022 at 19:25

1 Answer 1

1
$\begingroup$

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.)

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .