Therefore, our empirical study results imply that addressing the AI system’s
robustness issue should be an essential and urgent need. Robustness
enhancement techniques and software repair can be two possible solutions to
address this issue. However, even though there have been a lot of studies on a
single DNN model’s robustness enhancement or repair [93, 89, 82, 73, 95], yet
there is no study on the system-level robustness enhancement or repair. Hence,
both these two topics can be interesting research opportunities.

Another interesting finding from analyzing an AI system’s performance on
corrupted data is that different configurations can have different
performances when against the same corruption pattern for an identical system
and the same dataset. One example is the SpeechBrain system’s performance on
the corrupted LibriSpeech dataset, where configuration 1 has much better
performance against Background noise corruption than configuration 2. By
contrast, configuration 2 is much better at against Distortion corruption.
These results imply that different system configurations can bring different
robustness performances. As suggested by international standards ISO 26262
[74] and ISO/PAS 21448 (SOTIF) [75], modular redundancy is an important
strategy to improve system quality and reliability. Therefore, developers can
try to deploy differently configured systems in parallel to achieve better
system reliability.

###  5.2 OOD detection and uncertainty estimation techniques are still
effective when analyzing AI systems, while not every module’s OOD or
uncertainty information could reveal the system’s potential risks.

From our empirical study results in RQ1 and RQ2, we find that two general OOD
detection approaches and two general uncertainty estimation metrics are
effective when analyzing subject AI systems, while the choice of each can be
varied. Two general OOD detection techniques used in this study can
effectively measure the data distribution difference between clean data and
corrupted data. Surprisingly, we find that applying OOD detection on different
modules can reveal different OOD data for an identical system. For instance,
in the SpeechBrain system, we can observe that OOD data identified by two
modules are very different, whereas in most cases, the similarity between
top-1% OOD data identified by two modules is less than 5%. Hence the results
imply that when utilizing OOD detection techniques as an approach to risk
assessment, it’s needed to measure the data distribution difference based on
different AI modules in a system.

In the SpeechBrain system and ByteTrack system, two uncertainty estimation
metrics show promising performance in capturing the system’s unstable
performance on corrupted data. While in the ConvLab-2 system, we find that
these uncertainty estimation techniques are only effective in analyzing NLU
(Natural language understanding) module. Different from other systems, in
ConvLab-2, each module in the system will be executed several times in a
complete dialogue loop. Certain modules like DST (dialogue state tracking) and
POL (dialogue policy) will maintain complex internal state information
throughout the dialogue. Thus the system becomes much more complex compared
with other ones. These results imply that more advanced uncertainty estimation
techniques are needed when analyzing a more complex AI system.

Nevertheless, we also find that not every module’s OOD or uncertainty
information in an AI system can be used as indicators of the system’s
potential risks. For instance, in the SpeechBrain system, both the OOD and
uncertainty scores from the second module have a strong correlation with
system output error. However, there is almost no correlation when it comes to
the first module. These results reveal that, when using OOD or uncertainty to
do a risk assessment, it’s important for developers to decide which module in
the system should be measured. Besides, even though the second module shows
more obvious effects on the system output error, it’s still unclear how the
output from the previous module has affected it. When the system triggers a
prediction error, how to do fault localization inside the system can be
another promising research direction.

###  5.3 More advanced risk assessment techniques are needed to enhance the AI
system’s quality.

In RQ4, as an early attempt or AI system risk assessment, we use simple
strategies of combining OOD and uncertainty information to reject data
instances that the system might not be able to handle. Our results show that
even such a simple strategy can further improve the system’s performance
against common corruption patterns by up to 38% (Table 26 & 27). These
empirical study results show the potential of making risk assessments for
improving the AI system’s robustness, as well as the potential of using OOD
detection techniques and uncertainty estimation as risk assessment approaches.
While in future work, more complex strategies can be explicitly designed for
different AI systems and usage scenarios toward building a more reliable AI
system.

Besides, in this study, we try to focus on making the risk assessment from the
data distribution and uncertainty perspectives. However, the general AI system
quality assurance can also be achieved from some other perspectives, _e.g_.,
leveraging testing and debugging methods to identify an AI system’s ability
boundaries. In the unit-level analysis, DNN debugging and testing techniques
have been widely used to improve a single DNN model’s quality and reliability.
Yet there is neither an exploratory study on how these unit-level
testing/debugging approaches perform on system-level analysis nor methods
specifically designed for AI systems at multiple levels. Therefore, we believe
studying the quality assurance of AI systems from multiple angles and levels
could be another important research direction.

##  6 Threat to Validity

In this section, we summarize the threats to the validity of our study.

Internal validity. The internal threats to validity mainly lie in both the
choice of measurement metrics for OOD detection and uncertainty and the choice
of performance metrics for different AI systems. For the measurement metrics,
specifically Maha-d and KDE in OOD detection and VR and VRO in uncertainty,
their effectiveness has mostly been shown in image classification tasks but
not in other domains. They may not be suitable for tasks like audio and object
tracking. Besides, for modules that do not directly handle prediction tasks,
e.g., the first module of the SpeechBrain system, uncertainty calculation may
not be representative enough. To mitigate these threats, we first extend the
definition of uncertainty measurement in Section 3.2 and perform comprehensive
experiments for all the measurement metrics in RQ1 and RQ2, demonstrating that
they are indeed capable of capturing the differences between clean data and
corrupted data regardless of different kinds of modules. For the performance
metrics, if they are not chosen properly, there may be bias when evaluating
the AI systems. To reduce the bias, we report the standard metrics commonly
used in their own domains.

Construct validity. The construct threats mainly lie in the randomness
inherent in our experiments. Specifically, the measurement of uncertainty for
AI systems and creating corrupted datasets both incur randomness in our
evaluation. To reduce the influence of randomness on uncertainty measurement,
we repeat our evaluations on a large number of corrupted test instances for
multiple different configurations of each AI system. We additionally analyze
the results using two uncertainty metrics, namely VR and VRO. To reduce the
influence of randomness in the data corruption process, we select multiple
(_i.e_., three) different levels for each corruption and apply the corruption
to diverse datasets. We further perform a large-scale statistical analysis
using tests like Wilcoxon signed-rank test and Kendall’s $\tau$ test on the
results to reduce the threats.

External validity. The external threats to validity include the selection of
AI systems and the selection of data corruption patterns. Improper selection
may hinder the generalization of the results. To mitigate the threat, we first
perform a survey on available AI systems and select representative systems
across diverse input domains and tasks under the criteria described in Sec.
3.1. To reduce the bias of the AI systems towards one specific corruption
pattern, we select 5 common image corruption patterns, 5 text corruption
patterns and 5 audio corruption patterns with three different severity levels.
Overall, our large-scale evaluation takes 5000+ (NVIDIA A4000 16GB VRAM) GPU
hours, on which we tried our greatest extent based on the available
computational resource. We hope the results obtained through such a large-
scale evaluation could be helpful to inspire future research toward designing
more advanced quality assurance techniques for building trustworthy AI
systems. Finally, our research process and empirical results on selected AI
systems might not generalize to all AI systems. To mitigate this threat, we
tried our best to select AI systems from four different domains.

##  7 Conclusion

In this paper, we propose an exploratory study of AI system risk assessment.
We first investigate a wide range of AI-enabled systems and select four
representative ones from different domains and applications as our study
subjects: Udacity [3], ByteTrack [99], SpeechBrain [68], and ConvLab-2 [43].
Then we leverage the common corruption patterns from different data domains,
_i.e_., image, audio, and text, to create large-scale corrupted datasets to
evaluate these AI systems’ robustness. Overall, we find that these common
corruption patterns pose certain quality and reliability threats to AI
systems. As an early exploratory study of AI systems at the system level, we
then use two general OOD detection techniques and two uncertainty estimation
metrics to analyze the AI systems. We find that these two different aspects’
techniques are still effective in analyzing different modules in the system.
We further investigate the error propagation in AI systems. We also find that
not every module in the AI system strongly correlates to the system output
error. Based on these experimental results and findings, as an early attempt
at AI system risk assessment, we use several simple combination strategies of
OOD and uncertainty information to show their potential to improve the
reliability of AI systems. At the end of this paper, we also discuss several
implications from our exploratory study as well as future research
opportunities.

As an early step in studying the complex AI system, our work confirms that the
quality assurance of AI systems is an urgent demand. Yet the effectiveness of
unit-level analysis techniques can be either limited or hasn’t been
demonstrated. Considering AI systems have been widely used in real-world
environments, we call for software engineering researchers and practitioners’
attention to step forward to not only focusing on the unit (model) level but
also on a more advanced system level which is a largely explored realm. When
moving towards the AI system level, we believe experience and practice from
the software engineering community could bring much potential and played a key
important role towards building safe, reliable and trustworthy AI systems.

## References

  * [1] Who’s responsible when an autonomous car crashes?, http://money.cnn.com/2016/07/07/technology/tesla-liability-risk/index.html (2016)
  * [2] Tesla’s Self-Driving System Cleared in Deadly Crash, https://www.nytimes.com/2017/01/19/business/tesla-model-s-autopilot-fatal-crash.html (2016)
  * [3] Udacity self driving car challenge. https://github.com/udacity/self-driving-car (2016)
  * [4] Apollo (2022), https://apollo.auto/
  * [5] Pedalboard: a python library for manipulating audio. https://github.com/spotify/pedalboard (2022)
  * [6] Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P., Cao, X., Khosravi, A., Acharya, U.R., et al.: A review of uncertainty quantification in deep learning: Techniques, applications and challenges. Information Fusion 76, 243–297 (2021)
  * [7] Amodei, D., Ananthanarayanan, S., Anubhai, R., Bai, J., Battenberg, E., Case, C., Casper, J., Catanzaro, B., Cheng, Q., Chen, G., et al.: Deep speech 2: End-to-end speech recognition in english and mandarin. In: International conference on machine learning. pp. 173–182. PMLR (2016)
  * [8] Amodei, D., Olah, C., Steinhardt, J., Christiano, P., Schulman, J., Mané, D.: Concrete problems in ai safety. arXiv preprint arXiv:1606.06565 (2016)
  * [9] Ardila, R., Branson, M., Davis, K., Kohler, M., Meyer, J., Henretty, M., Morais, R., Saunders, L., Tyers, F., Weber, G.: Common voice: A massively-multilingual speech corpus. In: Proceedings of the 12th Language Resources and Evaluation Conference. pp. 4218–4222 (2020)
  * [10] Ayhan, M.S., Berens, P.: Test-time data augmentation for estimation of heteroscedastic aleatoric uncertainty in deep neural networks (2018)
  * [11] Baevski, A., Zhou, Y., Mohamed, A., Auli, M.: wav2vec 2.0: A framework for self-supervised learning of speech representations. Advances in Neural Information Processing Systems 33, 12449–12460 (2020)
  * [12] Bernardin, K., Elbs, A., Stiefelhagen, R.: Multiple object tracking performance metrics and evaluation in a smart room environment. In: Sixth IEEE International Workshop on Visual Surveillance, in conjunction with ECCV. vol. 90. Citeseer (2006)
  * [13] Bewley, A., Ge, Z., Ott, L., Ramos, F., Upcroft, B.: Simple online and realtime tracking. In: 2016 IEEE international conference on image processing (ICIP). pp. 3464–3468. IEEE (2016)
  * [14] Blundell, C., Cornebise, J., Kavukcuoglu, K., Wierstra, D.: Weight uncertainty in neural network. In: International conference on machine learning. pp. 1613–1622. PMLR (2015)
  * [15] Bojarski, M., Chen, C., Daw, J., Değirmenci, A., Deri, J., Firner, B., Flepp, B., Gogri, S., Hong, J., Jackel, L., et al.: The nvidia pilotnet experiments. arXiv preprint arXiv:2010.08776 (2020)
  * [16] Bojarski, M., Del Testa, D., Dworakowski, D., Firner, B., Flepp, B., Goyal, P., Jackel, L.D., Monfort, M., Muller, U., Zhang, J., et al.: End to end learning for self-driving cars. arXiv preprint arXiv:1604.07316 (2016)
  * [17] Chi, L., Mu, Y.: Deep steering: Learning end-to-end driving model from spatial and temporal visual cues. arXiv preprint arXiv:1708.03798 (2017)
  * [18] Dahl, G.E., Yu, D., Deng, L., Acero, A.: Context-dependent pre-trained deep neural networks for large-vocabulary speech recognition. IEEE Transactions on audio, speech, and language processing 20(1), 30–42 (2011)
  * [19] Deng, J., Guo, J., Xue, N., Zafeiriou, S.: Arcface: Additive angular margin loss for deep face recognition. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 4690–4699 (2019)
  * [20] Du, S., Guo, H., Simpson, A.: Self-driving car steering angle prediction based on image recognition. arXiv preprint arXiv:1912.05440 (2019)
  * [21] Feng, D., Rosenbaum, L., Dietmayer, K.: Towards safe autonomous driving: Capture uncertainty in the deep neural network for lidar 3d vehicle detection. In: 2018 21st International Conference on Intelligent Transportation Systems (ITSC). pp. 3266–3273. IEEE (2018)
  * [22] Gal, Y., Ghahramani, Z.: Bayesian convolutional neural networks with Bernoulli approximate variational inference. In: 4th International Conference on Learning Representations (ICLR) workshop track (2016)
  * [23] Gal, Y., Ghahramani, Z.: Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In: international conference on machine learning. pp. 1050–1059. PMLR (2016)
  * [24] Gawlikowski, J., Tassi, C.R.N., Ali, M., Lee, J., Humt, M., Feng, J., Kruspe, A., Triebel, R., Jung, P., Roscher, R., et al.: A survey of uncertainty in deep neural networks. arXiv preprint arXiv:2107.03342 (2021)
  * [25] Ge, Z., Liu, S., Wang, F., Li, Z., Sun, J.: Yolox: Exceeding yolo series in 2021\. arXiv preprint arXiv:2107.08430 (2021)
  * [26] Geifman, Y., El-Yaniv, R.: Selectivenet: A deep neural network with an integrated reject option. In: International Conference on Machine Learning. pp. 2151–2159. PMLR (2019)
  * [27] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial examples. In: 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (2015)
  * [28] Graves, A.: Sequence transduction with recurrent neural networks. arXiv preprint arXiv:1211.3711 (2012)
  * [29] Graves, A., Jaitly, N.: Towards end-to-end speech recognition with recurrent neural networks. In: International conference on machine learning. pp. 1764–1772. PMLR (2014)
  * [30] Gustafsson, F.K., Danelljan, M., Schon, T.B.: Evaluating scalable bayesian deep learning methods for robust computer vision. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops. pp. 318–319 (2020)
  * [31] Hannun, A., Case, C., Casper, J., Catanzaro, B., Diamos, G., Elsen, E., Prenger, R., Satheesh, S., Sengupta, S., Coates, A., et al.: Deep speech: Scaling up end-to-end speech recognition. arXiv preprint arXiv:1412.5567 (2014)
  * [32] Hendrycks, D., Dietterich, T.: Benchmarking neural network robustness to common corruptions and perturbations. Proceedings of the International Conference on Learning Representations (2019)
  * [33] Hendrycks, D., Gimpel, K.: A baseline for detecting misclassified and out-of-distribution examples in neural networks. In: 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings (2017), https://openreview.net/forum?id=Hkg4TI9xl
  * [34] Hendrycks, D., Mazeika, M., Dietterich, T.: Deep anomaly detection with outlier exposure. In: International Conference on Learning Representations (2019), https://openreview.net/forum?id=HyxCxhRcY7
  * [35] Hou, Y., Ma, Z., Liu, C., Loy, C.C.: Learning to steer by mimicking features from heterogeneous auxiliary networks. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 33, pp. 8433–8440 (2019)
  * [36] Hsu, Y.C., Shen, Y., Jin, H., Kira, Z.: Generalized odin: Detecting out-of-distribution image without learning from out-of-distribution data. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 10951–10960 (2020)
  * [37] Kahn, J., Lee, A., Hannun, A.: Self-training for end-to-end speech recognition. In: ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). pp. 7084–7088. IEEE (2020)
  * [38] Kim, J., Feldt, R., Yoo, S.: Guiding deep learning system testing using surprise adequacy. In: 2019 IEEE/ACM 41st International Conference on Software Engineering (ICSE). pp. 1039–1049. IEEE (2019)
  * [39] Koh, P.W., Sagawa, S., Marklund, H., Xie, S.M., Zhang, M., Balsubramani, A., Hu, W., Yasunaga, M., Phillips, R.L., Gao, I., et al.: Wilds: A benchmark of in-the-wild distribution shifts. In: International Conference on Machine Learning. pp. 5637–5664. PMLR (2021)
  * [40] Lakshminarayanan, B., Pritzel, A., Blundell, C.: Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in neural information processing systems 30 (2017)
  * [41] Lee, K., Lee, H., Lee, K., Shin, J.: Training confidence-calibrated classifiers for detecting out-of-distribution samples. In: International Conference on Learning Representations (2018), https://openreview.net/forum?id=ryiAv2xAZ
  * [42] Lee, K., Lee, K., Lee, H., Shin, J.: A simple unified framework for detecting out-of-distribution samples and adversarial attacks. Advances in neural information processing systems 31 (2018)
  * [43] Lee, S., Zhu, Q., Takanobu, R., Li, X., Zhang, Y., Zhang, Z., Li, J., Peng, B., Li, X., Huang, M., et al.: Convlab: Multi-domain end-to-end dialog system platform. arXiv preprint arXiv:1904.08637 (2019)
  * [44] Li, H., Li, Y., Porikli, F.: Deeptrack: Learning discriminative feature representations online for robust visual tracking. IEEE Transactions on Image Processing 25(4), 1834–1848 (2015)
  * [45] Li, X., Chen, Y.N., Li, L., Gao, J., Celikyilmaz, A.: End-to-end task-completion neural dialogue systems. arXiv preprint arXiv:1703.01008 (2017)
  * [46] Liang, C., Zhang, Z., Lu, Y., Zhou, X., Li, B., Ye, X., Zou, J.: Rethinking the competition between detection and reid in multi-object tracking. arXiv preprint arXiv:2010.12138 (2020)
  * [47] Liang, S., Li, Y., Srikant, R.: Enhancing the reliability of out-of-distribution image detection in neural networks. In: International Conference on Learning Representations (2018), https://openreview.net/forum?id=H1VGkIxRZ
  * [48] Lin, T.Y., Goyal, P., Girshick, R., He, K., Dollár, P.: Focal loss for dense object detection. In: Proceedings of the IEEE international conference on computer vision. pp. 2980–2988 (2017)
  * [49] Lipton, Z.C.: The mythos of model interpretability: In machine learning, the concept of interpretability is both important and slippery. Queue 16(3), 31–57 (2018)
  * [50] Liu, W., Wang, X., Owens, J., Li, Y.: Energy-based out-of-distribution detection. Advances in Neural Information Processing Systems 33, 21464–21475 (2020)
  * [51] Ma, L., Juefei-Xu, F., Xue, M., Li, B., Li, L., Liu, Y., Zhao, J.: Deepct: Tomographic combinatorial testing for deep learning systems. In: 2019 IEEE 26th International Conference on Software Analysis, Evolution and Reengineering (SANER). pp. 614–618. IEEE (2019)
  * [52] Ma, L., Juefei-Xu, F., Zhang, F., Sun, J., Xue, M., Li, B., Chen, C., Su, T., Li, L., Liu, Y., et al.: Deepgauge: Multi-granularity testing criteria for deep learning systems. In: Proceedings of the 33rd ACM/IEEE International Conference on Automated Software Engineering. pp. 120–131 (2018)
  * [53] Ma, L., Zhang, F., Sun, J., Xue, M., Li, B., Juefei-Xu, F., Xie, C., Li, L., Liu, Y., Zhao, J., et al.: Deepmutation: Mutation testing of deep learning systems. In: 2018 IEEE 29th International Symposium on Software Reliability Engineering (ISSRE). pp. 100–111. IEEE (2018)
  * [54] Ma, W., Papadakis, M., Tsakmalis, A., Cordy, M., Traon, Y.L.: Test selection for deep learning systems. ACM Transactions on Software Engineering and Methodology (TOSEM) 30(2), 1–22 (2021)
  * [55] Mairesse, F., Gasic, M., Jurcicek, F., Keizer, S., Thomson, B., Yu, K., Young, S.: Spoken language understanding from unaligned data using discriminative classification models. In: 2009 IEEE International Conference on Acoustics, Speech and Signal Processing. pp. 4749–4752. IEEE (2009)
  * [56] Malinin, A., Gales, M.: Predictive uncertainty estimation via prior networks. Advances in neural information processing systems 31 (2018)
  * [57] Milan, A., Leal-Taixé, L., Reid, I., Roth, S., Schindler, K.: Mot16: A benchmark for multi-object tracking. arXiv preprint arXiv:1603.00831 (2016)
  * [58] Morteza, P., Li, Y.: Provable guarantees for understanding out-of-distribution detection. Proceedings of the AAAI Conference on Artificial Intelligence (2022)
  * [59] Ovadia, Y., Fertig, E., Ren, J., Nado, Z., Sculley, D., Nowozin, S., Dillon, J., Lakshminarayanan, B., Snoek, J.: Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift. Advances in neural information processing systems 32 (2019)
  * [60] Panayotov, V., Chen, G., Povey, D., Khudanpur, S.: Librispeech: an asr corpus based on public domain audio books. In: 2015 IEEE international conference on acoustics, speech and signal processing (ICASSP). pp. 5206–5210. IEEE (2015)
  * [61] Papadopoulos, A.A., Rajati, M.R., Shaikh, N., Wang, J.: Outlier exposure with confidence control for out-of-distribution detection. Neurocomputing 441, 138–150 (2021)
  * [62] Papakipos, Z., Bitton, J.: Augly: Data augmentations for robustness (2022)
  * [63] Papineni, K., Roukos, S., Ward, T., Zhu, W.J.: Bleu: a method for automatic evaluation of machine translation. In: Proceedings of the 40th annual meeting of the Association for Computational Linguistics. pp. 311–318 (2002)
  * [64] Pei, K., Cao, Y., Yang, J., Jana, S.: Deepxplore: Automated whitebox testing of deep learning systems. In: proceedings of the 26th Symposium on Operating Systems Principles. pp. 1–18 (2017)
  * [65] Peng, Z., Yang, J., Chen, T.H., Ma, L.: A first look at the integration of machine learning models in complex autonomous driving systems: A case study on apollo. In: Proceedings of the 28th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering. pp. 1240–1250 (2020)
  * [66] Raghu, M., Blumer, K., Sayres, R., Obermeyer, Z., Kleinberg, B., Mullainathan, S., Kleinberg, J.: Direct uncertainty prediction for medical second opinions. In: International Conference on Machine Learning. pp. 5281–5290. PMLR (2019)
  * [67] Ramalho, T., Miranda, M.: Density estimation in representation space to predict model uncertainty. In: International Workshop on Engineering Dependable and Secure Machine Learning Systems. pp. 84–96. Springer (2020)
  * [68] Ravanelli, M., Parcollet, T., Plantinga, P., Rouhe, A., Cornell, S., Lugosch, L., Subakan, C., Dawalatabad, N., Heba, A., Zhong, J., et al.: Speechbrain: A general-purpose speech toolkit. arXiv preprint arXiv:2106.04624 (2021)
  * [69] Redmon, J., Farhadi, A.: Yolov3: An incremental improvement. arXiv preprint arXiv:1804.02767 (2018)
  * [70] Rychalska, B., Basaj, D., Gosiewska, A., Biecek, P.: Models in the wild: On corruption robustness of neural nlp systems. In: International Conference on Neural Information Processing. pp. 235–247. Springer (2019)
  * [71] Sensoy, M., Kaplan, L., Kandemir, M.: Evidential deep learning to quantify classification uncertainty. Advances in Neural Information Processing Systems 31 (2018)
  * [72] Shanmugam, D., Blalock, D., Balakrishnan, G., Guttag, J.: When and why test-time augmentation works. arXiv e-prints pp. arXiv–2011 (2020)
  * [73] Sotoudeh, M., Thakur, A.V.: Correcting deep neural networks with small, generalizing patches. In: Workshop on Safety and Robustness in Decision Making (2019)
  * [74] for Standardization, I.O.: ISO 26262:Road vehicles - Functional safety (2011)
  * [75] for Standardization, I.O.: ISO/PAS 21448: Road vehicles - Safety of the intended functionality (2019)
  * [76] Synnaeve, G., Xu, Q., Kahn, J., Likhomanenko, T., Grave, E., Pratap, V., Sriram, A., Liptchinsky, V., Collobert, R.: End-to-end asr: from supervised to semi-supervised learning with modern architectures. arXiv preprint arXiv:1911.08460 (2019)
  * [77] Takanobu, R., Zhu, H., Huang, M.: Guided dialog policy learning: Reward estimation for multi-domain task-oriented dialog. In: Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP). pp. 100–110. Association for Computational Linguistics, Hong Kong, China (Nov 2019), https://aclanthology.org/D19-1010
  * [78] Tian, Y., Pei, K., Jana, S., Ray, B.: Deeptest: Automated testing of deep-neural-network-driven autonomous cars. In: Proceedings of the 40th international conference on software engineering. pp. 303–314 (2018)
  * [79] Ulmer, D., Cinà, G.: Know your limits: Uncertainty estimation with relu classifiers fails at reliable ood detection. In: Uncertainty in Artificial Intelligence. pp. 1766–1776. PMLR (2021)
  * [80] Van Amersfoort, J., Smith, L., Teh, Y.W., Gal, Y.: Uncertainty estimation using a single deep deterministic neural network. In: International conference on machine learning. pp. 9690–9700. PMLR (2020)
  * [81] Wang, G., Li, W., Aertsen, M., Deprest, J., Ourselin, S., Vercauteren, T.: Aleatoric uncertainty estimation with test-time augmentation for medical image segmentation with convolutional neural networks. Neurocomputing 338, 34–45 (2019)
  * [82] Wang, H., Ustun, B., Calmon, F.: Repairing without retraining: Avoiding disparate impact with counterfactual distributions. In: International Conference on Machine Learning. pp. 6618–6627. PMLR (2019)
  * [83] Wang, H., Wang, Y., Zhou, Z., Ji, X., Gong, D., Zhou, J., Li, Z., Liu, W.: Cosface: Large margin cosine loss for deep face recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 5265–5274 (2018)
  * [84] Wang, L., Ouyang, W., Wang, X., Lu, H.: Visual tracking with fully convolutional networks. In: Proceedings of the IEEE international conference on computer vision. pp. 3119–3127 (2015)
  * [85] Watanabe, S., Hori, T., Kim, S., Hershey, J.R., Hayashi, T.: Hybrid ctc/attention architecture for end-to-end speech recognition. IEEE Journal of Selected Topics in Signal Processing 11(8), 1240–1253 (2017)
  * [86] Wen, T.H., Gasic, M., Mrksic, N., Su, P.H., Vandyke, D., Young, S.: Semantically conditioned lstm-based natural language generation for spoken dialogue systems. arXiv preprint arXiv:1508.01745 (2015)
  * [87] Wen, T.H., Vandyke, D., Mrksic, N., Gasic, M., Rojas-Barahona, L.M., Su, P.H., Ultes, S., Young, S.: A network-based end-to-end trainable task-oriented dialogue system. arXiv preprint arXiv:1604.04562 (2016)
  * [88] Wojke, N., Bewley, A., Paulus, D.: Simple online and realtime tracking with a deep association metric. In: 2017 IEEE international conference on image processing (ICIP). pp. 3645–3649. IEEE (2017)
  * [89] Xie, X., Guo, W., Ma, L., Le, W., Wang, J., Zhou, L., Liu, Y., Xing, X.: Rnnrepair: Automatic rnn repair via model-based analysis. In: International Conference on Machine Learning. pp. 11383–11392. PMLR (2021)
  * [90] Xie, X., Ma, L., Juefei-Xu, F., Xue, M., Chen, H., Liu, Y., Zhao, J., Li, B., Yin, J., See, S.: Deephunter: a coverage-guided fuzz testing framework for deep neural networks. In: Proceedings of the 28th ACM SIGSOFT International Symposium on Software Testing and Analysis. pp. 146–157 (2019)
  * [91] Yang, L., Liang, X., Wang, T., Xing, E.: Real-to-virtual domain unification for end-to-end autonomous driving. In: Proceedings of the European conference on computer vision (ECCV). pp. 530–545 (2018)
  * [92] Yang, Z., Zhang, Y., Yu, J., Cai, J., Luo, J.: End-to-end multi-modal multi-task vehicle control for self-driving cars with visual perceptions. In: 2018 24th International Conference on Pattern Recognition (ICPR). pp. 2289–2294. IEEE (2018)
  * [93] Yu, B., Qi, H., Qing, G., Juefei-Xu, F., Xie, X., Ma, L., Zhao, J.: DeepRepair: Style-guided repairing for deep neural networks in the real-world operational environment. IEEE Transactions on Reliability pp. 1–16 (2021). https://doi.org/10.1109/tr.2021.3096332, https://doi.org/10.1109/tr.2021.3096332
  * [94] Zang, X., Rastogi, A., Sunkara, S., Gupta, R., Zhang, J., Chen, J.: Multiwoz 2.2: A dialogue dataset with additional annotation corrections and state tracking baselines (2020)
  * [95] Zhang, H., Chan, W.: Apricot: A weight-adaptation approach to fixing deep learning models. In: 2019 34th IEEE/ACM International Conference on Automated Software Engineering (ASE). pp. 376–387. IEEE (2019)
  * [96] Zhang, M., Zhang, Y., Zhang, L., Liu, C., Khurshid, S.: Deeproad: Gan-based metamorphic testing and input validation framework for autonomous driving systems. In: 2018 33rd IEEE/ACM International Conference on Automated Software Engineering (ASE). pp. 132–142. IEEE (2018)
  * [97] Zhang, X.: Uncertainty-guided testing and robustness enhancement for deep learning systems. In: 2020 IEEE/ACM 42nd International Conference on Software Engineering: Companion Proceedings (ICSE-Companion). pp. 101–103. IEEE (2020)
  * [98] Zhang, X., Xie, X., Ma, L., Du, X., Hu, Q., Liu, Y., Zhao, J., Sun, M.: Towards characterizing adversarial defects of deep learning software from the lens of uncertainty. In: 2020 IEEE/ACM 42nd International Conference on Software Engineering (ICSE). pp. 739–751. IEEE (2020)
  * [99] Zhang, Y., Sun, P., Jiang, Y., Yu, D., Yuan, Z., Luo, P., Liu, W., Wang, X.: Bytetrack: multi-object tracking by associating every detection box. arXiv:2110.06864 [cs] (Oct 2021), arXiv: 2110.06864
  * [100] Zhou, H., Li, W., Kong, Z., Guo, J., Zhang, Y., Yu, B., Zhang, L., Liu, C.: Deepbillboard: Systematic physical-world testing of autonomous driving systems. In: 2020 IEEE/ACM 42nd International Conference on Software Engineering (ICSE). pp. 347–358. IEEE (2020)
  * [101] Zhou, X., Wang, D., Krähenbühl, P.: Objects as points. arXiv preprint arXiv:1904.07850 (2019)
  * [102] Zhu, Q., Zhang, Z., Fang, Y., Li, X., Takanobu, R., Li, J., Peng, B., Gao, J., Zhu, X., Huang, M.: ConvLab-2: An Open-Source Toolkit for Building, Evaluating, and Diagnosing Dialogue Systems. In: Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics: System Demonstrations. pp. 142–149. Association for Computational Linguistics (2020)


# Isothermal and adiabatic elastic constants from virial fluctuations

Andrey Pereverzev<EMAIL_ADDRESS>Department of Chemistry,
University of Missouri,
Columbia, Missouri 65211-7600, USA

###### Abstract

We derive expressions for classical isothermal and adiabatic elastic constants
for periodic systems with the boundary contributions included explicitly. The
potential-dependent part of these expressions is written in terms of potential
energies of atomic groups that make up the total potential energy. It is shown
that in the thermodynamic limit the Born term, which depends on the second
derivatives of potential energy, can be expressed exactly in terms of
equilibrium averages that involve two types of atomic-group virials. As a
result the new form of the Born term involves only first derivatives of either
atomic-group or total potential energies. The derived elastic constants
expressions involving the two forms of the Born terms are tested and compared
using molecular dynamics simulations of crystalline argon and silicon. For
both materials the elastic constants obtained using the two forms of the Born
term are in good agreement. In particular, the new form of the Born term
converges to the same value as the original Born term but at a slower rate.
The results for silicon also agree well with the results from the previous
molecular dynamics studies.

##  I Introduction

Second order elastic constants describe the response of solid material stress
tensor to imposed strain and, as such, they are important for characterizing
materials properties under elastic deformation. The microscopic theory of
elastic constants was first developed by Born (see Ref. [1] and references
therein) for the case of $T=0$. The theory was later extended to finite
temperatures [2, 3, 4, 5, 6].

The computational approaches for calculating elastic constants can be grouped
into two categories: equilibrium fluctuation methods and the direct method.

The equilibrium fluctuation methods can be further subdivided into the strain-
fluctuation and the stress-fluctuation approaches. In the strain-fluctuation
approach, the compliance tensor is calculated from the strain fluctuations in
the isothermal-isostress ensemble [3, 4, 5, 6]. The compliance tensor is then
inverted to obtain elastic constants. The main disadvantage of this method is
that it shows slow convergence and, as a result, requires very long simulation
trajectories.

In the stress-fluctuation method [2, 7, 8, 6, 9], the elastic constants are
calculated from the canonical or microcanonical ensemble averages and require
evaluation of the so-called Born, kinetic, and stress fluctuations terms. The
disadvantage of the strain-fluctuation approach is the necessity to calculate
the Born term which involves the second derivatives of potential energy and
can be quite complicated. The complexity of the Born term calculation for
arbitrary potentials makes it difficult to implement the strain-fluctuation
approach in standard software packages without great effort and risk of
introducing undetectable errors. Recently, a capability to numerically
calculate the original form of the Born term for all potentials was added to
LAMMPS simulation package [10] based on the approach of Zhen and Chu [11].

Numerical tests comparing the stress- and strain-fluctuation methods show that
the stress-fluctuation approach provides faster convergence of the elastic
constants compared to the strain-fluctuation approach [12, 13].

Because of the disadvantages of the stress- and strain-fluctuation methods
just discussed, the elastic constants are often calculated using the so-called
direct method, in which finite strains are applied to the system and the
resulting changes in stress tensor are computed [14]. This approach is fairly
straightforward to implement but it is computationally demanding, especially
for low-symmetry crystals.

Note that the formal elastic constants expressions in the case of the stress-
fluctuation method are derived for infinite systems. Numerical simulations are
performed with finite systems, typically under periodic boundary conditions.
The effect of periodic boundaries and how they are treated in simulations are
rarely discussed explicitly [2, 9].

In this work we revisit elastic constants equations for the stress-fluctuation
approach. Our goal is two-fold. First, we re-derive expressions for isothermal
and adiabatic elastic constants in the form that explicitly includes the
effect of periodic boundary conditions. Second, we show that the Born term can
be written in a form that involves only first derivatives of the potential
energy. Apart from fundamental theoretical interest, the new expression for
the Born term can be used as the foundation of a new method of numerical
calculation of isothertmal and adiabatic elastic constants.

##  II Isothermal and adiabatic elastic constants under periodic boundary
conditions

###  II.1 The Hamiltonian of a periodic system

A conventional way to treat periodic boundary conditions is to assume that the
system of interest which is placed in the central parallelepiped (or
triclinic) box is surrounded by an infinite number of translated identical
image boxes that fill up the space [15, 16, 17, 18]. It is convenient to
specify these boxes with vector ${\bf n}$ given by

${\bf n}=\xi_{a}{\bf a}+\xi_{b}{\bf b}+\xi_{c}{\bf c},$ (1)

where ${\bf a}$, ${\bf b}$, and ${\bf c}$ are the vectors specifying the three
edges of the central box and $\xi_{a}$, $\xi_{b}$, and $\xi_{c}$ can take any
integer values. The central box corresponds to ${\bf n}=0$. Then, the
classical Hamiltonian of a system of $N$ atoms subject to periodic boundary
conditions can be written as

$H=H_{0}+U\left(\\{{\bf x}_{{\bf n}i}\\}\right).$ (2)

Here

$H_{0}=\sum_{i}^{N}\sum_{\alpha}\frac{p^{2}_{i\alpha}}{2m_{i}}$ (3)

is the kinetic energy, in which $p_{i\alpha}$ is the $\alpha$th component of
the momentum of atom $i$ with mass $m_{i}$ and $U\left(\\{{\bf x}_{{\bf
n}i}\\}\right)$ is the system potential energy, which depends on the set of
atomic coordinates $\\{{\bf x}_{{\bf n}i}\\}$ both for the central box and the
image boxes. These coordinates are given by

${\bf x}_{i{\bf n}}={\bf x}_{i}+{\bf n},$ (4)

where ${\bf x}_{i}$ is the coordinate vector of atom $i$ in the central box.
For example, in the case of pair interactions, the potential energy can be
written as [15, 18]

$U\left(\\{{\bf x}_{{\bf n}i}\\}\right)=\sum_{\bf
n}\sum_{\begin{subarray}{c}i=1\\\ j>i\end{subarray}}^{N}u(|{\bf x}_{i}-{\bf
x}_{j{\bf n}}|).$ (5)

Here the sum over ${\bf n}$ extends over all image boxes. However, for
potentials with a cutoff and under the minimum image convention only the
central box and its 26 nearest-neighbor image boxes contribute to the sum over
${\bf n}$.

An alternative but equivalent form of the potential energy for a periodic
system that is useful for numerical applications can be written in terms of
the atomic-group potential energies as was done in Ref. [18]. The system
potential energy is usually given by a sum of two-atom, three-atom, four-atom
or, in general, few-atom interaction terms. Each of these few-atom terms is
commonly referred to as a group [18] and can be specified with an index $k$.
The number of groups, the number of atoms in each group, and the number of
groups in which an atom participates depends on the total potential energy but
otherwise are completely arbitrary. The potential energy of the group $k$ is
denoted by $u_{k}(\\{{\bf x}_{i}^{k}\\})$, where $\\{{\bf x}_{i}^{k}\\}$ is
the set of coordinates for group $k$ with the index $i$ running from 1 to
$N_{k}$, the total number of atom in the group $k$. (Note that the same ${\bf
x}_{i}$ for atom $i$ appears as ${\bf x}_{i}^{k}$ with a different $k$ for
every group in which atom $i$ participates.) The Hamiltonian of a periodic
system in the atomic-group form can be written as

$H=H_{0}+\sum_{k}u_{k}(\\{{\bf x}_{i}^{k}\\}),$ (6)

where the sum over $k$ runs over all groups associated with the central box.
Thompson et al. [18] discuss in detail how such groups are defined and
selected. Both Eq. (2) and Eq. (6) forms of the Hamiltonian are used below to
derive expressions for the elastic constants.

###  II.2 The choice of strain

Elastic constants describe system response to the imposed strain and, in
general, depend on a particular form of strain used. Consider a
parallelepiped-shaped sample of solid material specified by three basis
vectors starting from the common vertex. Any homogeneous transformation of
this sample can be represented by a real $3\times 3$ matrix M that transforms
vectors ${\bf x}$ in the unstrained system to vectors ${\bf x}^{\prime}$ in
the strained system

${\bf x}^{\prime}={\bf M}{\bf x}.$ (7)

The matrix ${\bf M}$ has nine independent components of which six represent
pure strain and three represent rotation or reflection. A convenient way to
separate strain from rotation or reflection is to use the polar decomposition
of the matrix ${\bf M}$ wherein ${\bf M}$ is uniquely written as [19]

${\bf M}={\bf O}{\bf S},$ (8)

where ${\bf O}$ is an orthogonal matrix and ${\bf S}$ is a symmetric matrix.
Thus, the action of matrix ${\bf M}$ can be viewed as inducing a pure strain
with the matrix ${\bf S}$ followed by a pure rotation or reflection, given by
the matrix ${\bf O}$, in which the strain state remains unchanged. Clearly, if
${\bf M}$ is symmetric then ${\bf M}={\bf S}$ and the matrix ${\bf M}$ itself
induces only pure strain with no rotation or reflection. The difference
between the matrix ${\bf S}$ and the identity matrix ${\bf I}$

$\bm{\mathcal{E}}={\bf S}-{\bf I}$ (9)

arguably represents the simplest definition of strain in that it is given by a
linear deviation of the matrix ${\bf S}$ from the identity matrix. From Eq.
(9) the matrix S is written in terms of $\bm{\mathcal{E}}$ as

${\bf{S}}={\bf I}+\bm{\mathcal{E}}.$ (10)

Using Eq. (8) the stress $\bm{\mathcal{E}}$ can be re-expressed through ${\bf
M}$ as

$\bm{\mathcal{E}}=\left({\bf M}^{\text{T}}{\bf M}\right)^{\frac{1}{2}}-{\bf
I},$ (11)

where superscript T denotes the matrix transpose. The strain
$\bm{\mathcal{E}}$ is sometimes referred to as the Biot strain [20].

A measure of strain most often used in the literature on elastic constants is
the Lagrangian strain $\bm{\mathcal{H}}$ given by

$\bm{\mathcal{H}}=\frac{1}{2}\left({\bf M}^{\text{T}}{\bf M}-{\bf
I}\right)=\frac{1}{2}\left({\bf S}^{2}-{\bf
I}\right)=\bm{\mathcal{E}}+\frac{\bm{\mathcal{E}}^{2}}{2},$ (12)

for which

${\bf S}=\left({\bf I}+2\bm{\mathcal{H}}\right)^{\frac{1}{2}}.$ (13)

In this work we choose to define elastic constants using $\bm{\mathcal{E}}$
rather than $\bm{\mathcal{H}}$ for the following reasons. First, it is more
natural to study system response to linear changes in geometry as given by Eq.
(10) rather than the nonlinear ones as given by Eq. (13). Second, because of
the linearity of Eq. (10) the use of $\bm{\mathcal{E}}$ leads to simpler
general expressions for the elastic constants. Third, there exists a one-to-
one correspondence between the elastic constants defined using
$\bm{\mathcal{E}}$ and those defined using $\bm{\mathcal{H}}$ (see Appendix
A). Thus one set of constants can always be converted to the other one if
necessary. (Note that for pair potentials the use of $\bm{\mathcal{H}}$ rather
than $\bm{\mathcal{E}}$ leads to slightly simpler microscopic expressions for
elastic constants. This is generally not true for more complex potentials.)

###  II.3 Isothermal elastic constants

The second-order isothermal elastic constants can be defined using Helmholtz
free energy $A$ as [21]

$C^{T}_{\alpha\beta\mu\nu}=\frac{1}{V}\left(\frac{\partial^{2}A}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right)_{T,\bm{\mathcal{E}}=0}.$
(14)

Here $V$ is the system volume and $\varepsilon_{\alpha\beta}$ denotes the
component of the symmetric strain tensor $\bm{\mathcal{E}}$; subscripts $T$
and $\bm{\mathcal{E}}=0$ mean that the derivatives are evaluated at constant
temperature and for $\bm{\mathcal{E}}=0$. Since the tensor $\bm{\mathcal{E}}$
is symmetric, the derivative with respect to $\varepsilon_{\alpha\beta}$ in
Eq. (14) means
$(\partial/\partial\varepsilon_{\alpha\beta}+\partial/\partial\varepsilon_{\beta\alpha})/2$.

The Helmholtz free energy in Eq. (14) is given by

$A=-k_{B}T\ln Z,$ (15)

where

$Z=\int d{{\bf p}}^{3N}d{{\bf x}}^{3N}e^{-\frac{H}{k_{B}T}}.$ (16)

Here the integration extends over the whole phase space. The dependence of $A$
on strain comes from the dependence of the integral in Eq. (16) on integration
limits in coordinate space. A standard technique to avoid the mathematically
inconvenient strain-dependence of the integration limits is to transfer this
dependence to the Hamiltonian explicitly while holding the integration limits
fixed. This can be achieved by using a suitable variable transformation. This
approach was first used by Bogoliubov to derive the microscopic expression for
pressure as a derivative of $A$ with respect to volume [22]. It was later
generalized to express the dependence of $A$ on the full strain tensor [2, 4,
6].

Here we use the following strain-dependent canonical transformation of
coordinates and momenta of the Hamiltonian (2)

$\widetilde{\bf x}_{{\bf n}i}=\left({\bf I}+\bm{\mathcal{E}}\right){\bf
x}_{{\bf n}i},\qquad\widetilde{\bf p}_{i}=\left({\bf
I}+\bm{\mathcal{E}}\right)^{-1}{\bf p}_{i}.$ (17)

Note that the transformation of coordinates affects both the atoms in the
central box and the image boxes. Inserting $\widetilde{\bf p}_{i}$,
$\widetilde{\bf x}_{{\bf n}i}$ in the Hamiltonian (2) leads to the strain-
dependent transformed Hamiltonian
$\widetilde{H}=\widetilde{H}_{0}+\widetilde{U}$, and the strain-dependent
Helmholtz free energy $A$ for which the derivatives with respect to the strain
components can now be evaluated. Using Eqs. (14, 15, 16) we obtains

$C_{\alpha\beta\mu\nu}^{T}=F_{\alpha\beta\mu\nu}^{T}+K_{\alpha\beta\mu\nu}^{T}+B_{\alpha\beta\mu\nu}^{T},$
(18)

where

$F_{\alpha\beta\mu\nu}^{T}=\frac{1}{k_{B}TV}\left(\left\langle\frac{\partial\widetilde{H}}{\partial\varepsilon_{\alpha\beta}}\right\rangle\left\langle\frac{\partial\widetilde{H}}{\partial\varepsilon_{\mu\nu}}\right\rangle-\left\langle\frac{\partial\widetilde{H}}{\partial\varepsilon_{\alpha\beta}}\frac{\partial\widetilde{H}}{\partial\varepsilon_{\mu\nu}}\right\rangle\right)$
(19)

is the so-called fluctuation term,

$K_{\alpha\beta\mu\nu}^{T}=\frac{1}{V}\left\langle\frac{\partial^{2}\widetilde{H}_{0}}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right\rangle$
(20)

is the kinetic term, and

$B_{\alpha\beta\mu\nu}^{T}=\frac{1}{V}\left\langle\frac{\partial^{2}\widetilde{U}}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right\rangle$
(21)

is the Born term. The derivatives in Eqs. (19, 20, 21) are evaluated at
constant $T$ and for $\bm{\mathcal{E}}=0$ as in Eq. (14) with the
corresponding subscripts dropped for brevity; the brackets denote averaging
over the canonical ensemble.

The derivatives of $\widetilde{H}_{0}$ and $\widetilde{U}$ with respect to
strain components in Eqs. (19, 20, 21) can be rewritten in terms of
derivatives with respect to momenta and coordinates using Eq. (17).
Calculation of these derivatives is rather straightforward for $\widetilde{U}$
because of the linear dependence of $\widetilde{\bf x}_{{\bf n}i}$ on
$\bm{\mathcal{E}}$. The derivatives of $\widetilde{H}_{0}$ can also be
calculated easily if one recalls that
$\sum_{\alpha}\widetilde{p}_{i\alpha}^{2}$ is a dot product of vector
$\widetilde{\bf p}_{i}$ with itself and

$\displaystyle\widetilde{\bf p}_{i}\cdot\widetilde{\bf p}_{i}$
$\displaystyle=$ $\displaystyle{\bf p}_{i}\cdot\left[\left(({\bf
I}+\bm{\mathcal{E}})^{-1}\right)^{2}{\bf p}_{i}\right]$ (22) $\displaystyle=$
$\displaystyle{\bf p}_{i}\cdot\left[({\bf
I}-2\bm{\mathcal{E}}+3\bm{\mathcal{E}}^{2}){\bf
p}_{i}\right]+O(\bm{\mathcal{E}}^{3}).$

Calculating the derivatives explicitly in Eq. (19) we obtain for the
fluctuation term

$F_{\alpha\beta\mu\nu}^{T}=-\frac{V}{k_{B}T}\big{(}\left\langle\sigma_{\alpha\beta}\sigma_{\mu\nu}\right\rangle-\left\langle\sigma_{\alpha\beta}\right\rangle\left\langle\sigma_{\mu\nu}\right\rangle\big{)},$
(23)

where

$\displaystyle\sigma_{\alpha\beta}=-\frac{1}{V}\sum_{i=1}^{N}\frac{p_{i\alpha}p_{i\beta}}{m_{i}}$
$\displaystyle+\frac{1}{2V}\sum_{{\bf n}}\sum_{i=1}^{N}\left(x_{{\bf
n}i\alpha}\frac{\partial U}{\partial x_{{\bf n}i\beta}}+x_{{\bf
n}i\beta}\frac{\partial U}{\partial x_{{\bf n}i\alpha}}\right)$ (24)

is the microscopic stress tensor. The second line of the last equation is the
symmetrized form of the virial tensor (divided by volume) for a system with
periodic boundaries which was derived by Thompson et al. [18] without using
strain derivatives explicitly. This form of the virial was referred to as the
atom form in Ref. [18].

Using Eqs. (20,22) the kinetic term is calculated to be

$K_{\alpha\beta\mu\nu}^{T}=\frac{3Nk_{B}T}{2V}(\delta_{\alpha\mu}\delta_{\beta\nu}+\delta_{\alpha\nu}\delta_{\beta\mu}).$
(25)

Note that the kinetic term (25) is $\frac{3}{4}$ of the kinetic term defined
using the Lagrangian strain [12]. We will use this fact in subsection III.2
when comparing the elastic constants defined using the Biot and Lagrangian
strains.

The Born term is

$\displaystyle B_{\alpha\beta\mu\nu}^{T}=$ $\displaystyle\frac{1}{V}\sum_{{\bf
n},{\bf m}}\sum_{i,j=1}^{N}\left\langle x_{{\bf n}i\alpha}x_{{\bf
m}j\mu}\frac{\partial^{2}U}{\partial x_{{\bf n}i\beta}\partial x_{{\bf
m}j\nu}}\right\rangle_{{\text{sym}}}.$ (26)

Here the subscript sym is introduced to shorten the corresponding expression.
It means that the expression in brackets is the average of four terms, namely
the term shown and the three terms obtained from it by the following subscript
changes $(\alpha\leftrightarrow\beta)$, $(\mu\leftrightarrow\nu)$, and
$(\alpha\leftrightarrow\beta,\mu\leftrightarrow\nu)$. The last equation
represents the atom form of the Born term for a system with periodic
boundaries.

For practical applications it is useful to derive the equivalent atomic-group
forms of Eqs. (24,26). These can be obtained from the atomic-group form of the
Hamiltonian (6) by using the canonical transformation

$\widetilde{\bf x}_{i}^{k}=\left({\bf I}+\bm{\mathcal{E}}\right){\bf
x}_{i}^{k},\qquad\widetilde{\bf p}_{i}=\left({\bf
I}+\bm{\mathcal{E}}\right)^{-1}{\bf p}_{i}.$ (27)

and repeating the steps that lead to Eqs. (24,26). This gives

$\sigma_{\alpha\beta}=-\frac{1}{V}\sum_{i=1}^{N}\frac{p_{i\alpha}p_{i\beta}}{m_{i}}-\frac{1}{V}\sum_{k}W_{\alpha\beta}^{k},$
(28)

where

$W_{\alpha\beta}^{k}=-\frac{1}{2}\sum_{i=1}^{N_{k}}\left(x_{i\alpha}^{k}\frac{\partial
u_{k}}{\partial x_{i\beta}^{k}}+x_{i\beta}^{k}\frac{\partial u_{k}}{\partial
x_{i\alpha}^{k}}\right)$ (29)

is the symmetrized form of the atomic-group virial tensor for a system with
periodic boundaries which was also derived in Ref. [18] without using strain
derivatives explicitly. Thus, the atomic-group virial $W_{\alpha\beta}^{k}$
given by Eq. (29) is expressed in terms of coordinate components of atoms in a
given group and the partial force components due to the group potential energy
$u_{k}$.

For the atomic-group form of the Born term we obtain

$B_{\alpha\beta\mu\nu}^{T}=\frac{1}{V}\sum_{k}\sum_{i,j=1}^{N_{k}}\left\langle
x_{i\alpha}^{k}x_{j\mu}^{k}\frac{\partial^{2}u_{k}}{\partial
x_{i\beta}^{k}\partial x_{j\nu}^{k}}\right\rangle_{\text{sym}}.$ (30)

Equation (30) can be viewed as an extension of Eq. (29) for the group-form of
the virial for periodic systems: the virial is obtained from the first
derivatives of $\widetilde{U}$ with respect to strain components whereas Eq.
(30) is derived using the second derivatives. Equation (30) can be used
directly to calculate the Born term contribution to the elastic constants:
once the potential energy in the group form is expressed in Cartesian
coordinates the second derivatives in Eq. (30) can be calculated. These
expressions can be complicated for potentials that involve few-particle
interactions.

Here we show that the right-hand side of Eq. (30) can be reduced to a form
that involves only first derivatives by using integration by parts with
respect to one of the derivatives. The basic steps of the derivation are
similar to those used in the derivation of the equipartition theorem [23, 24].
A similar approach was used in Ref. [25] to express an ensemble-averaged
matrix of second derivatives in terms of the force-force covariance matrix.
The details of the derivation are given in Appendix B. The final expression,
which is the central result of this work, is

$\displaystyle B_{\alpha\beta\mu\nu}^{T}$ $\displaystyle=$
$\displaystyle\frac{1}{2k_{B}TV}\sum_{k}\left\langle\left(\widehat{W}^{k}_{\alpha\beta}W^{k}_{\mu\nu}+W^{k}_{\alpha\beta}\widehat{W}^{k}_{\mu\nu}\right)\right\rangle$
(31) $\displaystyle+$
$\displaystyle\frac{1}{2V}\sum_{k}(N_{k}-1)\left(\delta_{\alpha\beta}\left\langle
W_{\mu\nu}^{k}\right\rangle+\delta_{\mu\nu}\left\langle
W_{\alpha\beta}^{k}\right\rangle\right)$ $\displaystyle+$
$\displaystyle\frac{1}{V}\sum_{k}\left(\delta_{\alpha\mu}\left\langle
W^{k}_{\beta\nu}\right\rangle\right)_{\text{sym}},$

where the virial $\widehat{W}_{\alpha\beta}^{k}$ is given by

$\displaystyle\widehat{W}_{\alpha\beta}^{k}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{N_{k}}\left(x_{i\alpha}^{k}\frac{\partial
U}{\partial x_{i\beta}^{k}}+x_{i\beta}^{k}\frac{\partial U}{\partial
x_{i\alpha}^{k}}\right)$ (32)
$\displaystyle-\frac{1}{2}\left(\overline{x}_{\alpha}^{k}\overline{f}_{\beta}^{k}+\overline{x}_{\beta}^{k}\overline{f}_{\alpha}^{k}\right),$

with

$\overline{x}_{\alpha}^{k}=\frac{1}{N_{k}}\sum_{i=1}^{N_{k}}x_{i\alpha}^{k},\qquad\overline{f}_{\alpha}^{k}=-\sum_{i=1}^{N_{k}}\frac{\partial
U}{\partial x_{i\alpha}^{k}}.$ (33)

Comparison of $\widehat{W}_{\alpha\beta}^{k}$ given by Eq. (32) and
$W_{\alpha\beta}^{k}$ given by Eq. (29) shows that
$\widehat{W}_{\alpha\beta}^{k}$ depends on the same atomic coordinate
components but the partial forces are replaced with the total force component
for a given atom. The second line of Eq. (32) involves products of the
arithmetic average of coordinate components of all atoms in the group and the
sum of force components for atoms in the group; these terms ensure that
$\widehat{W}_{\alpha\beta}^{k}$ is translationally invariant. Translational
invariance of $\widehat{W}_{\alpha\beta}^{k}$ in Eq. (32) can be verified by
displacing the system along any of the Cartesian directions and using Eqs.
(33) together with the translational invariance of the total force components.
Thus, the Born term given by Eq. (31) depends on the first derivatives of
either atomic-group or total potential energies but does not contain the
second derivatives. Because of this fact, Eq. (31) can be of interest for
numerical calculation of elastic constants, particularly for complex
potentials that involve few-particles interaction terms.

###  II.4 Adiabatic elastic constants

Adiabatic elastic constants are defined using system internal energy $E$ as
[21]

$C^{S}_{\alpha\beta\mu\nu}=\frac{1}{V}\left(\frac{\partial^{2}E}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right)_{S,\bm{\mathcal{E}}=0},$
(34)

where the derivatives are now evaluated at constant entropy $S$ and for
$\bm{\mathcal{E}}=0$. It was shown in Refs. [26, 7, 6] that the adiabatic
elastic constants $C_{\alpha\beta\mu\nu}^{S}$ can be calculated using the same
expressions as for the isothermal ones but with all canonical-ensemble
averages replaced with averages over the microcanonical ensemble. Thus, Eqs.
(18,19,20,21,23,25,26,30) remain valid for adiabatic constants with the
superscript $S$ replacing the superscript $T$ and the brackets now denoting
microcanonical averaging. We show in Appendix C that Eq. (31) also remains
valid for the adiabatic case (with averages performed over the microcanonical
ensemble). Thus, adiabatic elastic constants can also be calculated using
atomic-group virials.

##  III Numerical verification

The main goal of our numerical verification is to confirm that elastic
constants calculated using Eq. (18) with the original Born term given by Eq.
(30) and the Born term calculated using virials (31) are indeed numerically
identical. To this end, we calculated elastic constants for crystalline argon
and silicon using molecular dynamics (MD). All MD simulations were performed
using the LAMMPS package [10].

###  III.1 Isothermal elastic constants of argon

In the case of argon isothermal elastic constants were calculated for several
values of pressure and temperature. The results for adiabatic constants for
argon are qualitatively similar and will not be reported here. Solid argon
forms a face-centered cubic (fcc) lattice. There are three independent elastic
constants for cubic crystals [27], which are commonly chosen to be
$C_{1111}^{T},C_{1122}^{T}$, and $C_{1212}^{T}$ whose values we report below.

The atomic interactions in argon were modeled with the Lennard-Jones potential

$V(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(35)

where $r$ is the distance between the two atoms, the parameters $\sigma$ and
$\epsilon$ were set to be $3.4$ Å and $1.67\times 10^{-21}$ J, respectively
[28, 27]. The cubic simulation cell consisting of $12\times 12\times 12$
conventional four-atom fcc unit cells was used. The cutoff value for the
potential was set at the relatively high value of 28 Å to minimize any
possible cutoff effects. A time step of 0.5 fs was used. Atomic coordinates
were recorded every 100 fs. These coordinates were used in Eqs. (30) and (31)
to calculate two forms of the Born term as functions of time. The stress
tensor components for the fluctuation term were calculated using LAMMPS.
Elastic constants were calculated at zero pressure and at $P=1$ GPa. For each
of the two pressures three temperature values were considered. These
temperature values corresponded, approximately, to 2%, 40%, and 90% of the
argon melting temperature at a given pressure, which we will refer to as the
low, moderate and high temperature, respectively. The melting temperature of
argon is approximately 84 K at zero pressure [29] and 250 K at 1 GPa [30].
Thus, the temperatures of 2 K, 33 K, and 75 K for zero pressure and 5 K, 100
K, and 225 K for the 1 GPa pressure were considered. The fluctuation terms,
the kinetic term, and the two forms of the Born term were calculated from 1 ns
equilibrium NVT trajectory for each set of pressure and temperature values. As
an example, the behavior of these terms as functions of time for
$C_{1111}^{T}$ is shown in Fig. 1 for the case of zero pressure and $T=33$ K.

Figure 1:  Constitutive components of $C_{1111}^{T}$ at zero pressure and
$33\,\,{\text{K}}$ as functions of time: (a) $K_{1111}^{T}$ (blue) and
$F_{1111}^{T}$ (gray), (b) $B_{1111}^{T}$ calculated using Eqs. (30) (red) and
(31) (black). Note different scales for the $y$-axis for the two panels.

Tables 1 and 2 list the elastic constants along with the kinetic, fluctuation,
and Born terms for each set of pressure and temperature values. In the case of
the Born terms two values are listed: the first line shows the values obtained
from Eq. (30) and the second line gives the values obtained from Eq. (31). The
two values for the elastic constants in Tables 1 and 2 are calculated using
the two forms of the Born term with the same kinetic and fluctuation terms.
The uncertainties in Tables 1 and 2 correspond to one-sigma standard error.
They were calculated using statistical analysis for the time-dependent
correlated data [31].

Table 1: Isothermal elastic constants of argon at zero pressure. The constitutive kinetic, fluctuation, and Born terms are also shown. The two values for the Born term are obtained from Eqs. (30) and (31). The units are GPa.  | 2 K | 33 K | 75 K
---|---|---|---
$K_{1111}^{T}$ | $0.0022947\pm 0.0000003$ | $0.036483\pm 0.000005$ | $0.07742\pm 0.00001$
$F_{1111}^{T}$ | $-0.0440\pm 0.0004$ | $-0.695\pm 0.006$ | $-1.50\pm 0.01$
$B_{1111}^{T}$ | $4.4324\pm 0.0001$ | $3.9840\pm 0.0003$ | $3.3035\pm 0.0003$
$4.43\pm 0.01$ | $3.981\pm 0.003$ | $3.302\pm 0.001$
$C_{1111}^{T}$ | $4.3907\pm 0.0004$ | $3.326\pm 0.006$ | $1.88\pm 0.01$
$4.38\pm 0.01$ | $3.323\pm 0.007$ | $1.88\pm 0.01$
$K_{1122}^{T}$ | $0$ | $0$ | $0$
$F_{1122}^{T}$ | $-0.0268\pm 0.0003$ | $-0.381\pm 0.005$ | $-0.80\pm 0.01$
$B_{1122}^{T}$ | $2.5339\pm 0.0003$ | $2.2525\pm 0.0001$ | $1.8146\pm 0.0001$
$2.53\pm 0.03$ | $2.252\pm 0.001$ | $1.814\pm 0.001$
$C_{1122}^{T}$ | $2.5071\pm 0.0004$ | $1.872\pm 0.005$ | $1.02\pm 0.01$
$2.50\pm 0.03$ | $1.871\pm 0.005$ | $1.02\pm 0.01$
$K_{1212}^{T}$ | $0.0011473\pm 0.0000002$ | $0.018242\pm 0.000002$ | $0.038708\pm 0.000005$
$F_{1212}^{T}$ | $-0.0184\pm 0.0002$ | $-0.306\pm 0.002$ | $-0.625\pm 0.005$
$B_{1212}^{T}$ | $2.5343\pm 0.0003$ | $2.2586\pm 0.0001$ | $1.8276\pm 0.0001$
$2.53\pm 0.04$ | $2.258\pm 0.001$ | $1.827\pm 0.001$
$C_{1212}^{T}$ | $2.5170\pm 0.0004$ | $1.971\pm 0.002$ | $1.241\pm 0.005$
$2.51\pm 0.04$ | $1.970\pm 0.002$ | $1.241\pm 0.005$
Table 2: Isothermal elastic constants of argon at $P=1$ GPa. The constitutive kinetic, fluctuation, and Born terms are also shown. The two values for the Born term are obtained from Eqs. (30) and (31). The units are GPa.  | 5 K | 100 K | 225 K
---|---|---|---
$K_{1111}^{T}$ | $0.006753\pm 0.000001$ | $0.13024\pm 0.00001$ | $0.27713\pm 0.00004$
$F_{1111}^{T}$ | $-0.112\pm 0.001$ | $-2.08\pm 0.02$ | $-4.43\pm 0.04$
$B_{1111}^{T}$ | $13.2845\pm 0.0001$ | $12.5129\pm 0.0004$ | $11.5724\pm 0.0005$
$13.29\pm 0.02$ | $12.512\pm 0.002$ | $11.572\pm 0.001$
$C_{1111}^{T}$ | $13.180\pm 0.001$ | $10.56\pm 0.02$ | $7.42\pm 0.04$
$13.18\pm 0.02$ | $10.56\pm 0.02$ | $7.42\pm 0.04$
$K_{1122}^{T}$ | $0$ | $0$ | $0$
$F_{1122}^{T}$ | $-0.0618\pm 0.0007$ | $-1.09\pm 0.01$ | $-2.00\pm 0.03$
$B_{1122}^{T}$ | $7.6040\pm 0.0001$ | $7.0648\pm 0.0003$ | $6.3332\pm 0.0005$
$7.605\pm 0.008$ | $7.065\pm 0.001$ | $6.333\pm 0.002$
$C_{1122}^{T}$ | $7.5422\pm 0.0007$ | $5.98\pm 0.01$ | $4.33\pm 0.03$
$7.543\pm 0.008$ | $5.98\pm 0.01$ | $4.33\pm 0.03$
$K_{1212}^{T}$ | $0.0033765\pm 0.0000005$ | $0.065118\pm 0.000006$ | $0.13856\pm 0.00002$
$F_{1212}^{T}$ | $-0.0489\pm 0.0004$ | $-0.898\pm 0.007$ | $-1.85\pm 0.02$
$B_{1212}^{T}$ | $7.1051\pm 0.0001$ | $6.5865\pm 0.0002$ | $5.8795\pm 0.0005$
$7.11\pm 0.01$ | $6.587\pm 0.001$ | $5.880\pm 0.002$
$C_{1212}^{T}$ | $7.0595\pm 0.0004$ | $5.753\pm 0.007$ | $4.16\pm 0.02$
$7.06\pm 0.01$ | $5.754\pm 0.007$ | $4.16\pm 0.02$

One can see from Tables 1 and 2 that elastic constants calculated using two
forms of the Born term are indeed the same within uncertainties for all
pressure and temperature values considered. Note that in all cases the Born
term calculated using Eq. (31) converges slower than the one obtained from the
original expression (30), as reflected by the corresponding uncertainties in
Tables 1 and 2. (See also panel (b) of Fig. 1). However, for moderate and high
temperatures the Born term in the form of Eq. (31) converges faster than the
fluctuation term (cf. Fig. 1), whereas for low temperatures it converges
slower than the fluctuation term. Thus, for moderate and high temperatures the
overall convergence of the elastic constants calculated using the Born term
given by Eq. (31) is approximately the same as that for the elastic constants
obtained using the Born term given by Eq. (30). For low temperatures the
elastic constants calculated using Eq. (31) converge slower.

###  III.2 Adiabatic elastic constants of silicon

In the case of silicon, adiabatic rather than isothermal elastic constants
were calculated. This was done for the following two reasons: to show that the
new expression for the Born term can be also be applied to obtain adiabatic
constants and to compare the calculated values to earlier MD-based adiabatic
results for elastic constants of silicon [32, 11]. Similarly to fcc argon,
silicon is a cubic crystal. Thus, only $C^{S}_{1111}$, $C^{S}_{1122}$, and
$C^{S}_{1212}$ needed to be calculated. The system was studied at zero
pressure and the following four temperatures: 300 K, 888 K, 1164 K, and 1477
K. The latter three temperatures corresponded to the ones considered in Refs.
[32, 11].

The atomic interactions in silicon were modeled with the Stillinger-Weber
potential [33] which consists of two- and three-atom interaction terms. A
cubic simulation cell of 216 silicon atoms was used. All force-field and other
simulation parameters reported in Refs. [32, 11] were kept except for the
total number of simulation steps. References [32] and [11] used 150000 and 1
million steps, respectively. We used 2 million steps for the three higher
temperatures and 6 million steps for 300 K to obtain good convergence of the
Born terms obtained using Eq. (31).

Kinetic, fluctuation, and Born terms were calculated using microcanonical
averaging, i. e. NVE trajectories, for each of the four temperatures. Two
forms of the Born terms were calculated: the original one using Eq. (30) and
the new form based on Eq. (31). For each of these two forms both two-atom and
three-atom contributions to the total Born terms were calculated. These
results are summarized in Table 3. The two values for the elastic constants
reported in Table 3 were calculated using the two forms of the total Born term
with the same kinetic and fluctuation terms.

One can see that elastic constants calculated using two forms of the Born term
agree with each other within uncertainties for all temperatures considered.
This remains true for the two-atom, three-atom, and total Born terms
calculated using the two approaches. For all cases the Born terms calculated
using Eq. (31) converges slower than the ones obtained from the original
expression (30), as reflected by the corresponding uncertainties. Comparison
of uncertainties for the fluctuation terms and the Born terms calculated from
Eq. (31) shows that $B_{1111}^{S}$ and $B_{1122}^{S}$ converge slower and
$B_{1212}^{S}$ much faster than the corresponding fluctuation terms for all
four temperatures.

The results of Ref. [11] for the elastic constants of silicon along with the
constitutive Born, kinetic, and fluctuation terms for 888 K, 1164 K, and 1477
K are also listed in Table 3. Also shown are the elastic constants reported in
Ref. [32] for the same three temperatures. When comparing our results to the
previously published data we need to bear in mind that the elastic constants
reported in Refs. [32, 11] were calculated using the Lagrangian strain whereas
our results were based of the Biot strain. Fortunately, the transformation
between the Lagrangian-strain and Biot-strain elastic constants (and their
constitutive terms) is straightforward. It follows from Eq. (40) that at zero
pressure the adiabatic elastic constants obtained using the Lagrangian and
Biot strains are identical, thus they can be compared directly. However, this
is not true for the kinetic and Born terms for the two strains considered
separately. It was noted following Eq. (25) that the Biot kinetic term is
$\frac{3}{4}$ of the Lagrangian one. The fluctuation terms for the two strain
measures are identical because the stress tensors obtained using the two
strains are identical. This implies that the Biot Born term is equal the
Lagrangian Born term plus one fourth of the Lagrangian kinetic term. Thus, the
original Lagrangian results of Ref. [11] for the kinetic and total Born terms
were converted to the Biot form using the recipe just above when reporting
them in Table 3.

One can see that the total Born terms of Ref. [11] and both forms of the total
Born terms calculated by us are generally in good agreement: for most cases
they are identical within uncertainties. Elastic constants calculated by us
also compare well to the elastic constants from Refs. [32] and [11]. The
larger differences for $C_{1212}^{S}$ at 1164 K and 1477 K are still within
uncertainties and are due to the very large uncertainties of the corresponding
fluctuation terms.

Table 3: Adiabatic elastic constants of Si at zero pressure. The constitutive kinetic, fluctuation, and two-atom, three-atom and total Born terms are also shown. The two listed values for the Born terms are obtained from Eqs. (30) and (31). Results of Refs. [11] and [32] are also shown. The units are $10^{10}$ Pa. | 300 K | 888 K | 1164 K | 1477 K
---|---|---|---|---
$K_{1111}^{S}$ | $0.0619\pm 0.0006$ | $0.1819\pm 0.0005$ | $0.2378\pm 0.0007$ | $0.3015\pm 0.0008$
$K_{1111}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $0.18$ | $0.25$ | $0.30$
$F_{1111}^{S}$ | $-0.087\pm 0.002$ | $-0.345\pm 0.009$ | $-0.50\pm 0.01$ | $-0.70\pm 0.01$
$F_{1111}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $-0.34$ | $-0.53$ | $-0.71$
$B_{1111}^{S}{\text{(2-atom)}}$ | $10.1571\pm 0.0004$ | $9.687\pm 0.006$ | $9.17\pm 0.01$ | $8.44\pm 0.02$
$10.14\pm 0.05$ | $9.68\pm 0.06$ | $9.21\pm 0.05$ | $8.49\pm 0.06$
$B_{1111}^{S}{\text{(3-atom)}}$ | $4.781\pm 0.001$ | $4.669\pm 0.004$ | $4.883\pm 0.006$ | $5.31\pm 0.01$
$4.76\pm 0.04$ | $4.68\pm 0.02$ | $4.88\pm 0.01$ | $5.31\pm 0.02$
$B_{1111}^{S}{\text{(total)}}$ | $14.938\pm 0.001$ | $14.356\pm 0.007$ | $14.05\pm 0.01$ | $13.76\pm 0.02$
$14.90\pm 0.06$ | $14.36\pm 0.06$ | $14.09\pm 0.05$ | $13.81\pm 0.06$
$B_{1111}^{S}{\text{(total) Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $14.35$ | $14.01$ | $13.76$
$C_{1111}^{S}$ | $14.912\pm 0.002$ | $14.19\pm 0.01$ | $13.79\pm 0.01$ | $13.35\pm 0.03$
$14.88\pm 0.06$ | $14.20\pm 0.09$ | $13.82\pm 0.05$ | $13.42\pm 0.06$
$C_{1111}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $14.19\pm 0.02$ | $13.73\pm 0.02$ | $13.35\pm 0.02$
$C_{1111}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{Kluge}{}{}]}}}$ |  | $14.14\pm 0.01$ | $13.73\pm 0.03$ | $13.32\pm 0.01$
$K_{1122}^{S}$ | $0$ | $0$ | $0$ | $0$
$F_{1122}^{S}$ | $0.017\pm 0.002$ | $0.006\pm 0.004$ | $-0.01\pm 0.01$ | $-0.04\pm 0.01$
$F_{1122}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $0.01$ | $-0.02$ | $-0.05$
$B_{1122}^{S}{\text{(2-atom)}}$ | $9.980\pm 0.002$ | $9.407\pm 0.006$ | $9.17\pm 0.01$ | $8.41\pm 0.01$
$9.96\pm 0.05$ | $9.40\pm 0.06$ | $9.21\pm 0.05$ | $8.46\pm 0.06$
$B_{1122}^{S}{\text{(3-atom)}}$ | $-2.353\pm 0.001$ | $-1.891\pm 0.005$ | $-1.513\pm 0.007$ | $-0.94\pm 0.01$
$-2.34\pm 0.02$ | $-1.90\pm 0.01$ | $-1.51\pm 0.01$ | $-0.95\pm 0.02$
$B_{1122}^{S}{\text{(total)}}$ | $7.626\pm 0.002$ | $7.516\pm 0.008$ | $7.47\pm 0.01$ | $7.47\pm 0.01$
$7.62\pm 0.05$ | $7.50\pm 0.06$ | $7.50\pm 0.05$ | $7.51\pm 0.06$
$B_{1122}^{S}{\text{(total) Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $7.52$ | $7.47$ | $7.48$
$C_{1122}^{S}$ | $7.643\pm 0.003$ | $7.522\pm 0.009$ | $7.46\pm 0.01$ | $7.42\pm 0.01$
$7.64\pm 0.05$ | $7.51\pm 0.06$ | $7.49\pm 0.05$ | $7.47\pm 0.06$
$C_{1122}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $7.53\pm 0.01$ | $7.45\pm 0.01$ | $7.43\pm 0.01$
$C_{1122}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{Kluge}{}{}]}}}$ |  | $7.52\pm 0.00$ | $7.43\pm 0.01$ | $7.39\pm 0.04$
$K_{1212}^{S}$ | $0.0309\pm 0.0003$ | $0.0909\pm 0.0002$ | $0.1189\pm 0.0004$ | $0.1507\pm 0.0004$
$K_{1212}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $0.09$ | $0.12$ | $0.15$
$F_{1212}^{S}$ | $-5.6\pm 0.4$ | $-5.0\pm 0.3$ | $-5.4\pm 0.2$ | $-5.6\pm 0.3$
$F_{1212}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $-4.92$ | $-5.66$ | $-5.87$
$B_{1212}^{S}$(2-atom) | $10.003\pm 0.001$ | $9.483\pm 0.006$ | $9.089\pm 0.007$ | $8.56\pm 0.01$
$9.99\pm 0.05$ | $9.48\pm 0.06$ | $9.12\pm 0.05$ | $8.61\pm 0.06$
$B_{1212}^{S}$(3-atom) | $0.8285\pm 0.0002$ | $0.969\pm 0.002$ | $1.161\pm 0.004$ | $1.468\pm 0.007$
$0.82\pm 0.01$ | $0.97\pm 0.01$ | $1.16\pm 0.01$ | $1.47\pm 0.01$
$B_{1212}^{S}$(total) | $10.832\pm 0.001$ | $10.451\pm 0.006$ | $10.251\pm 0.007$ | $10.03\pm 0.01$
$10.81\pm 0.05$ | $10.45\pm 0.06$ | $10.28\pm 0.05$ | $10.07\pm 0.06$
$B_{1212}^{S}{\text{(total) Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $10.45$ | $10.22$ | $10.03$
$C_{1212}^{S}$ | $5.3\pm 0.4$ | $5.6\pm 0.3$ | $4.9\pm 0.2$ | $4.6\pm 0.3$
$5.3\pm 0.4$ | $5.6\pm 0.3$ | $5.0\pm 0.2$ | $4.7\pm 0.3$
$C_{1212}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{ZhenChu}{}{}]}}}$ |  | $5.62\pm 0.54$ | $4.68\pm 0.21$ | $4.31\pm 0.36$
$C_{1212}^{S}{\text{ Ref.\cite[cite]{[\@@bibref{Number}{Kluge}{}{}]}}}$ |  | $5.24\pm 0.84$ | $4.57\pm 1.14$ | $4.20\pm 0.83$

##  IV Conclusions

We derived expressions for isothermal and adiabatic elastic constants that
explicitly incorporate the effect of periodic boundaries. These expressions
can be used for numerical calculations of elastic constants and represent the
second-derivative generalizations of expressions for the virials of periodic
systems [18, 17].

We also showed that the Born term can be expressed in the form that involves
only the first derivatives of the atomic-group and total potential energies.
This fact is important from a fundamental standpoint because it means that
knowledge of atomic coordinates and suitably chosen partial and total forces
along with system volume and temperature is sufficient to define isothermal
and adiabatic elastic constants. Equation (31) is also of interest for
numerical calculation of elastic constants (as was done for argon and silicon
in this work): molecular dynamics simulation packages such as LAMMPS can
compute atomic virial tensors given by Eq. (29) and calculation of the Born
term using Eq. (31) can be done with modest code modifications.

Application of Eq. (31) to the Lennard-Jones argon and the Stillinger-Weber
silicon shows very good agreement with the original Born expression given by
Eq. (30) but with a slower convergence. Numerical validation of Eq. (31) for
the cases of more complex potentials along with the studies of the system
size, stress, and temperature dependence is intended by us in the near future.

###### Acknowledgements.

This research was funded by Air Force Office of Scientific Research, grant
numbers FA9550-19-1-0318 and FA9550-22-1-0212, and by AFOSR DURIP equipment
grant FA9550-20-1-0205. The author is grateful to Yuriy Pereverzev and Tommy
Sewell for valuable comments and suggestions.

## Appendix A

Here we derive the relationship between the elastic constants
$C_{\alpha\beta\mu\nu}^{T}$ defined using the Biot strain $\bm{\mathcal{E}}$
and the elastic constants $\widehat{C}_{\alpha\beta\mu\nu}^{T}$ defined using
the Lagrangian strain $\bm{\mathcal{H}}$. We treat tensor $\bm{\mathcal{H}}$
(whose components are denoted by $\eta_{\alpha\beta}$) as a function of
$\bm{\mathcal{E}}$ as given by the last equality of Eq. (12) and apply the
chain rule to express derivatives of $A$ with respect to
$\varepsilon_{\alpha\beta}$ in Eq. (14) in terms of derivatives with respect
to $\eta_{\alpha\beta}$. We have (using Einstein notation and dropping
subscript $T$)

$\displaystyle\left(\frac{\partial^{2}A}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right)_{\bm{\mathcal{E}}=0}=\left(\frac{\partial^{2}A}{\partial\eta_{\tau\upsilon}\eta_{\phi\omega}}\frac{\partial\eta_{\tau\upsilon}}{\partial\varepsilon_{\alpha\beta}}\frac{\partial\eta_{\phi\omega}}{\partial\varepsilon_{\mu\nu}}\right.$
$\displaystyle\left.+\frac{\partial
A}{\partial\eta_{\tau\upsilon}}\frac{\partial^{2}\eta_{\tau\upsilon}}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right)_{\bm{\mathcal{E}}=0}.$
(36)

Using Eq. (12) we can evaluate the following derivatives

$\displaystyle\left(\frac{\partial\eta_{\tau\upsilon}}{\partial\varepsilon_{\alpha\beta}}\right)_{\bm{\mathcal{E}}=0}=\frac{1}{2}(\delta_{\tau\alpha}\delta_{\upsilon\beta}+\delta_{\tau\beta}\delta_{\upsilon\alpha}),$
$\displaystyle\left(\frac{\partial^{2}\eta_{\tau\upsilon}}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right)_{\bm{\mathcal{E}}=0}=\frac{1}{8}\left(\delta_{\tau\alpha}\delta_{\beta\mu}\delta_{\upsilon\nu}+\delta_{\tau\beta}\delta_{\alpha\mu}\delta_{\upsilon\nu}\right.$
$\displaystyle\left.+\delta_{\tau\alpha}\delta_{\beta\nu}\delta_{\upsilon\mu}+\delta_{\tau\beta}\delta_{\alpha\nu}\delta_{\upsilon\mu}+\delta_{\alpha\nu}\delta_{\upsilon\beta}\delta_{\tau\mu}+\delta_{\nu\beta}\delta_{\upsilon\alpha}\delta_{\tau\mu}\right.$
$\displaystyle\left.+\delta_{\alpha\mu}\delta_{\upsilon\beta}\delta_{\tau\nu}+\delta_{\mu\beta}\delta_{\upsilon\alpha}\delta_{\tau\nu}\right).$
(37)

Inserting these into Eq. (36) and using the fact that $\bm{\mathcal{H}}=0$
when $\bm{\mathcal{E}}=0$ we obtain

$\displaystyle\left(\frac{\partial^{2}A}{\partial\varepsilon_{\alpha\beta}\partial\varepsilon_{\mu\nu}}\right)_{\bm{\mathcal{E}}=0}=\left(\frac{\partial^{2}A}{\partial\eta_{\alpha\beta}\partial\eta_{\mu\nu}}\right)_{\bm{\mathcal{H}}=0}$
$\displaystyle+\frac{1}{4}\left(\frac{\partial
A}{\partial\eta_{\alpha\nu}}\delta_{\beta\mu}+\frac{\partial
A}{\partial\eta_{\beta\nu}}\delta_{\alpha\mu}+\frac{\partial
A}{\partial\eta_{\alpha\mu}}\delta_{\beta\nu}+\frac{\partial
A}{\partial\eta_{\beta\mu}}\delta_{\alpha\nu}\right)_{\bm{\mathcal{H}}=0}.$

Dividing both sides of the last equation by volume and using the fact that

$\frac{1}{V}\left(\frac{\partial
A}{\partial\eta_{\alpha\nu}}\right)_{\bm{\mathcal{H}}=0}=\frac{1}{V}\left(\frac{\partial
A}{\partial\varepsilon_{\alpha\nu}}\right)_{\bm{\mathcal{E}}=0}=\langle\sigma_{\alpha\beta}\rangle$
(39)

we obtain the sought-after relationship between the two sets of isothermal
elastic constants,

$\displaystyle C_{\alpha\beta\mu\nu}^{T}$ $\displaystyle=$
$\displaystyle\widehat{C}_{\alpha\beta\mu\nu}^{T}+\frac{1}{4}\big{(}\langle\sigma_{\alpha\nu}\rangle\delta_{\beta\mu}+\langle\sigma_{\beta\nu}\rangle\delta_{\alpha\mu}+\langle\sigma_{\alpha\mu}\rangle\delta_{\beta\nu}$
(40) $\displaystyle+\langle\sigma_{\beta\mu}\rangle\delta_{\alpha\nu}\big{)}.$

Applying a similar argument to the system energy $E$, one can verify that this
relationship remains true for adiabatic constants as well.

## Appendix B

Here we outline the derivation of Eq. (31) from Eq. (30). Before proceeding
with the derivation let us note that the potential energy $u_{k}$ of a group
of $N_{k}$ atoms is translationally invariant and, therefore, depends on at
most $3(N_{k}-1)$ variables. We want to account for this translational
invariance when doing integration by parts. This can be achieved by doing an
orthogonal transformation from coordinates ${\bf x}_{i}$ to ${\bf y}_{i}$. For
each group $k$ the transformation involves $N_{k}$ atomic-group coordinates
${\bf x}_{i}^{k}$ and transforms them to new $N_{k}$ coordinates ${\bf
y}_{i}^{k}$. The transformation within this group of $N_{k}$ coordinates can
be arbitrary apart from the fact that one of the new vector variables (assumed
to be the last here) is expressed through the arithmetic average of ${\bf
x}_{i}^{k}$, i. e. ${\bf
y}_{N_{k}}^{k}=\frac{1}{\sqrt{N_{k}}}\sum_{i=1}^{N_{k}}{\bf x}_{i}^{k}$. For
coordinates not involved in the group $k$, we take ${\bf y}_{i}={\bf x}_{i}$.
The following relationship between the old and new variables within the group
is used in the analysis below,

$\sum_{i=1}^{N_{k}}x_{i\alpha}^{k}\frac{\partial g(\\{{\bf
x}_{i}^{k}\\})}{\partial
x_{i\beta}^{k}}=\sum_{i=1}^{N_{k}}y_{i\alpha}^{k}\frac{\partial g(\\{{\bf
y}_{i}^{k}\\})}{\partial y_{i\beta}^{k}}.$ (41)

Here $g(\\{{\bf x}_{i}^{k}\\})$ is an arbitrary function of the old variables
${\bf x}_{i}^{k}$ and the same symbol is used for the transformed function
$g(\\{{\bf y}_{i}^{k}\\})$ of the new variables ${\bf y}_{i}^{k}$.
Transforming to the new variables and using Eq. (41) the double sum over atoms
in one group in the Born term (30) can be written as

$\displaystyle\sum_{i,j=1}^{N_{k}}$ $\displaystyle\left\langle
x_{i\alpha}^{k}x_{j\mu}^{k}\frac{\partial^{2}u_{k}}{\partial
x_{i\beta}^{k}\partial x_{j\nu}^{k}}\right\rangle_{\text{sym}}$ (42)
$\displaystyle=\sum_{i,j=1}^{N_{k}-1}\left\langle
y_{i\alpha}^{k}y_{j\mu}^{k}\frac{\partial^{2}u_{k}}{\partial
y_{i\beta}^{k}\partial y_{j\nu}^{k}}\right\rangle_{\text{sym}},$

where summations in the last expression now run from $1$ to $N_{k}-1$ because
$u_{k}$ is translationally invariant in the space of $N_{k}$ atomic-group
coordinates and, therefore, does not depend on ${\bf y}_{N_{k}}^{k}$. Consider
one term in the summand of the last expression with integrations written out
explicitly and the superscript $k$ and the subscript sym dropped for brevity:

$\left\langle y_{i\alpha}y_{j\mu}\frac{\partial^{2}u_{k}}{\partial
y_{i\beta}\partial y_{j\nu}}\right\rangle=\frac{1}{Q}\int d{\bf
y}^{3N}e^{-\frac{U}{k_{B}T}}y_{i\alpha}y_{j\mu}\frac{\partial^{2}u_{k}}{\partial
y_{i\beta}\partial y_{j\nu}},$ (43)

where $Q=\int d{\bf y}^{3N}e^{-\frac{U}{k_{B}T}}$. Integrating by parts with
respect to $y_{i\beta}$ gives

$\displaystyle\left\langle
y_{i\alpha}y_{j\mu}\frac{\partial^{2}u_{k}}{\partial y_{i\beta}\partial
y_{j\nu}}\right\rangle$ (44) $\displaystyle=$ $\displaystyle\frac{1}{Q}\int
d{\bf y}^{3N-1}e^{-\frac{U}{k_{B}T}}y_{i\alpha}y_{j\mu}\frac{\partial
u_{k}}{\partial y_{j\nu}}\bigg{|}_{y_{i\beta}=l_{1}}^{y_{i\beta}=l_{2}}$
$\displaystyle-\frac{1}{Q}\int d{\bf y}^{3N}\frac{\partial u_{k}}{\partial
y_{j\nu}}\frac{\partial}{\partial
y_{i\beta}}\left(e^{-\frac{U}{k_{B}T}}y_{i\alpha}y_{j\mu}\right)$
$\displaystyle=$ $\displaystyle\frac{1}{Q}\int d{\bf
y}^{3N-1}e^{-\frac{U}{k_{B}T}}y_{i\alpha}y_{j\mu}\frac{\partial
u_{k}}{\partial y_{j\nu}}\bigg{|}_{y_{i\beta}=l_{1}}^{y_{i\beta}=l_{2}}$
$\displaystyle-\frac{1}{Q}\int d{\bf y}^{3N}\frac{\partial u_{k}}{\partial
y_{j\nu}}e^{-\frac{U}{k_{B}T}}\bigg{(}-\frac{1}{k_{B}T}\frac{\partial
U}{\partial y_{i\beta}}y_{i\alpha}y_{j\mu}$
$\displaystyle+\delta_{\alpha\beta}y_{j\mu}+\delta_{ij}\delta_{\beta\mu}y_{i\alpha}\bigg{)}.$

Here $l_{1}$ and $l_{2}$ in the second and forth lines are the integration
limits for $y_{i\beta}$.

Up to this point the derivation is exact. The key step in deriving Eq. (31) is
the omission of the integrated part, i. e. the fourth line of Eq. (44). The
integrated part will vanish exactly in the thermodynamic limit (as the
integration limits tend to $\pm\infty$) for typical atomic-group potential
functions found in solids. More specifically, for potentials that decrease
with inter-atomic distance (such as the Lennard-Jones or Coulombic potentials)
$u_{k}$ and its derivatives vanish at infinity which makes the integrated part
(the first line of Eq. (44)) vanish exactly. For potentials that grow with
inter-atomic distance, such as the harmonic bond potential, the exponential
factor in the integrated part of Eq. (44) will tend to zero as the potential
itself goes to infinity at the boundaries, which, again, makes the integrated
part vanish exactly. If the system is finite, the integrated part may not
vanish exactly but its contribution will become progressively smaller as the
system size increases. In the following derivation we assume that the system
is large enough to be treated as being in the thermodynamic limit and omit the
integrated part.

Clearly, integration by parts in (43) can also be done with respect to
$y_{j\nu}$ with the result

$\displaystyle\left\langle
y_{i\alpha}y_{j\mu}\frac{\partial^{2}u_{k}}{\partial y_{i\beta}\partial
y_{j\nu}}\right\rangle$ (45) $\displaystyle=$ $\displaystyle-\frac{1}{Q}\int
d{\bf y}^{3N}\frac{\partial u_{k}}{\partial
y_{i\beta}}e^{-\frac{U}{k_{B}T}}\bigg{(}-\frac{1}{k_{B}T}\frac{\partial
U}{\partial y_{j\nu}}y_{i\alpha}y_{j\mu}$
$\displaystyle+\delta_{\mu\nu}y_{i\alpha}+\delta_{ij}\delta_{\alpha\nu}y_{j\mu}\bigg{)}.$

To maintain proper symmetry of elastic constants the average of (44) and (45)
has to be taken. Substituting this average into (42), performing summations,
where it can be done explicitly, and restoring the superscript $k$ and the
subscript sym we obtain

$\displaystyle\sum_{i,j=1}^{N_{k}-1}\left\langle
y_{i\alpha}^{k}y_{j\mu}^{k}\frac{\partial^{2}u_{k}}{\partial
y_{i\beta}^{k}\partial y_{j\nu}^{k}}\right\rangle_{\text{sym}}$ (46)
$\displaystyle=$
$\displaystyle\frac{1}{2k_{B}T}\sum_{i,j=1}^{N_{k}-1}\left\langle
y_{i\alpha}^{k}y_{j\mu}^{k}\left(\frac{\partial u_{k}}{\partial
y_{j\nu}^{k}}\frac{\partial U}{\partial y_{i\beta}^{k}}+\frac{\partial
u_{k}}{\partial y_{i\beta}^{k}}\frac{\partial U}{\partial
y_{j\nu}^{k}}\right)\right\rangle_{\text{sym}}$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{N_{k}-1}\left\langle\big{(}(N_{k}-1)\delta_{\alpha\beta}y_{i\mu}^{k}+\delta_{\beta\mu}y_{i\alpha}^{k}\big{)}\frac{\partial
u_{k}}{\partial y_{i\nu}^{k}}\right\rangle_{\text{sym}}$
$\displaystyle-\frac{1}{2}\sum_{i=1}^{N_{k}-1}\left\langle\big{(}(N_{k}-1)\delta_{\mu\nu}y_{i\alpha}^{k}+\delta_{\alpha\nu}y_{i\mu}^{k}\big{)}\frac{\partial
u_{k}}{\partial y_{i\beta}^{k}}\right\rangle_{\text{sym}}.$

We now transform Eq. (46) back to the original coordinates ${\bf x}_{i}^{k}$
using Eq. (41). The sums involving $u_{k}$ transform as follows

$\sum_{i=1}^{N_{k}-1}y_{i\alpha}^{k}\frac{\partial u_{k}}{\partial
y_{i\beta}^{k}}=\sum_{i=1}^{N_{k}}y_{i\alpha}^{k}\frac{\partial
u_{k}}{\partial
y_{i\beta}^{k}}=\sum_{i=1}^{N_{k}}x_{i\alpha}^{k}\frac{\partial
u_{k}}{\partial x_{i\beta}^{k}}.$ (47)

Here the summation in the second expression is extended to $N_{k}$ because
$u_{k}$ does not depend on ${\bf y}_{N_{k}}^{k}$. The total potential energy
$U$ in Eq. (46) does depend on ${\bf y}_{N_{k}}^{k}$, in general. Therefore,
when transforming back to ${\bf x}_{i}^{k}$,

$\displaystyle\sum_{i=1}^{N_{k}-1}y_{i\alpha}^{k}\frac{\partial U}{\partial
y_{i\beta}^{k}}=\left(\sum_{i=1}^{N_{k}}y_{i\alpha}^{k}\frac{\partial
U}{\partial y_{i\beta}^{k}}\right)-y_{N_{k}\alpha}^{k}\frac{\partial
U}{\partial y_{N_{k}\beta}^{k}}$
$\displaystyle=\left(\sum_{i=1}^{N_{k}}x_{i\alpha}^{k}\frac{\partial
U}{\partial
x_{i\beta}^{k}}\right)-\frac{1}{N_{k}}\left(\sum_{i=1}^{N_{k}}x_{i\alpha}^{k}\right)\left(\sum_{i=1}^{N_{k}}\frac{\partial
U}{\partial x_{i\beta}^{k}}\right).$

Here the last term involving the product of two sums is obtained by
transforming $y_{N_{k}\alpha}^{k}\partial U/{\partial y_{N_{k}\beta}^{k}}$ to
coordinates ${\bf x}_{i}^{k}$. Substituting Eqs. (47) and (B) into Eq. (46),
performing symmetrizations, summing over all groups, and dividing by the
volume lead to Eq. (31).

## Appendix C

In this Appendix we show that the Born term (30) calculated using the
microcanonical ensemble,

$\rho=\frac{1}{\Omega}\delta(E-H),$ (49)

where $\Omega=\int d{\bf p}^{3N}d{\bf x}^{3N}\delta(E-H)$, can also be
rewritten in the form of Eq. (31), in which all averages are taken with the
microcanonical ensemble as well. The initial steps are the same as in Appendix
B but with the canonical averages replaced with the microcanonical ones. Thus
we have the following microcanonical analogue of (43)

$\displaystyle\left\langle
y_{i\alpha}y_{j\mu}\frac{\partial^{2}u_{k}}{\partial y_{i\beta}\partial
y_{j\nu}}\right\rangle$ $\displaystyle=\frac{1}{\Omega}\int d{\bf p}^{3N}d{\bf
y}^{3N}\delta(E-H)y_{i\alpha}y_{j\mu}\frac{\partial^{2}u_{k}}{\partial
y_{i\beta}\partial y_{j\nu}},$ (50)

Performing integration by parts as in (44) and omitting the integrated term we
obtain

$\displaystyle-\frac{1}{\Omega}\int d{\bf p}^{3N}d{\bf y}^{3N}\frac{\partial
u_{k}}{\partial y_{j\nu}}\frac{\partial}{\partial
y_{i\beta}}\big{(}\delta(E-H)y_{i\alpha}y_{j\mu}\big{)}$ (51) $\displaystyle=$
$\displaystyle-\frac{1}{\Omega}\int d{\bf p}^{3N}d{\bf y}^{3N}\frac{\partial
u_{k}}{\partial y_{j\nu}}\bigg{(}\frac{\partial\delta(E-H)}{\partial
y_{i\beta}}y_{i\alpha}y_{j\mu}$
$\displaystyle+\delta(E-H)(\delta_{\alpha\beta}y_{j\mu}+\delta_{ij}\delta_{\beta\mu}y_{i\alpha})\bigg{)}.$

The key step now is to show that

$\displaystyle-\frac{1}{\Omega}\int d{\bf p}^{3N}d{\bf y}^{3N}\frac{\partial
u_{k}}{\partial y_{j\nu}}\frac{\partial\delta(E-H)}{\partial
y_{i\beta}}y_{i\alpha}y_{j\mu}$
$\displaystyle=\frac{1}{k_{B}T}\left\langle\frac{\partial u_{k}}{\partial
y_{j\nu}}\frac{\partial U}{\partial
y_{i\beta}}y_{i\alpha}y_{j\mu}\right\rangle,$ (52)

where the brackets now denote averaging over the microcanonical ensemble. Here
we follow the approach used in Ref. [34]. Let us denote $(\partial
u_{k}/\partial y_{j\nu})y_{i\alpha}y_{j\mu}$ in the first line of (52) by $G$.
We have

$\displaystyle-\frac{1}{\Omega(E)}\int d{\bf p}^{3N}d{\bf
y}^{3N}\frac{\partial\delta(E-H)}{\partial
y_{i\beta}}G=\frac{1}{\Omega(E)}\int d{\bf p}^{3N}d{\bf
y}^{3N}\frac{\partial\delta(E-H)}{\partial E}\frac{\partial U}{\partial
y_{i\beta}}G$ $\displaystyle=\frac{1}{\Omega(E)}\frac{\partial}{\partial
E}\int d{\bf p}^{3N}d{\bf
y}^{3N}\frac{\delta(E-H)}{\Omega(E)}\Omega(E)\frac{\partial U}{\partial
y_{i\beta}}G=\int d{\bf p}^{3N}d{\bf
y}^{3N}\frac{\delta(E-H)}{\Omega(E)}\left(\frac{\partial\ln\Omega(E)}{\partial
E}\right)\frac{\partial U}{\partial y_{i\beta}}G$
$\displaystyle+\frac{\partial}{\partial E}\int d{\bf p}^{3N}d{\bf
y}^{3N}\frac{\delta(E-H)}{\Omega(E)}\frac{\partial U}{\partial
y_{i\beta}}G=\frac{1}{k_{B}T}\left\langle\frac{\partial U}{\partial
y_{i\beta}}G\right\rangle+\frac{\partial}{\partial
E}\left\langle\frac{\partial U}{\partial y_{i\beta}}G\right\rangle.$ (53)

Here the second expression follows because $\delta(E-H)$ depends on
$y_{i\beta}$ only through the potential energy $U$; the third expression is
obtained by moving the derivative with respect to $E$ in front of the integral
and inserting $\Omega(E)/\Omega(E)=1$ in the integrand; the fourth and fifth
expressions follow because

$\displaystyle\frac{\partial}{\partial
E}\left(\frac{\delta(E-H)\Omega(E)}{\Omega(E)}\right)=\frac{\delta(E-H)}{\Omega(E)}\frac{\partial\Omega(E)}{\partial
E}$ $\displaystyle+\Omega(E)\frac{\partial}{\partial
E}\left(\frac{\delta(E-H)}{\Omega(E)}\right);$ (54)

and the last equality in Eq. (53) is obtained using the fact that
$\ln\Omega(E)=S(E)/k_{B}$ and ${(\partial S/\partial E)}_{V}=1/T$. If
$\langle(\partial U/\partial y_{i\beta})G\rangle$ in (53) is finite then the
last term in (53) vanishes in the thermodynamic limit because it represents a
derivative of $\langle(\partial U/\partial y_{i\beta})G\rangle$ with respect
to an extensive variable. This proves the relationship (52). The rest of the
proof follows the same steps as in Appendix B.

## References

  * Born and Huang [1954] M. Born and K. Huang,  _Dynamical Theory of Crystal Lattices_ (Oxford University Press, Oxford, 1954) p. 129.
  * Squire  _et al._ [1969] D. R. Squire, A. C. Holt, and W. G. Hoover, Physica 42, 388 (1969).
  * Parrinello and Rahman [1981] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
  * Parrinello and Rahman [1982] M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982).
  * Sprik  _et al._ [1984] M. Sprik, R. W. Impey, and M. L. Klein, Phys. Rev. B 29, 4368 (1984).
  * Ray [1988] J. R. Ray, Comput. Phys. Rep. 8, 109 (1988).
  * Ray and Rahman [1984] J. R. Ray and A. Rahman, J. Chem. Phys. 80, 4423 (1984).
  * Ray and Rahman [1985] J. R. Ray and A. Rahman, J. Chem. Phys. 82, 4243 (1985).
  * Lutsko [1989] J. F. Lutsko, J. Appl. Phys. 65, 2991 (1989).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022), also see https://www.lammps.org/.
  * Zhen and Chu [2012] Y. Zhen and C. Chu, Comput. Phys. Commun. 183, 261 (2012).
  * Gao  _et al._ [2006] G. Gao, K. Van Workum, J. D. Schall, and J. A. Harrison, J. Phys.: Condens. Matter 18, S1737 (2006).
  * Clavier  _et al._ [2017] G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend, and B. Rousseau, Mol. Simul. 43, 1413 (2017).
  * Hooks  _et al._ [2015] D. E. Hooks, K. J. Ramos, C. A. Bolme, and M. J. Cawkwell, Propellants Explos. Pyrotech. 40, 333 (2015).
  * Erpenbeck and Wood [1977] J. J. Erpenbeck and W. W. Wood, in  _Modern Theoretical Chemistry, Statistical Mechanics, Part B: Time Dependent Processes_ , edited by B. J. Berne (Plenum, New York, 1977) p. 6.
  * Todd and Daivis [2017] B. D. Todd and P. J. Daivis,  _Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications_ (Cambridge University Press, Cambridge, 2017).
  * Bekker  _et al._ [1995] H. Bekker, E. J. Dijkstra, M. K. R. Renardus, and H. J. C. Berendsen, Mol. Simul. 14, 137 (1995).
  * Thompson  _et al._ [2009] A. P. Thompson, S. J. Plimpton, and W. Mattson, J. Chem. Phys. 131, 154107 (2009).
  * Korn and Korn [2000] G. A. Korn and T. M. Korn,  _Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review_ (Dover Publications, Mineola, New York, 2000) p. 411.
  * Başar and Weichert [2000] Y. Başar and D. Weichert,  _Nonlinear Continuum Mechanics of Solids_ (Springer, Berlin, 2000) p. 63.
  * Wallace [1972] D. C. Wallace,  _Thermodynamics of Crystals_ (John Wiley & Sons, Inc., New York, 1972) p. 19.
  * Bogoliubov [1960] N. N. Bogoliubov,  _Problems of Dynamic Theory in Statistical Physics_ (Oak Ridge, Tennessee, 1960).
  * Tolman [1938] R. C. Tolman,  _The Principles of Statistical Mechanics_ (Clarendon Press, Oxford, 1938) p. 95.
  * Pathria [1996] R. K. Pathria,  _Statistical Mechanics_ , 2nd ed. (Butterworth-Heinemann, Oxford, 1996) p. 63.
  * Pereverzev and Sewell [2015] A. Pereverzev and T. D. Sewell, J. Chem. Phys. 142, 134110 (2015).
  * Ray and Graben [1981] J. R. Ray and H. W. Graben, Mol. Phys. 43, 1293 (1981).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Harcourt College Publishers, Forth Worth, 1976).
  * Bernardes [1958] N. Bernardes, Phys. Rev. 112, 1534 (1958).
  * Haynes  _et al._ [2016] W. M. Haynes, D. R. Lide, and T. J. Bruno,  _CRC handbook of Chemistry and Physics : a Ready-reference Book of Chemical and Physical Data_ , 97th ed. (CRC Press, Boca Raton, Florida, 2016).
  * Stishov and Fedosimov [1971] S. M. Stishov and V. I. Fedosimov, Pis’ma Zh. Eksp. Teor. Fiz. 14, 326 (1971), [Sov. Phys. JETP Lett. 14, 217 (1971)].
  * Flyvbjerg and Petersen [1989] H. Flyvbjerg and H. G. Petersen, J. Chem. Phys. 91, 461 (1989).
  * Kluge  _et al._ [1986] M. D. Kluge, J. R. Ray, and A. Rahman, J. Chem. Phys. 85, 4028 (1986).
  * Stillinger and Weber [1985] F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  * Zubarev [1974] D. N. Zubarev,  _Nonequilibrium Statistical Thermodynamics_ (Consultants Bureau, New York, 1974) p. 47.