MaginnGroup / PyLAT

GNU General Public License v3.0
104 stars 54 forks source link

Question about viscosity expression #30

Closed Alexandrina-Chen closed 8 months ago

Alexandrina-Chen commented 9 months ago

Dear PyLAT developers,

Thank you for writing this awesome package! I have a question about the calculation of viscosity, to be more specific, the average of the pressure tensor autocorrelation.

The original paper (J. Chem. Theory Comput. 2015, 11, 3537−3546) doesn't describe how different components of the pressure tensor autocorrelation function are averaged ( $\lt P{\alpha\beta}(t)P{\alpha\beta}(0) \gt$ ). And according to another paper by Prof. Maginn, Living J. Comp. Mol. Sci. 2019, 1(1), 6324, the average over the autocorrelation function of the traceless symmetric part of the stress tensor is recommended, as shown below $$P{\alpha\beta} = (\sigma{\alpha\beta} + \sigma{\beta\alpha}) / 2 -\frac{\delta{\alpha\beta}}{3} (\sum{\gamma} \sigma{\gamma\gamma})$$

Then the viscosity is given by $$\eta=\frac{V}{10kB T}\sum{\sigma}\sum_{\beta}\int0^{\inf} \lt P{\alpha\beta}(t)P_{\alpha\beta}(0) \gt dt$$

The factor 10 in the denominator results from assigning weighting factors of 3/3 and 4/3 fir each of the six off-diagonal terms and the three diagonal terms, respectively.

But in visco.py line 156 to 164

        a1=self.llog['pxy'][cutoff:]
        a2=self.llog['pxz'][cutoff:]
        a3=self.llog['pyz'][cutoff:]
        a4=self.llog['pxx'][cutoff:]-self.llog['pyy'][cutoff:]
        a5=self.llog['pyy'][cutoff:]-self.llog['pzz'][cutoff:]
        a6=self.llog['pxx'][cutoff:]-self.llog['pzz'][cutoff:]
        array_array=[a1,a2,a3,a4,a5,a6]
        pv=p.map(autocorrelate,array_array)
        pcorr = (pv[0]+pv[1]+pv[2])/6+(pv[3]+pv[4]+pv[5])/24

I didn't quite get the reason why the average is like this (pv[0]+pv[1]+pv[2])/6+(pv[3]+pv[4]+pv[5])/24, especially the denominators 6 and 24. The off-diagonal terms in the equation above and the code share the same form and have nothing special. Then I have tried to transform the diagonal terms in the equation above into the expression of the code (namely, rewrite the autocorrelation of $P{\alpha\alpha} = \sigma{\alpha\alpha} -\frac{\delta{\alpha\beta}}{3} (\sum{\gamma} \sigma{\gamma\gamma})$ to the combination of the autocorrelations of $\sigma{xx} - \sigma{yy}$, $\sigma{xx} - \sigma{zz}$, and $\sigma{yy} - \sigma_{zz}$), but I don't think they are equal.

I was wondering if there is any reason behind the weighting factors in the code, or if I can be guided to any paper that explain this. Thank you in advance!

Best, Sijia

zhangyongND commented 9 months ago

Hi Sijia, in the code, we followed the convention from this paper: https://doi.org/10.1021/jp062885s.

Alexandrina-Chen commented 8 months ago

@zhangyongND Hi Yong, thank you a lot for guiding me to the paper!