Hello!
Sorry for the late reply, although I did some research on this topic a while ago it took me some time to read up on this again...
You are correct that the ML_FF file is read and processed via LAMMPS exactly the same way it is used inside VASP when using the VASPml C++ library (ML_LIB = .TRUE.). Hence, the predicted energies and forces should be identical up to numeric rounding errors. To be precise: the major part of computing the predicted force contributions (descriptors, kernel, per-atom and per-neighbor terms) is handled inside the VASPml library and is therefore processed by the same source code. However, since LAMMPS has a much more sophisticated MD parallelization (spatial domain decomposition) the final summation of force contributions over neighboring cells and the computation of stress is handled inside LAMMPS, VASPml just computes the contributions and stores them in arrays provided by LAMMPS.
Now, if we want to calculate the heat flux via per-atom stress contributions we need this virial term:
\(\sum_i \mathsf{W}_i \cdot \mathbf{v}_i
\qquad\text{with}\qquad
\mathsf{W}_i = \sum_{j \ne i} \mathbf{r}_{ij} \otimes \frac{\partial E_j}{\partial \mathbf{r}_i},
\qquad
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i
\)
Atom \(i\) collects, for every neighbor \(j\) whose environment contains \(i\), the vector to that neighbor times the force that the atomic energy of \(j\) exerts on \(i\). The whole term belongs to atom \(i\), because it needs to be multiplied with \(\mathbf{v}_i\) . Unfortunately, atoms \(i\) and \(j\) may not "live" on the same processor, hence they cannot fill their respective arrays without communication. Therefore, we would need to first compute the contributions for all local and ghost atoms and then let LAMMPS communicate/sum up the contributions of \(\mathsf{W}_i\) similarly to the forces (this time 9 components for each atom). This is certainly possible but at the time of writing the LAMMPS interface this was simply not implemented. Hence, the corresponding arrays which are needed by "compute stress/atom" and "compute centroid/stress/atom" are not filled and therefore the heat flux computation cannot work via this method for the VASP force fields inside LAMMPS. As far as I know the situation is similar for other manybody force fields. On the other hand, as you already mentioned, inside VASP the formula above is implemented and accessible via the ML_LHEAT tag.
I am sorry that this was not communicated properly so far... I will update the Wiki and put the implementation of per-atom stress within LAMMPS on our ToDo list. I think the infrastructure in LAMMPS and VASPml is available and we know in principle what needs to be done, so I personally would like to have this feature available. However, since the Müller-Plathe method is our recommended alternative it will not have top priority and hence I cannot give a timeline for this feature.
All the best,
Andreas Singraber
